IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 31, 2010, 5:00:42 PM (16 years ago)
Author:
eugene
Message:

updates relative to 20091201, fixes for all psphot variants

Location:
branches/eam_branches/psModules.stack.20100120
Files:
30 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/psModules.stack.20100120

  • branches/eam_branches/psModules.stack.20100120/src/camera/pmFPAMaskWeight.c

    r26076 r26747  
    370370
    371371    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts
    372 
    373372    int numCols = image->numCols, numRows = image->numRows; // Size of image
    374     int numPix = numCols * numRows;                         // Number of pixels
    375     int num = PS_MAX(sample, numPix);                       // Number we care about
     373
     374    int xMin, xMax, yMin, yMax;         // Bounds of image
     375    if (mask) {
     376        xMin = numCols;
     377        xMax = 0;
     378        yMin = numRows;
     379        yMax = 0;
     380        for (int y = 0; y < numRows; y++) {
     381            for (int x = 0; x < numCols; x++) {
     382                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
     383                    continue;
     384                }
     385                xMin = PS_MIN(xMin, x);
     386                xMax = PS_MAX(xMax, x);
     387                yMin = PS_MIN(yMin, y);
     388                yMax = PS_MAX(yMax, y);
     389            }
     390        }
     391    } else {
     392        xMin = 0;
     393        xMax = numCols;
     394        yMin = 0;
     395        yMax = numRows;
     396    }
     397
     398    int xNum = xMax - xMin, yNum = yMax - yMin; // Number of pixels
     399
     400    int numPix = xNum * yNum;                                  // Number of pixels
     401    int num = PS_MAX(sample, numPix);                          // Number we care about
    376402    psVector *signoise = psVectorAllocEmpty(num, PS_TYPE_F32);   // Signal-to-noise values
    377403
    … …  
    379405        // We have an image smaller than Nsubset, so just loop over the image pixels
    380406        int index = 0;                  // Index for vector
    381         for (int y = 0; y < numRows; y++) {
    382             for (int x = 0; x < numCols; x++) {
     407        for (int y = yMin; y < yMax; y++) {
     408            for (int x = xMin; x < xMax; x++) {
    383409                if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
    384410                    !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) {
    … …  
    397423            // Pixel coordinates
    398424            int pixel = numPix * psRandomUniform(rng);
    399             int x = pixel % numCols;
    400             int y = pixel / numCols;
     425            int x = xMin + pixel % xNum;
     426            int y = yMin + pixel / xNum;
    401427
    402428            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
    … …  
    433459    }
    434460
    435     psBinaryOp(variance, variance, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
     461    psImage *subImage = psImageSubset(variance, psRegionSet(xMin, xMax, yMin, yMax)); // Smaller image
     462    psBinaryOp(subImage, subImage, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
     463    psFree(subImage);
    436464
    437465    pmHDU *hdu = pmHDUFromReadout(readout); // HDU for readout

  • branches/eam_branches/psModules.stack.20100120/src/config/Makefile.am

    r23794 r26747  
    3232        pmConfigDump.c \
    3333        pmConfigRun.c \
     34        pmConfigRecipeValue.c \
    3435        pmVersion.c \
    3536        pmErrorCodes.c
    … …  
    4344        pmConfigDump.h \
    4445        pmConfigRun.h \
     46        pmConfigRecipeValue.h \
    4547        pmVersion.h \
    4648        pmErrorCodes.h

  • branches/eam_branches/psModules.stack.20100120/src/detrend/pmPattern.c

    r24905 r26747  
    6363    float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
    6464    float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data
     65    float background = stats->robustMedian;
    6566    psFree(stats);
    6667    psFree(rng);
    … …  
    113114
    114115        for (int x = 0; x < numCols; x++) {
    115             image->data.F32[y][x] -= solution->data.F32[x];
     116            image->data.F32[y][x] += (background - solution->data.F32[x]);
    116117        }
    117118        psFree(solution);

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmStackReject.c

    r25468 r26747  
    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);
    … …  
    150150    pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
    151151    inRO->image = image;
     152    convRO->image = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
    152153    for (int i = 0; i < numRegions; i++) {
    153154        psRegion *region = subRegions->data[i]; // Region of interest
    … …  
    165166
    166167        // 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;
     168        float xNorm, yNorm;             // Normalised coordinates
     169        p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, 0.5 * (region->x1 - region->x0),
     170                                            0.5 * (region->y1 - region->y0),
     171                                            kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax);
    171172        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
    172173        if (!kernel) {

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.c

    r26667 r26747  
    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
    … …  
    11581161        PM_ASSERT_READOUT_NON_NULL(ro1, false);
    11591162        PM_ASSERT_READOUT_IMAGE(ro1, false);
     1163        PM_ASSERT_READOUT_IMAGE(out1, false);
    11601164        numCols = ro1->image->numCols;
    11611165        numRows = ro1->image->numRows;
    … …  
    11671171        PM_ASSERT_READOUT_NON_NULL(ro2, false);
    11681172        PM_ASSERT_READOUT_IMAGE(ro2, false);
     1173        PM_ASSERT_READOUT_IMAGE(out2, false);
    11691174        if (numCols == 0 && numRows == 0) {
    11701175            numCols = ro2->image->numCols;
    … …  
    12001205    bool threaded = pmSubtractionThreaded(); // Running threaded?
    12011206
    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     }
    12251207    psImage *convMask = NULL;           // Convolved mask image (common to inputs 1 and 2)
    12261208    if (subMask) {
    12271209        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);
    12321210            convMask = out1->mask;
    12331211        }
    12341212        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);
    12391213            if (!convMask) {
    12401214                convMask = out2->mask;
    … …  
    12561230
    12571231    // Get region for convolution: [xMin:xMax,yMin:yMax]
    1258     int xMin = size, xMax = numCols - size;
    1259     int yMin = size, yMax = numRows - size;
     1232    int xMin = kernels->xMin + size, xMax = kernels->xMax - size;
     1233    int yMin = kernels->yMin + size, yMax = kernels->yMax - size;
    12601234    if (region) {
    12611235        xMin = PS_MAX(region->x0, xMin);

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtraction.h

    r25279 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionAnalysis.c

    r26593 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.c

    r26667 r26747  
    2828static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
    2929                                  psVector *vector, // Least-squares vector, updated
     30                                  double *norm,     // Normalisation, updated
    3031                                  const psKernel *input, // Input image (target)
    3132                                  const psKernel *reference, // Reference image (convolution source)
    … …  
    3637                                  const psImage *polyValues, // Spatial polynomial values
    3738                                  int footprint, // (Half-)Size of stamp
     39                                  int normWindow, // Window (half-)size for normalisation measurement
    3840                                  const pmSubtractionEquationCalculationMode mode
    3941                                  )
    … …  
    156158            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    157159                vector->data.F64[iIndex] = sumIC * poly[iTerm];
    158                 // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
    159                 // instead, calculate A - \sum(k x B), with full hermitians
    160                 if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
     160                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
    161161                    // subtract norm * sumRC * poly[iTerm]
    162162                    psAssert (kernels->solution1, "programming error: define solution first!");
    … …  
    174174    double sumR = 0.0;                  // Sum of the reference
    175175    double sumI = 0.0;                  // Sum of the input
     176    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
    176177    for (int y = - footprint; y <= footprint; y++) {
    177178        for (int x = - footprint; x <= footprint; x++) {
    … …  
    181182            double rr = PS_SQR(ref);
    182183            double one = 1.0;
     184
     185            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     186                normI1 += ref;
     187                normI2 += in;
     188            }
     189
    183190            if (weight) {
    184191                float wtVal = weight->kernel[y][x];
    … …  
    204211        }
    205212    }
     213
     214    *norm = normI2 / normI1;
     215
    206216    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    207217        matrix->data.F64[normIndex][normIndex] = sumRR;
    … …  
    217227        matrix->data.F64[bgIndex][normIndex] = sumR;
    218228    }
     229
     230    // check for any NAN values in the result, skip if found:
     231    for (int iy = 0; iy < matrix->numRows; iy++) {
     232        for (int ix = 0; ix < matrix->numCols; ix++) {
     233            if (!isfinite(matrix->data.F64[iy][ix])) {
     234                fprintf (stderr, "WARNING: NAN in matrix\n");
     235                return false;
     236            }
     237        }
     238    }
     239    for (int ix = 0; ix < vector->n; ix++) {
     240        if (!isfinite(vector->data.F64[ix])) {
     241            fprintf (stderr, "WARNING: NAN in vector\n");
     242            return false;
     243        }
     244    }
     245
    219246    return true;
    220247}
    … …  
    224251static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated
    225252                                      psVector *vector, // Least-squares vector, updated
     253                                      double *norm,     // Normalisation, updated
    226254                                      const psKernel *image1, // Image 1
    227255                                      const psKernel *image2, // Image 2
    … …  
    232260                                      const pmSubtractionKernels *kernels, // Kernels
    233261                                      const psImage *polyValues, // Spatial polynomial values
    234                                       int footprint // (Half-)Size of stamp
     262                                      int footprint, // (Half-)Size of stamp
     263                                      int normWindow, // Window (half-)size for normalisation measurement
     264                                      const pmSubtractionEquationCalculationMode mode
    235265                                      )
    236266{
    … …  
    272302    }
    273303
     304
     305    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
     306    // choose to calculate
     307    psImageInit(matrix, 0.0);
     308    psVectorInit(vector, 1.0);
     309    for (int i = 0; i < matrix->numCols; i++) {
     310        matrix->data.F64[i][i] = 1.0;
     311    }
    274312
    275313    for (int i = 0; i < numKernels; i++) {
    … …  
    306344            }
    307345
    308             // Spatial variation
    309             for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    310                 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    311                     double aa = sumAA * poly2[iTerm][jTerm];
    312                     double bb = sumBB * poly2[iTerm][jTerm];
    313                     double ab = sumAB * poly2[iTerm][jTerm];
    314 
    315                     matrix->data.F64[iIndex][jIndex] = aa;
    316                     matrix->data.F64[jIndex][iIndex] = aa;
    317 
    318                     matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
    319                     matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
    320 
    321                     matrix->data.F64[iIndex][jIndex + numParams] = ab;
    322                     matrix->data.F64[jIndex + numParams][iIndex] = ab;
     346            // Spatial variation of kernel coeffs
     347            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     348                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     349                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     350                        double aa = sumAA * poly2[iTerm][jTerm];
     351                        double bb = sumBB * poly2[iTerm][jTerm];
     352                        double ab = sumAB * poly2[iTerm][jTerm];
     353
     354                        matrix->data.F64[iIndex][jIndex] = aa;
     355                        matrix->data.F64[jIndex][iIndex] = aa;
     356
     357                        matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
     358                        matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
     359
     360                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
     361                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
     362                    }
    323363                }
    324364            }
    … …  
    340380            }
    341381
    342             // Spatial variation
    343             for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    344                 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    345                     double ab = sumAB * poly2[iTerm][jTerm];
    346                     matrix->data.F64[iIndex][jIndex + numParams] = ab;
    347                     matrix->data.F64[jIndex + numParams][iIndex] = ab;
     382            // Spatial variation of kernel coeffs
     383            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     384                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     385                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     386                        double ab = sumAB * poly2[iTerm][jTerm];
     387                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
     388                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
     389                    }
    348390                }
    349391            }
    … …  
    403445            double bi2 = sumBI2 * poly[iTerm];
    404446            double ai1 = sumAI1 * poly[iTerm];
    405             double a = sumA * poly[iTerm];
     447            double a   = sumA * poly[iTerm];
    406448            double bi1 = sumBI1 * poly[iTerm];
    407             double b = sumB * poly[iTerm];
    408 
    409             matrix->data.F64[iIndex][normIndex] = ai1;
    410             matrix->data.F64[normIndex][iIndex] = ai1;
    411             matrix->data.F64[iIndex][bgIndex] = a;
    412             matrix->data.F64[bgIndex][iIndex] = a;
    413             matrix->data.F64[iIndex + numParams][normIndex] = bi1;
    414             matrix->data.F64[normIndex][iIndex + numParams] = bi1;
    415             matrix->data.F64[iIndex + numParams][bgIndex] = b;
    416             matrix->data.F64[bgIndex][iIndex + numParams] = b;
    417             vector->data.F64[iIndex] = ai2;
    418             vector->data.F64[iIndex + numParams] = bi2;
     449            double b   = sumB * poly[iTerm];
     450
     451            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     452                matrix->data.F64[iIndex][normIndex] = ai1;
     453                matrix->data.F64[normIndex][iIndex] = ai1;
     454                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
     455                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
     456            }
     457            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     458                matrix->data.F64[iIndex][bgIndex] = a;
     459                matrix->data.F64[bgIndex][iIndex] = a;
     460                matrix->data.F64[iIndex + numParams][bgIndex] = b;
     461                matrix->data.F64[bgIndex][iIndex + numParams] = b;
     462            }
     463            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     464                vector->data.F64[iIndex] = ai2;
     465                vector->data.F64[iIndex + numParams] = bi2;
     466                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
     467                    // subtract norm * sumRC * poly[iTerm]
     468                    psAssert (kernels->solution1, "programming error: define solution first!");
     469                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     470                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
     471                    vector->data.F64[iIndex] -= norm * ai1;
     472                    vector->data.F64[iIndex + numParams] -= norm * bi1;
     473                }
     474            }
    419475        }
    420476    }
    … …  
    425481    double sumI2 = 0.0;                 // Sum of I_2 (for vector, background)
    426482    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
     483    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
    427484    for (int y = - footprint; y <= footprint; y++) {
    428485        for (int x = - footprint; x <= footprint; x++) {
    … …  
    433490            double one = 1.0;
    434491            double i1i2 = i1 * i2;
     492
     493            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     494                normI1 += i1;
     495                normI2 += i2;
     496            }
    435497
    436498            if (weight) {
    … …  
    457519        }
    458520    }
    459     matrix->data.F64[bgIndex][normIndex] = sumI1;
    460     matrix->data.F64[normIndex][bgIndex] = sumI1;
    461     matrix->data.F64[normIndex][normIndex] = sumI1I1;
    462     matrix->data.F64[bgIndex][bgIndex] = sum1;
    463     vector->data.F64[bgIndex] = sumI2;
    464     vector->data.F64[normIndex] = sumI1I2;
     521
     522    *norm = normI2 / normI1;
     523
     524    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     525        matrix->data.F64[normIndex][normIndex] = sumI1I1;
     526        vector->data.F64[normIndex] = sumI1I2;
     527    }
     528    if (mode & PM_SUBTRACTION_EQUATION_BG) {
     529        matrix->data.F64[bgIndex][bgIndex] = sum1;
     530        vector->data.F64[bgIndex] = sumI2;
     531    }
     532    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
     533        matrix->data.F64[bgIndex][normIndex] = sumI1;
     534        matrix->data.F64[normIndex][bgIndex] = sumI1;
     535    }
     536
     537    // check for any NAN values in the result, skip if found:
     538    for (int iy = 0; iy < matrix->numRows; iy++) {
     539        for (int ix = 0; ix < matrix->numCols; ix++) {
     540            if (!isfinite(matrix->data.F64[iy][ix])) {
     541                fprintf (stderr, "WARNING: NAN in matrix\n");
     542                return false;
     543            }
     544        }
     545    }
     546    for (int ix = 0; ix < vector->n; ix++) {
     547        if (!isfinite(vector->data.F64[ix])) {
     548            fprintf (stderr, "WARNING: NAN in vector\n");
     549            return false;
     550        }
     551    }
     552
    465553
    466554    return true;
    … …  
    469557#if 1
    470558// Add in penalty term to least-squares vector
    471 static bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
     559bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
    472560                             psVector *vector,                    // Vector to which to add in penalty term
    473561                             const pmSubtractionKernels *kernels, // Kernel parameters
    … …  
    482570    int spatialOrder = kernels->spatialOrder; // Order of spatial variations
    483571    int numKernels = kernels->num; // Number of kernel components
     572    int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations
     573    int numParams = numKernels * numSpatial;                 // Number of kernel parameters
     574
     575    // order is :
     576    // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0]
     577    // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0]
     578    // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1]
     579    // [norm]
     580    // [bg]
     581    // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0]
     582    // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0]
     583    // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1]
     584
    484585    for (int i = 0; i < numKernels; i++) {
    485586        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
    … …  
    488589                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
    489590                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
     591                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     592                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
     593                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
     594                    // penalties scale with second moments
     595                    //
     596                }
    490597            }
    491598        }
    … …  
    668775    switch (kernels->mode) {
    669776      case PM_SUBTRACTION_MODE_1:
    670         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1,
     777        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
    671778                                       weight, window, stamp->convolutions1, kernels,
    672                                        polyValues, footprint, mode);
     779                                       polyValues, footprint, stamps->normWindow, mode);
    673780        break;
    674781      case PM_SUBTRACTION_MODE_2:
    675         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2,
     782        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
    676783                                       weight, window, stamp->convolutions2, kernels,
    677                                        polyValues, footprint, mode);
     784                                       polyValues, footprint, stamps->normWindow, mode);
    678785        break;
    679786      case PM_SUBTRACTION_MODE_DUAL:
    680         status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2,
     787        status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm,
     788                                           stamp->image1, stamp->image2,
    681789                                           weight, window, stamp->convolutions1, stamp->convolutions2,
    682                                            kernels, polyValues, footprint);
     790                                           kernels, polyValues, footprint, stamps->normWindow, mode);
    683791        break;
    684792      default:
    … …  
    721829}
    722830
    723 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode)
     831bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     832                                    const pmSubtractionEquationCalculationMode mode)
    724833{
    725834    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    843952        psVectorInit(sumVector, 0.0);
    844953        psImageInit(sumMatrix, 0.0);
    845         int numStamps = 0;              // Number of good stamps
    846         for (int i = 0; i < stamps->num; i++) {
    847             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    848 
    849             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
    850 
    851 #ifdef TESTING
    852               // XXX double-check for NAN in data:
    853                 for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
    854                     for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
    855                         if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
    856                             fprintf (stderr, "WARNING: NAN in matrix\n");
    857                         }
    858                     }
    859                 }
    860                 for (int ix = 0; ix < stamp->vector->n; ix++) {
    861                     if (!isfinite(stamp->vector->data.F64[ix])) {
    862                         fprintf (stderr, "WARNING: NAN in vector\n");
    863                     }
    864                 }
    865 #endif
    866 
    867                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
    868                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
    869                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    870                 numStamps++;
    871             } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
    872                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
    873             }
    874         }
    875 
    876 #ifdef TESTING
    877         for (int ix = 0; ix < sumVector->n; ix++) {
    878             if (!isfinite(sumVector->data.F64[ix])) {
    879                 fprintf (stderr, "WARNING: NAN in vector\n");
    880             }
    881         }
    882 #endif
    883 
    884 #if 0
    885         int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    886         calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
    887 #endif
    888 
    889 #ifdef TESTING
    890         for (int ix = 0; ix < sumVector->n; ix++) {
    891             if (!isfinite(sumVector->data.F64[ix])) {
    892                 fprintf (stderr, "WARNING: NAN in vector\n");
    893             }
    894         }
    895         {
    896             psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
    897             psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
    898             psFree(inverse);
    899         }
    900         {
    901             psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
    902             psImage *Xt = psMatrixTranspose(NULL, X);
    903             psImage *XtX = psMatrixMultiply(NULL, Xt, X);
    904             psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
    905             psFree(X);
    906             psFree(Xt);
    907             psFree(XtX);
    908         }
    909 #endif
    910 
    911 # define SVD_ANALYSIS 0
    912 # define COEFF_SIG 0.0
    913 # define SVD_TOL 0.0
    914 
    915         psVector *solution = NULL;
    916         psVector *solutionErr = NULL;
    917 
    918 #ifdef TESTING
    919         psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
    920         psVectorWriteFile ("B.dat", sumVector);
    921 #endif
    922 
    923         // Use SVD to determine the kernel coeffs (and validate)
    924         if (SVD_ANALYSIS) {
    925 
    926             // We have sumVector and sumMatrix.  we are trying to solve the following equation:
    927             // sumMatrix * x = sumVector.
    928 
    929             // we can use any standard matrix inversion to solve this.  However, the basis
    930             // functions in general have substantial correlation, so that the solution may be
    931             // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
    932             // system of equations may be statistically ill-conditioned.  Noise in the image
    933             // will drive insignificant, but correlated, terms in the solution.  To avoid these
    934             // problems, we can use SVD to identify numerically unconstrained values and to
    935             // avoid statistically badly determined value.
    936 
    937             // A = sumMatrix, B = sumVector
    938             // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
    939             // x = V (1/w) (U^T B)
    940             // \sigma_x = sqrt(diag(A^{-1}))
    941             // solve for x and A^{-1} to get x & dx
    942             // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
    943             // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
    944 
    945             // If I use the SVD trick to re-condition the matrix, I need to break out the
    946             // kernel and normalization terms from the background term.
    947             // XXX is this true?  or was this due to an error in the analysis?
    948 
    949             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    950 
    951             // now pull out the kernel elements into their own square matrix
    952             psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
    953             psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
    954 
    955             for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
    956                 if (ix == bgIndex) continue;
    957                 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
    958                     if (iy == bgIndex) continue;
    959                     kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
    960                     ky++;
    961                 }
    962                 kernelVector->data.F64[kx] = sumVector->data.F64[ix];
    963                 kx++;
    964             }
    965 
    966             psImage *U = NULL;
    967             psImage *V = NULL;
    968             psVector *w = NULL;
    969             if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
    970                 psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
    971                 return NULL;
    972             }
    973 
    974             // calculate A_inverse:
    975             // Ainv = V * w * U^T
    976             psImage *wUt  = p_pmSubSolve_wUt (w, U);
    977             psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
    978             psImage *Xvar = NULL;
    979             psFree (wUt);
    980 
    981 # ifdef TESTING
    982             // kernel terms:
    983             for (int i = 0; i < w->n; i++) {
    984                 fprintf (stderr, "w: %f\n", w->data.F64[i]);
    985             }
    986 # endif
    987             // loop over w adding in more and more of the values until chisquare is no longer
    988             // dropping significantly.
    989             // XXX this does not seem to work very well: we seem to need all terms even for
    990             // simple cases...
    991 
    992             psVector *Xsvd = NULL;
    993             {
    994                 psVector *Ax = NULL;
    995                 psVector *UtB = NULL;
    996                 psVector *wUtB = NULL;
    997 
    998                 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
    999                 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
    1000                 psVectorInit (wMask, 1); // start by masking everything
    1001 
    1002                 double chiSquareLast = NAN;
    1003                 int maxWeight = 0;
    1004 
    1005                 double Axx, Bx, y2;
    1006 
    1007                 // XXX this is an attempt to exclude insignificant modes.
    1008                 // it was not successful with the ISIS kernel set: removing even
    1009                 // the least significant mode leaves additional ringing / noise
    1010                 // because the terms are so coupled.
    1011                 for (int k = 0; false && (k < w->n); k++) {
    1012 
    1013                     // unmask the k-th weight
    1014                     wMask->data.U8[k] = 0;
    1015                     p_pmSubSolve_SetWeights(wApply, w, wMask);
    1016 
    1017                     // solve for x:
    1018                     // x = V * w * (U^T * B)
    1019                     p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1020                     p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1021                     p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1022 
    1023                     // chi-square for this system of equations:
    1024                     // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1025                     // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1026                     p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1027                     p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1028                     p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1029                     p_pmSubSolve_y2 (&y2, kernels, stamps);
    1030 
    1031                     // apparently, this works (compare with the brute force value below
    1032                     double chiSquare = Axx - 2.0*Bx + y2;
    1033                     double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
    1034                     chiSquareLast = chiSquare;
    1035 
    1036                     // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
    1037                     if (k && !maxWeight && (deltaChi < 1.0)) {
    1038                         maxWeight = k;
    1039                     }
    1040                 }
    1041 
    1042                 // keep all terms or we get extra ringing
    1043                 maxWeight = w->n;
    1044                 psVectorInit (wMask, 1);
    1045                 for (int k = 0; k < maxWeight; k++) {
    1046                     wMask->data.U8[k] = 0;
    1047                 }
    1048                 p_pmSubSolve_SetWeights(wApply, w, wMask);
    1049 
    1050                 // solve for x:
    1051                 // x = V * w * (U^T * B)
    1052                 p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1053                 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1054                 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1055 
    1056                 // chi-square for this system of equations:
    1057                 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1058                 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1059                 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1060                 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1061                 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1062                 p_pmSubSolve_y2 (&y2, kernels, stamps);
    1063 
    1064                 // apparently, this works (compare with the brute force value below
    1065                 double chiSquare = Axx - 2.0*Bx + y2;
    1066                 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
    1067 
    1068                 // re-calculate A^{-1} to get new variances:
    1069                 // Ainv = V * w * U^T
    1070                 // XXX since we keep all terms, this is identical to Ainv
    1071                 psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
    1072                 Xvar = p_pmSubSolve_VwUt (V, wUt);
    1073                 psFree (wUt);
    1074 
    1075                 psFree (Ax);
    1076                 psFree (UtB);
    1077                 psFree (wUtB);
    1078                 psFree (wApply);
    1079                 psFree (wMask);
    1080             }
    1081 
    1082             // copy the kernel solutions to the full solution vector:
    1083             solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1084             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1085 
    1086             for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
    1087                 if (ix == bgIndex) {
    1088                     solution->data.F64[ix] = 0;
    1089                     solutionErr->data.F64[ix] = 0.001;
    1090                     continue;
    1091                 }
    1092                 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
    1093                 solution->data.F64[ix] = Xsvd->data.F64[kx];
    1094                 kx++;
    1095             }
    1096 
    1097             psFree (kernelMatrix);
    1098             psFree (kernelVector);
    1099 
    1100             psFree (U);
    1101             psFree (V);
    1102             psFree (w);
    1103 
    1104             psFree (Ainv);
    1105             psFree (Xsvd);
    1106         } else {
    1107             psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    1108             psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    1109             if (!luMatrix) {
    1110                 psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
    1111                 psFree(solution);
    1112                 psFree(sumVector);
    1113                 psFree(sumMatrix);
    1114                 psFree(luMatrix);
    1115                 psFree(permutation);
    1116                 return NULL;
    1117             }
    1118 
    1119             solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
    1120             psFree(luMatrix);
    1121             psFree(permutation);
    1122             if (!solution) {
    1123                 psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
    1124                 psFree(solution);
    1125                 psFree(sumVector);
    1126                 psFree(sumMatrix);
    1127                 return NULL;
    1128             }
    1129 
    1130             // XXX LUD does not provide A^{-1}?  fake the error for now
    1131             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1132             for (int ix = 0; ix < sumVector->n; ix++) {
    1133                 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
    1134             }
    1135         }
    1136 
    1137         if (!kernels->solution1) {
    1138             kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
    1139             psVectorInit (kernels->solution1, 0.0);
    1140         }
    1141 
    1142         // only update the solutions that we chose to calculate:
    1143         if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    1144             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1145             kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
    1146         }
    1147         if (mode & PM_SUBTRACTION_EQUATION_BG) {
    1148             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    1149             kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
    1150         }
    1151         if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    1152             int numKernels = kernels->num;
    1153             int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
    1154             int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    1155             for (int i = 0; i < numKernels * numPoly; i++) {
    1156                 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
    1157                 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
    1158                     // XXX fprintf (stderr, "drop\n");
    1159                     kernels->solution1->data.F64[i] = 0.0;
    1160                 } else {
    1161                     // XXX fprintf (stderr, "keep\n");
    1162                     kernels->solution1->data.F64[i] = solution->data.F64[i];
    1163                 }
    1164             }
    1165         }
    1166         // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
    1167         // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
    1168 
    1169         psFree(solution);
    1170         psFree(sumVector);
    1171         psFree(sumMatrix);
    1172 
    1173 #ifdef TESTING
    1174         // XXX double-check for NAN in data:
    1175         for (int ix = 0; ix < kernels->solution1->n; ix++) {
    1176             if (!isfinite(kernels->solution1->data.F64[ix])) {
    1177                 fprintf (stderr, "WARNING: NAN in vector\n");
    1178             }
    1179         }
    1180 #endif
    1181 
    1182     } else {
    1183         // Dual convolution solution
    1184 
    1185         // Accumulation of stamp matrices/vectors
    1186         psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    1187         psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
    1188         psImageInit(sumMatrix, 0.0);
    1189         psVectorInit(sumVector, 0.0);
     954
     955        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
    1190956
    1191957        int numStamps = 0;              // Number of good stamps
    … …  
    1195961                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
    1196962                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
     963                psVectorAppend(norms, stamp->norm);
     964                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
     965                numStamps++;
     966            } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
     967                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
     968            }
     969        }
     970
     971#if 0
     972        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     973        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
     974#endif
     975
     976        psVector *solution = NULL;                       // Solution to equation!
     977        solution = psVectorAlloc(numParams, PS_TYPE_F64);
     978        psVectorInit(solution, 0);
     979
     980#if 0
     981        // Regular, straight-forward solution
     982        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     983#else
     984        {
     985            // Solve normalisation and background separately
     986            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     987            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     988
     989            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
     990            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
     991                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
     992                psFree(stats);
     993                psFree(sumMatrix);
     994                psFree(sumVector);
     995                psFree(norms);
     996                return false;
     997            }
     998
     999            double normValue = stats->robustMedian;
     1000            // double bgValue = 0.0;
     1001
     1002            psFree(stats);
     1003
     1004            fprintf(stderr, "Norm: %lf\n", normValue);
     1005
     1006            // Solve kernel components
     1007            for (int i = 0; i < numSolution1; i++) {
     1008                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
     1009
     1010                sumMatrix->data.F64[i][normIndex] = 0.0;
     1011                sumMatrix->data.F64[normIndex][i] = 0.0;
     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;
     1016
     1017            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
     1018            sumVector->data.F64[normIndex] = 0.0;
     1019
     1020            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1021
     1022            solution->data.F64[normIndex] = normValue;
     1023        }
     1024# endif
     1025
     1026        if (!kernels->solution1) {
     1027            kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1028            psVectorInit(kernels->solution1, 0.0);
     1029        }
     1030
     1031        // only update the solutions that we chose to calculate:
     1032        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     1033            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1034            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     1035        }
     1036        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     1037            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1038            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     1039        }
     1040        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     1041            int numKernels = kernels->num;
     1042            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     1043            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     1044            for (int i = 0; i < numKernels * numPoly; i++) {
     1045                kernels->solution1->data.F64[i] = solution->data.F64[i];
     1046            }
     1047        }
     1048
     1049        psFree(solution);
     1050        psFree(sumVector);
     1051        psFree(sumMatrix);
     1052
     1053#ifdef TESTING
     1054        // XXX double-check for NAN in data:
     1055        for (int ix = 0; ix < kernels->solution1->n; ix++) {
     1056            if (!isfinite(kernels->solution1->data.F64[ix])) {
     1057                fprintf (stderr, "WARNING: NAN in vector\n");
     1058            }
     1059        }
     1060#endif
     1061
     1062    } else {
     1063        // Dual convolution solution
     1064
     1065        // Accumulation of stamp matrices/vectors
     1066        psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1067        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
     1068        psImageInit(sumMatrix, 0.0);
     1069        psVectorInit(sumVector, 0.0);
     1070
     1071        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
     1072
     1073        int numStamps = 0;              // Number of good stamps
     1074        for (int i = 0; i < stamps->num; i++) {
     1075            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     1076            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
     1077                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
     1078                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
     1079
     1080                psVectorAppend(norms, stamp->norm);
    11971081
    11981082                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    … …  
    12071091
    12081092#if 1
    1209         int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    1210         calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
     1093        // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     1094        // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
     1095
     1096        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1097        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0);
    12111098#endif
    12121099
    … …  
    12161103
    12171104#if 0
     1105        // Regular, straight-forward solution
     1106        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1107#else
    12181108        {
    1219             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    1220             if (!solution) {
    1221                 psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
     1109            // Solve normalisation and background separately
     1110            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1111            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     1112
     1113#if 0
     1114            psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64);
     1115            psVector *normVector = psVectorAlloc(2, PS_TYPE_F64);
     1116
     1117            normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex];
     1118            normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex];
     1119            normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex];
     1120
     1121            normVector->data.F64[0] = sumVector->data.F64[normIndex];
     1122            normVector->data.F64[1] = sumVector->data.F64[bgIndex];
     1123
     1124            psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN);
     1125
     1126            double normValue = normSolution->data.F64[0];
     1127            double bgValue = normSolution->data.F64[1];
     1128
     1129            psFree(normMatrix);
     1130            psFree(normVector);
     1131            psFree(normSolution);
     1132#endif
     1133
     1134            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
     1135            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
     1136                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
     1137                psFree(stats);
    12221138                psFree(sumMatrix);
    12231139                psFree(sumVector);
    1224                 return NULL;
    1225             }
    1226 
    1227 #ifdef TESTING
    1228             for (int i = 0; i < numParams; i++) {
    1229                 fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
    1230             }
    1231 #endif
    1232 
    1233             int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
    1234             int numKernels = kernels->num; // Number of kernel basis functions
    1235 
    1236 #define MASK_BASIS(INDEX)                                               \
    1237             {                                                           \
    1238                 for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
    1239                     for (int k = 0; k < numParams; k++) {               \
    1240                         sumMatrix->data.F64[index][k] = sumMatrix->data.F64[k][index] = 0.0; \
    1241                     }                                                   \
    1242                     sumMatrix->data.F64[index][index] = 1.0;            \
    1243                     sumVector->data.F64[index] = 0.0;                   \
    1244                 }                                                       \
    1245             }
    1246 
    1247             #define TOL 1.0e-3
    1248             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1249             double norm = fabs(solution->data.F64[normIndex]);  // Normalisation
    1250             double thresh = norm * TOL;                         // Threshold for low parameters
    1251             for (int i = 0; i < numKernels; i++) {
    1252                 // Getting 0th order parameter value.  In the presence of spatial variation, the actual value
    1253                 // of the parameter will vary over the image.  We are in effect getting the value in the
    1254                 // centre of the image.  If we use different polynomial functions (e.g., Chebyshev), we may
    1255                 // have to change this to properly determine the value of the parameter at the centre.
    1256                 double param1 = solution->data.F64[i],
    1257                     param2 = solution->data.F64[numSolution1 + i]; // Parameters of interest
    1258                 bool mask1 = false, mask2 = false;              // Masked the parameter?
    1259                 if (fabs(param1) < thresh) {
    1260                     psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i);
    1261                     MASK_BASIS(i);
    1262                     mask1 = true;
    1263                 }
    1264                 if (fabs(param2) < thresh) {
    1265                     psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i);
    1266                     MASK_BASIS(numSolution1 + i);
    1267                     mask2 = true;
    1268                 }
    1269 
    1270                 if (!mask1 && !mask2) {
    1271                     if (fabs(param1) > fabs(param2)) {
    1272                         psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i);
    1273                         MASK_BASIS(numSolution1 + i);
    1274                     } else {
    1275                         psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i);
    1276                         MASK_BASIS(i);
    1277                     }
    1278                 }
    1279             }
    1280 
    1281 #ifdef TESTING
    1282             psFitsWriteImageSimple ("sumMatrixFix.fits", sumMatrix, NULL);
    1283             psVectorWriteFile("sumVectorFix.dat", sumVector);
    1284 #endif
    1285         }
    1286 #endif /*** kernel-masking block ***/
    1287         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1140                psFree(norms);
     1141                return false;
     1142            }
     1143
     1144            double normValue = stats->robustMedian;
     1145
     1146            psFree(stats);
     1147
     1148            fprintf(stderr, "Norm: %lf\n", normValue);
     1149
     1150            // Solve kernel components
     1151            for (int i = 0; i < numSolution2; i++) {
     1152                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
     1153                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1];
     1154
     1155                sumMatrix->data.F64[i][normIndex] = 0.0;
     1156                sumMatrix->data.F64[normIndex][i] = 0.0;
     1157
     1158                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
     1159                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
     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;
     1164
     1165            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
     1166
     1167            sumVector->data.F64[normIndex] = 0.0;
     1168
     1169            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1170
     1171            solution->data.F64[normIndex] = normValue;
     1172        }
     1173#endif
     1174
    12881175
    12891176#ifdef TESTING
    … …  
    12961183        psFree(sumVector);
    12971184
     1185        psFree(norms);
     1186
    12981187        if (!kernels->solution1) {
    12991188            kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64);
     1189            psVectorInit (kernels->solution1, 0.0);
    13001190        }
    13011191        if (!kernels->solution2) {
    13021192            kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
    1303         }
     1193            psVectorInit (kernels->solution2, 0.0);
     1194        }
     1195
     1196        // only update the solutions that we chose to calculate:
     1197        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     1198            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1199            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     1200        }
     1201        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     1202            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1203            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     1204        }
     1205        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     1206            int numKernels = kernels->num;
     1207            for (int i = 0; i < numKernels * numSpatial; i++) {
     1208                // XXX fprintf (stderr, "keep\n");
     1209                kernels->solution1->data.F64[i] = solution->data.F64[i];
     1210                kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
     1211            }
     1212        }
     1213
    13041214
    13051215        memcpy(kernels->solution1->data.F64, solution->data.F64,
    … …  
    13641274    }
    13651275    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
     1276    if (!isfinite(sum))  return false;
     1277    if (!isfinite(dmax)) return false;
     1278    if (!isfinite(dmin)) return false;
     1279    if (!isfinite(peak)) return false;
     1280
    13661281    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
    13671282    psVectorAppend(fSigRes, sigma/sum);
    … …  
    19731888}
    19741889
     1890
     1891# if 0
     1892
     1893#ifdef TESTING
     1894        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
     1895        psVectorWriteFile ("B.dat", sumVector);
     1896#endif
     1897
     1898# define SVD_ANALYSIS 0
     1899# define COEFF_SIG 0.0
     1900# define SVD_TOL 0.0
     1901
     1902        // Use SVD to determine the kernel coeffs (and validate)
     1903        if (SVD_ANALYSIS) {
     1904
     1905            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
     1906            // sumMatrix * x = sumVector.
     1907
     1908            // we can use any standard matrix inversion to solve this.  However, the basis
     1909            // functions in general have substantial correlation, so that the solution may be
     1910            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
     1911            // system of equations may be statistically ill-conditioned.  Noise in the image
     1912            // will drive insignificant, but correlated, terms in the solution.  To avoid these
     1913            // problems, we can use SVD to identify numerically unconstrained values and to
     1914            // avoid statistically badly determined value.
     1915
     1916            // A = sumMatrix, B = sumVector
     1917            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
     1918            // x = V (1/w) (U^T B)
     1919            // \sigma_x = sqrt(diag(A^{-1}))
     1920            // solve for x and A^{-1} to get x & dx
     1921            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
     1922            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
     1923
     1924            // If I use the SVD trick to re-condition the matrix, I need to break out the
     1925            // kernel and normalization terms from the background term.
     1926            // XXX is this true?  or was this due to an error in the analysis?
     1927
     1928            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1929
     1930            // now pull out the kernel elements into their own square matrix
     1931            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
     1932            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
     1933
     1934            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
     1935                if (ix == bgIndex) continue;
     1936                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
     1937                    if (iy == bgIndex) continue;
     1938                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
     1939                    ky++;
     1940                }
     1941                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
     1942                kx++;
     1943            }
     1944
     1945            psImage *U = NULL;
     1946            psImage *V = NULL;
     1947            psVector *w = NULL;
     1948            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
     1949                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
     1950                return NULL;
     1951            }
     1952
     1953            // calculate A_inverse:
     1954            // Ainv = V * w * U^T
     1955            psImage *wUt  = p_pmSubSolve_wUt (w, U);
     1956            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
     1957            psImage *Xvar = NULL;
     1958            psFree (wUt);
     1959
     1960# ifdef TESTING
     1961            // kernel terms:
     1962            for (int i = 0; i < w->n; i++) {
     1963                fprintf (stderr, "w: %f\n", w->data.F64[i]);
     1964            }
     1965# endif
     1966            // loop over w adding in more and more of the values until chisquare is no longer
     1967            // dropping significantly.
     1968            // XXX this does not seem to work very well: we seem to need all terms even for
     1969            // simple cases...
     1970
     1971            psVector *Xsvd = NULL;
     1972            {
     1973                psVector *Ax = NULL;
     1974                psVector *UtB = NULL;
     1975                psVector *wUtB = NULL;
     1976
     1977                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
     1978                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
     1979                psVectorInit (wMask, 1); // start by masking everything
     1980
     1981                double chiSquareLast = NAN;
     1982                int maxWeight = 0;
     1983
     1984                double Axx, Bx, y2;
     1985
     1986                // XXX this is an attempt to exclude insignificant modes.
     1987                // it was not successful with the ISIS kernel set: removing even
     1988                // the least significant mode leaves additional ringing / noise
     1989                // because the terms are so coupled.
     1990                for (int k = 0; false && (k < w->n); k++) {
     1991
     1992                    // unmask the k-th weight
     1993                    wMask->data.U8[k] = 0;
     1994                    p_pmSubSolve_SetWeights(wApply, w, wMask);
     1995
     1996                    // solve for x:
     1997                    // x = V * w * (U^T * B)
     1998                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1999                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     2000                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     2001
     2002                    // chi-square for this system of equations:
     2003                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     2004                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     2005                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     2006                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     2007                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     2008                    p_pmSubSolve_y2 (&y2, kernels, stamps);
     2009
     2010                    // apparently, this works (compare with the brute force value below
     2011                    double chiSquare = Axx - 2.0*Bx + y2;
     2012                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
     2013                    chiSquareLast = chiSquare;
     2014
     2015                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
     2016                    if (k && !maxWeight && (deltaChi < 1.0)) {
     2017                        maxWeight = k;
     2018                    }
     2019                }
     2020
     2021                // keep all terms or we get extra ringing
     2022                maxWeight = w->n;
     2023                psVectorInit (wMask, 1);
     2024                for (int k = 0; k < maxWeight; k++) {
     2025                    wMask->data.U8[k] = 0;
     2026                }
     2027                p_pmSubSolve_SetWeights(wApply, w, wMask);
     2028
     2029                // solve for x:
     2030                // x = V * w * (U^T * B)
     2031                p_pmSubSolve_UtB (&UtB, U, kernelVector);
     2032                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     2033                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     2034
     2035                // chi-square for this system of equations:
     2036                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     2037                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     2038                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     2039                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     2040                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     2041                p_pmSubSolve_y2 (&y2, kernels, stamps);
     2042
     2043                // apparently, this works (compare with the brute force value below
     2044                double chiSquare = Axx - 2.0*Bx + y2;
     2045                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
     2046
     2047                // re-calculate A^{-1} to get new variances:
     2048                // Ainv = V * w * U^T
     2049                // XXX since we keep all terms, this is identical to Ainv
     2050                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
     2051                Xvar = p_pmSubSolve_VwUt (V, wUt);
     2052                psFree (wUt);
     2053
     2054                psFree (Ax);
     2055                psFree (UtB);
     2056                psFree (wUtB);
     2057                psFree (wApply);
     2058                psFree (wMask);
     2059            }
     2060
     2061            // copy the kernel solutions to the full solution vector:
     2062            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2063            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2064
     2065            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
     2066                if (ix == bgIndex) {
     2067                    solution->data.F64[ix] = 0;
     2068                    solutionErr->data.F64[ix] = 0.001;
     2069                    continue;
     2070                }
     2071                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
     2072                solution->data.F64[ix] = Xsvd->data.F64[kx];
     2073                kx++;
     2074            }
     2075
     2076            psFree (kernelMatrix);
     2077            psFree (kernelVector);
     2078
     2079            psFree (U);
     2080            psFree (V);
     2081            psFree (w);
     2082
     2083            psFree (Ainv);
     2084            psFree (Xsvd);
     2085        } else {
     2086            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
     2087            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
     2088            if (!luMatrix) {
     2089                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
     2090                psFree(solution);
     2091                psFree(sumVector);
     2092                psFree(sumMatrix);
     2093                psFree(luMatrix);
     2094                psFree(permutation);
     2095                return NULL;
     2096            }
     2097
     2098            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
     2099            psFree(luMatrix);
     2100            psFree(permutation);
     2101            if (!solution) {
     2102                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     2103                psFree(solution);
     2104                psFree(sumVector);
     2105                psFree(sumMatrix);
     2106                return NULL;
     2107            }
     2108
     2109            // XXX LUD does not provide A^{-1}?  fake the error for now
     2110            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2111            for (int ix = 0; ix < sumVector->n; ix++) {
     2112                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
     2113            }
     2114        }
     2115
     2116        if (!kernels->solution1) {
     2117            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
     2118            psVectorInit (kernels->solution1, 0.0);
     2119        }
     2120
     2121        // only update the solutions that we chose to calculate:
     2122        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     2123            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     2124            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     2125        }
     2126        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     2127            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     2128            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     2129        }
     2130        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     2131            int numKernels = kernels->num;
     2132            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     2133            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     2134            for (int i = 0; i < numKernels * numPoly; i++) {
     2135                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
     2136                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
     2137                    // XXX fprintf (stderr, "drop\n");
     2138                    kernels->solution1->data.F64[i] = 0.0;
     2139                } else {
     2140                    // XXX fprintf (stderr, "keep\n");
     2141                    kernels->solution1->data.F64[i] = solution->data.F64[i];
     2142                }
     2143            }
     2144        }
     2145        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
     2146        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
     2147
     2148        psFree(solution);
     2149        psFree(sumVector);
     2150        psFree(sumMatrix);
     2151# endif
     2152
     2153#ifdef TESTING
     2154              // XXX double-check for NAN in data:
     2155                for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
     2156                    for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
     2157                        if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
     2158                            fprintf (stderr, "WARNING: NAN in matrix\n");
     2159                        }
     2160                    }
     2161                }
     2162                for (int ix = 0; ix < stamp->vector->n; ix++) {
     2163                    if (!isfinite(stamp->vector->data.F64[ix])) {
     2164                        fprintf (stderr, "WARNING: NAN in vector\n");
     2165                    }
     2166                }
     2167#endif
     2168
     2169#ifdef TESTING
     2170        for (int ix = 0; ix < sumVector->n; ix++) {
     2171            if (!isfinite(sumVector->data.F64[ix])) {
     2172                fprintf (stderr, "WARNING: NAN in vector\n");
     2173            }
     2174        }
     2175#endif
     2176
     2177#ifdef TESTING
     2178        for (int ix = 0; ix < sumVector->n; ix++) {
     2179            if (!isfinite(sumVector->data.F64[ix])) {
     2180                fprintf (stderr, "WARNING: NAN in vector\n");
     2181            }
     2182        }
     2183        {
     2184            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
     2185            psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
     2186            psFree(inverse);
     2187        }
     2188        {
     2189            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
     2190            psImage *Xt = psMatrixTranspose(NULL, X);
     2191            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
     2192            psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
     2193            psFree(X);
     2194            psFree(Xt);
     2195            psFree(XtX);
     2196        }
     2197#endif
     2198

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.h

    r26577 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionIO.c

    r23378 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.c

    r26657 r26747  
    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);
    … …  
    306306
    307307    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
    308                                                               spatialOrder, penalty, mode); // The kernels
     308                                                              spatialOrder, penalty, bounds, mode); // Kernels
    309309    psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    310310
    … …  
    332332pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder,
    333333                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    334                                                       float penalty, pmSubtractionMode mode)
     334                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
    335335{
    336336    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    362362    }
    363363
    364     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
    365366    psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    366367
    … …  
    389390pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
    390391                                               const psVector *fwhmsIN, const psVector *ordersIN,
    391                                                float penalty, pmSubtractionMode mode)
     392                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    392393{
    393394    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    418419    }
    419420
    420     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
    421423    psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    422424
    423     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);
    424427    psFree(params);
    425428
    … …  
    440443
    441444pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
    442                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    443                                                      float penalty, pmSubtractionMode mode)
     445                                                      const psVector *fwhmsIN, const psVector *ordersIN,
     446                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
    444447{
    445448    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    470473    }
    471474
    472     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
    473477    psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    474478
    … …  
    545549
    546550pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
    547                                                 int size, int spatialOrder, float penalty,
     551                                                int size, int spatialOrder, float penalty, psRegion bounds,
    548552                                                pmSubtractionMode mode)
    549553{
    … …  
    559563    kernels->uStop = NULL;
    560564    kernels->vStop = NULL;
     565    kernels->xMin = bounds.x0;
     566    kernels->xMax = bounds.x1;
     567    kernels->yMin = bounds.y0;
     568    kernels->yMax = bounds.y1;
    561569    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    562570    kernels->penalty = penalty;
    … …  
    567575    kernels->bgOrder = 0;
    568576    kernels->mode = mode;
    569     kernels->numCols = 0;
    570     kernels->numRows = 0;
    571577    kernels->solution1 = NULL;
    572578    kernels->solution2 = NULL;
    … …  
    641647}
    642648
    643 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
     649pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, psRegion bounds,
    644650                                               pmSubtractionMode mode)
    645651{
    … …  
    650656
    651657    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
    652                                                               spatialOrder, penalty, mode); // The kernels
     658                                                              spatialOrder, penalty, bounds, mode); // Kernels
    653659    psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
    654660    psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
    … …  
    665671pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
    666672                                               const psVector *fwhms, const psVector *orders,
    667                                                float penalty, pmSubtractionMode mode)
     673                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    668674{
    669675    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    670                                                                   penalty, mode); // Kernels
     676                                                                  penalty, bounds, mode); // Kernels
    671677    if (!kernels) {
    672678        return NULL;
    … …  
    677683
    678684pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
    679                                                float penalty, pmSubtractionMode mode)
     685                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    680686{
    681687    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    698704
    699705    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
    700                                                               spatialOrder, penalty, mode); // The kernels
     706                                                              spatialOrder, penalty, bounds, mode); // Kernels
    701707    kernels->inner = inner;
    702708    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
    … …  
    769775
    770776pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
    771                                                 pmSubtractionMode mode)
     777                                                psRegion bounds, pmSubtractionMode mode)
    772778{
    773779    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    796802
    797803    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
    798                                                               spatialOrder, penalty, mode); // The kernels
     804                                                              spatialOrder, penalty, bounds, mode); // Kernels
    799805    kernels->inner = inner;
    800806    psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
    … …  
    865871pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
    866872                                               const psVector *orders, int inner, float penalty,
    867                                                pmSubtractionMode mode)
     873                                               psRegion bounds, pmSubtractionMode mode)
    868874{
    869875    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    878884
    879885    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    880                                                                   penalty, mode); // Kernels
     886                                                                  penalty, bounds, mode); // Kernels
    881887    kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
    882888    psStringPrepend(&kernels->description, "GUNK=");
    … …  
    894900// RINGS --- just what it says
    895901pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
    896                                                 float penalty, pmSubtractionMode mode)
     902                                                float penalty, psRegion bounds, pmSubtractionMode mode)
    897903{
    898904    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    925931
    926932    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
    927                                                               spatialOrder, penalty, mode); // The kernels
     933                                                              spatialOrder, penalty, bounds, mode); // Kernels
    928934    kernels->inner = inner;
    929935    psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
    … …  
    10471053pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
    10481054                                                   const psVector *fwhms, const psVector *orders, int inner,
    1049                                                    int binning, int ringsOrder, float penalty,
     1055                                                   int binning, int ringsOrder, float penalty, psRegion bounds,
    10501056                                                   pmSubtractionMode mode)
    10511057{
    10521058    switch (type) {
    10531059      case PM_SUBTRACTION_KERNEL_POIS:
    1054         return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode);
     1060        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, bounds, mode);
    10551061      case PM_SUBTRACTION_KERNEL_ISIS:
    1056         return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
     1062        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10571063      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
    1058         return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, mode);
     1064        return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10591065      case PM_SUBTRACTION_KERNEL_HERM:
    1060         return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode);
     1066        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10611067      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
    1062         return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode);
     1068        return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10631069      case PM_SUBTRACTION_KERNEL_SPAM:
    1064         return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
     1070        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, bounds, mode);
    10651071      case PM_SUBTRACTION_KERNEL_FRIES:
    1066         return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode);
     1072        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, bounds, mode);
    10671073      case PM_SUBTRACTION_KERNEL_GUNK:
    1068         return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode);
     1074        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, bounds, mode);
    10691075      case PM_SUBTRACTION_KERNEL_RINGS:
    1070         return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode);
     1076        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode);
    10711077      default:
    10721078        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
    … …  
    11041110
    11051111pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder,
    1106                                                           pmSubtractionMode mode)
     1112                                                          psRegion bounds, pmSubtractionMode mode)
    11071113{
    11081114    PS_ASSERT_STRING_NON_EMPTY(description, NULL);
    … …  
    11291135
    11301136    type = pmSubtractionKernelsTypeFromString (description);
    1131     char *ptr = strchr(description, '(');
     1137    char *ptr = strchr(description, '(') + 1;
    11321138    psAssert (ptr, "description is missing kernel parameters");
    11331139
    … …  
    11571163        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    11581164        penalty = parseStringFloat(ptr);
    1159 
    1160         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
    1161 
     1165        break;
    11621166      case PM_SUBTRACTION_KERNEL_RINGS:
    11631167        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
    … …  
    11661170        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    11671171        PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
    1168         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
     1172        break;
    11691173      default:
    11701174        psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
    11711175    }
    1172     return NULL;
     1176
     1177    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning,
     1178                                        ringsOrder, penalty, bounds, mode);
    11731179}
    11741180
    … …  
    12461252    out->bgOrder = in->bgOrder;
    12471253    out->mode = in->mode;
    1248     out->numCols = in->numCols;
    1249     out->numRows = in->numRows;
     1254    out->xMin = in->xMin;
     1255    out->xMax = in->xMax;
     1256    out->yMin = in->yMin;
     1257    out->yMax = in->yMax;
    12501258    out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
    12511259    out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionKernels.h

    r26593 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionMask.c

    r24257 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionMask.h

    r23851 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.c

    r26667 r26747  
    3333# else
    3434# define SUBMODE PM_SUBTRACTION_EQUATION_ALL
    35 // # define SUBMODE PM_SUBTRACTION_EQUATION_KERNELS
    3635# endif
    3736
    … …  
    7170                                 const psImage *subMask, // Mask for subtraction, or NULL
    7271                                 psImage *variance,  // Variance map
    73                                  const psRegion *region, // Region of interest, or NULL
     72                                 const psRegion *region, // Region of interest
    7473                                 float thresh1,  // Threshold for stamp finding on readout 1
    7574                                 float thresh2,  // Threshold for stamp finding on readout 2
    7675                                 float stampSpacing, // Spacing between stamps
     76                                 float normFrac,     // Fraction of flux in window for normalisation window
    7777                                 float sysError,     // Relative systematic error in images
    7878                                 float skyError,     // Relative systematic error in images
    … …  
    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
    … …  
    8794
    8895    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
    89                                       size, footprint, stampSpacing, sysError, skyError, mode);
     96                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
    9097    if (!*stamps) {
    9198        psError(psErrorCodeLast(), false, "Unable to find stamps.");
    … …  
    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;
    … …  
    110117                                  const pmReadout *ro1, const pmReadout *ro2, // Input images
    111118                                  int stride, // Size for convolution patches
     119                                  float normFrac,           // Fraction of window for normalisation window
    112120                                  float sysError,           // Systematic error in images
    113121                                  float skyError,           // Systematic error in images
    … …  
    158166    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
    159167    PS_ASSERT_INT_NONNEGATIVE(stride, false);
     168    if (isfinite(normFrac)) {
     169        PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false);
     170        PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false);
     171    }
    160172    if (isfinite(sysError)) {
    161173        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
    … …  
    210222}
    211223
     224bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {
     225
     226    if (!readout) return true;
     227
     228    psImage *image = readout->image;
     229    psImage *mask  = readout->mask;
     230    psImage *variance = readout->variance;
     231    for (int y = 0; y < image->numRows; y++) {
     232        for (int x = 0; x < image->numCols; x++) {
     233            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;
     234            bool valid = false;
     235            valid = isfinite(image->data.F32[y][x]);
     236            if (variance) {
     237                valid &= isfinite(variance->data.F32[y][x]);
     238            }
     239            if (valid) continue;
     240            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;
     241        }
     242    }
     243
     244    return true;
     245}
    212246
    213247bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
    … …  
    282316    }
    283317
    284     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, kernelError,
     318    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError,
    285319                               maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
    286320        psFree(kernels);
    … …  
    289323    }
    290324
    291     psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
     325    pmSubtractionMaskInvalid(ro1, maskVal);
     326    pmSubtractionMaskInvalid(ro2, maskVal);
     327
     328    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     329
     330    psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0,
    292331                                         badFrac, mode); // Subtraction mask
    293332    if (!subMask) {
    … …  
    348387                        int inner, int ringsOrder, int binning, float penalty,
    349388                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    350                         int iter, float rej, float sysError, float skyError, float kernelError, psImageMaskType maskVal,
    351                         psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac,
    352                         float badFrac, pmSubtractionMode subMode)
     389                        int iter, float rej, float normFrac, float sysError, float skyError,
     390                        float kernelError, psImageMaskType maskVal, psImageMaskType maskBad,
     391                        psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
    353392{
    354     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, skyError, kernelError,
     393    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError,
    355394                               maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
    356395        return false;
    … …  
    411450    // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
    412451    psImage *subMask = NULL;            // Mask for subtraction
    413     psRegion *region = NULL;            // Iso-kernel region
     452    psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region
    414453    psString regionString = NULL;       // String for region
    415454    pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
    … …  
    424463    memCheck("start");
    425464
    426     subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
    427                                 badFrac, subMode);
     465    pmSubtractionMaskInvalid(ro1, maskVal);
     466    pmSubtractionMaskInvalid(ro2, maskVal);
     467
     468    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
     469
     470    subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode);
    428471    if (!subMask) {
    429472        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    … …  
    435478    // Get region of interest
    436479    int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
    437     float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
     480    float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions
    438481    if (isfinite(regionSize) && regionSize != 0.0) {
    439         xRegions = numCols / regionSize + 1;
    440         yRegions = numRows / regionSize + 1;
    441         xRegionSize = (float)numCols / (float)xRegions;
    442         yRegionSize = (float)numRows / (float)yRegions;
    443         region = psRegionAlloc(NAN, NAN, NAN, NAN);
    444     }
    445 
     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;
     489    }
     490
     491    // General background subtraction and measurement of stamp threshold
    446492    float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
    447493    {
    448         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
    449495        if (ro1) {
     496            psStatsInit(bg);
    450497            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
    451498                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
    … …  
    453500                goto MATCH_ERROR;
    454501            }
    455             stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
     502            stampThresh1 = threshold * bg->robustStdev;
     503            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
    456504        }
    457505        if (ro2) {
     506            psStatsInit(bg);
    458507            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
    459508                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
    … …  
    461510                goto MATCH_ERROR;
    462511            }
    463             stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
    464         }
     512            stampThresh2 = threshold * bg->robustStdev;
     513            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
     514       }
    465515        psFree(bg);
    466516    }
     517
     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_2 || 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
    467558
    468559    // Iterate over iso-kernel regions
    … …  
    471562            psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
    472563                    j * xRegions + i + 1, xRegions * yRegions);
    473             if (region) {
    474                 *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
    475                                       (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
    476                 psFree(regionString);
    477                 regionString = psRegionToString(*region);
    478                 psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
    479                         regionString, numCols, numRows);
    480             }
     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);
    481572
    482573            if (stampsName && strlen(stampsName) > 0) {
    483574                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
    484                                                         footprint, stampSpacing, sysError, skyError, subMode);
     575                                                        footprint, stampSpacing, normFrac,
     576                                                        sysError, skyError, subMode);
    485577            } else if (sources) {
    486578                stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
    487                                                            footprint, stampSpacing, sysError, skyError, subMode);
     579                                                           footprint, stampSpacing, normFrac,
     580                                                           sysError, skyError, subMode);
    488581            }
    489582
    490583            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    491584            // doesn't matter.
    492             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
    493                                       stampSpacing, sysError, skyError, size, footprint, subMode)) {
     585            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
     586                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
    494587                goto MATCH_ERROR;
    495588            }
    … …  
    498591            // generate the window function from the set of stamps
    499592            if (!pmSubtractionStampsGetWindow(stamps, size)) {
     593                psError(PS_ERR_UNKNOWN, false, "Unable to get stamp window.");
    500594                goto MATCH_ERROR;
    501595            }
    … …  
    503597            // Define kernel basis functions
    504598            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    505                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    506                                                           stamps, footprint, optThreshold, penalty, subMode);
     599                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
     600                                                          optFWHMs, optOrder, stamps, footprint,
     601                                                          optThreshold, penalty, bounds, subMode);
    507602                if (!kernels) {
    508603                    psErrorClear();
    … …  
    513608                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    514609                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    515                                                        inner, binning, ringsOrder, penalty, subMode);
     610                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
    516611                // pmSubtractionVisualShowKernels(kernels);
    517612            }
    … …  
    563658
    564659                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    565                                           stampThresh1, stampThresh2, stampSpacing, sysError, skyError,
    566                                           size, footprint, subMode)) {
     660                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
     661                                          sysError, skyError, size, footprint, subMode)) {
    567662                    goto MATCH_ERROR;
    568663                }
    … …  
    570665                // generate the window function from the set of stamps
    571666                if (!pmSubtractionStampsGetWindow(stamps, size)) {
     667                    psError(PS_ERR_UNKNOWN, false, "Unable to get stamps window.");
    572668                    goto MATCH_ERROR;
    573669                }

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionMatch.h

    r26667 r26747  
    3939                        int iter,       ///< Rejection iterations
    4040                        float rej,      ///< Rejection threshold
     41                        float normFrac, ///< Fraction of flux in window for normalisation window
    4142                        float sysError, ///< Relative systematic error in images
    4243                        float skyError, ///< Relative systematic error in images

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionParams.c

    r26491 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionParams.h

    r18287 r26747  
    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/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.c

    r26562 r26747  
    176176
    177177pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
    178                                                     int footprint, float spacing, float sysErr, float skyErr)
     178                                                    int footprint, float spacing, float normFrac,
     179                                                    float sysErr, float skyErr)
    179180{
    180181    pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
    … …  
    217218    list->flux = NULL;
    218219    list->window = NULL;
     220    list->normFrac = normFrac;
     221    list->normWindow = 0;
    219222    list->footprint = footprint;
    220223    list->sysErr = sysErr;
    … …  
    239242    out->window = psMemIncrRefCounter(in->window);
    240243    out->footprint = in->footprint;
     244    out->normWindow = in->normWindow;
    241245
    242246    for (int i = 0; i < num; i++) {
    … …  
    308312    stamp->matrix = NULL;
    309313    stamp->vector = NULL;
     314    stamp->norm = NAN;
    310315
    311316    return stamp;
    … …  
    316321                                                const psImage *image2, const psImage *subMask,
    317322                                                const psRegion *region, float thresh1, float thresh2,
    318                                                 int size, int footprint, float spacing, float sysErr, float skyErr,
    319                                                 pmSubtractionMode mode)
     323                                                int size, int footprint, float spacing, float normFrac,
     324                                                float sysErr, float skyErr, pmSubtractionMode mode)
    320325{
    321326    if (!image1 && !image2) {
    … …  
    373378
    374379    if (!stamps) {
    375         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr);
     380        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
     381                                             normFrac, sysErr, skyErr);
    376382    }
    377383
    … …  
    407413                // Take stamps off the top of the (sorted) list
    408414                for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
    409                     int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
    410 
    411415                    // Chop off the top of the list
    412416                    xList->n = j;
    … …  
    414418                    fluxList->n = j;
    415419
     420#if 0
    416421                    // Fish around a bit to see if we can find a pixel that isn't masked
     422                    // This is not a good idea if we're using the window feature
    417423                    psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n",
    418424                            i, xCentre, yCentre);
    419425
    420426                    // Search bounds
     427                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
    421428                    int search = footprint - size; // Search radius
    422429                    int xMin = PS_MAX(border, xCentre - search);
    … …  
    427434                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
    428435                                            subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
    429 
    430                     // XXX reset for a test:
    431                     xStamp = xList->data.F32[j];
    432                     yStamp = yList->data.F32[j];
     436#else
     437                    // Only search the exact centre pixel
     438                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
     439                                            subMask, xList->data.F32[j], xList->data.F32[j],
     440                                            yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
     441#endif
    433442                    // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
    434443                }
    … …  
    486495                                               const psImage *image, const psImage *subMask,
    487496                                               const psRegion *region, int size, int footprint,
    488                                                float spacing, float sysErr, float skyErr, pmSubtractionMode mode)
     497                                               float spacing, float normFrac, float sysErr, float skyErr,
     498                                               pmSubtractionMode mode)
    489499
    490500{
    … …  
    507517    pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
    508518                                                                 region, footprint, spacing,
    509                                                                  sysErr, skyErr); // Stamp list
     519                                                                 normFrac, sysErr, skyErr); // Stamp list
    510520    int numStamps = stamps->num;        // Number of stamps
    511521
    … …  
    612622    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    613623
    614     int size = kernelSize + stamps->footprint; // Size of postage stamps
     624    int size = stamps->footprint; // Size of postage stamps
    615625
    616626    psFree (stamps->window);
    … …  
    641651    // storage vector for flux data
    642652    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    643     psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
    644 
    645     double maxValue = 0.0;
     653    psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
     654    psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
    646655
    647656    // generate the window pixels
     657    double sum = 0.0;                   // Sum inside the window
     658    float maxValue = 0.0;               // Maximum value, for normalisation
    648659    for (int y = -size; y <= size; y++) {
    649660        for (int x = -size; x <= size; x++) {
    650661
    651             flux->n = 0;
     662            flux1->n = 0;
     663            flux2->n = 0;
    652664            for (int i = 0; i < stamps->num; i++) {
    653665                pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    … …  
    656668                if (!stamp->image2) continue;
    657669
    658                 psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
    659                 psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
     670                psVectorAppend(flux1, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
     671                psVectorAppend(flux2, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
    660672            }
    661673
    662674            psStatsInit (stats);
    663             if (!psVectorStats (stats, flux, NULL, NULL, 0)) {
     675            if (!psVectorStats (stats, flux1, NULL, NULL, 0)) {
    664676                psAbort ("failed to generate stats");
    665677            }
    666             stamps->window->kernel[y][x] = stats->robustMedian;
    667             if (maxValue < stats->robustMedian) {
    668                 maxValue = stats->robustMedian;
    669             }
    670         }
     678            float f1 = stats->robustMedian;
     679            psStatsInit (stats);
     680            if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
     681                psAbort ("failed to generate stats");
     682            }
     683            float f2 = stats->robustMedian;
     684
     685            stamps->window->kernel[y][x] = f1 + f2;
     686            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
     687                sum += stamps->window->kernel[y][x];
     688            }
     689            maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
     690        }
     691    }
     692
     693    psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
     694            sum, (1.0 - stamps->normFrac) * sum);
     695    bool done = false;
     696    for (int radius = 1; radius <= size && !done; radius++) {
     697        double within = 0.0;
     698        for (int y = -radius; y <= radius; y++) {
     699            for (int x = -radius; x <= radius; x++) {
     700                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
     701                    within += stamps->window->kernel[y][x];
     702                }
     703            }
     704        }
     705        psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
     706        if (within > (1.0 - stamps->normFrac) * sum) {
     707            stamps->normWindow = radius;
     708            done = true;
     709        }
     710    }
     711
     712    psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
     713    if (stamps->normWindow == 0 || stamps->normWindow >= size) {
     714        psError(PS_ERR_UNKNOWN, true, "Unable to determine normalisation window size.");
     715        return false;
    671716    }
    672717
    … …  
    678723    }
    679724
    680 # ifdef TESTING
     725#if 0
    681726    {
    682727        psFits *fits = psFitsOpen ("window.fits", "w");
    … …  
    684729        psFitsClose (fits);
    685730    }
    686 # endif
     731#endif
    687732
    688733    psFree (stats);
    689     psFree (flux);
     734    psFree (flux1);
     735    psFree(flux2);
    690736    psFree (norm1);
    691737    psFree (norm2);
    … …  
    694740
    695741bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
    696                                 psImage *variance, int kernelSize)
     742                                psImage *variance, int kernelSize, psRegion bounds)
    697743{
    698744    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    711757    }
    712758
    713     int numCols = image1->numCols, numRows = image1->numRows; // Size of images
    714759    int size = kernelSize + stamps->footprint; // Size of postage stamps
    715760
    … …  
    720765        }
    721766
    722         if (isnan(stamp->xNorm)) {
    723             stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
    724         }
    725         if (isnan(stamp->yNorm)) {
    726             stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
    727         }
     767        p_pmSubtractionPolynomialNormCoords(&stamp->xNorm, &stamp->yNorm, stamp->x, stamp->y,
     768                                            bounds.x0, bounds.x1, bounds.y0, bounds.y1);
    728769
    729770        int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
    730         if (x < size || x > numCols - size || y < size || y > numRows - size) {
    731             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);
    732773            return false;
    733774        }
    … …  
    788829pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
    789830                                                          const psImage *subMask, const psRegion *region,
    790                                                           int size, int footprint, float spacing, float sysErr, float skyErr,
     831                                                          int size, int footprint, float spacing,
     832                                                          float normFrac, float sysErr, float skyErr,
    791833                                                          pmSubtractionMode mode)
    792834{
    … …  
    819861
    820862    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
    821                                                             footprint, spacing, sysErr, skyErr, mode); // Stamps
     863                                                            footprint, spacing, normFrac,
     864                                                            sysErr, skyErr, mode); // Stamps
    822865    psFree(x);
    823866    psFree(y);
    … …  
    833876pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
    834877                                                       const psImage *subMask, const psRegion *region,
    835                                                        int size, int footprint, float spacing, float sysErr, float skyErr,
    836                                                        pmSubtractionMode mode)
     878                                                       int size, int footprint, float spacing, float normFrac,
     879                                                       float sysErr, float skyErr, pmSubtractionMode mode)
    837880{
    838881    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
    … …  
    851894
    852895    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
    853                                                             spacing, sysErr, skyErr, mode);
     896                                                            spacing, normFrac, sysErr, skyErr, mode);
    854897    psFree(data);
    855898

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionStamps.h

    r26578 r26747  
    2525    int footprint;                      ///< Half-size of stamps
    2626    psKernel *window;                   ///< window function generated from ensemble of stamps
     27    float normFrac;                     ///< Fraction of flux in window for normalisation window
     28    int normWindow;                     ///< Size of window for measuring normalisation
    2729    float sysErr;                       ///< Systematic error
    2830    float skyErr;                       ///< increase effective readnoise
    … …  
    3032
    3133/// Allocate a list of stamps
    32 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image
    33                                                     int numRows, // Number of rows in image
    34                                                     const psRegion *region, // Region for stamps, or NULL
    35                                                     int footprint, // Half-size of stamps
    36                                                     float spacing, // Rough average spacing between stamps
    37                                                     float sysErr,  // Relative systematic error or NAN
    38                                                     float skyErr  // Relative systematic error or NAN
     34pmSubtractionStampList *pmSubtractionStampListAlloc(
     35    int numCols, // Number of columns in image
     36    int numRows, // Number of rows in image
     37    const psRegion *region, // Region for stamps, or NULL
     38    int footprint, // Half-size of stamps
     39    float spacing, // Rough average spacing between stamps
     40    float normFrac, // Fraction of flux in window for normalisation window
     41    float sysErr,  // Relative systematic error or NAN
     42    float skyErr  // Relative systematic error or NAN
    3943    );
    4044
    … …  
    7882    psImage *matrix;                    ///< Least-squares matrix, or NULL
    7983    psVector *vector;                   ///< Least-squares vector, or NULL
     84    double norm;                        ///< Normalisation difference
    8085    pmSubtractionStampStatus status;    ///< Status of stamp
    8186} pmSubtractionStamp;
    … …  
    8590
    8691/// Find stamps on an image
    87 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL
    88                                                 const psImage *image1, ///< Image for which to find stamps
    89                                                 const psImage *image2, ///< Image for which to find stamps
    90                                                 const psImage *mask, ///< Mask, or NULL
    91                                                 const psRegion *region, ///< Region to search, or NULL
    92                                                 float thresh1, ///< Threshold for stamps in image 1
    93                                                 float thresh2, ///< Threshold for stamps in image 2
    94                                                 int size, ///< Kernel half-size
    95                                                 int footprint, ///< Half-size for stamps
    96                                                 float spacing, ///< Rough spacing for stamps
    97                                                 float sysErr,  ///< Relative systematic error in images
    98                                                 float skyErr,  ///< Relative systematic error in images
    99                                                 pmSubtractionMode mode ///< Mode for subtraction
     92pmSubtractionStampList *pmSubtractionStampsFind(
     93    pmSubtractionStampList *stamps, ///< Output stamps, or NULL
     94    const psImage *image1, ///< Image for which to find stamps
     95    const psImage *image2, ///< Image for which to find stamps
     96    const psImage *mask, ///< Mask, or NULL
     97    const psRegion *region, ///< Region to search, or NULL
     98    float thresh1, ///< Threshold for stamps in image 1
     99    float thresh2, ///< Threshold for stamps in image 2
     100    int size, ///< Kernel half-size
     101    int footprint, ///< Half-size for stamps
     102    float spacing, ///< Rough spacing for stamps
     103    float normFrac, // Fraction of flux in window for normalisation window
     104    float sysErr,  ///< Relative systematic error in images
     105    float skyErr,  ///< Relative systematic error in images
     106    pmSubtractionMode mode ///< Mode for subtraction
    100107    );
    101108
    102109/// Set stamps based on a list of x,y
    103 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp
    104                                                const psVector *y, ///< y coordinates for each stamp
    105                                                const psImage *image, ///< Image for flux of stamp
    106                                                const psImage *mask, ///< Mask, or NULL
    107                                                const psRegion *region, ///< Region to search, or NULL
    108                                                int size, ///< Kernel half-size
    109                                                int footprint, ///< Half-size for stamps
    110                                                float spacing, ///< Rough spacing for stamps
    111                                                float sysErr,  ///< Systematic error in images
    112                                                float skyErr,  ///< Systematic error in images
    113                                                pmSubtractionMode mode ///< Mode for subtraction
     110pmSubtractionStampList *pmSubtractionStampsSet(
     111    const psVector *x, ///< x coordinates for each stamp
     112    const psVector *y, ///< y coordinates for each stamp
     113    const psImage *image, ///< Image for flux of stamp
     114    const psImage *mask, ///< Mask, or NULL
     115    const psRegion *region, ///< Region to search, or NULL
     116    int size, ///< Kernel half-size
     117    int footprint, ///< Half-size for stamps
     118    float spacing, ///< Rough spacing for stamps
     119    float normFrac, // Fraction of flux in window for normalisation window
     120    float sysErr,  ///< Systematic error in images
     121    float skyErr,  ///< Systematic error in images
     122    pmSubtractionMode mode ///< Mode for subtraction
    114123    );
    115124
    … …  
    123132    int footprint,                      ///< Half-size for stamps
    124133    float spacing,                      ///< Rough spacing for stamps
     134    float normFrac, // Fraction of flux in window for normalisation window
    125135    float sysErr,                       ///< Systematic error in images
    126136    float skyErr,                       ///< Systematic error in images
    … …  
    137147    int footprint,                      ///< Half-size for stamps
    138148    float spacing,                      ///< Rough spacing for stamps
     149    float normFrac, // Fraction of flux in window for normalisation window
    139150    float sysErr,                       ///< Systematic error in images
    140151    float skyErr,                       ///< Systematic error in images
    … …  
    142153    );
    143154
    144 bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize);
     155/// Calculate the window and normalisation window from the stamps
     156bool pmSubtractionStampsGetWindow(
     157    pmSubtractionStampList *stamps,     ///< List of stamps
     158    int kernelSize                      ///< Half-size of kernel
     159    );
    145160
    146161/// Extract stamps from the images
    … …  
    149164                                psImage *image2, ///< Input image (or NULL)
    150165                                psImage *variance, ///< Variance map
    151                                 int kernelSize ///< Kernel half-size
     166                                int kernelSize, ///< Kernel half-size
     167                                psRegion bounds ///< Bounds of validity
    152168    );
    153169

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.c

    r26632 r26747  
    2727  psFree (detections->oldPeaks);
    2828  psFree (detections->oldFootprints);
     29
     30  psFree (detections->newSources);
     31  psFree (detections->allSources);
    2932  return;
    3033}
    … …  
    4043    detections->oldPeaks      = NULL;
    4144    detections->oldFootprints = NULL;
     45    detections->newSources    = NULL;
     46    detections->allSources    = NULL;
    4247    detections->last          = 0;
    4348

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmDetections.h

    r26631 r26747  
    2121typedef struct {
    2222  psArray *footprints;        // collection of footprints in the image
     23  psArray *oldFootprints;     // collection of footprints previously found
    2324  psArray *peaks;             // collection of all peaks contained by the footprints
    2425  psArray *oldPeaks;          // collection of all peaks previously found
    25   psArray *oldFootprints;     // collection of footprints previously found
     26  psArray *newSources;        // collection of sources
     27  psArray *allSources;        // collection of sources
    2628  int last;
    2729} pmDetections;

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.c

    r26623 r26747  
    275275// psphot-specific function which applies the recipe values
    276276// only apply selection to sources within specified region
    277 pmPSFClump pmSourcePSFClump(psRegion *region, psArray *sources, psMetadata *recipe)
     277pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX)
    278278{
    279279    psTrace("psModules.objects", 10, "---- begin ----\n");
    … …  
    285285
    286286    PS_ASSERT_PTR_NON_NULL(sources, errorClump);
    287     PS_ASSERT_PTR_NON_NULL(recipe, errorClump);
    288 
    289     bool status = true;                 // Status of MD lookup
    290     float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    291     if (!status) {
    292         PSF_SN_LIM = 0;
    293     }
    294     float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
    295     if (!status) {
    296         PSF_CLUMP_GRID_SCALE = 0.1;
    297     }
    298287
    299288    // find the sigmaX, sigmaY clump
    300289    {
    301         psF32 SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
    302         if (!status) {
    303             psWarning("MOMENTS_SX_MAX not set in recipe");
    304             SX_MAX = 10.0;
    305         }
    306         psF32 SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
    307         if (!status) {
    308             psWarning("MOMENTS_SY_MAX not set in recipe");
    309             SY_MAX = 10.0;
    310         }
    311         psF32 AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX");
    312         if (!status) {
    313             psWarning("MOMENTS_AR_MAX not set in recipe");
    314             AR_MAX =  3.0;
    315         }
    316290        psF32 AR_MIN = 1.0 / AR_MAX;
    317291
    … …  
    399373        psfClump.nSigma = stats->sampleStdev;
    400374
    401         const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
    402         if (keep_psf_clump)
    403         {
    404             psMetadataAdd(recipe, PS_LIST_TAIL,
    405                           "PSF_CLUMP", PS_DATA_IMAGE, "Image of PSF coefficients", splane);
     375        if (savedImage) {
     376            *savedImage = psMemIncrRefCounter(splane);
    406377        }
    407378        psFree (splane);
    … …  
    530501*****************************************************************************/
    531502
    532 bool pmSourceRoughClass(psRegion *region, psArray *sources, psMetadata *recipe, pmPSFClump clump, psImageMaskType maskSat)
     503bool pmSourceRoughClass(psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_NSIGMA, pmPSFClump clump, psImageMaskType maskSat)
    533504{
    534505    psTrace("psModules.objects", 10, "---- begin ----");
    535506
    536507    PS_ASSERT_PTR_NON_NULL(sources, false);
    537     PS_ASSERT_PTR_NON_NULL(recipe, false);
    538508
    539509    int Nsat     = 0;
    … …  
    548518    psVector *starsn_peaks = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    549519    psVector *starsn_moments = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    550 
    551     // get basic parameters, or set defaults
    552     bool status;
    553     float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
    554     if (!status) PSF_SN_LIM = 20.0;
    555     float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    556     if (!status) PSF_CLUMP_NSIGMA = 1.5;
    557 
    558     // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    559520
    560521    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmSource.h

    r26594 r26747  
    176176 *
    177177 * The return value indicates the success (TRUE) of the operation.
    178  *
    179  * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
    180  * XXX: Save the clump parameters on the Metadata (TBD)
    181  *
    182  */
     178 */
     179
    183180pmPSFClump pmSourcePSFClump(
     181    psImage **savedImage,
    184182    psRegion *region,                   ///< restrict measurement to specified region
    185183    psArray *source,                    ///< The input pmSource
    186     psMetadata *metadata                ///< Contains classification parameters
     184    float PSF_SN_LIM,
     185    float PSF_CLUMP_GRID_SCALE,
     186    psF32 SX_MAX,
     187    psF32 SY_MAX,
     188    psF32 AR_MAX
    187189);
    188190
    … …  
    200202    psRegion *region,                   ///< restrict measurement to specified region
    201203    psArray *sources,                    ///< The input pmSources
    202     psMetadata *metadata,               ///< Contains classification parameters
     204    float PSF_SN_LIM,                    ///< min S/N for source to be used for PSF model
     205    float PSF_CLUMP_NSIGMA,              ///< size of region around peak of clump for PSF stars
    203206    pmPSFClump clump,                   ///< Statistics about the PSF clump
    204207    psImageMaskType maskSat             ///< Mask value for saturated pixels

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmSourceIO.c

    r26260 r26747  
    4040#include "pmPSF.h"
    4141#include "pmModel.h"
     42#include "pmDetections.h"
    4243#include "pmSource.h"
    4344#include "pmModelClass.h"
    … …  
    344345
    345346    // if sources is NULL, write out an empty table
    346     // input / output sources are stored on the readout->analysis as "PSPHOT.SOURCES" -- a better name might be something like PM_SOURCE_DATA
    347     psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     347    // input / output sources are stored on the readout->analysis as "PSPHOT.DETECTIONS"
     348
     349    psArray *sources = NULL;
     350    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     351    if (detections) {
     352        sources = detections->allSources;
     353    }
    348354    if (!sources) {
     355        detections = pmDetectionsAlloc();
    349356        sources = psArrayAlloc(0);
    350         psMetadataAddArray(readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "Blank array of sources", sources);
    351         psFree(sources); // Held onto by the metadata, so we can continue to use
     357        detections->allSources = sources;
     358        psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "Blank array of sources", detections);
     359        psFree(detections); // Held onto by the metadata, so we can continue to use
    352360    }
    353361
    … …  
    10311039    }
    10321040    readout->data_exists = true;
    1033     status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "input sources", sources);
    1034     psFree (sources);
     1041
     1042    pmDetections *detections = pmDetectionsAlloc();
     1043    detections->allSources = sources;
     1044    status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "input sources", detections);
     1045    psFree (detections);
    10351046    return true;
    10361047}
    … …  
    11241135    bool status;
    11251136
    1126     // select the psf of interest
    1127     pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
    1128     if (!psf) return false;
    1129     return true;
    1130 }
    1131 
    1132 
     1137    // select the detections of interest
     1138    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     1139    if (!detections) return false;
     1140    if (!detections->allSources) return false;
     1141    return true;
     1142}
     1143
     1144

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotApResid.c

    r20937 r26747  
    3434#include "pmPSF.h"
    3535#include "pmModel.h"
     36#include "pmDetections.h"
    3637#include "pmSource.h"
    3738#include "pmSourcePlots.h"
    … …  
    5354    PS_ASSERT_PTR_NON_NULL(layout, false);
    5455
     56    bool status;
    5557    Graphdata graphdata;
    5658    KapaSection section;
    … …  
    6163    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    6264
    63     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    64     if (sources == NULL)
    65         return false;
     65    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     66    if (detections == NULL) return false;
     67
     68    psArray *sources = detections->allSources;
     69    if (sources == NULL) return false;
    6670
    6771    int kapa = pmKapaOpen (false);

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotMoments.c

    r20937 r26747  
    3737#include "pmPSF.h"
    3838#include "pmModel.h"
     39#include "pmDetections.h"
    3940#include "pmSource.h"
    4041#include "pmSourcePlots.h"
    … …  
    5455    PS_ASSERT_PTR_NON_NULL(layout, false);
    5556
     57    bool status;
    5658    Graphdata graphdata;
    5759    KapaSection section;
    … …  
    6264    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    6365
    64     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    65     if (sources == NULL)
    66         return false;
     66    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     67    if (detections == NULL) return false;
     68
     69    psArray *sources = detections->allSources;
     70    if (sources == NULL) return false;
    6771
    6872    int kapa = pmKapaOpen (false);

  • branches/eam_branches/psModules.stack.20100120/src/objects/pmSourcePlotPSFModel.c

    r20937 r26747  
    3737#include "pmPSF.h"
    3838#include "pmModel.h"
     39#include "pmDetections.h"
    3940#include "pmSource.h"
    4041#include "pmSourcePlots.h"
    … …  
    5657    PS_ASSERT_PTR_NON_NULL(layout, false);
    5758
     59    bool status;
    5860    Graphdata graphdata;
    5961    KapaSection section;
    … …  
    6466    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    6567
    66     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    67     if (sources == NULL)
    68         return false;
     68    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     69    if (detections == NULL) return false;
     70
     71    psArray *sources = detections->allSources;
     72    if (sources == NULL) return false;
    6973
    7074    int kapa = pmKapaOpen (false);

  • branches/eam_branches/psModules.stack.20100120/src/psmodules.h

    r26486 r26747  
    2929#include <pmConfigDump.h>
    3030#include <pmConfigRun.h>
     31#include <pmConfigRecipeValue.h>
    3132#include <pmVersion.h>
    3233
Note: See TracChangeset for help on using the changeset viewer.