IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26703


Ignore:
Timestamp:
Jan 28, 2010, 4:16:19 AM (16 years ago)
Author:
Paul Price
Message:

Solve for the normalisation of the convolution kernel completely separately from that of the other kernel terms. We use a 'normalisation window', basically an aperture in which we measure the flux of the stamps to determine the normalisation value. This introduces a new recipe parameter, NORM.FRAC, which is the fraction of the total in the window at which we set the normalisation window/aperture. The quality of the result seems to be rather sensitive to this value. Also found a bug in the imposition of the penalty function in dual convolution mode: it wasn't being set for the second solution part. This seems to produce decent dual convolution subtractions.

Location:
branches/eam_branches/20091201
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/ippconfig/recipes/ppSub.config

    r26689 r26703  
    22
    33KERNEL.TYPE     STR     ISIS_RADIAL     # Kernel type to use (POIS|ISIS|SPAM|FRIES|GUNK|RINGS)
    4 KERNEL.SIZE     S32     30              # Kernel half-size (pixels)
     4KERNEL.SIZE     S32     35              # Kernel half-size (pixels)
    55SPATIAL.ORDER   S32     1               # Spatial polynomial order
    66REGION.SIZE     F32     0               # Iso-kernel region size (pixels)
    77SOURCE.RADIUS   F32     3.0             # Source matching radius (pixels)
    88STAMP.SPACING   F32     300             # Typical spacing between stamps (pixels)
    9 STAMP.FOOTPRINT S32     35              # Size of stamps (pixels)
     9STAMP.FOOTPRINT S32     40              # Size of stamps (pixels)
    1010STAMP.THRESHOLD F32     5               # Flux threshold for stamps (stdev above background)
    1111STRIDE          S32     100             # Size of convolution patches (pixels)
    12 ITER            S32     10              # Number of rejection iterations
     12ITER            S32     2               # Number of rejection iterations
    1313REJ             F32     2.5             # Rejection level (std dev)
    1414KERNEL.ERR      F32     0.05            # Relative systematic error in kernel
    1515SYS.ERR         F32     0.1             # Relative systematic error in images
    1616SKY.ERR         F32     0.0             # Relative systematic error in images
     17NORM.FRAC       F32     0.03            # Fraction of window for normalisation window
    1718
    1819MASK.VAL        STR     MASK.VALUE,CONV.BAD     # Mask value for input
    … …  
    7374STACK   METADATA
    7475        KERNEL.TYPE     STR     ISIS    # Kernel type
     76        SCALE.REF       F32     9.0     # Reference for kernel parameter scaling
    7577END
    7678

  • branches/eam_branches/20091201/ppStack/src/ppStackMatch.c

    r26697 r26703  
    298298            float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
    299299            float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel
     300            float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw
    300301            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images
    301302            float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky
    … …  
    444445                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
    445446                                        orders, inner, ringsOrder, binning, penalty,
    446                                         optimum, optWidths, optOrder, optThresh, iter, rej, sysError, skyErr,
    447                                         kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac,
    448                                         PM_SUBTRACTION_MODE_2)) {
     447                                        optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
     448                                        sysError, skyErr, kernelError, maskVal, maskBad, maskPoor,
     449                                        poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
    449450                    psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
    450451                    psFree(fake);

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

    r26667 r26703  
    2222#include "ppSub.h"
    2323
     24#define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
     25
    2426// Normalise a region on an image
    2527static void normaliseRegion(psImage *image, // Image to normalise
    … …  
    3840}
    3941
     42#if 0
    4043// Measure the PSF for an image
    4144static float subImagePSF(ppSubData *data, // Processing data
    … …  
    8588    return fwhm;
    8689}
     90#endif
    8791
    8892// Scale the kernel parameters according to the PSFs
    … …  
    114118    }
    115119
     120#if 0
    116121    // Input images
    117122    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
    … …  
    133138    float inFWHM = subImagePSF(data, inRO, inSources); // FWHM for input
    134139    float refFWHM = subImagePSF(data, refRO, refSources); // FWHM for reference
     140#else
     141    float inFWHM = 7.631371;
     142    float refFWHM = 10.005879;
     143#endif
     144
    135145    psLogMsg("ppSub", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM);
    136146    if (!isfinite(inFWHM) || !isfinite(refFWHM)) {
    … …  
    252262    float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
    253263    float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel
     264    float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw
    254265    float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images
    255266    float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky
    … …  
    306317    }
    307318
     319    if (inRO->covariance) {
     320        psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC);
     321    }
     322    if (refRO->covariance) {
     323        psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC);
     324    }
     325
    308326    // Match the PSFs
    309327    bool success = false;               // Operation was successful?
    … …  
    315333                                     spacing, threshold, sources, data->stamps, type, size, order,
    316334                                     widths, orders, inner, ringsOrder, binning, penalty, optimum,
    317                                      optWidths, optOrder, optThresh, iter, rej, sysErr, skyErr, kernelErr, maskVal,
    318                                      maskBad, maskPoor, poorFrac, badFrac, subMode);
     335                                     optWidths, optOrder, optThresh, iter, rej, normFrac,
     336                                     sysErr, skyErr, kernelErr, maskVal, maskBad, maskPoor,
     337                                     poorFrac, badFrac, subMode);
    319338    }
    320339
    … …  
    407426    pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true);
    408427
     428    if (inConv->covariance) {
     429        psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC);
     430    }
     431    if (refConv->covariance) {
     432        psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC);
     433    }
     434
    409435    if (inConv->variance) {
    410436        psImageCovarianceTransfer(inConv->variance, inConv->covariance);

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

    r26702 r26703  
    222222static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated
    223223                                      psVector *vector, // Least-squares vector, updated
     224                                      double *norm,     // Normalisation, updated
    224225                                      const psKernel *image1, // Image 1
    225226                                      const psKernel *image2, // Image 2
    … …  
    231232                                      const psImage *polyValues, // Spatial polynomial values
    232233                                      int footprint, // (Half-)Size of stamp
    233                                       const pmSubtractionEquationCalculationMode mode
     234                                      int normWindow, // Window (half-)size for normalisation measurement
     235                                      const pmSubtractionEquationCalculationMode mode
    234236                                      )
    235237{
    … …  
    315317            // Spatial variation of kernel coeffs
    316318            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    317                 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    318                     for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    319                         double aa = sumAA * poly2[iTerm][jTerm];
    320                         double bb = sumBB * poly2[iTerm][jTerm];
    321                         double ab = sumAB * poly2[iTerm][jTerm];
    322 
    323                         matrix->data.F64[iIndex][jIndex] = aa;
    324                         matrix->data.F64[jIndex][iIndex] = aa;
    325 
    326                         matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
    327                         matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
    328 
    329                         matrix->data.F64[iIndex][jIndex + numParams] = ab;
    330                         matrix->data.F64[jIndex + numParams][iIndex] = ab;
    331                     }
    332                 }
    333             }
     319                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     320                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     321                        double aa = sumAA * poly2[iTerm][jTerm];
     322                        double bb = sumBB * poly2[iTerm][jTerm];
     323                        double ab = sumAB * poly2[iTerm][jTerm];
     324
     325                        matrix->data.F64[iIndex][jIndex] = aa;
     326                        matrix->data.F64[jIndex][iIndex] = aa;
     327
     328                        matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
     329                        matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
     330
     331                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
     332                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
     333                    }
     334                }
     335            }
    334336        }
    335337        for (int j = 0; j < i; j++) {
    … …  
    351353            // Spatial variation of kernel coeffs
    352354            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    353                 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    354                     for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    355                         double ab = sumAB * poly2[iTerm][jTerm];
    356                         matrix->data.F64[iIndex][jIndex + numParams] = ab;
    357                         matrix->data.F64[jIndex + numParams][iIndex] = ab;
    358                     }
    359                 }
    360             }
     355                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     356                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     357                        double ab = sumAB * poly2[iTerm][jTerm];
     358                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
     359                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
     360                    }
     361                }
     362            }
    361363        }
    362364
    … …  
    419421
    420422            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    421                 matrix->data.F64[iIndex][normIndex] = ai1;
    422                 matrix->data.F64[normIndex][iIndex] = ai1;
    423                 matrix->data.F64[iIndex + numParams][normIndex] = bi1;
    424                 matrix->data.F64[normIndex][iIndex + numParams] = bi1;
    425             }
     423                matrix->data.F64[iIndex][normIndex] = ai1;
     424                matrix->data.F64[normIndex][iIndex] = ai1;
     425                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
     426                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
     427            }
    426428            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    427                 matrix->data.F64[iIndex][bgIndex] = a;
    428                 matrix->data.F64[bgIndex][iIndex] = a;
    429                 matrix->data.F64[iIndex + numParams][bgIndex] = b;
    430                 matrix->data.F64[bgIndex][iIndex + numParams] = b;
    431             }
     429                matrix->data.F64[iIndex][bgIndex] = a;
     430                matrix->data.F64[bgIndex][iIndex] = a;
     431                matrix->data.F64[iIndex + numParams][bgIndex] = b;
     432                matrix->data.F64[bgIndex][iIndex + numParams] = b;
     433            }
    432434            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    433                 vector->data.F64[iIndex] = ai2;
    434                 vector->data.F64[iIndex + numParams] = bi2;
     435                vector->data.F64[iIndex] = ai2;
     436                vector->data.F64[iIndex + numParams] = bi2;
    435437                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
    436438                    // subtract norm * sumRC * poly[iTerm]
    … …  
    441443                    vector->data.F64[iIndex + numParams] -= norm * bi1;
    442444                }
    443             }
     445            }
    444446        }
    445447    }
    … …  
    450452    double sumI2 = 0.0;                 // Sum of I_2 (for vector, background)
    451453    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
     454    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
    452455    for (int y = - footprint; y <= footprint; y++) {
    453456        for (int x = - footprint; x <= footprint; x++) {
    … …  
    458461            double one = 1.0;
    459462            double i1i2 = i1 * i2;
     463
     464            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     465                normI1 += i1;
     466                normI2 += i2;
     467            }
    460468
    461469            if (weight) {
    … …  
    482490        }
    483491    }
     492
     493    *norm = normI2 / normI1;
     494
     495    fprintf(stderr, "Sums: %f %f %f\n", normI1, normI2, *norm);
     496
    484497    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    485         matrix->data.F64[normIndex][normIndex] = sumI1I1;
    486         vector->data.F64[normIndex] = sumI1I2;
     498        matrix->data.F64[normIndex][normIndex] = sumI1I1;
     499        vector->data.F64[normIndex] = sumI1I2;
    487500    }
    488501    if (mode & PM_SUBTRACTION_EQUATION_BG) {
    489         matrix->data.F64[bgIndex][bgIndex] = sum1;
    490         vector->data.F64[bgIndex] = sumI2;
     502        matrix->data.F64[bgIndex][bgIndex] = sum1;
     503        vector->data.F64[bgIndex] = sumI2;
    491504    }
    492505    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
    493         matrix->data.F64[bgIndex][normIndex] = sumI1;
    494         matrix->data.F64[normIndex][bgIndex] = sumI1;
     506        matrix->data.F64[bgIndex][normIndex] = sumI1;
     507        matrix->data.F64[normIndex][bgIndex] = sumI1;
    495508    }
    496509    return true;
    … …  
    512525    int spatialOrder = kernels->spatialOrder; // Order of spatial variations
    513526    int numKernels = kernels->num; // Number of kernel components
    514     for (int i = 0; i < numKernels; i++) {
     527    int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations
     528    int numParams = numKernels * numSpatial;                 // Number of kernel parameters
     529    for (int i = 0; i < numParams; i++) {
    515530        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
    516531            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    … …  
    518533                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
    519534                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
     535                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     536                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
     537                }
    520538            }
    521539        }
    … …  
    708726        break;
    709727      case PM_SUBTRACTION_MODE_DUAL:
    710         status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2,
     728        status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm,
     729                                           stamp->image1, stamp->image2,
    711730                                           weight, window, stamp->convolutions1, stamp->convolutions2,
    712                                            kernels, polyValues, footprint, mode);
     731                                           kernels, polyValues, footprint, stamps->normWindow, mode);
    713732        break;
    714733      default:
    … …  
    751770}
    752771
    753 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode)
     772bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     773                                    const pmSubtractionEquationCalculationMode mode)
    754774{
    755775    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    12191239        psVectorInit(sumVector, 0.0);
    12201240
     1241        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
     1242
    12211243        int numStamps = 0;              // Number of good stamps
    12221244        for (int i = 0; i < stamps->num; i++) {
    … …  
    12261248                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
    12271249
     1250                psVectorAppend(norms, stamp->norm);
     1251
    12281252                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    12291253                numStamps++;
    … …  
    12461270
    12471271#if 0
     1272        // Regular, straight-forward solution
     1273        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1274#else
    12481275        {
    1249             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    1250             if (!solution) {
    1251                 psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
     1276            // Solve normalisation and background separately
     1277            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1278
     1279#if 0
     1280            psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64);
     1281            psVector *normVector = psVectorAlloc(2, PS_TYPE_F64);
     1282
     1283            normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex];
     1284            normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex];
     1285            normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex];
     1286
     1287            normVector->data.F64[0] = sumVector->data.F64[normIndex];
     1288            normVector->data.F64[1] = sumVector->data.F64[bgIndex];
     1289
     1290            psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN);
     1291
     1292            double normValue = normSolution->data.F64[0];
     1293            double bgValue = normSolution->data.F64[1];
     1294
     1295            psFree(normMatrix);
     1296            psFree(normVector);
     1297            psFree(normSolution);
     1298#endif
     1299
     1300            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
     1301            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
     1302                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
     1303                psFree(stats);
    12521304                psFree(sumMatrix);
    12531305                psFree(sumVector);
    1254                 return NULL;
    1255             }
    1256 
    1257 #ifdef TESTING
    1258             for (int i = 0; i < numParams; i++) {
    1259                 fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
    1260             }
    1261 #endif
    1262 
    1263             int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
    1264             int numKernels = kernels->num; // Number of kernel basis functions
    1265 
    1266 #define MASK_BASIS(INDEX)                                               \
    1267             {                                                           \
    1268                 for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
    1269                     for (int k = 0; k < numParams; k++) {               \
    1270                         sumMatrix->data.F64[index][k] = sumMatrix->data.F64[k][index] = 0.0; \
    1271                     }                                                   \
    1272                     sumMatrix->data.F64[index][index] = 1.0;            \
    1273                     sumVector->data.F64[index] = 0.0;                   \
    1274                 }                                                       \
    1275             }
    1276 
    1277             #define TOL 1.0e-3
    1278             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1279             double norm = fabs(solution->data.F64[normIndex]);  // Normalisation
    1280             double thresh = norm * TOL;                         // Threshold for low parameters
    1281             for (int i = 0; i < numKernels; i++) {
    1282                 // Getting 0th order parameter value.  In the presence of spatial variation, the actual value
    1283                 // of the parameter will vary over the image.  We are in effect getting the value in the
    1284                 // centre of the image.  If we use different polynomial functions (e.g., Chebyshev), we may
    1285                 // have to change this to properly determine the value of the parameter at the centre.
    1286                 double param1 = solution->data.F64[i],
    1287                     param2 = solution->data.F64[numSolution1 + i]; // Parameters of interest
    1288                 bool mask1 = false, mask2 = false;              // Masked the parameter?
    1289                 if (fabs(param1) < thresh) {
    1290                     psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i);
    1291                     MASK_BASIS(i);
    1292                     mask1 = true;
    1293                 }
    1294                 if (fabs(param2) < thresh) {
    1295                     psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i);
    1296                     MASK_BASIS(numSolution1 + i);
    1297                     mask2 = true;
    1298                 }
    1299 
    1300                 if (!mask1 && !mask2) {
    1301                     if (fabs(param1) > fabs(param2)) {
    1302                         psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i);
    1303                         MASK_BASIS(numSolution1 + i);
    1304                     } else {
    1305                         psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i);
    1306                         MASK_BASIS(i);
    1307                     }
    1308                 }
    1309             }
    1310 
    1311 #ifdef TESTING
    1312             psFitsWriteImageSimple ("sumMatrixFix.fits", sumMatrix, NULL);
    1313             psVectorWriteFile("sumVectorFix.dat", sumVector);
    1314 #endif
    1315         }
    1316 #endif /*** kernel-masking block ***/
    1317         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1306                psFree(norms);
     1307                return false;
     1308            }
     1309
     1310            double normValue = stats->robustMedian;
     1311//            double bgValue = 0.0;
     1312
     1313            psFree(stats);
     1314
     1315            fprintf(stderr, "Norm: %lf\n", normValue);
     1316//            fprintf(stderr, "BG: %lf\n", bgValue);
     1317
     1318            // Solve kernel components
     1319            for (int i = 0; i < numSolution2; i++) {
     1320                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
     1321                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; // + bgValue * sumMatrix->data.F64[bgIndex][i + numSolution1];
     1322
     1323                sumMatrix->data.F64[i][normIndex] = 0.0;
     1324                sumMatrix->data.F64[normIndex][i] = 0.0;
     1325
     1326//                sumMatrix->data.F64[i][bgIndex] = 0.0;
     1327//                sumMatrix->data.F64[bgIndex][i] = 0.0;
     1328                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
     1329                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
     1330//                sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0;
     1331//                sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0;
     1332            }
     1333
     1334            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
     1335//            sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
     1336//            sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1337
     1338            sumVector->data.F64[normIndex] = 0.0;
     1339//            sumVector->data.F64[bgIndex] = 0.0;
     1340
     1341            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1342
     1343            solution->data.F64[normIndex] = normValue;
     1344//            solution->data.F64[bgIndex] = bgValue;
     1345        }
     1346#endif
     1347
    13181348
    13191349#ifdef TESTING
    … …  
    13251355        psFree(sumMatrix);
    13261356        psFree(sumVector);
     1357
     1358        psFree(norms);
    13271359
    13281360        if (!kernels->solution1) {
    … …  
    13471379            int numKernels = kernels->num;
    13481380            for (int i = 0; i < numKernels * numSpatial; i++) {
    1349                 // XXX fprintf (stderr, "keep\n");
    1350                 kernels->solution1->data.F64[i] = solution->data.F64[i];
    1351                 kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
     1381                // XXX fprintf (stderr, "keep\n");
     1382                kernels->solution1->data.F64[i] = solution->data.F64[i];
     1383                kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
    13521384            }
    13531385        }

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

    r26702 r26703  
    2828static bool useFFT = true;              // Do convolutions using FFT
    2929
    30 # define SEPARATE 1
     30# define SEPARATE 0
    3131# if (SEPARATE)
    3232# define SUBMODE PM_SUBTRACTION_EQUATION_NORM
    … …  
    7474                                 float thresh2,  // Threshold for stamp finding on readout 2
    7575                                 float stampSpacing, // Spacing between stamps
     76                                 float normFrac,     // Fraction of flux in window for normalisation window
    7677                                 float sysError,     // Relative systematic error in images
    7778                                 float skyError,     // Relative systematic error in images
    … …  
    8687
    8788    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
    88                                       size, footprint, stampSpacing, sysError, skyError, mode);
     89                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
    8990    if (!*stamps) {
    9091        psError(psErrorCodeLast(), false, "Unable to find stamps.");
    … …  
    109110                                  const pmReadout *ro1, const pmReadout *ro2, // Input images
    110111                                  int stride, // Size for convolution patches
     112                                  float normFrac,           // Fraction of window for normalisation window
    111113                                  float sysError,           // Systematic error in images
    112114                                  float skyError,           // Systematic error in images
    … …  
    157159    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
    158160    PS_ASSERT_INT_NONNEGATIVE(stride, false);
     161    if (isfinite(normFrac)) {
     162        PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false);
     163        PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false);
     164    }
    159165    if (isfinite(sysError)) {
    160166        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
    … …  
    281287    }
    282288
    283     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, kernelError,
     289    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError,
    284290                               maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
    285291        psFree(kernels);
    … …  
    347353                        int inner, int ringsOrder, int binning, float penalty,
    348354                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    349                         int iter, float rej, float sysError, float skyError, float kernelError, psImageMaskType maskVal,
    350                         psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac,
    351                         float badFrac, pmSubtractionMode subMode)
    352 {
    353     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, skyError, kernelError,
     355                        int iter, float rej, float normFrac, float sysError, float skyError,
     356                        float kernelError, psImageMaskType maskVal, psImageMaskType maskBad,
     357                        psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
     358{
     359    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError,
    354360                               maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
    355361        return false;
    … …  
    443449    }
    444450
     451    // General background subtraction and measurement of stamp threshold
    445452    float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
    446453    {
    447454        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
    448455        if (ro1) {
     456            psStatsInit(bg);
    449457            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
    450458                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
    … …  
    452460                goto MATCH_ERROR;
    453461            }
    454             stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
     462            stampThresh1 = threshold * bg->robustStdev;
     463            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
    455464        }
    456465        if (ro2) {
     466            psStatsInit(bg);
    457467            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
    458468                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
    … …  
    460470                goto MATCH_ERROR;
    461471            }
    462             stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
    463         }
     472            stampThresh2 = threshold * bg->robustStdev;
     473            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
     474       }
    464475        psFree(bg);
    465476    }
    … …  
    481492            if (stampsName && strlen(stampsName) > 0) {
    482493                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
    483                                                         footprint, stampSpacing, sysError, skyError, subMode);
     494                                                        footprint, stampSpacing, normFrac,
     495                                                        sysError, skyError, subMode);
    484496            } else if (sources) {
    485497                stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
    486                                                            footprint, stampSpacing, sysError, skyError, subMode);
     498                                                           footprint, stampSpacing, normFrac,
     499                                                           sysError, skyError, subMode);
    487500            }
    488501
    … …  
    490503            // doesn't matter.
    491504            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
    492                                       stampSpacing, sysError, skyError, size, footprint, subMode)) {
     505                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
    493506                goto MATCH_ERROR;
    494507            }
    … …  
    562575                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    563576
     577#if 0
    564578                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    565                                           stampThresh1, stampThresh2, stampSpacing, sysError, skyError,
    566                                           size, footprint, subMode)) {
    567                     goto MATCH_ERROR;
    568                 }
     579                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
     580                                          sysError, skyError, size, footprint, subMode)) {
     581                    goto MATCH_ERROR;
     582                }
     583#endif
    569584
    570585                // generate the window function from the set of stamps

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

    r26667 r26703  
    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/20091201/psModules/src/imcombine/pmSubtractionStamps.c

    r26562 r26703  
    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;
    … …  
    308311    stamp->matrix = NULL;
    309312    stamp->vector = NULL;
     313    stamp->norm = NAN;
    310314
    311315    return stamp;
    … …  
    316320                                                const psImage *image2, const psImage *subMask,
    317321                                                const psRegion *region, float thresh1, float thresh2,
    318                                                 int size, int footprint, float spacing, float sysErr, float skyErr,
    319                                                 pmSubtractionMode mode)
     322                                                int size, int footprint, float spacing, float normFrac,
     323                                                float sysErr, float skyErr, pmSubtractionMode mode)
    320324{
    321325    if (!image1 && !image2) {
    … …  
    373377
    374378    if (!stamps) {
    375         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr);
     379        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
     380                                             normFrac, sysErr, skyErr);
    376381    }
    377382
    … …  
    407412                // Take stamps off the top of the (sorted) list
    408413                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 
    411414                    // Chop off the top of the list
    412415                    xList->n = j;
    … …  
    414417                    fluxList->n = j;
    415418
     419#if 0
    416420                    // Fish around a bit to see if we can find a pixel that isn't masked
     421                    // This is not a good idea if we're using the window feature
    417422                    psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n",
    418423                            i, xCentre, yCentre);
    419424
    420425                    // Search bounds
     426                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
    421427                    int search = footprint - size; // Search radius
    422428                    int xMin = PS_MAX(border, xCentre - search);
    … …  
    427433                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
    428434                                            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];
     435#else
     436                    // Only search the exact centre pixel
     437                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
     438                                            subMask, xList->data.F32[j], xList->data.F32[j],
     439                                            yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
     440#endif
    433441                    // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
    434442                }
    … …  
    486494                                               const psImage *image, const psImage *subMask,
    487495                                               const psRegion *region, int size, int footprint,
    488                                                float spacing, float sysErr, float skyErr, pmSubtractionMode mode)
     496                                               float spacing, float normFrac, float sysErr, float skyErr,
     497                                               pmSubtractionMode mode)
    489498
    490499{
    … …  
    507516    pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
    508517                                                                 region, footprint, spacing,
    509                                                                  sysErr, skyErr); // Stamp list
     518                                                                 normFrac, sysErr, skyErr); // Stamp list
    510519    int numStamps = stamps->num;        // Number of stamps
    511520
    … …  
    612621    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    613622
    614     int size = kernelSize + stamps->footprint; // Size of postage stamps
     623    int size = stamps->footprint; // Size of postage stamps
    615624
    616625    psFree (stamps->window);
    … …  
    641650    // storage vector for flux data
    642651    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;
     652    psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
     653    psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
    646654
    647655    // generate the window pixels
     656    double sum = 0.0;                   // Sum inside the window
     657    float maxValue = 0.0;               // Maximum value, for normalisation
    648658    for (int y = -size; y <= size; y++) {
    649659        for (int x = -size; x <= size; x++) {
    650660
    651             flux->n = 0;
     661            flux1->n = 0;
     662            flux2->n = 0;
    652663            for (int i = 0; i < stamps->num; i++) {
    653664                pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    … …  
    656667                if (!stamp->image2) continue;
    657668
    658                 psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
    659                 psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
     669                psVectorAppend(flux1, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
     670                psVectorAppend(flux2, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
    660671            }
    661672
    662673            psStatsInit (stats);
    663             if (!psVectorStats (stats, flux, NULL, NULL, 0)) {
     674            if (!psVectorStats (stats, flux1, NULL, NULL, 0)) {
    664675                psAbort ("failed to generate stats");
    665676            }
    666             stamps->window->kernel[y][x] = stats->robustMedian;
    667             if (maxValue < stats->robustMedian) {
    668                 maxValue = stats->robustMedian;
    669             }
    670         }
     677            float f1 = stats->robustMedian;
     678            psStatsInit (stats);
     679            if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
     680                psAbort ("failed to generate stats");
     681            }
     682            float f2 = stats->robustMedian;
     683
     684            stamps->window->kernel[y][x] = f1 + f2;
     685            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
     686                sum += stamps->window->kernel[y][x];
     687            }
     688            maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
     689        }
     690    }
     691
     692    psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
     693            sum, (1.0 - stamps->normFrac) * sum);
     694    bool done = false;
     695    for (int radius = 1; radius <= size && !done; radius++) {
     696        double within = 0.0;
     697        for (int y = -radius; y <= radius; y++) {
     698            for (int x = -radius; x <= radius; x++) {
     699                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
     700                    within += stamps->window->kernel[y][x];
     701                }
     702            }
     703        }
     704        psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
     705        if (within > (1.0 - stamps->normFrac) * sum) {
     706            stamps->normWindow = radius;
     707            done = true;
     708        }
     709    }
     710    psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
     711    if (stamps->normWindow == 0 || stamps->normWindow >= size) {
     712        psError(PS_ERR_UNKNOWN, true, "Unable to determine normalisation window size.");
     713        return false;
    671714    }
    672715
    … …  
    678721    }
    679722
    680 # ifdef TESTING
     723#if 1
    681724    {
    682725        psFits *fits = psFitsOpen ("window.fits", "w");
    … …  
    684727        psFitsClose (fits);
    685728    }
    686 # endif
     729#endif
    687730
    688731    psFree (stats);
    689     psFree (flux);
     732    psFree (flux1);
     733    psFree(flux2);
    690734    psFree (norm1);
    691735    psFree (norm2);
    … …  
    788832pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
    789833                                                          const psImage *subMask, const psRegion *region,
    790                                                           int size, int footprint, float spacing, float sysErr, float skyErr,
     834                                                          int size, int footprint, float spacing,
     835                                                          float normFrac, float sysErr, float skyErr,
    791836                                                          pmSubtractionMode mode)
    792837{
    … …  
    819864
    820865    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
    821                                                             footprint, spacing, sysErr, skyErr, mode); // Stamps
     866                                                            footprint, spacing, normFrac,
     867                                                            sysErr, skyErr, mode); // Stamps
    822868    psFree(x);
    823869    psFree(y);
    … …  
    833879pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
    834880                                                       const psImage *subMask, const psRegion *region,
    835                                                        int size, int footprint, float spacing, float sysErr, float skyErr,
    836                                                        pmSubtractionMode mode)
     881                                                       int size, int footprint, float spacing, float normFrac,
     882                                                       float sysErr, float skyErr, pmSubtractionMode mode)
    837883{
    838884    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
    … …  
    851897
    852898    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
    853                                                             spacing, sysErr, skyErr, mode);
     899                                                            spacing, normFrac, sysErr, skyErr, mode);
    854900    psFree(data);
    855901

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

    r26578 r26703  
    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
Note: See TracChangeset for help on using the changeset viewer.