Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 14606)
+++ trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 14671)
@@ -9,4 +9,5 @@
 #include "pmHDU.h"
 #include "pmFPA.h"
+#include "pmSubtractionParams.h"
 #include "pmSubtractionKernels.h"
 #include "pmSubtractionStamps.h"
@@ -44,9 +45,77 @@
 
 
+static bool getStamps(psArray **stamps, // Stamps to read
+                      const psArray *stampsData, // Stamp data from a file
+                      const pmReadout *reference, // Reference readout
+                      const pmReadout *input, // Input readout, or NULL to generate fake stamps
+                      const psImage *subMask, // Mask for subtraction, or NULL
+                      psImage *weight,  // Weight map
+                      const psRegion *region, // Region of interest, or NULL
+                      float threshold,  // Threshold for stamp finding
+                      float stampSpacing, // Spacing between stamps
+                      float targetWidth,// Target width for fake stamps
+                      int size,         // Kernel half-size
+                      int footprint     // Convolution footprint for stamps
+    )
+{
+
+    psImage *inImage = NULL; // Input image
+    if (input) {
+        inImage = input->image;
+    }
+
+    if (stampsData) {
+        psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
+        if (input) {
+            // We have x, y because the target is provided by the input image
+            xStamp = stampsData->data[0];
+            yStamp = stampsData->data[1];
+        } else {
+            // We have x, y and flux in order to generate a target
+            xStamp = stampsData->data[0];
+            yStamp = stampsData->data[1];
+            fluxStamp = stampsData->data[2];
+        }
+
+        psFree(*stamps);
+        *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region);
+    } else {
+        psTrace("psModules.imcombine", 3, "Finding stamps...\n");
+        *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region,
+                                          threshold, stampSpacing);
+    }
+    if (!stamps) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
+        return false;
+    }
+
+    memCheck("  find stamps");
+
+    if (!input && !pmSubtractionGenerateStamps(*stamps, targetWidth, footprint, size)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
+    if (!pmSubtractionExtractStamps(*stamps, reference->image, inImage, weight, footprint, size)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
+        return false;
+    }
+
+    memCheck("   extract stamps");
+
+    return true;
+}
+
+
+
+
+
 bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input,
                         int footprint, float regionSize, float stampSpacing, float threshold,
                         const char *stampsName, float targetWidth, pmSubtractionKernelsType type,
-                        int size, int order, const psVector *isisWidths, const psVector *isisOrders,
-                        int inner, int ringsOrder, int binning, int iter, float rej, psMaskType maskBad,
+                        int size, int spatialOrder, const psVector *isisWidths, const psVector *isisOrders,
+                        int inner, int ringsOrder, int binning, bool optimum, psVector *optFWHMs,
+                        int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
                         psMaskType maskBlank)
 {
@@ -88,5 +157,5 @@
     // We'll check kernel type when we allocate the kernels
     PS_ASSERT_INT_POSITIVE(size, false);
-    PS_ASSERT_INT_NONNEGATIVE(order, false);
+    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, false);
     if (isisWidths || isisOrders) {
         PS_ASSERT_VECTOR_NON_NULL(isisWidths, false);
@@ -99,4 +168,10 @@
     PS_ASSERT_INT_NONNEGATIVE(ringsOrder, false);
     PS_ASSERT_INT_POSITIVE(binning, false);
+    if (optimum) {
+        PS_ASSERT_VECTOR_NON_NULL(optFWHMs, false);
+        PS_ASSERT_INT_NONNEGATIVE(optOrder, false);
+        PS_ASSERT_FLOAT_LARGER_THAN(optThreshold, 0.0, false);
+        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false);
+    }
     PS_ASSERT_INT_POSITIVE(iter, false);
     PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
@@ -136,4 +211,24 @@
     }
 
+    // Read stamps from file
+    psArray *stampsData = NULL;         // Stamps data read from file
+    if (stampsName && strlen(stampsName) > 0) {
+        psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
+        psVector *xStamp = NULL, *yStamp = NULL; // Stamp positions and fluxes
+        if (input) {
+            // We have x, y because the target is provided by the input image
+            stampsData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions
+        } else {
+            // We have x, y and flux in order to generate a target
+            stampsData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions
+        }
+
+        // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
+        xStamp = stampsData->data[0];
+        yStamp = stampsData->data[1];
+        psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+        psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+    }
+
     int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
 
@@ -145,10 +240,36 @@
     memCheck("mask");
 
-    // Kernel basis functions
-    pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, isisWidths, isisOrders,
-                                                                 inner, binning, ringsOrder);
+
+    psArray *stamps = NULL;             // Stamps for matching PSF
+
+    pmSubtractionKernels *kernels = NULL; // Kernel basis functions
+    if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
+        if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, NULL,
+                       threshold, stampSpacing, targetWidth, size, footprint)) {
+            goto ERROR;
+        }
+        kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
+                                                  stamps, footprint, optThreshold);
+        // XXX This is not quite the best thing to do here--- we'll find the same set of stamps again later,
+        // but since we may not be interested in the entire image on the first pass through, we have to blow
+        // the whole lot away.  The alternative is to mark stamps that are outside the region of interest, and
+        // ignore them when generating sums and rejecting.
+        psFree(stamps);
+        stamps = NULL;
+
+        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);
+    }
+
     psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", PS_DATA_UNKNOWN,
                      "Subtraction kernels", kernels);
-    psArray *stamps = NULL;             // Stamps for matching PSF
     psVector *solution = NULL;          // Solution to match PSF
 
@@ -186,51 +307,8 @@
                 psTrace("psModules.imcombine", 2, "Iteration %d...\n", k);
 
-                if (stampsName && strlen(stampsName) > 0) {
-                    psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
-                    iter = 1;           // There is no iterating because we use all the stamps we have
-                    psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
-                    if (input) {
-                        // We have x, y because the target is provided by the input image
-                        psArray *stampData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions
-                        xStamp = stampData->data[0];
-                        yStamp = stampData->data[1];
-                    } else {
-                        // We have x, y and flux in order to generate a target
-                        psArray *stampData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions
-                        xStamp = stampData->data[0];
-                        yStamp = stampData->data[1];
-                        fluxStamp = stampData->data[2];
-                    }
-
-                    // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
-                    psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-                    psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-
-                    stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region);
-                } else {
-                    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
-                    stamps = pmSubtractionFindStamps(stamps, reference->image, subMask, region,
-                                                     threshold, stampSpacing);
-                }
-                if (!stamps) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
+                if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, region,
+                               threshold, stampSpacing, targetWidth, size, footprint)) {
                     goto ERROR;
                 }
-
-                memCheck("  find stamps");
-
-                if (!input && !pmSubtractionGenerateStamps(stamps, targetWidth, footprint, kernels)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
-                    goto ERROR;
-                }
-
-                psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-                if (!pmSubtractionExtractStamps(stamps, reference->image, inImage, weight,
-                                                footprint, kernels)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
-                    goto ERROR;
-                }
-
-                memCheck("   extract stamps");
 
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
@@ -376,4 +454,6 @@
     psFree(subMask);
     subMask = NULL;
+    psFree(stampsData);
+    stampsData = NULL;
 
     if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, kernels, maskBlank)) {
@@ -395,4 +475,5 @@
     psFree(stamps);
     psFree(solution);
+    psFree(stampsData);
     return false;
 }