Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 29004) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 29543) @@ -28,10 +28,5 @@ static bool useFFT = true; // Do convolutions using FFT -# define SEPARATE 0 -# if (SEPARATE) -# define SUBMODE PM_SUBTRACTION_EQUATION_NORM -# else # define SUBMODE PM_SUBTRACTION_EQUATION_ALL -# endif //#define TESTING @@ -280,26 +275,26 @@ } -bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) { - - if (!readout) return true; - - psImage *image = readout->image; - psImage *mask = readout->mask; - psImage *variance = readout->variance; - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue; - bool valid = false; - valid = isfinite(image->data.F32[y][x]); - if (variance) { - valid &= isfinite(variance->data.F32[y][x]); - } - if (valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal; - } - } - - return true; -} +// bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) { +// +// if (!readout) return true; +// +// psImage *image = readout->image; +// psImage *mask = readout->mask; +// psImage *variance = readout->variance; +// for (int y = 0; y < image->numRows; y++) { +// for (int x = 0; x < image->numCols; x++) { +// if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue; +// bool valid = false; +// valid = isfinite(image->data.F32[y][x]); +// if (variance) { +// valid &= isfinite(variance->data.F32[y][x]); +// } +// if (valid) continue; +// mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal; +// } +// } +// +// return true; +// } bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, @@ -392,6 +387,7 @@ } - pmSubtractionMaskInvalid(ro1, maskVal); - pmSubtractionMaskInvalid(ro2, maskVal); + // XXX this is done before calling this function + // pmSubtractionMaskInvalid(ro1, maskVal); + // pmSubtractionMaskInvalid(ro2, maskVal); // General background subtraction, since this is done in pmSubtractionMatch @@ -565,6 +561,6 @@ memCheck("start"); - pmSubtractionMaskInvalid(ro1, maskVal); - pmSubtractionMaskInvalid(ro2, maskVal); + // pmSubtractionMaskInvalid(ro1, maskVal); + // pmSubtractionMaskInvalid(ro2, maskVal); psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels @@ -672,5 +668,4 @@ kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, inner, binning, ringsOrder, penalty, bounds, subMode); - // pmSubtractionVisualShowKernels(kernels); } @@ -678,27 +673,4 @@ if (subMode == PM_SUBTRACTION_MODE_UNSURE) { -#if 0 - // Get backgrounds - psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background - psVector *buffer = NULL;// Buffer for stats - if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) { - psError(PM_ERR_DATA, false, "Unable to measure background of image 1."); - psFree(bgStats); - psFree(buffer); - goto MATCH_ERROR; - } - float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1 - if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) { - psError(PM_ERR_DATA, false, "Unable to measure background of image 2."); - psFree(bgStats); - psFree(buffer); - goto MATCH_ERROR; - } - float bg2 = psStatsGetValue(bgStats, BG_STAT); // Background for image 2 - psFree(bgStats); - psFree(buffer); - - pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use -#endif pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); switch (newMode) { @@ -732,6 +704,21 @@ } - // XXX step 1: calculate normalization - psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); + // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue + psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); + if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + // step 0b : calculate the moments, pass along to the next steps via stamps->normValue + // XXX this step is now skipped -- delete + psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); + if (!pmSubtractionCalculateMoments(kernels, stamps)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + // step 1: generate the elements of the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); @@ -739,5 +726,6 @@ } - psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n"); + // step 2: solve the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Solving kernel equations...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); @@ -746,23 +734,4 @@ memCheck(" solve equation"); -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters - psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" calculate equation"); - - psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" solve equation"); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -770,5 +739,4 @@ goto MATCH_ERROR; } - memCheck(" calculate deviations"); @@ -787,5 +755,6 @@ // if we hit the max number of iterations and we have rejected stamps, re-solve if (numRejected > 0) { - // XXX step 1: calculate normalization + + // step 1: generate the elements of the matrix equation Ax = B psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { @@ -794,5 +763,5 @@ } - // solve normalization + // step 2: solve the matrix equation Ax = B psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { @@ -801,23 +770,4 @@ } -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters - psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - // solve kernel parameters - psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" solve equation"); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -837,5 +787,4 @@ goto MATCH_ERROR; } - memCheck("diag outputs"); @@ -1134,21 +1083,4 @@ } -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } -# endif - psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations @@ -1181,22 +1113,4 @@ } psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); - -# if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); -# endif psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations @@ -1288,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); @@ -1294,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); @@ -1301,9 +1213,13 @@ PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); -// float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference + 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) {