Changeset 25407 for branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Sep 15, 2009, 4:03:13 PM (17 years ago)

Author:

eugene

Message:

updates from trunk

File:

• branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c (modified) (17 diffs)

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c

@@ -10 +10 @@
 #include "pmHDU.h"
 #include "pmFPA.h"
+#include "pmHDUUtils.h"
 #include "pmSubtractionParams.h"
 #include "pmSubtractionKernels.h"
@@ -58 +59 @@
 
 
-static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read
-                      const pmReadout *ro1, // Readout 1
-                      const pmReadout *ro2, // Readout 2
-                      const psImage *subMask, // Mask for subtraction, or NULL
-                      psImage *variance,  // Variance map
-                      const psRegion *region, // Region of interest, or NULL
-                      float thresh1,  // Threshold for stamp finding on readout 1
-                      float thresh2,  // Threshold for stamp finding on readout 2
-                      float stampSpacing, // Spacing between stamps
-                      int size,         // Kernel half-size
-                      int footprint,     // Convolution footprint for stamps
-                      pmSubtractionMode mode // Mode for subtraction
+static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read
+                                 const pmReadout *ro1, // Readout 1
+                                 const pmReadout *ro2, // Readout 2
+                                 const psImage *subMask, // Mask for subtraction, or NULL
+                                 psImage *variance,  // Variance map
+                                 const psRegion *region, // Region of interest, or NULL
+                                 float thresh1,  // Threshold for stamp finding on readout 1
+                                 float thresh2,  // Threshold for stamp finding on readout 2
+                                 float stampSpacing, // Spacing between stamps
+                                 int size,         // Kernel half-size
+                                 int footprint,     // Convolution footprint for stamps
+                                 pmSubtractionMode mode // Mode for subtraction
     )
 {
@@ -163 +164 @@
 }
 
+static void subtractionAnalysisUpdate(pmReadout *conv1, pmReadout *conv2, // Convolved images
+                                      const psMetadata *analysis, // Analysis metadata
+                                      const psMetadata *header // Header metadata
+    )
+{
+    if (conv1) {
+        conv1->analysis = psMetadataCopy(conv1->analysis, analysis);
+    }
+    if (conv2) {
+        conv2->analysis = psMetadataCopy(conv2->analysis, analysis);
+    }
+
+    if (conv1 && conv1->parent) {
+        pmHDU *hdu = pmHDUFromCell(conv1->parent);
+        if (hdu) {
+            hdu->header = psMetadataCopy(hdu->header, header);
+        }
+    }
+    if (conv2 && conv2->parent) {
+        pmHDU *hdu = pmHDUFromCell(conv2->parent);
+        if (hdu) {
+            hdu->header = psMetadataCopy(hdu->header, header);
+        }
+    }
+
+    return;
+}
 
 
@@ -253 +281 @@
 
     psMetadata *outAnalysis = psMetadataAlloc(); // Output analysis values
+    psMetadata *outHeader = psMetadataAlloc(); // Output header values
 
     psTrace("psModules.imcombine", 2, "Convolving...\n");
@@ -259 +288 @@
         psRegion *region = regions->data[i]; // Region of interest
 
-        if (!pmSubtractionAnalysis(outAnalysis, kernel, region, ro1->image->numCols, ro1->image->numRows)) {
+        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region,
+                                   ro1->image->numCols, ro1->image->numRows)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
             psFree(outAnalysis);
+            psFree(outHeader);
             psFree(subMask);
             psFree(kernels);
@@ -272 +303 @@
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
             psFree(outAnalysis);
+            psFree(outHeader);
             psFree(subMask);
             psFree(kernels);
@@ -283 +315 @@
     psFree(regions);
 
-    if (conv1) {
-        psMetadataCopy(conv1->analysis, outAnalysis);
-    }
-    if (conv2) {
-        psMetadataCopy(conv2->analysis, outAnalysis);
-    }
+    subtractionAnalysisUpdate(conv1, conv2, outAnalysis, outHeader);
     psFree(outAnalysis);
+    psFree(outHeader);
 
     return true;
@@ -369 +397 @@
     pmSubtractionKernels *kernels = NULL; // Kernel basis functions
     psMetadata *analysis = psMetadataAlloc(); // QA data
+    psMetadata *header = psMetadataAlloc(); // QA data for header
 
     int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
@@ -442 +471 @@
             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
             // doesn't matter.
-            if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
+            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
                            stampSpacing, size, footprint, subMode)) {
                 goto MATCH_ERROR;
             }
 
+
+            // 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 +517 @@
 
                 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 +533 @@
             }
 
-            // 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++) {
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-                if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
-                               stampSpacing, size, footprint, subMode)) {
+                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
+                                          stampThresh1, stampThresh2, stampSpacing,
+                                          size, footprint, subMode)) {
                     goto MATCH_ERROR;
                 }
@@ -535 +569 @@
 
                 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 +591 @@
                     goto MATCH_ERROR;
                 }
-                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint);
+                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
                 psFree(deviations);
             }
@@ -565 +599 @@
             memCheck("solution");
 
-            if (!pmSubtractionAnalysis(analysis, kernels, region, numCols, numRows)) {
+            if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
                 goto MATCH_ERROR;
@@ -601 +635 @@
     memCheck("convolution");
 
-    if (conv1) {
-        psMetadataCopy(conv1->analysis, analysis);
-    }
-    if (conv2) {
-        psMetadataCopy(conv2->analysis, analysis);
-    }
+    subtractionAnalysisUpdate(conv1, conv2, analysis, header);
     psFree(analysis);
+    psFree(header);
 
 #ifdef TESTING
@@ -629 +659 @@
 MATCH_ERROR:
     psFree(analysis);
+    psFree(header);
     psFree(region);
     psFree(regionString);
@@ -842 +873 @@
     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, "convolve 1", subMask1, rej)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 1");
+        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, "convolve 2", subMask2, rej)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 2");
+        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,
+             "Image 1: %f +/- %f from %d stamps\nImage 2: %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;
+}
+