Changeset 25148 for branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Aug 20, 2009, 11:47:38 AM (17 years ago)

Author:

watersc1

Message:

 

File:

• branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.c (modified) (6 diffs)

branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.c

@@ -447 +447 @@
             }
 
+
+            // 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
@@ -469 +488 @@
 
                 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:
@@ -483 +504 @@
             }
 
-            // 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++) {
@@ -535 +539 @@
 
                 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.");
@@ -557 +561 @@
                     goto MATCH_ERROR;
                 }
-                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint);
+                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
                 psFree(deviations);
             }
@@ -842 +846 @@
     return mode;
 }
+
+
+// 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
+                                psImage *subMask,  // Subtraction mask
+                                float rej               // Rejection threshold
+                                )
+{
+    assert(stamps);
+    assert(kernels);
+
+    psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Solving %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, "Calculate %s deviations...\n", description);
+    psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
+    if (!deviations) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
+    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
+    if (numRejected < 0) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+        psFree(deviations);
+        return false;
+    }
+    psFree(deviations);
+
+    if (numRejected > 0) {
+        // 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) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+            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);
+    }
+
+    return true;
+}
+
+
+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;
+
+    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;
+}
+