Index: branches/eam_branches/20091201/ippconfig/recipes/ppSub.config =================================================================== --- branches/eam_branches/20091201/ippconfig/recipes/ppSub.config (revision 26702) +++ branches/eam_branches/20091201/ippconfig/recipes/ppSub.config (revision 26703) @@ -2,17 +2,18 @@ KERNEL.TYPE STR ISIS_RADIAL # Kernel type to use (POIS|ISIS|SPAM|FRIES|GUNK|RINGS) -KERNEL.SIZE S32 30 # Kernel half-size (pixels) +KERNEL.SIZE S32 35 # Kernel half-size (pixels) SPATIAL.ORDER S32 1 # Spatial polynomial order REGION.SIZE F32 0 # Iso-kernel region size (pixels) SOURCE.RADIUS F32 3.0 # Source matching radius (pixels) STAMP.SPACING F32 300 # Typical spacing between stamps (pixels) -STAMP.FOOTPRINT S32 35 # Size of stamps (pixels) +STAMP.FOOTPRINT S32 40 # Size of stamps (pixels) STAMP.THRESHOLD F32 5 # Flux threshold for stamps (stdev above background) STRIDE S32 100 # Size of convolution patches (pixels) -ITER S32 10 # Number of rejection iterations +ITER S32 2 # Number of rejection iterations REJ F32 2.5 # Rejection level (std dev) KERNEL.ERR F32 0.05 # Relative systematic error in kernel SYS.ERR F32 0.1 # Relative systematic error in images SKY.ERR F32 0.0 # Relative systematic error in images +NORM.FRAC F32 0.03 # Fraction of window for normalisation window MASK.VAL STR MASK.VALUE,CONV.BAD # Mask value for input @@ -73,4 +74,5 @@ STACK METADATA KERNEL.TYPE STR ISIS # Kernel type + SCALE.REF F32 9.0 # Reference for kernel parameter scaling END Index: branches/eam_branches/20091201/ppStack/src/ppStackMatch.c =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackMatch.c (revision 26702) +++ branches/eam_branches/20091201/ppStack/src/ppStackMatch.c (revision 26703) @@ -298,4 +298,5 @@ float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky @@ -444,7 +445,7 @@ threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, - optimum, optWidths, optOrder, optThresh, iter, rej, sysError, skyErr, - kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, - PM_SUBTRACTION_MODE_2)) { + optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, + sysError, skyErr, kernelError, maskVal, maskBad, maskPoor, + poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { psError(PS_ERR_UNKNOWN, false, "Unable to match images."); psFree(fake); Index: branches/eam_branches/20091201/ppSub/src/ppSubMatchPSFs.c =================================================================== --- branches/eam_branches/20091201/ppSub/src/ppSubMatchPSFs.c (revision 26702) +++ branches/eam_branches/20091201/ppSub/src/ppSubMatchPSFs.c (revision 26703) @@ -22,4 +22,6 @@ #include "ppSub.h" +#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix + // Normalise a region on an image static void normaliseRegion(psImage *image, // Image to normalise @@ -38,4 +40,5 @@ } +#if 0 // Measure the PSF for an image static float subImagePSF(ppSubData *data, // Processing data @@ -85,4 +88,5 @@ return fwhm; } +#endif // Scale the kernel parameters according to the PSFs @@ -114,4 +118,5 @@ } +#if 0 // Input images pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout @@ -133,4 +138,9 @@ float inFWHM = subImagePSF(data, inRO, inSources); // FWHM for input float refFWHM = subImagePSF(data, refRO, refSources); // FWHM for reference +#else + float inFWHM = 7.631371; + float refFWHM = 10.005879; +#endif + psLogMsg("ppSub", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM); if (!isfinite(inFWHM) || !isfinite(refFWHM)) { @@ -252,4 +262,5 @@ float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky @@ -306,4 +317,11 @@ } + if (inRO->covariance) { + psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC); + } + if (refRO->covariance) { + psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC); + } + // Match the PSFs bool success = false; // Operation was successful? @@ -315,6 +333,7 @@ spacing, threshold, sources, data->stamps, type, size, order, widths, orders, inner, ringsOrder, binning, penalty, optimum, - optWidths, optOrder, optThresh, iter, rej, sysErr, skyErr, kernelErr, maskVal, - maskBad, maskPoor, poorFrac, badFrac, subMode); + optWidths, optOrder, optThresh, iter, rej, normFrac, + sysErr, skyErr, kernelErr, maskVal, maskBad, maskPoor, + poorFrac, badFrac, subMode); } @@ -407,4 +426,11 @@ pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true); + if (inConv->covariance) { + psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC); + } + if (refConv->covariance) { + psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC); + } + if (inConv->variance) { psImageCovarianceTransfer(inConv->variance, inConv->covariance); Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26702) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26703) @@ -222,4 +222,5 @@ static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated psVector *vector, // Least-squares vector, updated + double *norm, // Normalisation, updated const psKernel *image1, // Image 1 const psKernel *image2, // Image 2 @@ -231,5 +232,6 @@ const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp - const pmSubtractionEquationCalculationMode mode + int normWindow, // Window (half-)size for normalisation measurement + const pmSubtractionEquationCalculationMode mode ) { @@ -315,21 +317,21 @@ // Spatial variation of kernel coeffs if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double aa = sumAA * poly2[iTerm][jTerm]; - double bb = sumBB * poly2[iTerm][jTerm]; - double ab = sumAB * poly2[iTerm][jTerm]; - - matrix->data.F64[iIndex][jIndex] = aa; - matrix->data.F64[jIndex][iIndex] = aa; - - matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; - matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; - - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double aa = sumAA * poly2[iTerm][jTerm]; + double bb = sumBB * poly2[iTerm][jTerm]; + double ab = sumAB * poly2[iTerm][jTerm]; + + matrix->data.F64[iIndex][jIndex] = aa; + matrix->data.F64[jIndex][iIndex] = aa; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; + matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; + + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } } for (int j = 0; j < i; j++) { @@ -351,12 +353,12 @@ // Spatial variation of kernel coeffs if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double ab = sumAB * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double ab = sumAB * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } } @@ -419,18 +421,18 @@ if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = ai1; - matrix->data.F64[normIndex][iIndex] = ai1; - matrix->data.F64[iIndex + numParams][normIndex] = bi1; - matrix->data.F64[normIndex][iIndex + numParams] = bi1; - } + matrix->data.F64[iIndex][normIndex] = ai1; + matrix->data.F64[normIndex][iIndex] = ai1; + matrix->data.F64[iIndex + numParams][normIndex] = bi1; + matrix->data.F64[normIndex][iIndex + numParams] = bi1; + } if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = a; - matrix->data.F64[bgIndex][iIndex] = a; - matrix->data.F64[iIndex + numParams][bgIndex] = b; - matrix->data.F64[bgIndex][iIndex + numParams] = b; - } + matrix->data.F64[iIndex][bgIndex] = a; + matrix->data.F64[bgIndex][iIndex] = a; + matrix->data.F64[iIndex + numParams][bgIndex] = b; + matrix->data.F64[bgIndex][iIndex + numParams] = b; + } if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = ai2; - vector->data.F64[iIndex + numParams] = bi2; + vector->data.F64[iIndex] = ai2; + vector->data.F64[iIndex + numParams] = bi2; if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { // subtract norm * sumRC * poly[iTerm] @@ -441,5 +443,5 @@ vector->data.F64[iIndex + numParams] -= norm * bi1; } - } + } } } @@ -450,4 +452,5 @@ double sumI2 = 0.0; // Sum of I_2 (for vector, background) double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) + double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -458,4 +461,9 @@ double one = 1.0; double i1i2 = i1 * i2; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { + normI1 += i1; + normI2 += i2; + } if (weight) { @@ -482,15 +490,20 @@ } } + + *norm = normI2 / normI1; + + fprintf(stderr, "Sums: %f %f %f\n", normI1, normI2, *norm); + if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumI1I1; - vector->data.F64[normIndex] = sumI1I2; + matrix->data.F64[normIndex][normIndex] = sumI1I1; + vector->data.F64[normIndex] = sumI1I2; } if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI2; + matrix->data.F64[bgIndex][bgIndex] = sum1; + vector->data.F64[bgIndex] = sumI2; } if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[bgIndex][normIndex] = sumI1; - matrix->data.F64[normIndex][bgIndex] = sumI1; + matrix->data.F64[bgIndex][normIndex] = sumI1; + matrix->data.F64[normIndex][bgIndex] = sumI1; } return true; @@ -512,5 +525,7 @@ int spatialOrder = kernels->spatialOrder; // Order of spatial variations int numKernels = kernels->num; // Number of kernel components - for (int i = 0; i < numKernels; i++) { + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations + int numParams = numKernels * numSpatial; // Number of kernel parameters + for (int i = 0; i < numParams; i++) { for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { @@ -518,4 +533,7 @@ psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty"); matrix->data.F64[index][index] += norm * penalties->data.F32[i]; + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i]; + } } } @@ -708,7 +726,8 @@ break; case PM_SUBTRACTION_MODE_DUAL: - status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2, + status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, + stamp->image1, stamp->image2, weight, window, stamp->convolutions1, stamp->convolutions2, - kernels, polyValues, footprint, mode); + kernels, polyValues, footprint, stamps->normWindow, mode); break; default: @@ -751,5 +770,6 @@ } -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode) +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, + const pmSubtractionEquationCalculationMode mode) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -1219,4 +1239,6 @@ psVectorInit(sumVector, 0.0); + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + int numStamps = 0; // Number of good stamps for (int i = 0; i < stamps->num; i++) { @@ -1226,4 +1248,6 @@ (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + psVectorAppend(norms, stamp->norm); + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; @@ -1246,74 +1270,80 @@ #if 0 + // Regular, straight-forward solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); +#else { - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - if (!solution) { - psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n"); + // Solve normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + +#if 0 + psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); + psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); + + normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; + normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; + normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; + + normVector->data.F64[0] = sumVector->data.F64[normIndex]; + normVector->data.F64[1] = sumVector->data.F64[bgIndex]; + + psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); + + double normValue = normSolution->data.F64[0]; + double bgValue = normSolution->data.F64[1]; + + psFree(normMatrix); + psFree(normVector); + psFree(normSolution); +#endif + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation"); + psFree(stats); psFree(sumMatrix); psFree(sumVector); - return NULL; - } - -#ifdef TESTING - for (int i = 0; i < numParams; i++) { - fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]); - } -#endif - - int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations - int numKernels = kernels->num; // Number of kernel basis functions - -#define MASK_BASIS(INDEX) \ - { \ - for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \ - for (int k = 0; k < numParams; k++) { \ - sumMatrix->data.F64[index][k] = sumMatrix->data.F64[k][index] = 0.0; \ - } \ - sumMatrix->data.F64[index][index] = 1.0; \ - sumVector->data.F64[index] = 0.0; \ - } \ - } - - #define TOL 1.0e-3 - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(solution->data.F64[normIndex]); // Normalisation - double thresh = norm * TOL; // Threshold for low parameters - for (int i = 0; i < numKernels; i++) { - // Getting 0th order parameter value. In the presence of spatial variation, the actual value - // of the parameter will vary over the image. We are in effect getting the value in the - // centre of the image. If we use different polynomial functions (e.g., Chebyshev), we may - // have to change this to properly determine the value of the parameter at the centre. - double param1 = solution->data.F64[i], - param2 = solution->data.F64[numSolution1 + i]; // Parameters of interest - bool mask1 = false, mask2 = false; // Masked the parameter? - if (fabs(param1) < thresh) { - psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i); - MASK_BASIS(i); - mask1 = true; - } - if (fabs(param2) < thresh) { - psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i); - MASK_BASIS(numSolution1 + i); - mask2 = true; - } - - if (!mask1 && !mask2) { - if (fabs(param1) > fabs(param2)) { - psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i); - MASK_BASIS(numSolution1 + i); - } else { - psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i); - MASK_BASIS(i); - } - } - } - -#ifdef TESTING - psFitsWriteImageSimple ("sumMatrixFix.fits", sumMatrix, NULL); - psVectorWriteFile("sumVectorFix.dat", sumVector); -#endif - } -#endif /*** kernel-masking block ***/ - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + psFree(norms); + return false; + } + + double normValue = stats->robustMedian; +// double bgValue = 0.0; + + psFree(stats); + + fprintf(stderr, "Norm: %lf\n", normValue); +// fprintf(stderr, "BG: %lf\n", bgValue); + + // Solve kernel components + for (int i = 0; i < numSolution2; i++) { + sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i]; + sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; // + bgValue * sumMatrix->data.F64[bgIndex][i + numSolution1]; + + sumMatrix->data.F64[i][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i] = 0.0; + +// sumMatrix->data.F64[i][bgIndex] = 0.0; +// sumMatrix->data.F64[bgIndex][i] = 0.0; + sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; +// sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0; +// sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0; + } + + sumMatrix->data.F64[normIndex][normIndex] = 1.0; +// sumMatrix->data.F64[bgIndex][bgIndex] = 1.0; +// sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0; + + sumVector->data.F64[normIndex] = 0.0; +// sumVector->data.F64[bgIndex] = 0.0; + + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + + solution->data.F64[normIndex] = normValue; +// solution->data.F64[bgIndex] = bgValue; + } +#endif + #ifdef TESTING @@ -1325,4 +1355,6 @@ psFree(sumMatrix); psFree(sumVector); + + psFree(norms); if (!kernels->solution1) { @@ -1347,7 +1379,7 @@ int numKernels = kernels->num; for (int i = 0; i < numKernels * numSpatial; i++) { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; } } Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26702) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26703) @@ -28,5 +28,5 @@ static bool useFFT = true; // Do convolutions using FFT -# define SEPARATE 1 +# define SEPARATE 0 # if (SEPARATE) # define SUBMODE PM_SUBTRACTION_EQUATION_NORM @@ -74,4 +74,5 @@ float thresh2, // Threshold for stamp finding on readout 2 float stampSpacing, // Spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window float sysError, // Relative systematic error in images float skyError, // Relative systematic error in images @@ -86,5 +87,5 @@ *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, - size, footprint, stampSpacing, sysError, skyError, mode); + size, footprint, stampSpacing, normFrac, sysError, skyError, mode); if (!*stamps) { psError(psErrorCodeLast(), false, "Unable to find stamps."); @@ -109,4 +110,5 @@ const pmReadout *ro1, const pmReadout *ro2, // Input images int stride, // Size for convolution patches + float normFrac, // Fraction of window for normalisation window float sysError, // Systematic error in images float skyError, // Systematic error in images @@ -157,4 +159,8 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false); PS_ASSERT_INT_NONNEGATIVE(stride, false); + if (isfinite(normFrac)) { + PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false); + } if (isfinite(sysError)) { PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false); @@ -281,5 +287,5 @@ } - if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, kernelError, + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) { psFree(kernels); @@ -347,9 +353,9 @@ int inner, int ringsOrder, int binning, float penalty, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, - int iter, float rej, float sysError, float skyError, float kernelError, psImageMaskType maskVal, - psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac, - float badFrac, pmSubtractionMode subMode) -{ - if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, skyError, kernelError, + int iter, float rej, float normFrac, float sysError, float skyError, + float kernelError, psImageMaskType maskVal, psImageMaskType maskBad, + psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode) +{ + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) { return false; @@ -443,8 +449,10 @@ } + // General background subtraction and measurement of stamp threshold float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images { psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun if (ro1) { + psStatsInit(bg); if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics."); @@ -452,7 +460,9 @@ goto MATCH_ERROR; } - stampThresh1 = bg->robustMedian + threshold * bg->robustStdev; + stampThresh1 = threshold * bg->robustStdev; + psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); } if (ro2) { + psStatsInit(bg); if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics."); @@ -460,6 +470,7 @@ goto MATCH_ERROR; } - stampThresh2 = bg->robustMedian + threshold * bg->robustStdev; - } + stampThresh2 = threshold * bg->robustStdev; + psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); + } psFree(bg); } @@ -481,8 +492,10 @@ if (stampsName && strlen(stampsName) > 0) { stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size, - footprint, stampSpacing, sysError, skyError, subMode); + footprint, stampSpacing, normFrac, + sysError, skyError, subMode); } else if (sources) { stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size, - footprint, stampSpacing, sysError, skyError, subMode); + footprint, stampSpacing, normFrac, + sysError, skyError, subMode); } @@ -490,5 +503,5 @@ // doesn't matter. if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2, - stampSpacing, sysError, skyError, size, footprint, subMode)) { + stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) { goto MATCH_ERROR; } @@ -562,9 +575,11 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); +#if 0 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, - stampThresh1, stampThresh2, stampSpacing, sysError, skyError, - size, footprint, subMode)) { - goto MATCH_ERROR; - } + stampThresh1, stampThresh2, stampSpacing, normFrac, + sysError, skyError, size, footprint, subMode)) { + goto MATCH_ERROR; + } +#endif // generate the window function from the set of stamps Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h (revision 26702) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h (revision 26703) @@ -39,4 +39,5 @@ int iter, ///< Rejection iterations float rej, ///< Rejection threshold + float normFrac, ///< Fraction of flux in window for normalisation window float sysError, ///< Relative systematic error in images float skyError, ///< Relative systematic error in images Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26702) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26703) @@ -176,5 +176,6 @@ pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, - int footprint, float spacing, float sysErr, float skyErr) + int footprint, float spacing, float normFrac, + float sysErr, float skyErr) { pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return @@ -217,4 +218,6 @@ list->flux = NULL; list->window = NULL; + list->normFrac = normFrac; + list->normWindow = 0; list->footprint = footprint; list->sysErr = sysErr; @@ -308,4 +311,5 @@ stamp->matrix = NULL; stamp->vector = NULL; + stamp->norm = NAN; return stamp; @@ -316,6 +320,6 @@ const psImage *image2, const psImage *subMask, const psRegion *region, float thresh1, float thresh2, - int size, int footprint, float spacing, float sysErr, float skyErr, - pmSubtractionMode mode) + int size, int footprint, float spacing, float normFrac, + float sysErr, float skyErr, pmSubtractionMode mode) { if (!image1 && !image2) { @@ -373,5 +377,6 @@ if (!stamps) { - stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr); + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, + normFrac, sysErr, skyErr); } @@ -407,6 +412,4 @@ // Take stamps off the top of the (sorted) list for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) { - int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre - // Chop off the top of the list xList->n = j; @@ -414,9 +417,12 @@ fluxList->n = j; +#if 0 // Fish around a bit to see if we can find a pixel that isn't masked + // This is not a good idea if we're using the window feature psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n", i, xCentre, yCentre); // Search bounds + int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre int search = footprint - size; // Search radius int xMin = PS_MAX(border, xCentre - search); @@ -427,8 +433,10 @@ goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, subMask, xMin, xMax, yMin, yMax, numCols, numRows, border); - - // XXX reset for a test: - xStamp = xList->data.F32[j]; - yStamp = yList->data.F32[j]; +#else + // Only search the exact centre pixel + goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, + subMask, xList->data.F32[j], xList->data.F32[j], + yList->data.F32[j], yList->data.F32[j], numCols, numRows, border); +#endif // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp); } @@ -486,5 +494,6 @@ const psImage *image, const psImage *subMask, const psRegion *region, int size, int footprint, - float spacing, float sysErr, float skyErr, pmSubtractionMode mode) + float spacing, float normFrac, float sysErr, float skyErr, + pmSubtractionMode mode) { @@ -507,5 +516,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, region, footprint, spacing, - sysErr, skyErr); // Stamp list + normFrac, sysErr, skyErr); // Stamp list int numStamps = stamps->num; // Number of stamps @@ -612,5 +621,5 @@ PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - int size = kernelSize + stamps->footprint; // Size of postage stamps + int size = stamps->footprint; // Size of postage stamps psFree (stamps->window); @@ -641,13 +650,15 @@ // storage vector for flux data psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); - - double maxValue = 0.0; + psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); + psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32); // generate the window pixels + double sum = 0.0; // Sum inside the window + float maxValue = 0.0; // Maximum value, for normalisation for (int y = -size; y <= size; y++) { for (int x = -size; x <= size; x++) { - flux->n = 0; + flux1->n = 0; + flux2->n = 0; for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -656,17 +667,49 @@ if (!stamp->image2) continue; - psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]); - psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]); + psVectorAppend(flux1, stamp->image1->kernel[y][x] / norm1->data.F32[i]); + psVectorAppend(flux2, stamp->image2->kernel[y][x] / norm2->data.F32[i]); } psStatsInit (stats); - if (!psVectorStats (stats, flux, NULL, NULL, 0)) { + if (!psVectorStats (stats, flux1, NULL, NULL, 0)) { psAbort ("failed to generate stats"); } - stamps->window->kernel[y][x] = stats->robustMedian; - if (maxValue < stats->robustMedian) { - maxValue = stats->robustMedian; - } - } + float f1 = stats->robustMedian; + psStatsInit (stats); + if (!psVectorStats (stats, flux2, NULL, NULL, 0)) { + psAbort ("failed to generate stats"); + } + float f2 = stats->robustMedian; + + stamps->window->kernel[y][x] = f1 + f2; + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) { + sum += stamps->window->kernel[y][x]; + } + maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]); + } + } + + psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n", + sum, (1.0 - stamps->normFrac) * sum); + bool done = false; + for (int radius = 1; radius <= size && !done; radius++) { + double within = 0.0; + for (int y = -radius; y <= radius; y++) { + for (int x = -radius; x <= radius; x++) { + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { + within += stamps->window->kernel[y][x]; + } + } + } + psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within); + if (within > (1.0 - stamps->normFrac) * sum) { + stamps->normWindow = radius; + done = true; + } + } + psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow); + if (stamps->normWindow == 0 || stamps->normWindow >= size) { + psError(PS_ERR_UNKNOWN, true, "Unable to determine normalisation window size."); + return false; } @@ -678,5 +721,5 @@ } -# ifdef TESTING +#if 1 { psFits *fits = psFitsOpen ("window.fits", "w"); @@ -684,8 +727,9 @@ psFitsClose (fits); } -# endif +#endif psFree (stats); - psFree (flux); + psFree (flux1); + psFree(flux2); psFree (norm1); psFree (norm2); @@ -788,5 +832,6 @@ pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, float sysErr, float skyErr, + int size, int footprint, float spacing, + float normFrac, float sysErr, float skyErr, pmSubtractionMode mode) { @@ -819,5 +864,6 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, - footprint, spacing, sysErr, skyErr, mode); // Stamps + footprint, spacing, normFrac, + sysErr, skyErr, mode); // Stamps psFree(x); psFree(y); @@ -833,6 +879,6 @@ pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, float sysErr, float skyErr, - pmSubtractionMode mode) + int size, int footprint, float spacing, float normFrac, + float sysErr, float skyErr, pmSubtractionMode mode) { PS_ASSERT_STRING_NON_EMPTY(filename, NULL); @@ -851,5 +897,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint, - spacing, sysErr, skyErr, mode); + spacing, normFrac, sysErr, skyErr, mode); psFree(data); Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h (revision 26702) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h (revision 26703) @@ -25,4 +25,6 @@ int footprint; ///< Half-size of stamps psKernel *window; ///< window function generated from ensemble of stamps + float normFrac; ///< Fraction of flux in window for normalisation window + int normWindow; ///< Size of window for measuring normalisation float sysErr; ///< Systematic error float skyErr; ///< increase effective readnoise @@ -30,11 +32,13 @@ /// Allocate a list of stamps -pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image - int numRows, // Number of rows in image - const psRegion *region, // Region for stamps, or NULL - int footprint, // Half-size of stamps - float spacing, // Rough average spacing between stamps - float sysErr, // Relative systematic error or NAN - float skyErr // Relative systematic error or NAN +pmSubtractionStampList *pmSubtractionStampListAlloc( + int numCols, // Number of columns in image + int numRows, // Number of rows in image + const psRegion *region, // Region for stamps, or NULL + int footprint, // Half-size of stamps + float spacing, // Rough average spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysErr, // Relative systematic error or NAN + float skyErr // Relative systematic error or NAN ); @@ -78,4 +82,5 @@ psImage *matrix; ///< Least-squares matrix, or NULL psVector *vector; ///< Least-squares vector, or NULL + double norm; ///< Normalisation difference pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp; @@ -85,31 +90,35 @@ /// Find stamps on an image -pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL - const psImage *image1, ///< Image for which to find stamps - const psImage *image2, ///< Image for which to find stamps - const psImage *mask, ///< Mask, or NULL - const psRegion *region, ///< Region to search, or NULL - float thresh1, ///< Threshold for stamps in image 1 - float thresh2, ///< Threshold for stamps in image 2 - int size, ///< Kernel half-size - int footprint, ///< Half-size for stamps - float spacing, ///< Rough spacing for stamps - float sysErr, ///< Relative systematic error in images - float skyErr, ///< Relative systematic error in images - pmSubtractionMode mode ///< Mode for subtraction +pmSubtractionStampList *pmSubtractionStampsFind( + pmSubtractionStampList *stamps, ///< Output stamps, or NULL + const psImage *image1, ///< Image for which to find stamps + const psImage *image2, ///< Image for which to find stamps + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float thresh1, ///< Threshold for stamps in image 1 + float thresh2, ///< Threshold for stamps in image 2 + int size, ///< Kernel half-size + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysErr, ///< Relative systematic error in images + float skyErr, ///< Relative systematic error in images + pmSubtractionMode mode ///< Mode for subtraction ); /// Set stamps based on a list of x,y -pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp - const psVector *y, ///< y coordinates for each stamp - const psImage *image, ///< Image for flux of stamp - const psImage *mask, ///< Mask, or NULL - const psRegion *region, ///< Region to search, or NULL - int size, ///< Kernel half-size - int footprint, ///< Half-size for stamps - float spacing, ///< Rough spacing for stamps - float sysErr, ///< Systematic error in images - float skyErr, ///< Systematic error in images - pmSubtractionMode mode ///< Mode for subtraction +pmSubtractionStampList *pmSubtractionStampsSet( + const psVector *x, ///< x coordinates for each stamp + const psVector *y, ///< y coordinates for each stamp + const psImage *image, ///< Image for flux of stamp + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + int size, ///< Kernel half-size + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysErr, ///< Systematic error in images + float skyErr, ///< Systematic error in images + pmSubtractionMode mode ///< Mode for subtraction ); @@ -123,4 +132,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window float sysErr, ///< Systematic error in images float skyErr, ///< Systematic error in images @@ -137,4 +147,5 @@ int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps + float normFrac, // Fraction of flux in window for normalisation window float sysErr, ///< Systematic error in images float skyErr, ///< Systematic error in images @@ -142,5 +153,9 @@ ); -bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize); +/// Calculate the window and normalisation window from the stamps +bool pmSubtractionStampsGetWindow( + pmSubtractionStampList *stamps, ///< List of stamps + int kernelSize ///< Half-size of kernel + ); /// Extract stamps from the images