Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26657) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26667) @@ -66,7 +66,7 @@ // Contribute to an image of the solved kernel component using the preCalculated image static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated - const pmSubtractionKernels *kernels, // Kernel basis functions - float value, // Normalisation value for basis function - int index // Index of basis function of interest + const pmSubtractionKernels *kernels, // Kernel basis functions + float value, // Normalisation value for basis function + int index // Index of basis function of interest ) { @@ -162,10 +162,10 @@ break; } - case PM_SUBTRACTION_KERNEL_ISIS: - case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: - case PM_SUBTRACTION_KERNEL_HERM: + case PM_SUBTRACTION_KERNEL_ISIS: + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + case PM_SUBTRACTION_KERNEL_HERM: case PM_SUBTRACTION_KERNEL_DECONV_HERM: { - solvedKernelPreCalc(kernel, kernels, value, i); - break; + solvedKernelPreCalc(kernel, kernels, value, i); + break; } case PM_SUBTRACTION_KERNEL_RINGS: { @@ -460,7 +460,7 @@ // Convolve a stamp using a pre-calculated kernel basis function static psKernel *convolveStampPreCalc(const psKernel *image, // Image to convolve - const pmSubtractionKernels *kernels, // Kernel basis functions - int index, // Index of basis function of interest - int footprint // Half-size of stamp + const pmSubtractionKernels *kernels, // Kernel basis functions + int index, // Index of basis function of interest + int footprint // Half-size of stamp ) { @@ -678,10 +678,10 @@ return convolved; } - case PM_SUBTRACTION_KERNEL_ISIS: - case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + case PM_SUBTRACTION_KERNEL_ISIS: + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: case PM_SUBTRACTION_KERNEL_HERM: case PM_SUBTRACTION_KERNEL_DECONV_HERM: { - return convolveStampPreCalc(image, kernels, index, footprint); - } + return convolveStampPreCalc(image, kernels, index, footprint); + } case PM_SUBTRACTION_KERNEL_RINGS: { psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image @@ -864,4 +864,6 @@ int numGood = 0; // Number of good stamps double newMean = 0.0; // New mean + psString log = NULL; // Log message + psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit); for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -877,4 +879,7 @@ psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i, + (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), + fabsf(deviations->data.F32[i] - mean)); numRejected++; for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { @@ -911,4 +916,10 @@ } newMean /= numGood; + + if (numRejected == 0) { + psStringAppend(&log, "\n"); + } + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log); + psFree(log); if (ds9) { Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26657) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26667) @@ -1425,4 +1425,5 @@ } + psString log = psStringCopy("Deviations:\n"); // Log message with deviations for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest @@ -1556,4 +1557,6 @@ psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); + psStringAppend(&log, "Stamp %d (%d,%d): %f\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); if (!isfinite(deviations->data.F32[i])) { stamp->status = PM_SUBTRACTION_STAMP_REJECTED; @@ -1600,4 +1603,7 @@ } + + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log); + psFree(log); // calculate and report the normalization and background for the image center Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26657) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c (revision 26667) @@ -496,8 +496,8 @@ - // generate the window function from the set of stamps - if (!pmSubtractionStampsGetWindow(stamps, size)) { - goto MATCH_ERROR; - } + // generate the window function from the set of stamps + if (!pmSubtractionStampsGetWindow(stamps, size)) { + goto MATCH_ERROR; + } // Define kernel basis functions @@ -514,5 +514,5 @@ kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, inner, binning, ringsOrder, penalty, subMode); - // pmSubtractionVisualShowKernels(kernels); + // pmSubtractionVisualShowKernels(kernels); } @@ -568,10 +568,10 @@ } - // generate the window function from the set of stamps - if (!pmSubtractionStampsGetWindow(stamps, size)) { - goto MATCH_ERROR; - } - - // XXX step 1: calculate normalization + // generate the window function from the set of stamps + if (!pmSubtractionStampsGetWindow(stamps, size)) { + goto MATCH_ERROR; + } + + // XXX step 1: calculate normalization psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { @@ -588,8 +588,8 @@ # if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters + // 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)) { @@ -626,7 +626,7 @@ } - // if we hit the max number of iterations and we have rejected stamps, re-solve + // if we hit the max number of iterations and we have rejected stamps, re-solve if (numRejected > 0) { - // XXX step 1: calculate normalization + // XXX step 1: calculate normalization psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { @@ -635,5 +635,5 @@ } - // solve normalization + // solve normalization psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { @@ -643,8 +643,8 @@ # if (SEPARATE) - // set USED -> CALCULATE - pmSubtractionStampsResetStatus (stamps); - - // XXX step 2: calculate kernel parameters + // 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)) { @@ -653,5 +653,5 @@ } - // solve kernel parameters + // solve kernel parameters psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) { @@ -992,5 +992,5 @@ return false; } -# endif +# endif psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); @@ -1012,9 +1012,9 @@ if (numRejected > 0) { // Allow re-fit with reduced stamps set - psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - return false; - } + psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); + if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); @@ -1026,12 +1026,12 @@ # 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(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - return false; - } + // set USED -> CALCULATE + pmSubtractionStampsResetStatus (stamps); + + psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); + if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); @@ -1129,2 +1129,38 @@ } + +bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, + float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax) +{ + 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(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); + PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false); + 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 scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor + + if (isfinite(scaleMin) && scale < scaleMin) { + scale = scaleMin; + } + if (isfinite(scaleMax) && scale > scaleMax) { + scale = scaleMax; + } + + for (int i = 0; i < widths->n; i++) { + 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); + + return true; +} Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h (revision 26657) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h (revision 26667) @@ -95,3 +95,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 fwhm1, float fwhm2, ///< FWHMs for inputs + float scaleRef, ///< Reference width for scaling + float scaleMin, ///< Minimum scaling ratio, or NAN + float scaleMax ///< Maximum scaling ratio, or NAN + ); + + #endif