IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16479


Ignore:
Timestamp:
Feb 14, 2008, 1:33:09 PM (18 years ago)
Author:
Paul Price
Message:

Rewinding state of psModules/src/imcombine to the last merge (pap_merge_080122), so that ppStack on the mainline will build (I had branched ppStack but made changes to psModules/src/imcombine on the mainline; this check in is restoring to the state of the mainline while the development proceeds on branch pap_branch_080214)

Location:
trunk/psModules/src/imcombine
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmPSFEnvelope.c

    r16398 r16479  
    3131
    3232
    33 //#define TESTING                         // Enable test output
     33#define TESTING                         // Enable test output
    3434#define PEAK_FLUX 1.0e3                 // Peak flux for each source
    3535#define SKY_VALUE 0.0e0                 // Sky value for fake image
    … …  
    9494            fake->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE);
    9595            fake->type = PM_SOURCE_TYPE_STAR;
     96
     97#ifdef TESTING
     98            // Required to generate model image
    9699            fake->psfMag = -2.5 * log10(PEAK_FLUX);
     100#endif
    97101
    98102            psTrace("psModules.imcombine", 5, "Source %d: %.2f,%.2f\n",

  • trunk/psModules/src/imcombine/pmStack.c

    r16420 r16479  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2008-02-13 02:55:33 $
     10 *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2008-02-14 23:33:09 $
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1313 *
    … …  
    486486        for (int j = 0; j < pixels->n; j++) {
    487487            int x = pixels->data[j].x, y = pixels->data[j].y; // Coordinates of interest
    488             if (x < 0 || x >= numCols || y < 0 || y >= numRows) {
    489                 psWarning("Bad pixel coordinate: %d,%d --- ignored.", x, y);
    490                 continue;
    491             }
    492488            psArray *columns = map->data[y]; // The columns for that row
    493489            psVector *images = columns->data[x]; // The images for that column

  • trunk/psModules/src/imcombine/pmStackReject.c

    r16476 r16479  
    1111#define PIXEL_LIST_BUFFER 100           // Number of pixels to add to list at a time
    1212
    13 psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold,
    14                         const psArray *subRegions, const psArray *kernels)
     13psPixels *pmStackReject(const psPixels *in, float threshold, const psArray *regions,
     14                        const psArray *solutions, const pmSubtractionKernels *kernels)
    1515{
    1616    PS_ASSERT_PIXELS_NON_NULL(in, NULL);
    1717    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(threshold, 0.0, NULL);
    1818    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL);
    19     PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL);
    20     PS_ASSERT_ARRAY_NON_NULL(kernels, NULL);
    21     PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL);
     19    PS_ASSERT_ARRAY_NON_NULL(regions, NULL);
     20    PS_ASSERT_ARRAY_NON_NULL(solutions, NULL);
     21    PS_ASSERT_ARRAYS_SIZE_EQUAL(regions, solutions, NULL);
     22    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
    2223
    2324    // Get the original image size
    24     int numRegions = subRegions->n;        // Number of regions
     25    int numRegions = regions->n;        // Number of regions
    2526    int numCols = 0, numRows = 0;       // Size of original image
    2627    int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image --- should be 0,0
    27     int size = 0;                       // Size of kernel
    2828    for (int i = 0; i < numRegions; i++) {
    29         psRegion *subRegion = subRegions->data[i]; // Region of interest
    30         if (subRegion->x0 < minCols) {
    31             minCols = subRegion->x0;
     29        psRegion *region = regions->data[i]; // Region of interest
     30        if (region->x0 < minCols) {
     31            minCols = region->x0;
    3232        }
    33         if (subRegion->y0 < minRows) {
    34             minRows = subRegion->y0;
     33        if (region->y0 < minRows) {
     34            minRows = region->y0;
    3535        }
    36         if (subRegion->x1 > numCols) {
    37             numCols = subRegion->x1;
     36        if (region->x1 > numCols) {
     37            numCols = region->x1;
    3838        }
    39         if (subRegion->y1 > numRows) {
    40             numRows = subRegion->y1;
    41         }
    42 
    43         pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
    44         if (size == 0) {
    45             size = kernel->size;
    46         } else if (kernel->size != size) {
    47             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d",
    48                     size, kernel->size);
    49             return NULL;
     39        if (region->y1 > numRows) {
     40            numRows = region->y1;
    5041        }
    5142    }
    52 
    53     // Adjust the size for the size of the subimage
    54     if (valid) {
    55         minCols = PS_MAX(valid->x0, minCols);
    56         minRows = PS_MAX(valid->y0, minRows);
    57         numCols = PS_MIN(valid->x1, numCols);
    58         numRows = PS_MIN(valid->y1, numRows);
     43    if (minCols != 0 || minRows != 0) {
     44        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     45                "Some error with image regions --- minimum coordinate is not 0,0");
     46        return NULL;
    5947    }
    6048
    61     psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1),
    62                                    0x01); // Mask
     49    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image
    6350    psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
    6451    psFree(mask);
    … …  
    6855    pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
    6956    inRO->image = image;
    70     inRO->col0 = minCols;
    71     inRO->row0 = minRows;
    7257    for (int i = 0; i < numRegions; i++) {
    73         psRegion *region = subRegions->data[i]; // Region of interest
    74         if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 ||
    75                       region->y0 > valid->y1 || region->y1 < valid->y0)) {
    76             // Region is outside of our sub-image
    77             continue;
    78         }
    79         pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
    80         if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, true)) {
     58        psRegion *region = regions->data[i]; // Region of interest
     59        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernels, true)) {
    8160            psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
    8261            psFree(convRO);
    … …  
    8968
    9069        // Image of the kernel at the centre of the region
    91         float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) /
    92             (float)kernel->numCols;
    93         float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) /
    94             (float)kernel->numRows;
    95         psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false);
    96         if (!image) {
     70        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - numCols/2.0) / (float)numCols;
     71        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - numRows/2.0) / (float)numRows;
     72        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
     73        if (!kernel) {
    9774            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
    9875            psFree(convRO);
    … …  
    10178        }
    10279        float sum = 0.0;
    103         for (int y = 0; y < image->numRows; y++) {
    104             for (int x = 0; x < image->numCols; x++) {
    105                 sum += image->data.F32[y][x];
     80        for (int y = 0; y < kernel->numRows; y++) {
     81            for (int x = 0; x < kernel->numCols; x++) {
     82                sum += kernel->data.F32[y][x];
    10683            }
    10784        }
    108         psFree(image);
     85        psFree(kernel);
    10986
    110         // Range for normalisation
    111         int yMin = PS_MAX(minRows, region->y0) - inRO->row0;
    112         int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0;
    113         int xMin = PS_MAX(minCols, region->x0) - inRO->col0;
    114         int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0;
    115         psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n",
    116                 sum, xMin, xMax, yMin, yMax);
    117         for (int y = yMin; y <= yMax; y++) {
    118             for (int x = xMin; x <= xMax; x++) {
    119                 convRO->image->data.F32[y][x] /= sum;
    120             }
    121         }
     87        psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image
     88        psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
    12289    }
    12390    psFree(inRO);
    … …  
    12794    // Threshold the convolved image
    12895    psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked
    129     for (int y = 0; y < convolved->numRows; y++) {
    130         for (int x = 0; x < convolved->numCols; x++) {
     96    for (int y = 0; y < numRows; y++) {
     97        for (int x = 0; x < numCols; x++) {
    13198            if (convolved->data.F32[y][x] > threshold) {
    13299                bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x, y);
    … …  
    136103    psFree(convolved);
    137104
    138     // Now, grow the mask to include everything that touches a bad pixel in the convolution
    139     mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0xff);
    140     assert(mask->numCols == numCols && mask->numRows == numRows);
     105    // Now, we want to convolve the original pixels properly
     106    mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols, 0, numRows), 0xff);
     107    int size = kernels->size;           // Size of kernels
    141108    for (int i = 0; i < bad->n; i++) {
    142109        int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest
    … …  
    148115        for (int y = yMin; y <= yMax; y++) {
    149116            for (int x = xMin; x <= xMax; x++) {
    150                 assert(x < mask->numCols && y < mask->numRows);
    151117                mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff;
    152118            }

  • trunk/psModules/src/imcombine/pmStackReject.h

    r16423 r16479  
    1010/// We apply a matched filter to the corresponding mask image, and threshold to find the original pixels
    1111psPixels *pmStackReject(const psPixels *in, ///< List of pixels in the convolved image
    12                         const psRegion *valid, ///< Valid region to consider
    1312                        float threshold, ///< Threshold on convolved image, 0..1
    1413                        const psArray *regions, ///< Array of image regions for image
    15                         const psArray *kernels ///< Array of kernel parameters for each region
     14                        const psArray *solutions, ///< Array of solution vectors for image
     15                        const pmSubtractionKernels *kernels ///< Kernel parameters
    1616    );
    1717

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r16476 r16479  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2008-02-14 23:18:34 $
     6 *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2008-02-14 23:33:09 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    5050
    5151    // Take the square of the normal kernel
     52    double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels
    5253    for (int v = yMin; v <= yMax; v++) {
    5354        for (int u = xMin; u <= xMax; u++) {
    54             out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
    55         }
    56     }
     55            float value = normalKernel->kernel[v][u]; // Value of interest
     56            float value2 = PS_SQR(value); // Value squared
     57            sumNormal += value;
     58            sumVariance += value2;
     59            out->kernel[v][u] = value2;
     60        }
     61    }
     62
     63    // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
     64    // This is required to keep the relative scaling between the image and the weight map
     65    psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
    5766
    5867    return out;
    … …  
    681690        PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false);
    682691    }
    683     if (region && psRegionIsNaN(*region)) {
    684         psString string = psRegionToString(*region);
    685         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
    686         psFree(string);
    687         return false;
     692    if (region) {
     693        if (psRegionIsNaN(*region)) {
     694            psString string = psRegionToString(*region);
     695            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
     696            psFree(string);
     697            return false;
     698        }
     699        if (region->x0 < 0 || region->x1 > ro1->image->numCols ||
     700            region->y0 < 0 || region->y1 > ro1->image->numRows) {
     701            psString string = psRegionToString(*region);
     702            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
     703                    string, ro1->image->numCols, ro1->image->numRows);
     704            psFree(string);
     705            return false;
     706        }
    688707    }
    689708
    … …  
    705724
    706725    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
    707     int x0 = ro1->col0, y0 = ro1->row0; // Image offset
    708726
    709727    psImage *convImage1 = out1->image;   // Convolved image
    … …  
    762780        yMax = PS_MIN(region->y1, yMax);
    763781    }
    764 
    765     // Size to use when calculating normalised coordinates (different from actual size when convolving
    766     // subimage)
    767     int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
    768     int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
    769782
    770783    psMaskType maskSource;              // Mask these pixels when convolving
    … …  
    794807    for (int j = yMin; j < yMax; j += fullSize) {
    795808        int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
    796         float yNorm = 2.0 * (float)(j + y0 + size + 1 - yNormSize/2.0) /
    797             (float)yNormSize; // Normalised coordinate
     809        float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // Normalised coordinate
    798810        for (int i = xMin; i < xMax; i += fullSize) {
    799811            int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
    800             float xNorm = 2.0 * (float)(i + x0 + size + 1 - xNormSize/2.0) /
    801                 (float)xNormSize; // Normalised coordinate
     812            float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate
    802813
    803814            // Only generate polynomial values every kernel footprint, since we have already assumed

  • trunk/psModules/src/imcombine/pmSubtractionEquation.c

    r16384 r16479  
    1414
    1515
    16 //#define TESTING
     16#define TESTING
    1717
    1818//////////////////////////////////////////////////////////////////////////////////////////////////////////////

  • trunk/psModules/src/imcombine/pmSubtractionKernels.c

    r16423 r16479  
    180180    kernels->bgOrder = 0;
    181181    kernels->mode = mode;
    182     kernels->numCols = 0;
    183     kernels->numRows = 0;
    184182    kernels->solution1 = NULL;
    185183    kernels->solution2 = NULL;

  • trunk/psModules/src/imcombine/pmSubtractionKernels.h

    r16423 r16479  
    3939    int bgOrder;                        ///< The order for the background fitting
    4040    pmSubtractionMode mode;             ///< Mode for subtraction
    41     int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
    4241    psVector *solution1, *solution2;    ///< Solution for the PSF matching
    4342} pmSubtractionKernels;

  • trunk/psModules/src/imcombine/pmSubtractionMatch.c

    r16406 r16479  
    306306                                                       inner, binning, ringsOrder, mode);
    307307            }
    308 
    309             // Add analysis metadata
    310             {
    311                 psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
    312                                  PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
    313                 if (conv2) {
    314                     psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
    315                                      PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
    316                 }
    317                 psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
    318                                  PS_META_DUPLICATE_OK, "Subtraction kernels", mode);
    319                 if (conv2) {
    320                     psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
    321                                      PS_META_DUPLICATE_OK, "Subtraction kernels", mode);
    322                 }
    323                 psRegion *subRegion;
    324                 if (region) {
    325                     subRegion = psMemIncrRefCounter(region);
    326                 } else {
    327                     subRegion = psRegionAlloc(0, numCols, 0, numRows);
    328                 }
    329                 psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    330                                   PS_DATA_REGION | PS_META_DUPLICATE_OK,
    331                                  "Region over which subtraction was performed", subRegion);
    332                 if (conv2) {
    333                     psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    334                                      PS_DATA_REGION | PS_META_DUPLICATE_OK,
    335                                      "Region over which subtraction was performed", subRegion);
    336                 }
    337                 psFree(subRegion);
    338             }
     308            psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL",
     309                             PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
    339310
    340311            memCheck("kernels");
    … …  
    470441            kernels = NULL;
    471442
     443#if 0
     444            // Put the solution on the metadata
     445            {
     446                psString comment = NULL; // Comment for metadata
     447                psStringAppend(&comment, "Subtraction solution for region %s", regionString);
     448                psMetadataAddVector(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION",
     449                                    PS_META_DUPLICATE_OK, comment, solution);
     450                psFree(comment);
     451                if (region) {
     452                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
     453                                     PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
     454                } else {
     455                    region = psRegionAlloc(0, numCols, 0, numRows);
     456                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
     457                                     PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
     458                    psFree(region);
     459                    region = NULL;
     460                }
     461            }
     462#endif
     463
    472464            // There is data in the readout now
    473465            conv1->data_exists = true;

  • trunk/psModules/src/imcombine/pmSubtractionMatch.h

    r16401 r16479  
    88#include <pmSubtractionKernels.h>
    99#include <pmSubtractionStamps.h>
    10 
    11 #define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Name of kernel in the analysis metadata
    12 #define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Name of subtraction mode in the analysis metadata
    13 #define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Name of subtraction region in the analysis MD
    14 
    1510
    1611/// Match two images
Note: See TracChangeset for help on using the changeset viewer.