Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26505) @@ -1695,8 +1695,8 @@ double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background); - } - - pmSubtractionVisualShowFit(); - pmSubtractionVisualPlotFit(kernels); + + pmSubtractionVisualShowFit(norm); + pmSubtractionVisualPlotFit(kernels); + } psFree(residual); Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26505) @@ -177,5 +177,5 @@ } } -#if 1 +#if 0 fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, moment); #endif @@ -211,5 +211,5 @@ } -#if 1 +#if 0 sum = 0.0; // Sum of kernel component min = FLT_MAX; @@ -409,5 +409,6 @@ preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false); + // XXX do we use Alard-Lupton normalization (last param true) or not? + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); // XXXX test demo that deconvolved kernel is valid Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26505) @@ -75,4 +75,5 @@ float stampSpacing, // Spacing between stamps float sysError, // Relative systematic error in images + float skyError, // Relative systematic error in images int size, // Kernel half-size int footprint, // Convolution footprint for stamps @@ -85,5 +86,5 @@ *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, - size, footprint, stampSpacing, sysError, mode); + size, footprint, stampSpacing, sysError, skyError, mode); if (!*stamps) { psError(psErrorCodeLast(), false, "Unable to find stamps."); @@ -109,4 +110,5 @@ int stride, // Size for convolution patches float sysError, // Systematic error in images + float skyError, // Systematic error in images float kernelError, // Systematic error in kernel psImageMaskType maskVal, // Value to mask for input @@ -159,4 +161,7 @@ PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false); } + if (isfinite(sysError)) { + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.0, false); + } if (isfinite(kernelError)) { PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false); @@ -276,5 +281,5 @@ } - if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, kernelError, + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) { psFree(kernels); @@ -342,9 +347,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 kernelError, psImageMaskType maskVal, + 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, kernelError, + if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, skyError, kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) { return false; @@ -476,8 +481,8 @@ if (stampsName && strlen(stampsName) > 0) { stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size, - footprint, stampSpacing, sysError, subMode); + footprint, stampSpacing, sysError, skyError, subMode); } else if (sources) { stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size, - footprint, stampSpacing, sysError, subMode); + footprint, stampSpacing, sysError, skyError, subMode); } @@ -485,5 +490,5 @@ // doesn't matter. if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2, - stampSpacing, sysError, size, footprint, subMode)) { + stampSpacing, sysError, skyError, size, footprint, subMode)) { goto MATCH_ERROR; } @@ -557,5 +562,5 @@ if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, - stampThresh1, stampThresh2, stampSpacing, sysError, + stampThresh1, stampThresh2, stampSpacing, sysError, skyError, size, footprint, subMode)) { goto MATCH_ERROR; Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h (revision 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h (revision 26505) @@ -40,4 +40,5 @@ float rej, ///< Rejection threshold float sysError, ///< Relative systematic error in images + float skyError, ///< Relative systematic error in images float kernelError, ///< Relative systematic error in kernel psImageMaskType maskVal, ///< Value to mask for input Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (revision 26505) @@ -181,5 +181,5 @@ pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, - int footprint, float spacing, float sysErr) + int footprint, float spacing, float sysErr, float skyErr) { pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return @@ -224,4 +224,5 @@ list->footprint = footprint; list->sysErr = sysErr; + list->skyErr = skyErr; return list; @@ -326,5 +327,5 @@ const psImage *image2, const psImage *subMask, const psRegion *region, float thresh1, float thresh2, - int size, int footprint, float spacing, float sysErr, + int size, int footprint, float spacing, float sysErr, float skyErr, pmSubtractionMode mode) { @@ -383,5 +384,5 @@ if (!stamps) { - stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr); + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr); } @@ -441,5 +442,5 @@ xStamp = xList->data.F32[j]; yStamp = yList->data.F32[j]; - fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp); + // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp); } } else { @@ -496,5 +497,5 @@ const psImage *image, const psImage *subMask, const psRegion *region, int size, int footprint, - float spacing, float sysErr, pmSubtractionMode mode) + float spacing, float sysErr, float skyErr, pmSubtractionMode mode) { @@ -517,5 +518,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, region, footprint, spacing, - sysErr); // Stamp list + sysErr, skyErr); // Stamp list int numStamps = stamps->num; // Number of stamps @@ -548,5 +549,5 @@ } - fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix); + // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix); bool found = false; @@ -742,5 +743,5 @@ return false; } - fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y); + // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y); // Catch memory leaks --- these should have been freed and NULLed before @@ -765,11 +766,14 @@ psImage *varSub = psImageSubset(variance, region); // Subimage with stamp psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp - if (isfinite(stamps->sysErr) && stamps->sysErr > 0) { + + if ((isfinite(stamps->skyErr) && (stamps->skyErr > 0)) || + (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) { float sysErr = 0.25 * PS_SQR(stamps->sysErr); // Systematic error + float skyErr = stamps->skyErr; psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images for (int y = -size; y <= size; y++) { for (int x = -size; x <= size; x++) { float additional = image1->kernel[y][x] + image2->kernel[y][x]; - weight->kernel[y][x] = 1.0 / (var->kernel[y][x] + sysErr * PS_SQR(additional)); + weight->kernel[y][x] = 1.0 / (skyErr + var->kernel[y][x] + sysErr * PS_SQR(additional)); } } @@ -795,5 +799,5 @@ pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, float sysErr, + int size, int footprint, float spacing, float sysErr, float skyErr, pmSubtractionMode mode) { @@ -819,5 +823,5 @@ y->data.F32[numOut] = source->peak->yf; } - fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]); + // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]); numOut++; } @@ -826,5 +830,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, - footprint, spacing, sysErr, mode); // Stamps + footprint, spacing, sysErr, skyErr, mode); // Stamps psFree(x); psFree(y); @@ -840,5 +844,5 @@ pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image, const psImage *subMask, const psRegion *region, - int size, int footprint, float spacing, float sysErr, + int size, int footprint, float spacing, float sysErr, float skyErr, pmSubtractionMode mode) { @@ -858,5 +862,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint, - spacing, sysErr, mode); + spacing, 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 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h (revision 26505) @@ -26,4 +26,5 @@ psKernel *window; ///< window function generated from ensemble of stamps float sysErr; ///< Systematic error + float skyErr; ///< Systematic error } pmSubtractionStampList; @@ -34,5 +35,6 @@ int footprint, // Half-size of stamps float spacing, // Rough average spacing between stamps - float sysErr // Relative systematic error or NAN + float sysErr, // Relative systematic error or NAN + float skyErr // Relative systematic error or NAN ); @@ -95,4 +97,5 @@ 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 ); @@ -108,4 +111,5 @@ float spacing, ///< Rough spacing for stamps float sysErr, ///< Systematic error in images + float skyErr, ///< Systematic error in images pmSubtractionMode mode ///< Mode for subtraction ); @@ -121,4 +125,5 @@ float spacing, ///< Rough spacing for stamps float sysErr, ///< Systematic error in images + float skyErr, ///< Systematic error in images pmSubtractionMode mode ///< Mode for subtraction ); @@ -134,4 +139,5 @@ float spacing, ///< Rough spacing for stamps float sysErr, ///< Systematic error in images + float skyErr, ///< Systematic error in images pmSubtractionMode mode ///< Mode for subtraction ); Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26505) @@ -535,5 +535,15 @@ } -bool pmSubtractionVisualShowFit() { +bool pmSubtractionVisualShowFit(double norm) { + + // XXX a dumb test : dump the residual image and exit + { + psMetadata *header = psMetadataAlloc(); + psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm); + psFits *fits = psFitsOpen("resid.fits", "w"); + psFitsWriteImage(fits, header, residualImage, 0, NULL); + psFitsClose(fits); + exit (0); + } // if (!pmVisualIsVisual()) return true; Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.h (revision 26504) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.h (revision 26505) @@ -9,5 +9,5 @@ bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps); bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index); -bool pmSubtractionVisualShowFit(); +bool pmSubtractionVisualShowFit(double norm); bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels); bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels);