Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmPSFEnvelope.c (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmPSFEnvelope.c (revision 29449) @@ -380,6 +380,5 @@ // measure the source moments: tophat windowing, no pixel S/N cutoff - // XXX probably should be passing the maskVal to this function so we can pass it along here... - if (!pmSourceMoments(source, maxRadius, 0.0, 1.0, maskVal)) { + if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, maskVal)) { // Can't do anything about it; limp along as best we can psErrorClear(); Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c (revision 29449) @@ -19,5 +19,7 @@ //#define USE_WEIGHT // Include weight (1/variance) in equation? -//#define USE_WINDOW // Include weight (1/variance) in equation? + +// XXX TEST: +# define USE_WINDOW // window to avoid neighbor contamination # define PENALTY false @@ -752,5 +754,5 @@ stamps->normValue2 = stats->robustMedian; - fprintf (stderr, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2); + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2); psFree(stats); @@ -797,5 +799,5 @@ stamp->normSquare2 = normSquare2; - fprintf (stderr, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); return true; @@ -1089,8 +1091,17 @@ } +#if 1 + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); + psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); +#endif + psVector *solution = NULL; // Solution to equation! solution = psVectorAlloc(numParams, PS_TYPE_F64); psVectorInit(solution, 0); + // XXX TEST: try some constraint on the svd solution + // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); @@ -1149,5 +1160,5 @@ } -#if 0 +#if 1 psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); @@ -1264,5 +1275,5 @@ if (!isfinite(dTotal)) return false; - fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum); + // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum); psVectorAppend(fResSigma, sigma/sum); psVectorAppend(fResOuter, dOuter/sum); Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c (revision 29449) @@ -258,5 +258,5 @@ } -#if 1 +#if 0 { double Sum = 0.0; // Sum of kernel component Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c (revision 29449) @@ -1202,5 +1202,5 @@ bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, - float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax) + float scaleRef, float scaleMin, float scaleMax) { PS_ASSERT_PTR_NON_NULL(kernelSize, false); @@ -1208,6 +1208,4 @@ PS_ASSERT_VECTOR_NON_NULL(widths, false); PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false); - PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false); PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false); PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false); @@ -1215,10 +1213,14 @@ 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"); + // 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 - // XXX save these values in a static for later use - pmSubtractionSetFWHMs(fwhm1, fwhm2); - if (isfinite(scaleMin) && scale < scaleMin) { scale = scaleMin; Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.h (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.h (revision 29449) @@ -104,5 +104,4 @@ int *stampSize, ///< Half-size of the stamp (footprint) psVector *widths, ///< ISIS widths - float fwhm1, float fwhm2, ///< FWHMs for inputs float scaleRef, ///< Reference width for scaling float scaleMin, ///< Minimum scaling ratio, or NAN Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c (revision 29449) @@ -51,4 +51,5 @@ psFree(list->y); psFree(list->flux); + psFree(list->window); psFree(list->window1); psFree(list->window2); @@ -96,6 +97,6 @@ // fprintf (stderr, "xMin, xMax: %d, %d -> ", xMin, xMax); - float xRaw = xMin; - float yRaw = yMin; + // float xRaw = xMin; + // float yRaw = yMin; // Ensure we're not going to go outside the bounds of the image @@ -126,5 +127,5 @@ if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) { - fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw); + // fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw); return false; } @@ -142,7 +143,8 @@ } - if (!found) { - fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw); - } + // if (!found) { + // fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw); + // } + return found; } @@ -245,4 +247,5 @@ list->normFrac = normFrac; list->normValue = NAN; + list->window = NULL; list->window1 = NULL; list->window2 = NULL; @@ -269,4 +272,5 @@ out->y = NULL; out->flux = NULL; + out->window = psMemIncrRefCounter(in->window); out->window1 = psMemIncrRefCounter(in->window1); out->window2 = psMemIncrRefCounter(in->window2); @@ -487,9 +491,4 @@ yList->data.F32[j], yList->data.F32[j], numCols, numRows, border); #endif - if (0 && goodStamp) { - fprintf (stderr, "find: %6.1f < %6.1f < %6.1f, %6.1f < %6.1f %6.1f\n", - region->x0 + size, xStamp, region->x1 - size, - region->y0 + size, yStamp, region->y1 - size); - } } } else { @@ -684,4 +683,24 @@ psImageInit(stamps->window2->image, 0.0); + // Generate a weighting window based on the fwhms (20% larger than the largest) + { + float fwhm1, fwhm2; + + // XXX this is annoyingly hack-ish + pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + + float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35; + + psFree (stamps->window); + stamps->window = psKernelAlloc(-size, size, -size, size); + psImageInit(stamps->window->image, 0.0); + + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma)); + } + } + } + // generate normalizations for each stamp psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32); @@ -759,8 +778,8 @@ { psFits *fits = NULL; - fits = psFitsOpen ("window1.fits", "w"); + fits = psFitsOpen ("window1.raw.fits", "w"); psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL); psFitsClose (fits); - fits = psFitsOpen ("window2.fits", "w"); + fits = psFitsOpen ("window2.raw.fits", "w"); psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL); psFitsClose (fits); @@ -818,5 +837,5 @@ // interpolate to the radius at which delta2 is normFrac: stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1); - fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1); + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1); done1 = true; } @@ -825,5 +844,5 @@ // interpolate to the radius at which delta2 is normFrac: stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2); - fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2); + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2); done2 = true; } @@ -873,8 +892,8 @@ { psFits *fits = NULL; - fits = psFitsOpen ("window1.fits", "w"); + fits = psFitsOpen ("window1.norm.fits", "w"); psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL); psFitsClose (fits); - fits = psFitsOpen ("window2.fits", "w"); + fits = psFitsOpen ("window2.norm.fits", "w"); psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL); psFitsClose (fits); @@ -962,5 +981,5 @@ psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty"); - fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2); + // fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2); return true; @@ -989,6 +1008,6 @@ penalty *= 1.0 / sum2; - if (1) { - fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2); + if (0) { + // fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2); // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty); } Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h (revision 29449) @@ -29,4 +29,5 @@ psKernel *window1; ///< window function generated from ensemble of stamps (input 1) psKernel *window2; ///< window function generated from ensemble of stamps (input 2) + psKernel *window; ///< weighting window function (sigma = 1.1 * MAX(fwhm)) float normWindow1; ///< Size of window for measuring normalisation float normWindow2; ///< Size of window for measuring normalisation Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c (revision 29448) +++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c (revision 29449) @@ -146,4 +146,5 @@ } pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true); + psFree(canvas); pmVisualAskUser(&plotStamps); @@ -281,8 +282,9 @@ } } - fprintf (stderr, "kernel %d, sum %f\n", i, sum); + // fprintf (stderr, "kernel %d, sum %f\n", i, sum); } pmVisualScaleImage(kapa1, output, "Image", 0, true); + psFree(output); pmVisualAskUser(&plotImage); return true; @@ -309,12 +311,12 @@ if (!isfinite(stamp->flux)) continue; if (!stamp->convolutions1 && !stamp->convolutions2) continue; - fprintf (stderr, "flux: %f, maxFlux: %f ", stamp->flux, maxFlux); + // fprintf (stderr, "flux: %f, maxFlux: %f ", stamp->flux, maxFlux); if (!maxStamp) { maxFlux = stamp->flux; maxStamp = stamp; - fprintf (stderr, "maxStamp %d\n", i); + // fprintf (stderr, "maxStamp %d\n", i); continue; } else { - fprintf (stderr, "\n"); + // fprintf (stderr, "\n"); } if (stamp->flux > maxFlux) { @@ -361,5 +363,5 @@ } } - fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); + // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); } pmVisualScaleImage(kapa2, output, "Image", 0, true); @@ -403,5 +405,5 @@ } } - fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); + // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); } pmVisualScaleImage(kapa2, output, "Image", 1, true); @@ -412,4 +414,5 @@ psFitsClose(fits); } + psFree(output); } @@ -708,4 +711,5 @@ psFree (x); psFree (y); + psFree (polyValues); pmVisualAskUser(NULL);