Changeset 25407 for branches/eam_branches/20090715/psModules/src/imcombine

Timestamp:

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

Author:

eugene

Message:

updates from trunk

Location:

branches/eam_branches/20090715/psModules/src/imcombine

Files:

• pmPSFEnvelope.c (modified) (2 diffs)
• pmStack.c (modified) (1 diff)
• pmSubtraction.c (modified) (3 diffs)
• pmSubtraction.h (modified) (2 diffs)
• pmSubtractionAnalysis.c (modified) (9 diffs)
• pmSubtractionAnalysis.h (modified) (1 diff)
• pmSubtractionKernels.c (modified) (1 diff)
• pmSubtractionKernels.h (modified) (1 diff)
• pmSubtractionMatch.c (modified) (17 diffs)
• pmSubtractionMatch.h (modified) (1 diff)
• pmSubtractionStamps.c (modified) (1 diff)
• pmSubtractionStamps.h (modified) (1 diff)

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

@@ -124 +124 @@
         pmResiduals *resid = psf->residuals;// PSF residuals
         psf->residuals = NULL;
-        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, xOffset, yOffset, psf,
+        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf,
                                       NAN, radius, true, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout.");
@@ -298 +298 @@
         }
 
-        // measure the source moments: tophat windowing, no pixel S/N cutoff
+        // measure the source moments: tophat windowing, no pixel S/N cutoff
         if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) {
             // Can't do anything about it; limp along as best we can

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

@@ -30 +30 @@
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
 #define PIXEL_MAP_BUFFER 2              // Number of entries to add to pixel map at a time
-#define ADD_VARIANCE                    // Allow additional variance (besides variance factor)?
+//#define ADD_VARIANCE                    // Allow additional variance (besides variance factor)?
 #define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 

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

@@ -733 +733 @@
 
 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
-                              const psVector *deviations, psImage *subMask, float sigmaRej, int footprint)
+                              const psVector *deviations, psImage *subMask, float sigmaRej)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
@@ -821 +821 @@
     ds9num++;
 
+    int footprint = stamps->footprint;  // Half-size of stamp region of interest
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
@@ -956 +957 @@
     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
 
-    psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return
-    psVector *fakeSolution = psVectorAlloc(solution->n, PS_TYPE_F64); // Fake solution vector
+    psArray *images = psArrayAlloc(kernels->solution1->n - 1); // Images of each kernel to return
+    psVector *fakeSolution = psVectorAlloc(kernels->solution1->n, PS_TYPE_F64); // Fake solution vector
     psVectorInit(fakeSolution, 0.0);
 
-    for (int i = 0; i < solution->n - 1; i++) {
-        fakeSolution->data.F64[i] = solution->data.F64[i];
+    for (int i = 0; i < kernels->solution1->n - 1; i++) {
+        fakeSolution->data.F64[i] = kernels->solution1->data.F64[i];
         images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual);
         fakeSolution->data.F64[i] = 0.0;

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.h

@@ -68 +68 @@
                               const psVector *deviations, ///< Deviations for each stamp
                               psImage *subMask, ///< Subtraction mask
-                              float sigmaRej, ///< Number of RMS deviations above zero at which to reject
-                              int footprint ///< Half-size of stamp
+                              float sigmaRej ///< Number of RMS deviations above zero at which to reject
     );
 
@@ -94 +93 @@
 
 /// Generate images of the convolution kernel elements
-psArray *pmSubtractionKernelSolutions(const psVector *solution, ///< Solution vector
-                                      const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                      float x, float y ///< Normalised position [-1,1] for images
+psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                      float x, float y, ///< Normalised position [-1,1] for images
+                                      bool wantDual ///< Calculate for the dual kernel?
     );
+
 
 /// Execute a thread job to convolve a patch of the image

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

@@ -17 +17 @@
 
 
-bool pmSubtractionAnalysis(psMetadata *analysis, pmSubtractionKernels *kernels, psRegion *region,
+bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header,
+                           pmSubtractionKernels *kernels, psRegion *region,
                            int numCols, int numRows)
 {
-    if (analysis) {
-        PS_ASSERT_METADATA_NON_NULL(analysis, false);
-    }
+    PS_ASSERT_METADATA_NON_NULL(analysis, false);
+    PS_ASSERT_METADATA_NON_NULL(header, false);
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false);
@@ -38 +38 @@
                          PS_DATA_REGION | PS_META_DUPLICATE_OK,
                          "Region over which subtraction was performed", subRegion);
+
+        psString string = psRegionToString(*subRegion);
         psFree(subRegion);
+
+        psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK,
+                         "Region over which subtraction was performed", string);
+        psFree(string);
     }
 
@@ -45 +51 @@
                      PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
     psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
-                     PS_META_DUPLICATE_OK, "Subtraction kernels", kernels->mode);
+                     PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
+    psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
+                     PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
 
     // Realisations of kernel
@@ -113 +121 @@
     {
         psMetadata *header = psMetadataAlloc(); // Header
-        for (int i = 0; i < solution->n; i++) {
+        for (int i = 0; i < kernels->solution1->n; i++) {
             psString name = NULL;       // Header keyword
             psStringAppend(&name, "SOLN%04d", i);
-            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
+            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, kernels->solution1->data.F64[i]);
             psFree(name);
         }
-        psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
+        psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, false);
         psFits *kernelFile = psFitsOpen("kernels.fits", "w");
         (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
@@ -128 +136 @@
 #endif
 
-
-    // Set the variance factors
-    float vf1 = 1.0, vf2 = 1.0;         // Variance factors for each image
-    switch (kernels->mode) {
-      case PM_SUBTRACTION_MODE_1:
-        vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
-        break;
-      case PM_SUBTRACTION_MODE_2:
-        vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
-        break;
-      case PM_SUBTRACTION_MODE_DUAL:
-        vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
-        vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, true);
-        break;
-      default:
-        psAbort("Invalid subtraction mode: %x", kernels->mode);
-    }
-
-    // Weight by the area
-    if (region) {
-        float norm = (region->x1 - region->x0 + 1) * (region->y1 - region->y0 + 1) / (numCols * numRows);
-        vf1 *= norm;
-        vf2 *= norm;
-    }
-
-    // Update the variance factor
-#define UPDATE_VARFACTOR(VF, ANALYSIS) { \
-    psMetadataItem *vfItem = psMetadataLookup(analysis, ANALYSIS); \
-    if (vfItem) { \
-        psAssert(vfItem->type == PS_TYPE_F32, "Should be the type we said."); \
-        vfItem->data.F32 += VF; \
-    } else { \
-        psMetadataAddF32(analysis, PS_LIST_TAIL, ANALYSIS, 0, "Variance factor weighted by the area", VF); \
-    } \
-}
-
-    UPDATE_VARFACTOR(vf1, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
-    UPDATE_VARFACTOR(vf2, PM_SUBTRACTION_ANALYSIS_VARFACTOR_2);
 
     // Kernel shape
@@ -201 +171 @@
         psFree(image);
 
-        psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
-        if (item) {
-            item->data.F32 = PS_MAX(item->data.F32, max);
-        } else {
-            psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
-                             "Maximum deconvolution fraction", max);
+        {
+            psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
+            if (item) {
+                max = item->data.F32 = PS_MAX(item->data.F32, max);
+            } else {
+                psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
+                                 "Maximum deconvolution fraction", max);
+            }
+        }
+
+        {
+            psMetadataItem *item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
+            if (item) {
+                item->data.F32 = max;
+            } else {
+                psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
+                                 "Maximum deconvolution fraction", max);
+            }
         }
     }
@@ -254 +236 @@
         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY,
                          PS_META_DUPLICATE_OK, "Moment in yy", m02);
+
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM,
+                         PS_META_DUPLICATE_OK, "Normalisation", m00);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX,
+                         PS_META_DUPLICATE_OK, "Moment in x", m10);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MY,
+                         PS_META_DUPLICATE_OK, "Moment in y", m01);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXX,
+                         PS_META_DUPLICATE_OK, "Moment in xx", m20);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXY,
+                         PS_META_DUPLICATE_OK, "Moment in xy", m11);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY,
+                         PS_META_DUPLICATE_OK, "Moment in yy", m02);
     }
 
@@ -263 +258 @@
 
         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF,
+                         PS_META_DUPLICATE_OK, "Background difference", bg);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF,
                          PS_META_DUPLICATE_OK, "Background difference", bg);
         psFree(polyValues);
@@ -275 +272 @@
         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
                          kernels->rms);
+
+        psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps",
+                         kernels->numStamps);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation",
+                         kernels->mean);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
+                         kernels->rms);
     }
 

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.h

@@ -27 +27 @@
 bool pmSubtractionAnalysis(
     psMetadata *analysis,               ///< Metadata container for QA information
+    psMetadata *header,                 ///< Metadata container for QA information to put in header
     pmSubtractionKernels *kernels,      ///< Kernels
     psRegion *region,                   ///< Region for subtraction

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

@@ -765 +765 @@
     return PM_SUBTRACTION_KERNEL_NONE;
 }
+
+pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in)
+{
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(in, NULL);
+
+    pmSubtractionKernels *out = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
+    psMemSetDeallocator(out, (psFreeFunc)subtractionKernelsFree);
+
+    out->type = in->type;
+    out->description = psMemIncrRefCounter(in->description);
+    out->num = in->num;
+    out->u = psMemIncrRefCounter(in->u);
+    out->v = psMemIncrRefCounter(in->v);
+    out->widths = psMemIncrRefCounter(in->widths);
+    out->preCalc = psMemIncrRefCounter(in->preCalc);
+    out->penalty = in->penalty;
+    out->penalties = psMemIncrRefCounter(in->penalties);
+    out->uStop = psMemIncrRefCounter(in->uStop);
+    out->vStop = psMemIncrRefCounter(in->vStop);
+    out->size = in->size;
+    out->inner = in->inner;
+    out->spatialOrder = in->spatialOrder;
+    out->bgOrder = in->bgOrder;
+    out->mode = in->mode;
+    out->numCols = in->numCols;
+    out->numRows = in->numRows;
+    out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
+    out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
+
+    return out;
+}

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.h

@@ -203 +203 @@
     );
 
+/// Copy kernels
+///
+/// A deep copy is performed on the solution only; the other components are merely pointers.
+pmSubtractionKernels *pmSubtractionKernelsCopy(
+    const pmSubtractionKernels *in      // Kernels to copy
+    );
+
 
 #endif

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;
+}
+

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

@@ -83 +83 @@
     );
 
+/// Determine best subtraction mode to use
+///
+/// Subtractions are attempted each way, and the mode with the lower residual is taken to be the best
+pmSubtractionMode pmSubtractionBestMode(
+    pmSubtractionStampList **stamps,    ///< Stamps to use for solution
+    pmSubtractionKernels **kernels,     ///< Kernels to use for solution
+    const psImage *subMask,             ///< Subtraction mask
+    float rej                           ///< Rejection threshold for stamps
+    );
 
 #endif

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

@@ -225 +225 @@
 }
 
+pmSubtractionStampList *pmSubtractionStampListCopy(const pmSubtractionStampList *in)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(in, NULL);
+
+    pmSubtractionStampList *out = psAlloc(sizeof(pmSubtractionStampList)); // Copied stamp list to return
+    psMemSetDeallocator(out, (psFreeFunc)subtractionStampListFree);
+
+    int num = out->num = in->num;       // Number of stamps
+    out->stamps = psArrayAlloc(num);
+    out->regions = psArrayAlloc(num);
+    out->x = NULL;
+    out->y = NULL;
+    out->flux = NULL;
+    out->footprint = in->footprint;
+
+    for (int i = 0; i < num; i++) {
+        psRegion *inRegion = in->regions->data[i]; // Input region
+        out->regions->data[i] = psRegionAlloc(inRegion->x0, inRegion->x1, inRegion->y0, inRegion->y1);
+
+        pmSubtractionStamp *inStamp = in->stamps->data[i]; // Input stamp
+        pmSubtractionStamp *outStamp = psAlloc(sizeof(pmSubtractionStamp));
+        psMemSetDeallocator(outStamp, (psFreeFunc)subtractionStampFree);
+        outStamp->x = inStamp->x;
+        outStamp->y = inStamp->y;
+        outStamp->flux = inStamp->flux;
+        outStamp->xNorm = inStamp->xNorm;
+        outStamp->yNorm = inStamp->yNorm;
+        outStamp->status = inStamp->status;
+
+        outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL;
+        outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL;
+        outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL;
+
+        if (inStamp->convolutions1) {
+            int size = inStamp->convolutions1->n; // Size of array
+            outStamp->convolutions1 = psArrayAlloc(size);
+            for (int j = 0; j < size; j++) {
+                psKernel *conv = inStamp->convolutions1->data[j]; // Convolution
+                outStamp->convolutions1->data[j] = conv ? psKernelCopy(conv) : NULL;
+            }
+        } else {
+            outStamp->convolutions1 = NULL;
+        }
+        if (inStamp->convolutions2) {
+            int size = inStamp->convolutions2->n; // Size of array
+            outStamp->convolutions2 = psArrayAlloc(size);
+            for (int j = 0; j < size; j++) {
+                psKernel *conv = inStamp->convolutions2->data[j]; // Convolution
+                outStamp->convolutions2->data[j] = conv ? psKernelCopy(conv) : NULL;
+            }
+        } else {
+            outStamp->convolutions2 = NULL;
+        }
+
+        outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL;
+        outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL;
+        outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL;
+        outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL;
+        outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL;
+
+        out->stamps->data[i] = outStamp;
+    }
+
+    return out;
+}
+
 pmSubtractionStamp *pmSubtractionStampAlloc(void)
 {

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.h

@@ -52 +52 @@
     } \
 }
+
+/// Copy a list of stamps
+///
+/// A deep copy is performed of the stamp list and the component stamps
+pmSubtractionStampList *pmSubtractionStampListCopy(
+    const pmSubtractionStampList *in    // Stamp list to copy
+    );
+
 
 /// A stamp for image subtraction