Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (revision 30382) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (revision 30383) @@ -1072,5 +1072,5 @@ // SINGLE solution # if (1) - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10); invMatrix = psMatrixInvert(NULL, sumMatrix, NULL); # endif @@ -1176,5 +1176,5 @@ // DUAL solution # if (1) - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10); invMatrix = psMatrixInvert(NULL, sumMatrix, NULL); # endif @@ -1306,5 +1306,5 @@ // given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq -bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) { +bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window) { # ifndef USE_WEIGHT @@ -1316,22 +1316,27 @@ int npix = 0; - float chisq = 0; + float chisqR = 0; + float chisqD = 0; // get the chisq for (int y = residual->yMin; y <= residual->yMax; y++) { for (int x = residual->xMin; x <= residual->xMax; x++) { - float value = PS_SQR(residual->kernel[y][x]); + float valueR = PS_SQR(residual->kernel[y][x]); if (weight) { - value *= weight->kernel[y][x]; + valueR *= weight->kernel[y][x]; } - // XXX NOTE: do NOT apply the window to the chisq portions of the calculation - if (false && window) { - value *= window->kernel[y][x]; + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation (that would bias the chisq) + chisqR += valueR; + + float valueD = PS_SQR(difference->kernel[y][x]); + if (weight) { + valueD *= weight->kernel[y][x]; } - chisq += value; + chisqD += valueD; npix ++; } } - psVectorAppend(chisqVector, chisq / npix); + psVectorAppend(chisqRVector, chisqR / npix); + psVectorAppend(chisqDVector, chisqD / npix); float value1 = 0; @@ -1430,5 +1435,5 @@ int Nelem = fluxesVector->n - 1; bool valid = true; - valid &= isfinite(chisqVector->data.F32[Nelem]); + valid &= isfinite(chisqRVector->data.F32[Nelem]); valid &= isfinite(fluxesVector->data.F32[Nelem]); valid &= isfinite(momentVector->data.F32[Nelem]); @@ -1453,5 +1458,6 @@ // XXX need to save these somewhere psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); - psVector *chisq = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *chisqD = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *chisqR = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK); @@ -1464,4 +1470,5 @@ // storage for the image (convolved2 is not used in SINGLE mode) psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + psKernel *difference = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image @@ -1474,5 +1481,6 @@ psVectorAppend(moments, NAN); psVectorAppend(fluxes, NAN); - psVectorAppend(chisq, NAN); + psVectorAppend(chisqD, NAN); + psVectorAppend(chisqR, NAN); psVectorAppend(stampMask, 0x01); continue; @@ -1487,4 +1495,5 @@ // Calculate residuals psImageInit(residual->image, 0.0); + psImageInit(difference->image, 0.0); psKernel *weight = NULL; @@ -1537,13 +1546,17 @@ } - // generate the residual image + // Generate the difference, residual, and convolved source images. Note the we + // accumulate the convolution of (A-B), so we need to replace it to generate the + // images of the convolved source image. for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { + difference->kernel[y][x] = target->kernel[y][x] - source->kernel[y][x] * norm - background; + residual->kernel[y][x] = difference->kernel[y][x] - convolved1->kernel[y][x]; convolved1->kernel[y][x] += source->kernel[y][x] * norm; - residual->kernel[y][x] = target->kernel[y][x] - convolved1->kernel[y][x] - background; } } + // XXX if we want to have a weight and window, we'll need to pass through to here - pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, NULL, residual, weight, window); + pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, NULL, difference, residual, weight, window); } else { @@ -1574,9 +1587,13 @@ } + // Generate the difference, residual, and convolved source images. Note the we + // accumulate the convolutions of (A-B), so we need to replace (A or B) to generate + // the images of the convolved source images. for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { + difference->kernel[y][x] = image2->kernel[y][x] - image1->kernel[y][x] * norm - background; + residual->kernel[y][x] = difference->kernel[y][x] + convolved2->kernel[y][x] - convolved1->kernel[y][x]; convolved1->kernel[y][x] += image1->kernel[y][x] * norm; convolved2->kernel[y][x] += image2->kernel[y][x]; - residual->kernel[y][x] = convolved2->kernel[y][x] - convolved1->kernel[y][x] - background; } } @@ -1589,5 +1606,5 @@ } - pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, convolved2, residual, weight, window); + pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, convolved2, difference, residual, weight, window); } } @@ -1595,6 +1612,10 @@ // find the mean chisq and mean moment psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - psVectorStats (stats, chisq, NULL, stampMask, 0xff); - float chisqValue = stats->sampleMean; + psVectorStats (stats, chisqD, NULL, stampMask, 0xff); + float chisqDValue = stats->sampleMean; + + psStatsInit(stats); + psVectorStats (stats, chisqR, NULL, stampMask, 0xff); + float chisqRValue = stats->sampleMean; psStatsInit(stats); @@ -1637,6 +1658,6 @@ // penalized by increasing the score somewhat. the 0.01 value is not well-chosen. float orderFactor = 0.01 * kernels->spatialOrder; - float score = 2.0 * chisqValue / (sumKernel1 + sumKernel2) + orderFactor; - psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqValue, momentValue, sumKernel1, sumKernel2, score); + float score = 2.0 * chisqRValue / (sumKernel1 + sumKernel2) + orderFactor; + psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, chisqD: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqRValue, chisqDValue, momentValue, sumKernel1, sumKernel2, score); // save this result if it is the first or the best (skip if bestMatch is NULL) @@ -1663,5 +1684,5 @@ match->nGood = nGood; match->fluxes = psMemIncrRefCounter(fluxes); - match->chisq = psMemIncrRefCounter(chisq); + match->chisq = psMemIncrRefCounter(chisqR); match->moments = psMemIncrRefCounter(moments); match->stampMask = psMemIncrRefCounter(stampMask); @@ -1669,8 +1690,9 @@ } - pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments); + pmSubtractionVisualPlotChisqAndMoments(fluxes, chisqR, moments); psFree(stats); - psFree(chisq); + psFree(chisqR); + psFree(chisqD); psFree(fluxes); psFree(moments); Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h (revision 30382) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h (revision 30383) @@ -82,5 +82,5 @@ bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window); -bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window); +bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window); bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, pmSubtractionStampList *stamps, pmSubtractionKernels *kernels); Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c (revision 30382) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c (revision 30383) @@ -476,5 +476,5 @@ const psArray *sources, const char *stampsName, pmSubtractionKernelsType type, int size, int spatialOrder, - const psVector *isisWidths, const psVector *isisOrders, + psVector *isisWidths, const psVector *isisOrders, int inner, int ringsOrder, int binning, float penalty, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, @@ -695,9 +695,9 @@ // check on the kernel scaling -- if the kron-based radial moments are very different, adjust to match them { - float fwhm1; - float fwhm2; - pmSubtractionGetFWHMs(&fwhm1, &fwhm2); - psAssert(isfinite(fwhm1), "fwhm 1 not set"); - psAssert(isfinite(fwhm2), "fwhm 2 not set"); + // float fwhm1; + // float fwhm2; + // pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + // psAssert(isfinite(fwhm1), "fwhm 1 not set"); + // psAssert(isfinite(fwhm2), "fwhm 2 not set"); // XXX this is BAD: depends on the relationship below: @@ -706,18 +706,19 @@ float radMoment1 = stamps->normWindow1 / 2.75; float radMoment2 = stamps->normWindow2 / 2.75; - - float maxFWHM = PS_MAX(fwhm1, fwhm2); - float maxRadial = PS_MAX(radMoment1, radMoment2); + pmSubtractionParamsScale(NULL, NULL, isisWidths, radMoment1, radMoment2); + + // float maxFWHM = PS_MAX(fwhm1, fwhm2); + // float maxRadial = PS_MAX(radMoment1, radMoment2); // if (fabs(2.0*(maxFWHM - maxRadial)/(maxFWHM + maxRadial)) > 0.25) { - if (1) { - - float scale = maxRadial / maxFWHM; - psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale); - - for (int i = 0; i < isisWidths->n; i++) { - isisWidths->data.F32[i] *= scale; - } - } + // if (1) { + // + // float scale = maxRadial / maxFWHM; + // psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale); + // + // for (int i = 0; i < isisWidths->n; i++) { + // isisWidths->data.F32[i] *= scale; + // } + // } } @@ -1252,32 +1253,47 @@ } - -bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, - float scaleRef, float scaleMin, float scaleMax) +static float scaleRefOption = NAN; +static float scaleMinOption = NAN; +static float scaleMaxOption = NAN; +static bool scaleOption = false; + +bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax) { + + if (scale) { + PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); + } + + scaleRefOption = scaleRef; + scaleMinOption = scaleMin; + scaleMaxOption = scaleMax; + scaleOption = scale; + + return true; +} + +bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2) { - PS_ASSERT_PTR_NON_NULL(kernelSize, false); - PS_ASSERT_PTR_NON_NULL(stampSize, false); + // PS_ASSERT_PTR_NON_NULL(kernelSize, false); + // PS_ASSERT_PTR_NON_NULL(stampSize, false); PS_ASSERT_VECTOR_NON_NULL(widths, false); PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); - - float fwhm1; - float fwhm2; - - pmSubtractionGetFWHMs(&fwhm1, &fwhm2); - psAssert(isfinite(fwhm1), "fwhm 1 not set"); - psAssert(isfinite(fwhm2), "fwhm 2 not set"); + + if (!scaleOption) return true; + + // pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + // psAssert(isfinite(fwhm1), "fwhm 1 not set"); + // psAssert(isfinite(fwhm2), "fwhm 2 not set"); // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference - float scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor - - if (isfinite(scaleMin) && scale < scaleMin) { - scale = scaleMin; - } - if (isfinite(scaleMax) && scale > scaleMax) { - scale = scaleMax; + float scale = PS_MAX(fwhm1, fwhm2) / scaleRefOption; // Scaling factor + + if (isfinite(scaleMinOption) && scale < scaleMinOption) { + scale = scaleMinOption; + } + if (isfinite(scaleMaxOption) && scale > scaleMaxOption) { + scale = scaleMaxOption; } @@ -1285,9 +1301,14 @@ widths->data.F32[i] *= scale; } - *kernelSize = *kernelSize * scale + 0.5; - *stampSize = *stampSize * scale + 0.5; - - psLogMsg("psModules.imcombine", PS_LOG_INFO, - "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize); + if (kernelSize) { + *kernelSize = *kernelSize * scale + 0.5; + } + if (stampSize) { + *stampSize = *stampSize * scale + 0.5; + } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Scaling kernel parameters by %f", scale); + if (kernelSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified kernel size %d", *kernelSize); + if (stampSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified stamp size %d", *stampSize); return true; Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h (revision 30382) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h (revision 30383) @@ -27,5 +27,5 @@ int size, ///< Kernel half-size int order, ///< Spatial polynomial order - const psVector *widths, ///< ISIS Gaussian widths + psVector *widths, ///< ISIS Gaussian widths const psVector *orders, ///< ISIS Polynomial orders int inner, ///< Inner radius for various kernel types @@ -101,12 +101,16 @@ /// Scale subtraction parameters according to the FWHMs of the inputs -bool pmSubtractionParamsScale( - int *kernelSize, ///< Half-size of the kernel - int *stampSize, ///< Half-size of the stamp (footprint) - psVector *widths, ///< ISIS widths - float scaleRef, ///< Reference width for scaling - float scaleMin, ///< Minimum scaling ratio, or NAN - float scaleMax ///< Maximum scaling ratio, or NAN - ); +// bool pmSubtractionParamsScale( +// int *kernelSize, ///< Half-size of the kernel +// int *stampSize, ///< Half-size of the stamp (footprint) +// psVector *widths, ///< ISIS widths +// float scaleRef, ///< Reference width for scaling +// float scaleMin, ///< Minimum scaling ratio, or NAN +// float scaleMax ///< Maximum scaling ratio, or NAN +// ); + +bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2); + +bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax); bool pmSubtractionMatchAttempt(