IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26739


Ignore:
Timestamp:
Jan 29, 2010, 6:02:38 PM (16 years ago)
Author:
Paul Price
Message:

Fairly extensive change to make subtractions more stable on images with limited lit area in the presence of spatial variation. Took the polynomials which were scaled over the entire image and scale them only over the lit area. Removed some fprintf statements. Fixed the problem causing the background term to always come out 1.0.

Location:
branches/eam_branches/20091201
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/ppSub/src/ppSubThreshold.c

    r24862 r26739  
    2727    psImageMaskType maskIgnore,         // Ignore pixels with this mask
    2828    psImageMaskType maskThresh,         // Give pixels this mask if below threshold
     29    psRegion *region,                   // Region of interest
    2930    const char *description             // Description of image
    3031    )
    … …  
    3738    psImage *image = ro->image;         // Image
    3839    psImage *mask = ro->mask;           // Mask
    39     int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image
     40
     41    psImage *subImage = psImageSubset(image, *region); // Image with region of interest
     42    psImage *subMask = psImageSubset(mask, *region);  // Maks with region of interest
    4043
    4144    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
    4245    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);                               // Random number generator
    43     if (!psImageBackground(stats, NULL, image, mask, maskIgnore, rng)) {
     46    if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) {
    4447        psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold.");
    4548        psFree(rng);
    … …  
    4952    psFree(rng);
    5053
     54    psFree(subImage);
     55    psFree(subMask);
     56
    5157    float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip
    5258    psFree(stats);
    5359    psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description);
    5460
    55     for (int y = 0; y < numRows; y++) {
    56         for (int x = 0; x < numCols; x++) {
     61    for (int y = region->y0; y < region->y1; y++) {
     62        for (int x = region->x0; x < region->x1; x++) {
    5763            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) {
    5864                continue;
    … …  
    94100        return false;
    95101    }
    96     if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) {
    97         psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
    98         return false;
    99     }
    100102
    101103    pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image
    … …  
    104106        return false;
    105107    }
    106     if (!lowThreshold(ref, thresh, maskVal, maskThresh, "reference convolved image")) {
    107         psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
    108         return false;
     108
     109    psMetadataIterator *regIter = psMetadataIteratorAlloc(in->analysis, PS_LIST_HEAD,
     110                                                          "^" PM_SUBTRACTION_ANALYSIS_REGION "$");
     111    psMetadataItem *regItem;        // Item with region
     112    while ((regItem = psMetadataGetAndIncrement(regIter))) {
     113        psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type");
     114        psRegion *region = regItem->data.V; // Region of interest
     115        if (!lowThreshold(in, thresh, maskVal, maskThresh, region, "input convolved image")) {
     116            psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
     117            return false;
     118        }
     119        if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) {
     120            psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
     121            return false;
     122        }
    109123    }
     124    psFree(regIter);
    110125
    111126    psFree(view);

  • branches/eam_branches/20091201/psModules/src/imcombine/pmStackReject.c

    r25468 r26739  
    3535{
    3636    int size = kernels->size;           // Half-size of convolution kernel
    37     psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
    38                                                               xMin + size + 1, yMin + size + 1); // Polynomial
     37    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + size + 1,
     38                                                              yMin + size + 1); // Polynomial
    3939    int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues, false, poorFrac); // Radius of bad box
    4040    psTrace("psModules.imcombine", 10, "Growing by %d", box);
    … …  
    165165
    166166        // Image of the kernel at the centre of the region
    167         float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) /
    168             (float)kernels->numCols;
    169         float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) /
    170             (float)kernels->numRows;
     167        float xNorm, yNorm;             // Normalised coordinates
     168        p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, 0.5 * (region->x1 - region->x0),
     169                                            0.5 * (region->y1 - region->y0),
     170                                            kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax);
    171171        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
    172172        if (!kernel) {

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c

    r26667 r26739  
    424424                } else if (*source & subConvPoor) {
    425425                    *target |= maskPoor;
     426                } else {
     427                    *target &= ~maskBad & ~maskPoor;
    426428                }
    427429            }
    … …  
    574576}
    575577
     578void p_pmSubtractionPolynomialNormCoords(float *xOut, float *yOut, float xIn, float yIn,
     579                                         int xMin, int xMax, int yMin, int yMax)
     580{
     581    float xNormSize = xMax - xMin, yNormSize = yMax - yMin; // Size to use for normalisation
     582    *xOut = 2.0 * (float)(xIn - xMin - xNormSize/2.0) / xNormSize;
     583    *yOut = 2.0 * (float)(yIn - yMin - yNormSize/2.0) / yNormSize;
     584    return;
     585}
     586
    576587psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels,
    577                                              int numCols, int numRows, int x, int y)
     588                                             int x, int y)
    578589{
    579590    assert(kernels);
    580     assert(numCols > 0 && numRows > 0);
    581 
    582     // Size to use when calculating normalised coordinates (different from actual size when convolving
    583     // subimage)
    584     int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
    585     int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
    586 
    587     // Normalised coordinates
    588     float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize;
    589     float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize;
    590 
     591
     592    float xNorm, yNorm;                 // Normalised coordinates
     593    p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, x, y,
     594                                        kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax);
    591595    return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm);
    592596}
    … …  
    10721076    // Only generate polynomial values every kernel footprint, since we have already assumed
    10731077    // (with the stamps) that it does not vary rapidly on this scale.
    1074     psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
    1075                                                               xMin + x0 + size + 1,
    1076                                                               yMin + y0 + size + 1);
     1078    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + x0 + size + 1,
     1079                                                              yMin + y0 + size + 1);        // Polynomial
    10771080    float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term
    10781081
    … …  
    12001203    bool threaded = pmSubtractionThreaded(); // Running threaded?
    12011204
    1202     // Outputs
    1203     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    1204         if (!out1->image) {
    1205             out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    1206         }
    1207         if (ro1->variance) {
    1208             if (!out1->variance) {
    1209                 out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    1210             }
    1211             psImageInit(out1->variance, 0.0);
    1212         }
    1213     }
    1214     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    1215         if (!out2->image) {
    1216             out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    1217         }
    1218         if (ro2->variance) {
    1219             if (!out2->variance) {
    1220                 out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    1221             }
    1222             psImageInit(out2->variance, 0.0);
    1223         }
    1224     }
    12251205    psImage *convMask = NULL;           // Convolved mask image (common to inputs 1 and 2)
    12261206    if (subMask) {
    12271207        if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    1228             if (!out1->mask) {
    1229                 out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    1230             }
    1231             psImageInit(out1->mask, 0);
    12321208            convMask = out1->mask;
    12331209        }
    12341210        if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    1235             if (!out2->mask) {
    1236                 out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    1237             }
    1238             psImageInit(out2->mask, 0);
    12391211            if (!convMask) {
    12401212                convMask = out2->mask;
    … …  
    12561228
    12571229    // Get region for convolution: [xMin:xMax,yMin:yMax]
    1258     int xMin = size, xMax = numCols - size;
    1259     int yMin = size, yMax = numRows - size;
     1230    int xMin = kernels->xMin + size, xMax = kernels->xMax - size;
     1231    int yMin = kernels->yMin + size, yMax = kernels->yMax - size;
    12601232    if (region) {
    12611233        xMin = PS_MAX(region->x0, xMin);

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.h

    r25279 r26739  
    127127    );
    128128
     129/// Return normalised coordinates
     130void p_pmSubtractionPolynomialNormCoords(
     131    float *xOut, float *yOut,           ///< Normalised coordinates, returned
     132    float xIn, float yIn,               ///< Input coordinates
     133    int xMin, int xMax, int yMin, int yMax ///< Bounds of validity
     134    );
     135
    129136/// Given (normalised) coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
    130137psImage *p_pmSubtractionPolynomial(psImage *output, ///< Output matrix, or NULL
    … …  
    138145psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, ///< Output matrix, or NULL
    139146                                             const pmSubtractionKernels *kernels, ///< Kernel parameters
    140                                              int numCols, int numRows, ///< Size of image of interest
    141147                                             int x, int y ///< Position of interest
    142148    );

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c

    r26593 r26739  
    118118
    119119    // sample difference images
    120     { 
     120    {
    121121        psMetadataAddArray(analysis, PS_LIST_TAIL, "SUBTRACTION.SAMPLE.STAMP.SET", PS_META_DUPLICATE_OK, "Sample Difference Stamps", kernels->sampleStamps);
    122122    }
    … …  
    273273    // Difference in background
    274274    {
    275         psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
    276                                                                   numCols / 2.0, numRows / 2.0); // Polynomial
     275        psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0); // Polynomial
    277276        float bg = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background difference
    278277

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c

    r26737 r26739  
    213213
    214214    *norm = normI2 / normI1;
    215     fprintf(stderr, "Sums: %f %f %f, normWindow: %d\n", normI1, normI2, *norm, normWindow);
    216215
    217216    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    … …  
    522521
    523522    *norm = normI2 / normI1;
    524     fprintf(stderr, "Sums: %f %f %f - I1I1: %f\n", normI1, normI2, *norm, sumI1I1);
    525523
    526524    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    … …  
    590588                // Contribution to chi^2: a_i^2 P_i
    591589                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
    592                 fprintf (stderr, "penalty main: %d %e vs %e x %f : %f\n",
    593                          index,
    594                          matrix->data.F64[index][index],
    595                          norm, penalties->data.F32[i],
    596                          (norm * penalties->data.F32[i]) / matrix->data.F64[index][index] );
    597590                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
    598591                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    599                     fprintf (stderr, "penalty dual: %d %e vs %e x %f : %f\n",
    600                              index + numParams + 2,
    601                              matrix->data.F64[index + numParams + 2][index + numParams + 2],
    602                              norm, penalties->data.F32[i],
    603                              (norm * penalties->data.F32[i]) / matrix->data.F64[index + numParams + 2][index + numParams + 2] );
    604592                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
    605593                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
    … …  
    997985            // Solve normalisation and background separately
    998986            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     987            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    999988
    1000989            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
    … …  
    10141003
    10151004            fprintf(stderr, "Norm: %lf\n", normValue);
    1016             // fprintf(stderr, "BG: %lf\n", bgValue);
    10171005
    10181006            // Solve kernel components
    10191007            for (int i = 0; i < numSolution1; i++) {
    1020                 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
     1008                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
    10211009
    10221010                sumMatrix->data.F64[i][normIndex] = 0.0;
    10231011                sumMatrix->data.F64[normIndex][i] = 0.0;
    1024 
    1025                 // sumMatrix->data.F64[i][bgIndex] = 0.0;
    1026                 // sumMatrix->data.F64[bgIndex][i] = 0.0;
    1027             }
     1012            }
     1013            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
     1014            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1015            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
    10281016
    10291017            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
    1030             // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
    1031             // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
    1032 
    10331018            sumVector->data.F64[normIndex] = 0.0;
    1034             // sumVector->data.F64[bgIndex] = 0.0;
    10351019
    10361020            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    10371021
    10381022            solution->data.F64[normIndex] = normValue;
    1039             // solution->data.F64[bgIndex] = bgValue;
    10401023        }
    10411024# endif
    10421025
    10431026        if (!kernels->solution1) {
    1044             kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
    1045             psVectorInit (kernels->solution1, 0.0);
     1027            kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1028            psVectorInit(kernels->solution1, 0.0);
    10461029        }
    10471030
    … …  
    11261109            // Solve normalisation and background separately
    11271110            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1111            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    11281112
    11291113#if 0
    … …  
    11591143
    11601144            double normValue = stats->robustMedian;
    1161             // double bgValue = 0.0;
    11621145
    11631146            psFree(stats);
    11641147
    11651148            fprintf(stderr, "Norm: %lf\n", normValue);
    1166             // fprintf(stderr, "BG: %lf\n", bgValue);
    11671149
    11681150            // Solve kernel components
    11691151            for (int i = 0; i < numSolution2; i++) {
    1170                 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
    1171                 sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; // + bgValue * sumMatrix->data.F64[bgIndex][i + numSolution1];
     1152                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
     1153                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1];
    11721154
    11731155                sumMatrix->data.F64[i][normIndex] = 0.0;
    11741156                sumMatrix->data.F64[normIndex][i] = 0.0;
    11751157
    1176                 // sumMatrix->data.F64[i][bgIndex] = 0.0;
    1177                 // sumMatrix->data.F64[bgIndex][i] = 0.0;
    11781158                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
    11791159                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
    1180                 // sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0;
    1181                 // sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0;
    1182             }
     1160            }
     1161            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
     1162            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1163            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
    11831164
    11841165            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
    1185             // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
    1186             // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
    11871166
    11881167            sumVector->data.F64[normIndex] = 0.0;
    1189             // sumVector->data.F64[bgIndex] = 0.0;
    11901168
    11911169            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    11921170
    11931171            solution->data.F64[normIndex] = normValue;
    1194             // solution->data.F64[bgIndex] = bgValue;
    11951172        }
    11961173#endif

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h

    r26577 r26739  
    1010    PM_SUBTRACTION_EQUATION_BG      = 0x02,
    1111    PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
    12     PM_SUBTRACTION_EQUATION_ALL     = 0x05, // value should be NORM | BG | KERNELS
     12    PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
    1313} pmSubtractionEquationCalculationMode;
    1414
    … …  
    2121                                         const pmSubtractionKernels *kernels, ///< Kernel parameters
    2222                                         int index, ///< Index of stamp
    23                                         const pmSubtractionEquationCalculationMode mode
     23                                        const pmSubtractionEquationCalculationMode mode
    2424    );
    2525
    … …  
    2727bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
    2828                                    const pmSubtractionKernels *kernels, ///< Kernel parameters
    29                                     const pmSubtractionEquationCalculationMode mode
     29                                    const pmSubtractionEquationCalculationMode mode
    3030    );
    3131
    … …  
    3333bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
    3434                                const pmSubtractionStampList *stamps, ///< Stamps
    35                                 const pmSubtractionEquationCalculationMode mode
     35                                const pmSubtractionEquationCalculationMode mode
    3636    );
    3737

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionIO.c

    r23378 r26739  
    8080        while ((item = psMetadataGetAndIncrement(iter))) {
    8181            assert(item->type == PS_DATA_UNKNOWN);
    82             // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
    83             // images by scans.
    8482            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
    85             kernel->numCols = ro->image->numCols;
    86             kernel->numRows = ro->image->numRows;
    87 
    8883            kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
    8984        }
    … …  
    126121                         kernel->bgOrder);
    127122        psMetadataAddS32(row, PS_LIST_TAIL, NAME_MODE,  0, "Matching mode (enum)", kernel->mode);
    128         psMetadataAddS32(row, PS_LIST_TAIL, NAME_COLS,  0, "Number of columns", kernel->numCols);
    129         psMetadataAddS32(row, PS_LIST_TAIL, NAME_ROWS,  0, "Number of rows", kernel->numRows);
    130123        if (kernel->mode == PM_SUBTRACTION_MODE_1 || kernel->mode == PM_SUBTRACTION_MODE_2) {
    131124            psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL1, 0, "Solution vector 1", kernel->solution1);
    … …  
    318311
    319312        TABLE_LOOKUP(int, S32, bg,      NAME_BG);
    320         TABLE_LOOKUP(int, S32, numCols, NAME_COLS);
    321         TABLE_LOOKUP(int, S32, numRows, NAME_ROWS);
    322313
    323314        TABLE_LOOKUP(float, F32, mean,      NAME_MEAN);
    … …  
    325316        TABLE_LOOKUP(int,   S32, numStamps, NAME_NUMSTAMPS);
    326317
    327         pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode);
    328         kernels->numCols = numCols;
    329         kernels->numRows = numRows;
     318        pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, *region, mode);
    330319        kernels->mean = mean;
    331320        kernels->rms = rms;

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

    r26734 r26739  
    197197    }
    198198
    199 #if 1
     199#if 0
    200200    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
    201201#endif
    … …  
    245245    }
    246246
    247 #if 1
     247#if 0
    248248    {
    249249        double sum = 0.0;   // Sum of kernel component
    … …  
    276276pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
    277277                                                    const psVector *fwhmsIN, const psVector *ordersIN,
    278                                                     float penalty, pmSubtractionMode mode)
     278                                                    float penalty, psRegion bounds, pmSubtractionMode mode)
    279279{
    280280    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    305305    }
    306306
    307     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, spatialOrder, penalty, mode); // The kernels
     307    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
     308                                                              spatialOrder, penalty, bounds, mode); // Kernels
    308309    psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    309310
    … …  
    331332pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder,
    332333                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    333                                                       float penalty, pmSubtractionMode mode)
     334                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
    334335{
    335336    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    361362    }
    362363
    363     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, spatialOrder, penalty, mode); // The kernels
     364    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size,
     365                                                              spatialOrder, penalty, bounds, mode); // Kernels
    364366    psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    365367
    … …  
    388390pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
    389391                                               const psVector *fwhmsIN, const psVector *ordersIN,
    390                                                float penalty, pmSubtractionMode mode)
     392                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    391393{
    392394    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    417419    }
    418420
    419     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, spatialOrder, penalty, mode); // The kernels
     421    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size,
     422                                                              spatialOrder, penalty, bounds, mode); // Kernels
    420423    psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    421424
    422     psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
     425    psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements",
     426             params, spatialOrder, num);
    423427    psFree(params);
    424428
    … …  
    439443
    440444pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
    441                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    442                                                      float penalty, pmSubtractionMode mode)
     445                                                      const psVector *fwhmsIN, const psVector *ordersIN,
     446                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
    443447{
    444448    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    469473    }
    470474
    471     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, spatialOrder, penalty, mode); // The kernels
     475    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size,
     476                                                              spatialOrder, penalty, bounds, mode); // Kernels
    472477    psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    473478
    … …  
    544549
    545550pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
    546                                                 int size, int spatialOrder, float penalty,
     551                                                int size, int spatialOrder, float penalty, psRegion bounds,
    547552                                                pmSubtractionMode mode)
    548553{
    … …  
    558563    kernels->uStop = NULL;
    559564    kernels->vStop = NULL;
     565    kernels->xMin = bounds.x0;
     566    kernels->xMax = bounds.x1;
     567    kernels->yMin = bounds.y0;
     568    kernels->yMax = bounds.y1;
    560569    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    561570    kernels->penalty = penalty;
    … …  
    566575    kernels->bgOrder = 0;
    567576    kernels->mode = mode;
    568     kernels->numCols = 0;
    569     kernels->numRows = 0;
    570577    kernels->solution1 = NULL;
    571578    kernels->solution2 = NULL;
    … …  
    640647}
    641648
    642 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
     649pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, psRegion bounds,
    643650                                               pmSubtractionMode mode)
    644651{
    … …  
    649656
    650657    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
    651                                                               spatialOrder, penalty, mode); // The kernels
     658                                                              spatialOrder, penalty, bounds, mode); // Kernels
    652659    psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
    653660    psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
    … …  
    664671pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
    665672                                               const psVector *fwhms, const psVector *orders,
    666                                                float penalty, pmSubtractionMode mode)
     673                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    667674{
    668675    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    669                                                                   penalty, mode); // Kernels
     676                                                                  penalty, bounds, mode); // Kernels
    670677    if (!kernels) {
    671678        return NULL;
    … …  
    676683
    677684pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
    678                                                float penalty, pmSubtractionMode mode)
     685                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    679686{
    680687    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    697704
    698705    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
    699                                                               spatialOrder, penalty, mode); // The kernels
     706                                                              spatialOrder, penalty, bounds, mode); // Kernels
    700707    kernels->inner = inner;
    701708    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
    … …  
    768775
    769776pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
    770                                                 pmSubtractionMode mode)
     777                                                psRegion bounds, pmSubtractionMode mode)
    771778{
    772779    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    795802
    796803    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
    797                                                               spatialOrder, penalty, mode); // The kernels
     804                                                              spatialOrder, penalty, bounds, mode); // Kernels
    798805    kernels->inner = inner;
    799806    psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
    … …  
    864871pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
    865872                                               const psVector *orders, int inner, float penalty,
    866                                                pmSubtractionMode mode)
     873                                               psRegion bounds, pmSubtractionMode mode)
    867874{
    868875    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    877884
    878885    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    879                                                                   penalty, mode); // Kernels
     886                                                                  penalty, bounds, mode); // Kernels
    880887    kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
    881888    psStringPrepend(&kernels->description, "GUNK=");
    … …  
    893900// RINGS --- just what it says
    894901pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
    895                                                 float penalty, pmSubtractionMode mode)
     902                                                float penalty, psRegion bounds, pmSubtractionMode mode)
    896903{
    897904    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    924931
    925932    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
    926                                                               spatialOrder, penalty, mode); // The kernels
     933                                                              spatialOrder, penalty, bounds, mode); // Kernels
    927934    kernels->inner = inner;
    928935    psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
    … …  
    10461053pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
    10471054                                                   const psVector *fwhms, const psVector *orders, int inner,
    1048                                                    int binning, int ringsOrder, float penalty,
     1055                                                   int binning, int ringsOrder, float penalty, psRegion bounds,
    10491056                                                   pmSubtractionMode mode)
    10501057{
    10511058    switch (type) {
    10521059      case PM_SUBTRACTION_KERNEL_POIS:
    1053         return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode);
     1060        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, bounds, mode);
    10541061      case PM_SUBTRACTION_KERNEL_ISIS:
    1055         return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
     1062        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10561063      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
    1057         return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, mode);
     1064        return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10581065      case PM_SUBTRACTION_KERNEL_HERM:
    1059         return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode);
     1066        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10601067      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
    1061         return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode);
     1068        return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10621069      case PM_SUBTRACTION_KERNEL_SPAM:
    1063         return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
     1070        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, bounds, mode);
    10641071      case PM_SUBTRACTION_KERNEL_FRIES:
    1065         return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode);
     1072        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, bounds, mode);
    10661073      case PM_SUBTRACTION_KERNEL_GUNK:
    1067         return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode);
     1074        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, bounds, mode);
    10681075      case PM_SUBTRACTION_KERNEL_RINGS:
    1069         return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode);
     1076        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode);
    10701077      default:
    10711078        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
    … …  
    11031110
    11041111pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder,
    1105                                                           pmSubtractionMode mode)
     1112                                                          psRegion bounds, pmSubtractionMode mode)
    11061113{
    11071114    PS_ASSERT_STRING_NON_EMPTY(description, NULL);
    … …  
    11561163        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    11571164        penalty = parseStringFloat(ptr);
    1158 
    1159         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
    1160 
     1165        break;
    11611166      case PM_SUBTRACTION_KERNEL_RINGS:
    11621167        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
    … …  
    11651170        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    11661171        PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
    1167         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
     1172        break;
    11681173      default:
    11691174        psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
    11701175    }
    1171     return NULL;
     1176
     1177    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning,
     1178                                        ringsOrder, penalty, bounds, mode);
    11721179}
    11731180
    … …  
    12451252    out->bgOrder = in->bgOrder;
    12461253    out->mode = in->mode;
    1247     out->numCols = in->numCols;
    1248     out->numRows = in->numRows;
     1254    out->xMin = in->xMin;
     1255    out->xMax = in->xMax;
     1256    out->yMin = in->yMin;
     1257    out->yMax = in->yMax;
    12491258    out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
    12501259    out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h

    r26593 r26739  
    1212    PM_SUBTRACTION_KERNEL_ISIS_RADIAL,  ///< ISIS + higher-order radial Hermitians
    1313    PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
    14     PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
     14    PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
    1515    PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
    1616    PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
    … …  
    3232    pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
    3333    psString description;               ///< Description of the kernel parameters
     34    int xMin, xMax, yMin, yMax;         ///< Bounds of image (for normalisation)
    3435    long num;                           ///< Number of kernel components (not including the spatial ones)
    3536    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM)
    … …  
    4445    int bgOrder;                        ///< The order for the background fitting
    4546    pmSubtractionMode mode;             ///< Mode for subtraction
    46     int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
    4747    psVector *solution1, *solution2;    ///< Solution for the PSF matching
    4848    // Quality information
    4949    float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
    5050    int numStamps;                      ///< Number of good stamps
    51     float fSigResMean;                  ///< mean fractional stdev of residuals
    52     float fSigResStdev;                 ///< stdev of fractional stdev of residuals
    53     float fMaxResMean;                  ///< mean fractional positive swing in residuals
    54     float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
    55     float fMinResMean;                  ///< mean fractional negative swing in residuals
    56     float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
    57     psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
     51    float fSigResMean;                  ///< mean fractional stdev of residuals
     52    float fSigResStdev;                 ///< stdev of fractional stdev of residuals
     53    float fMaxResMean;                  ///< mean fractional positive swing in residuals
     54    float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
     55    float fMinResMean;                  ///< mean fractional negative swing in residuals
     56    float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
     57    psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
    5858} pmSubtractionKernels;
    5959
    6060// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures
    6161typedef struct {
    62     psVector *uCoords;                  // used by RINGS
    63     psVector *vCoords;                  // used by RINGS
    64     psVector *poly;                     // used by RINGS
    65 
    66     psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
    67     psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
    68     psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
     62    psVector *uCoords;                  // used by RINGS
     63    psVector *vCoords;                  // used by RINGS
     64    psVector *poly;                     // used by RINGS
     65
     66    psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
     67    psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
     68    psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
    6969} pmSubtractionKernelPreCalc;
    7070
    … …  
    162162                                                int spatialOrder, ///< Order of spatial variations
    163163                                                float penalty, ///< Penalty for wideness
     164                                                psRegion bounds,       ///< Bounds for validity
    164165                                                pmSubtractionMode mode ///< Mode for subtraction
    165166    );
    … …  
    168169pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(
    169170    pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated)
    170     int uOrder,                    ///< order in x-direction
    171     int vOrder,                    ///< order in x-direction
    172     int size,                      ///< Half-size of the kernel
    173     float sigma                    ///< sigma of gaussian kernel
     171    int uOrder,                    ///< order in x-direction
     172    int vOrder,                    ///< order in x-direction
     173    int size,                      ///< Half-size of the kernel
     174    float sigma                    ///< sigma of gaussian kernel
    174175    );
    175176
    … …  
    179180                                               int spatialOrder, ///< Order of spatial variations
    180181                                               float penalty, ///< Penalty for wideness
     182                                               psRegion bounds,       ///< Bounds for validity
    181183                                               pmSubtractionMode mode ///< Mode for subtraction
    182184    );
    … …  
    188190                                                    const psVector *orders, ///< Polynomial order of gaussians
    189191                                                    float penalty, ///< Penalty for wideness
     192                                                    psRegion bounds,       ///< Bounds for validity
    190193                                                    pmSubtractionMode mode ///< Mode for subtraction
    191194    );
    … …  
    197200                                               const psVector *orders, ///< Polynomial order of gaussians
    198201                                               float penalty, ///< Penalty for wideness
     202                                               psRegion bounds,       ///< Bounds for validity
    199203                                               pmSubtractionMode mode ///< Mode for subtraction
    200204                                               );
    … …  
    202206/// Generate ISIS + RADIAL_HERM kernels
    203207pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, ///< Half-size of the kernel
    204                                                       int spatialOrder, ///< Order of spatial variations
    205                                                       const psVector *fwhms, ///< Gaussian FWHMs
    206                                                       const psVector *orders, ///< Polynomial order of gaussians
    207                                                       float penalty, ///< Penalty for wideness
    208                                                       pmSubtractionMode mode ///< Mode for subtraction
     208                                                      int spatialOrder, ///< Order of spatial variations
     209                                                      const psVector *fwhms, ///< Gaussian FWHMs
     210                                                      const psVector *orders, ///< Polynomial order of gaussians
     211                                                      float penalty, ///< Penalty for wideness
     212                                                      psRegion bounds,       ///< Bounds for validity
     213                                                      pmSubtractionMode mode ///< Mode for subtraction
    209214                                               );
    210215
    … …  
    215220                                               const psVector *orders, ///< order of hermitian polynomials
    216221                                               float penalty, ///< Penalty for wideness
     222                                               psRegion bounds,       ///< Bounds for validity
    217223                                               pmSubtractionMode mode ///< Mode for subtraction
    218224                                               );
    … …  
    220226/// Generate DECONV_HERM kernels
    221227pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, ///< Half-size of the kernel
    222                                                      int spatialOrder, ///< Order of spatial variations
    223                                                      const psVector *fwhms, ///< Gaussian FWHMs
    224                                                      const psVector *orders, ///< order of hermitian polynomials
    225                                                      float penalty, ///< Penalty for wideness
    226                                                      pmSubtractionMode mode ///< Mode for subtraction
     228                                                      int spatialOrder, ///< Order of spatial variations
     229                                                      const psVector *fwhms, ///< Gaussian FWHMs
     230                                                      const psVector *orders, ///< order of hermitian polynomials
     231                                                      float penalty, ///< Penalty for wideness
     232                                                      psRegion bounds,       ///< Bounds for validity
     233                                                      pmSubtractionMode mode ///< Mode for subtraction
    227234    );
    228235
    … …  
    233240                                               int binning, ///< Kernel binning factor
    234241                                               float penalty, ///< Penalty for wideness
     242                                               psRegion bounds,       ///< Bounds for validity
    235243                                               pmSubtractionMode mode ///< Mode for subtraction
    236244    );
    … …  
    241249                                                int inner, ///< Inner radius to preserve unbinned
    242250                                                float penalty, ///< Penalty for wideness
     251                                                psRegion bounds,       ///< Bounds for validity
    243252                                                pmSubtractionMode mode ///< Mode for subtraction
    244253    );
    … …  
    251260                                               int inner, ///< Inner radius containing grid of delta functions
    252261                                               float penalty, ///< Penalty for wideness
     262                                               psRegion bounds,       ///< Bounds for validity
    253263                                               pmSubtractionMode mode ///< Mode for subtraction
    254264    );
    … …  
    260270                                                int ringsOrder, ///< Polynomial order
    261271                                                float penalty, ///< Penalty for wideness
     272                                                psRegion bounds,       ///< Bounds for validity
    262273                                                pmSubtractionMode mode ///< Mode for subtraction
    263274    );
    … …  
    274285                                                   int ringsOrder, ///< Polynomial order for RINGS
    275286                                                   float penalty, ///< Penalty for wideness
     287                                                   psRegion bounds,       ///< Bounds for validity
    276288                                                   pmSubtractionMode mode ///< Mode for subtraction
    277289    );
    … …  
    281293    const char *description,            ///< Description of kernel
    282294    int bgOrder,                        ///< Polynomial order for background fitting
     295    psRegion bounds,                    ///< Bounds for validity
    283296    pmSubtractionMode mode              ///< Mode for subtraction
    284297    );

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.c

    r24257 r26739  
    3838//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    3939
    40 psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psImageMaskType maskVal,
    41                            int size, int footprint, float badFrac, pmSubtractionMode mode)
    42 {
    43     PS_ASSERT_IMAGE_NON_NULL(mask1, NULL);
    44     PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_IMAGE_MASK, NULL);
    45     if (mask2) {
    46         PS_ASSERT_IMAGE_NON_NULL(mask2, NULL);
    47         PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_IMAGE_MASK, NULL);
    48         PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL);
    49     }
     40psImage *pmSubtractionMask(psRegion *bounds, const pmReadout *ro1, const pmReadout *ro2,
     41                           psImageMaskType maskVal, int size, int footprint, float badFrac,
     42                           pmSubtractionMode mode)
     43{
     44    PM_ASSERT_READOUT_NON_NULL(ro1, NULL);
     45    PM_ASSERT_READOUT_IMAGE(ro1, NULL);
     46    PM_ASSERT_READOUT_MASK(ro1, NULL);
     47    PM_ASSERT_READOUT_NON_NULL(ro2, NULL);
     48    PM_ASSERT_READOUT_IMAGE(ro2, NULL);
     49    PM_ASSERT_READOUT_MASK(ro2, NULL);
     50    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL);
    5051    PS_ASSERT_INT_NONNEGATIVE(size, NULL);
    5152    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
    … …  
    5556    }
    5657
    57     int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images
     58    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of the images
    5859
    5960    // Dereference inputs for convenience
    60     psImageMaskType **data1 = mask1->data.PS_TYPE_IMAGE_MASK_DATA;
    61     psImageMaskType **data2 = NULL;
    62     if (mask2) {
    63         data2 = mask2->data.PS_TYPE_IMAGE_MASK_DATA;
    64     }
     61    psF32 **imageData1 = ro1->image->data.F32, **imageData2 = ro2->image->data.F32;
     62    psImageMaskType **maskData1 = ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA,
     63        **maskData2 = ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA;
    6564
    6665    // First, a pass through to determine the fraction of bad pixels
    67     if (isfinite(badFrac) && badFrac != 1.0) {
     66    if (bounds || (isfinite(badFrac) && badFrac != 1.0)) {
     67        int xMin = numCols, xMax = 0, yMin = numRows, yMax = 0; // Bounds of good pixels
    6868        int numBad = 0;                 // Number of bad pixels
    6969        for (int y = 0; y < numRows; y++) {
    7070            for (int x = 0; x < numCols; x++) {
    71                 if (data1[y][x] & maskVal) {
     71                if ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x])) {
    7272                    numBad++;
    7373                    continue;
    7474                }
    75                 if (data2 && data2[y][x] & maskVal) {
     75                if ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x])) {
    7676                    numBad++;
     77                    continue;
    7778                }
    78             }
    79         }
    80         if (numBad > badFrac * numCols * numRows) {
     79                xMin = PS_MIN(xMin, x);
     80                xMax = PS_MAX(xMax, x);
     81                yMin = PS_MIN(yMin, y);
     82                yMax = PS_MAX(yMax, y);
     83            }
     84        }
     85        if (bounds) {
     86            bounds->x0 = xMin;
     87            bounds->x1 = xMax;
     88            bounds->y0 = yMin;
     89            bounds->y1 = yMax;
     90        }
     91        if (isfinite(badFrac) && badFrac != 1.0 && numBad > badFrac * numCols * numRows) {
    8192            psError(PM_ERR_SMALL_AREA, true,
    8293                    "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n",
    … …  
    117128    for (int y = 0; y < numRows; y++) {
    118129        for (int x = 0; x < numCols; x++) {
    119             if (data1[y][x] & maskVal) {
     130            if (maskData1[y][x] & maskVal) {
    120131                maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1;
    121132            }
    122             if (data2 && data2[y][x] & maskVal) {
     133            if (maskData2[y][x] & maskVal) {
    123134                maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2;
    124135            }

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.h

    r23851 r26739  
    55
    66/// Generate a mask for use in the subtraction process
    7 psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved)
    8                            const psImage *inMask, ///< Mask for the input image, or NULL
    9                            psImageMaskType maskVal, ///< Value to mask out
    10                            int size, ///< Half-size of the kernel (pmSubtractionKernels.size)
    11                            int footprint, ///< Half-size of the kernel footprint
    12                            float badFrac, ///< Maximum fraction of bad input pixels to accept
    13                            pmSubtractionMode mode  ///< Subtraction mode
     7psImage *pmSubtractionMask(
     8    psRegion *bounds,                   ///< Bounds of valid pixels (or NULL), returned
     9    const pmReadout *ro1,               ///< Readout 1
     10    const pmReadout *ro2,               ///< Readout 2
     11    psImageMaskType maskVal,            ///< Value to mask out
     12    int size,                           ///< Half-size of the kernel (pmSubtractionKernels.size)
     13    int footprint,                      ///< Half-size of the kernel footprint
     14    float badFrac,                      ///< Maximum fraction of bad input pixels to accept
     15    pmSubtractionMode mode              ///< Subtraction mode
    1416    );
    1517

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c

    r26733 r26739  
    7070                                 const psImage *subMask, // Mask for subtraction, or NULL
    7171                                 psImage *variance,  // Variance map
    72                                  const psRegion *region, // Region of interest, or NULL
     72                                 const psRegion *region, // Region of interest
    7373                                 float thresh1,  // Threshold for stamp finding on readout 1
    7474                                 float thresh2,  // Threshold for stamp finding on readout 2
    … …  
    8282    )
    8383{
     84    PS_ASSERT_PTR_NON_NULL(stamps, false);
     85    PM_ASSERT_READOUT_NON_NULL(ro1, false);
     86    PM_ASSERT_READOUT_NON_NULL(ro2, false);
     87    PS_ASSERT_IMAGE_NON_NULL(subMask, false);
     88    PS_ASSERT_IMAGE_NON_NULL(variance, false);
     89    PS_ASSERT_PTR_NON_NULL(region, false);
     90
    8491    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    8592
    … …  
    96103
    97104    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    98     if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) {
     105    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
    99106        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
    100107        return false;
    … …  
    319326    pmSubtractionMaskInvalid(ro2, maskVal);
    320327
    321     psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
     328    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     329
     330    psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0,
    322331                                         badFrac, mode); // Subtraction mask
    323332    if (!subMask) {
    … …  
    441450    // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
    442451    psImage *subMask = NULL;            // Mask for subtraction
    443     psRegion *region = NULL;            // Iso-kernel region
     452    psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region
    444453    psString regionString = NULL;       // String for region
    445454    pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
    … …  
    457466    pmSubtractionMaskInvalid(ro2, maskVal);
    458467
    459     subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
    460                                 badFrac, subMode);
     468    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     469
     470    subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode);
    461471    if (!subMask) {
    462472        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    … …  
    468478    // Get region of interest
    469479    int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
    470     float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
     480    float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions
    471481    if (isfinite(regionSize) && regionSize != 0.0) {
    472         xRegions = numCols / regionSize + 1;
    473         yRegions = numRows / regionSize + 1;
    474         xRegionSize = (float)numCols / (float)xRegions;
    475         yRegionSize = (float)numRows / (float)yRegions;
    476         region = psRegionAlloc(NAN, NAN, NAN, NAN);
     482        xRegions = (bounds.x1 - bounds.x0) / regionSize + 1;
     483        yRegions = (bounds.y1 - bounds.y0) / regionSize + 1;
     484        xRegionSize = (float)(bounds.x1 - bounds.x0) / (float)xRegions;
     485        yRegionSize = (float)(bounds.y1 - bounds.y0) / (float)yRegions;
     486    } else {
     487        xRegionSize = bounds.x1 - bounds.x0;
     488        yRegionSize = bounds.y1 - bounds.y0;
    477489    }
    478490
    … …  
    480492    float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
    481493    {
    482         psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
     494        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    483495        if (ro1) {
    484496            psStatsInit(bg);
    … …  
    504516    }
    505517
     518    // Just in case the iso-kernel region doesn't cover the entire image...
     519    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
     520        subMode == PM_SUBTRACTION_MODE_DUAL) {
     521        if (!conv1->image) {
     522            conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     523        }
     524        psImageInit(conv1->image, NAN);
     525        if (ro1->variance) {
     526            if (!conv1->variance) {
     527                conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     528            }
     529            psImageInit(conv1->variance, NAN);
     530        }
     531        if (subMask) {
     532            if (!conv1->mask) {
     533                conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     534            }
     535            psImageInit(conv1->mask, maskBad);
     536        }
     537    }
     538    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
     539        subMode == PM_SUBTRACTION_MODE_DUAL) {
     540        if (!conv2->image) {
     541            conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     542        }
     543        psImageInit(conv2->image, NAN);
     544        if (ro2->variance) {
     545            if (!conv2->variance) {
     546                conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     547            }
     548            psImageInit(conv2->variance, NAN);
     549        }
     550        if (subMask) {
     551            if (!conv2->mask) {
     552                conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     553            }
     554            psImageInit(conv2->mask, maskBad);
     555        }
     556    }
     557
     558
    506559    // Iterate over iso-kernel regions
    507560    for (int j = 0; j < yRegions; j++) {
    … …  
    509562            psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
    510563                    j * xRegions + i + 1, xRegions * yRegions);
    511             if (region) {
    512                 *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
    513                                       (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
    514                 psFree(regionString);
    515                 regionString = psRegionToString(*region);
    516                 psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
    517                         regionString, numCols, numRows);
    518             }
     564            *region = psRegionSet(bounds.x0 + (int)(i * xRegionSize),
     565                                  bounds.x0 + (int)((i + 1) * xRegionSize),
     566                                  bounds.y0 + (int)(j * yRegionSize),
     567                                  bounds.y0 + (int)((j + 1) * yRegionSize));
     568            psFree(regionString);
     569            regionString = psRegionToString(*region);
     570            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
     571                    regionString, numCols, numRows);
    519572
    520573            if (stampsName && strlen(stampsName) > 0) {
    … …  
    530583            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    531584            // doesn't matter.
    532             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
     585            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
    533586                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
    534587                goto MATCH_ERROR;
    … …  
    544597            // Define kernel basis functions
    545598            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    546                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    547                                                           stamps, footprint, optThreshold, penalty, subMode);
     599                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
     600                                                          optFWHMs, optOrder, stamps, footprint,
     601                                                          optThreshold, penalty, bounds, subMode);
    548602                if (!kernels) {
    549603                    psErrorClear();
    … …  
    554608                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    555609                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    556                                                        inner, binning, ringsOrder, penalty, subMode);
     610                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
    557611                // pmSubtractionVisualShowKernels(kernels);
    558612            }
    … …  
    603657                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    604658
    605 #if 0
    606659                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    607660                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
    … …  
    609662                    goto MATCH_ERROR;
    610663                }
    611 #endif
    612664
    613665                // generate the window function from the set of stamps

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.c

    r26491 r26739  
    204204                                                      int spatialOrder, const psVector *fwhms, int maxOrder,
    205205                                                      const pmSubtractionStampList *stamps, int footprint,
    206                                                       float tolerance, float penalty, pmSubtractionMode mode)
     206                                                      float tolerance, float penalty, psRegion bounds,
     207                                                      pmSubtractionMode mode)
    207208{
    208209    if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
    … …  
    232233    psVectorInit(orders, maxOrder);
    233234    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    234                                                                   penalty, mode); // Kernels
     235                                                                  penalty, bounds, mode); // Kernels
    235236    psFree(orders);
    236237    psFree(kernels->description);
    … …  
    483484    if (type == PM_SUBTRACTION_KERNEL_ISIS) {
    484485
    485         // XXX in r26035, this code was just wrong.  we had:
    486 
    487         // psKernel *subtract = kernels->preCalc->data[0]
    488 
    489         // but, kernels->preCalc was an array of psArray, not an array of kernels.  It is now
    490         // an array of pmSubtractionKernelPreCalc.
     486        // XXX in r26035, this code was just wrong.  we had:
     487
     488        // psKernel *subtract = kernels->preCalc->data[0]
     489
     490        // but, kernels->preCalc was an array of psArray, not an array of kernels.  It is now
     491        // an array of pmSubtractionKernelPreCalc.
    491492
    492493        pmSubtractionKernelPreCalc *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.h

    r18287 r26739  
    1717                                                      float tolerance, ///< Maximum difference in chi^2
    1818                                                      float penalty, ///< Penalty for wideness
    19                                                       pmSubtractionMode mode // Mode for subtraction
     19                                                      psRegion bounds,       ///< Bounds of validity
     20                                                      pmSubtractionMode mode ///< Mode for subtraction
    2021    );
    2122

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c

    r26735 r26739  
    709709        }
    710710    }
    711     fprintf(stderr, "Normalization window radius: %d\n", stamps->normWindow);
    712711
    713712    psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
    … …  
    741740
    742741bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
    743                                 psImage *variance, int kernelSize)
     742                                psImage *variance, int kernelSize, psRegion bounds)
    744743{
    745744    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    758757    }
    759758
    760     int numCols = image1->numCols, numRows = image1->numRows; // Size of images
    761759    int size = kernelSize + stamps->footprint; // Size of postage stamps
    762760
    … …  
    767765        }
    768766
    769         if (isnan(stamp->xNorm)) {
    770             stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
    771         }
    772         if (isnan(stamp->yNorm)) {
    773             stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
    774         }
     767        p_pmSubtractionPolynomialNormCoords(&stamp->xNorm, &stamp->yNorm, stamp->x, stamp->y,
     768                                            bounds.x0, bounds.x1, bounds.y0, bounds.y1);
    775769
    776770        int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
    777         if (x < size || x > numCols - size || y < size || y > numRows - size) {
    778             psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
     771        if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) {
     772            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
    779773            return false;
    780774        }

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h

    r26703 r26739  
    164164                                psImage *image2, ///< Input image (or NULL)
    165165                                psImage *variance, ///< Variance map
    166                                 int kernelSize ///< Kernel half-size
     166                                int kernelSize, ///< Kernel half-size
     167                                psRegion bounds ///< Bounds of validity
    167168    );
    168169
Note: See TracChangeset for help on using the changeset viewer.