Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 25101) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 25120) @@ -476,5 +476,24 @@ } + + // Define kernel basis functions + if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { + kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, + stamps, footprint, optThreshold, penalty, subMode); + if (!kernels) { + psErrorClear(); + psWarning("Unable to derive optimum ISIS kernel --- switching to default."); + } + } + if (kernels == NULL) { + // Not an ISIS/GUNK kernel, or the optimum kernel search failed + kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, + inner, binning, ringsOrder, penalty, subMode); + } + + memCheck("kernels"); + if (subMode == PM_SUBTRACTION_MODE_UNSURE) { +#if 0 // Get backgrounds psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background @@ -498,4 +517,6 @@ pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use +#endif + pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); switch (newMode) { case PM_SUBTRACTION_MODE_1: @@ -512,21 +533,4 @@ } - // Define kernel basis functions - if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold, penalty, subMode); - if (!kernels) { - psErrorClear(); - psWarning("Unable to derive optimum ISIS kernel --- switching to default."); - } - } - if (kernels == NULL) { - // Not an ISIS/GUNK kernel, or the optimum kernel search failed - kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder, penalty, subMode); - } - - memCheck("kernels"); - int numRejected = -1; // Number of rejected stamps in each iteration for (int k = 0; k < iter && numRejected != 0; k++) { @@ -565,5 +569,5 @@ psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint); + numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); @@ -587,5 +591,5 @@ goto MATCH_ERROR; } - pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint); + pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); psFree(deviations); } @@ -871,58 +875,10 @@ -#if 0 -/// A list of stamps -typedef struct { - long num; ///< Number of stamps - psArray *stamps; ///< The stamps - psArray *regions; ///< Regions for each stamp - psArray *x, *y; ///< Coordinates for possible stamps (or NULL) - psArray *flux; ///< Fluxes for possible stamps (or NULL) - int footprint; ///< Half-size of stamps -} pmSubtractionStampList; - -/// A stamp for image subtraction -typedef struct { - float x, y; ///< Position - float flux; ///< Flux - float xNorm, yNorm; ///< Normalised position - psKernel *image1; ///< Reference image postage stamp - psKernel *image2; ///< Input image postage stamp - psKernel *variance; ///< Variance image postage stamp, or NULL - psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL - psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL - psImage *matrix1, *matrix2; ///< Least-squares matrices for each image, or NULL - psImage *matrixX; ///< Cross-matrix (for mode DUAL), or NULL - psVector *vector1, *vector2; ///< Least-squares vectors for each image, or NULL - pmSubtractionStampStatus status; ///< Status of stamp -} pmSubtractionStamp; - -/// Kernels specification -typedef struct { - pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels - psString description; ///< Description of the kernel parameters - long num; ///< Number of kernel components (not including the spatial ones) - psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) - psVector *widths; ///< Gaussian FWHMs (ISIS) - psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) - float penalty; ///< Penalty for wideness - psVector *penalties; ///< Penalty for each kernel component - int size; ///< The half-size of the kernel - int inner; ///< The size of an inner region - int spatialOrder; ///< The spatial order of the kernels - int bgOrder; ///< The order for the background fitting - pmSubtractionMode mode; ///< Mode for subtraction - int numCols, numRows; ///< Size of image (for normalisation), or zero to use image provided - psVector *solution1, *solution2; ///< Solution for the PSF matching - // Quality information - float mean, rms; ///< Mean and RMS of chi^2 from stamps - int numStamps; ///< Number of good stamps -} pmSubtractionKernels; - // Test a subtraction mode by performing a single iteration static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode pmSubtractionKernels *kernels, // Kernel description - const char *description // Description for trace + const char *description, // Description for trace + psImage *subMask, // Subtraction mask + float rej // Rejection threshold ) { @@ -939,5 +895,4 @@ if (!pmSubtractionSolveEquation(kernels, stamps)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - kernels->mode = oldMode; return false; } @@ -951,5 +906,5 @@ psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description); - long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint); + long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); @@ -960,9 +915,11 @@ if (numRejected > 0) { - psTrace("psModules.imcombine", 3, "Solving equation...\n"); + // Allow re-fit with reduced stamps set + psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); if (!pmSubtractionSolveEquation(kernels, stamps)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); return false; } + psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { @@ -970,4 +927,11 @@ return false; } + psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description); + long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); + if (numRejected < 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); + psFree(deviations); + return false; + } psFree(deviations); } @@ -977,25 +941,69 @@ -pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) -{ - PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, PM_SUBTRACTION_MODE_ERR); - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, PM_SUBTRACTION_MODE_ERR); - - // Copies of the inputs so we can try each way - pmSubtractionStampList *stamps1 = pmSubtractionStampsListCopy(stamps); - pmSubtractionStampList *stamps2 = pmSubtractionStampsListCopy(stamps); - - pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(kernels); - pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(kernels); - +pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels, + const psImage *subMask, float rej) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(*kernels, PM_SUBTRACTION_MODE_ERR); + + // Copies of the inputs + pmSubtractionStampList *stamps1 = pmSubtractionStampListCopy(*stamps); + pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(*kernels); + psImage *subMask1 = psImageCopy(NULL, subMask, subMask->type.type); kernels1->mode = PM_SUBTRACTION_MODE_1; + + if (!subtractionModeTest(stamps1, kernels1, "forward", subMask1, rej)) { + psError(PS_ERR_UNKNOWN, false, "Unable to test forward subtraction"); + psFree(stamps1); + psFree(kernels1); + psFree(subMask1); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(subMask1); + + // Copies of the inputs + pmSubtractionStampList *stamps2 = pmSubtractionStampListCopy(*stamps); + pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(*kernels); + psImage *subMask2 = psImageCopy(NULL, subMask, subMask->type.type); kernels2->mode = PM_SUBTRACTION_MODE_2; - - subtractionModeTest(stamps1, kernels1, "forward"); - subtractionModeTest(stamps2, kernels2, "backward"); - - // XXX Compare kernels1->mean, kernels2->mean -} - -#endif + if (!subtractionModeTest(stamps2, kernels2, "backward", subMask2, rej)) { + psError(PS_ERR_UNKNOWN, false, "Unable to test backward subtraction"); + psFree(stamps2); + psFree(kernels2); + psFree(subMask2); + psFree(stamps1); + psFree(kernels1); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(subMask2); + + + pmSubtractionStampList *bestStamps = NULL; // Best choice for stamps + pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels + psLogMsg("psModules.imcombine", PS_LOG_INFO, + "Forward: %f +/- %f from %d stamps\nBackward: %f +/- %f from %d stamps\n", + kernels1->mean, kernels1->rms, kernels1->numStamps, + kernels2->mean, kernels2->rms, kernels2->numStamps); + + if (kernels1->mean < kernels2->mean) { + bestStamps = stamps1; + bestKernels = kernels1; + } else { + bestStamps = stamps2; + bestKernels = kernels2; + } + + psFree(*stamps); + psFree(*kernels); + *stamps = psMemIncrRefCounter(bestStamps); + *kernels = psMemIncrRefCounter(bestKernels); + + psFree(stamps1); + psFree(stamps2); + psFree(kernels1); + psFree(kernels2); + + return bestKernels->mode; +} +