Index: trunk/psModules/src/imcombine/Makefile.am
===================================================================
--- trunk/psModules/src/imcombine/Makefile.am	(revision 14552)
+++ trunk/psModules/src/imcombine/Makefile.am	(revision 14580)
@@ -9,5 +9,6 @@
 	pmSubtraction.c		\
 	pmSubtractionKernels.c	\
-	pmSubtractionStamps.c
+	pmSubtractionStamps.c	\
+	pmSubtractionMatch.c
 
 pkginclude_HEADERS = \
@@ -16,5 +17,6 @@
 	pmSubtraction.h		\
 	pmSubtractionKernels.h	\
-	pmSubtractionStamps.h
+	pmSubtractionStamps.h	\
+	pmSubtractionMatch.h
 
 CLEANFILES = *~
Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 14580)
+++ trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 14580)
@@ -0,0 +1,381 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmSubtractionKernels.h"
+#include "pmSubtractionStamps.h"
+#include "pmSubtraction.h"
+#include "pmSubtractionMatch.h"
+
+
+//#define TESTING
+
+// Output memory usage information
+static void memCheck(const char *where)
+{
+#ifdef TESTING
+    psMemBlock **leaks = NULL;
+    int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
+    size_t largestSize = 0;
+    psMemId largest = 0;
+    size_t totalSize = 0;
+    for (int i = 0; i < numLeaks; i++) {
+        psMemBlock *mb = leaks[i];
+        totalSize += mb->userMemorySize;
+        if (mb->userMemorySize > largestSize) {
+            largestSize = mb->userMemorySize;
+            largest = mb->id;
+        }
+    }
+    psFree(leaks);
+    fprintf(stderr, "%s:\n", where);
+    fprintf(stderr, "    Memory in use: %zd\n", totalSize);
+    fprintf(stderr, "    Largest block: %ld\n", largest);
+    fprintf(stderr, "    sbrk(): %zd\n", (size_t)sbrk(0));
+#endif
+    return;
+}
+
+
+bool pmSubtractionMatch(pmReadout *convolved,
+                        const pmReadout *reference,
+                        const pmReadout *input,
+                        int footprint,
+                        float regionSize,
+                        float stampSpacing,
+                        float threshold,
+                        psString 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,
+                        psMaskType maskBlank
+    )
+{
+    PS_ASSERT_PTR_NON_NULL(convolved, false);
+    PS_ASSERT_PTR_NON_NULL(reference, false);
+    PS_ASSERT_IMAGE_NON_NULL(reference->image, false);
+    PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false);
+    if (reference->mask) {
+        PS_ASSERT_IMAGE_NON_NULL(reference->mask, false);
+        PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false);
+    }
+    if (reference->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(reference->weight, false);
+        PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false);
+    }
+    if (input) {
+        PS_ASSERT_IMAGE_NON_NULL(input->image, false);
+        PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false);
+        if (input->mask) {
+            PS_ASSERT_IMAGE_NON_NULL(input->mask, false);
+            PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false);
+        }
+        if (input->weight) {
+            PS_ASSERT_IMAGE_NON_NULL(input->weight, false);
+            PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, false);
+        }
+    }
+    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
+    // regionSize can be just about anything (except maybe negative, but it can be NAN)
+    PS_ASSERT_FLOAT_LARGER_THAN(stampSpacing, 0.0, false);
+    // Don't care what threshold is
+    // stampsName may be anything
+    // targetWidth can be just about anything (except maybe negative, but it can be NAN)
+    // We'll check kernel type when we allocate the kernels
+    PS_ASSERT_INT_POSITIVE(size, false);
+    PS_ASSERT_INT_NONNEGATIVE(order, false);
+    if (isisWidths || isisOrders) {
+        PS_ASSERT_VECTOR_NON_NULL(isisWidths, false);
+        PS_ASSERT_VECTOR_TYPE(isisWidths, PS_TYPE_F32, false);
+        PS_ASSERT_VECTOR_NON_NULL(isisOrders, false);
+        PS_ASSERT_VECTOR_TYPE(isisOrders, PS_TYPE_S32, false);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(isisWidths, isisOrders, false);
+    }
+    PS_ASSERT_INT_NONNEGATIVE(inner, false);
+    PS_ASSERT_INT_NONNEGATIVE(ringsOrder, false);
+    PS_ASSERT_INT_POSITIVE(binning, false);
+    PS_ASSERT_INT_POSITIVE(iter, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
+    // Don't care about maskBad
+    // Don't care about maskBlank
+
+    // Reset the output readout, just in case
+    if (convolved->image) {
+        psFree(convolved->image);
+        convolved->image = NULL;
+    }
+    if (convolved->mask) {
+        psFree(convolved->mask);
+        convolved->mask = NULL;
+    }
+    if (convolved->weight) {
+        psFree(convolved->weight);
+        convolved->weight = NULL;
+    }
+
+    psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight
+    if (input) {
+        inImage = input->image;
+        inMask = input->mask;
+    }
+
+    // Where does our weight map come from?
+    psImage *weight = NULL;             // Weight image to use
+    if (input && input->weight) {
+        weight = input->weight;
+    } else if (reference->weight) {
+        weight = reference->weight;
+    } else if (input) {
+        weight = input->image;
+    } else {
+        weight = reference->image;
+    }
+
+    int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
+
+    memCheck("start");
+
+    // Mask for subtraction
+    psImage *subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint);
+
+    memCheck("mask");
+
+    // Kernel basis functions
+    pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, isisWidths, isisOrders,
+                                                                 inner, binning, ringsOrder);
+    psMetadataAddStr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.DESCRIPTION", 0,
+                     "Subtraction kernel description", kernels->description);
+    psArray *stamps = NULL;             // Stamps for matching PSF
+    psVector *solution = NULL;          // Solution to match PSF
+
+
+    memCheck("kernels");
+
+    int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
+    float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
+    psRegion *region = NULL;            // Iso-kernel region
+    if (isfinite(regionSize) || regionSize == 0.0) {
+        xRegions = numCols / regionSize + 1;
+        yRegions = numRows / regionSize + 1;
+        xRegionSize = (float)numCols / (float)xRegions;
+        yRegionSize = (float)numRows / (float)yRegions;
+        region = psRegionAlloc(NAN, NAN, NAN, NAN);
+    }
+
+    // Iterate over iso-kernel regions
+    for (int j = 0; j < yRegions; j++) {
+        for (int i = 0; i < xRegions; i++) {
+            psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
+                    j * xRegions + i + 1, xRegions * yRegions);
+            if (region) {
+                *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
+                                      (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
+                psString string = psRegionToString(*region);
+                psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
+                        string, numCols, numRows);
+                psFree(string);
+            }
+
+            int numRejected = -1;               // Number of rejected stamps in each iteration
+            for (int k = 0; k < iter && numRejected != 0; k++) {
+                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.");
+                    goto ERROR;
+                }
+
+                memCheck("  find stamps");
+
+                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");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto ERROR;
+                }
+
+                memCheck("  calculate equation");
+
+                psTrace("psModules.imcombine", 3, "Solving equation...\n");
+                solution = pmSubtractionSolveEquation(solution, stamps);
+                if (!solution) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto ERROR;
+                }
+
+                memCheck("  solve equation");
+
+                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
+                numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels);
+                if (numRejected < 0) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+                    goto ERROR;
+                }
+                psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, k);
+
+                memCheck("  reject stamps");
+            }
+
+            if (numRejected > 0) {
+                psTrace("psModules.imcombine", 3, "Solving equation...\n");
+                solution = pmSubtractionSolveEquation(solution, stamps);
+                if (!solution) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto ERROR;
+                }
+            }
+            psFree(stamps);
+            stamps = NULL;
+
+            memCheck("solution");
+
+            {
+                psTrace("psModules.imcombine", 2, "Generating diagnostics...\n");
+                // Generate image with convolution kernels
+                int fullSize = 2 * size + 1 + 1;    // Full size of kernel
+                psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
+                psImageInit(convKernels, NAN);
+                for (int j = -2; j <= 2; j++) {
+                    for (int i = -2; i <= 2; i++) {
+                        psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
+                                                                   (float)j / 2.0); // Image of the kernel
+                        if (!kernel) {
+                            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+                            psFree(convKernels);
+                            goto ERROR;
+                        }
+
+                        if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
+                                                  (j + 2) * fullSize, "=") == 0) {
+                            psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
+                            psFree(kernel);
+                            psFree(convKernels);
+                            goto ERROR;
+                        }
+                        psFree(kernel);
+                    }
+                }
+
+                psString comment = NULL;// Comment for metadata
+                psString regionString = psRegionToString(*region); // String for region
+                psStringAppend(&comment, "Subtraction kernel for region %s", regionString);
+                psFree(regionString);
+                psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL",
+                                 PS_META_DUPLICATE_OK, comment, convKernels);
+                psFree(comment);
+                psFree(convKernels);
+            }
+
+#ifdef TESTING
+            {
+                // Generate images of the kernel components
+                psMetadata *header = psMetadataAlloc(); // Header
+                for (int i = 0; i < solution->n; i++) {
+                    psString name = NULL;       // Header keyword
+                    psStringAppend(&name, "SOLN%04d", i);
+                    psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
+                    psFree(name);
+                }
+                psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
+                psFits *kernelFile = psFitsOpen("kernels.fits", "w");
+                (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
+                psFitsClose(kernelFile);
+                psFree(kernelImages);
+                psFree(header);
+            }
+#endif
+
+            memCheck("diag outputs");
+
+            psTrace("psModules.imcombine", 2, "Convolving...\n");
+            if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask,
+                                       reference->image, reference->weight, subMask, maskBlank, region,
+                                       solution, kernels)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
+                goto ERROR;
+            }
+
+            psFree(solution);
+            solution = NULL;
+        }
+    }
+    psFree(kernels);
+
+    psFree(region);
+    region = NULL;
+    psFree(subMask);
+    subMask = NULL;
+
+    if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, kernels, maskBlank)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to set border region of convolved image.");
+        goto ERROR;
+    }
+
+    memCheck("convolution");
+
+
+    return true;
+
+ERROR:
+    psFree(region);
+    psFree(subMask);
+    psFree(kernels);
+    psFree(stamps);
+    psFree(solution);
+    return false;
+}
Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionMatch.h	(revision 14580)
+++ trunk/psModules/src/imcombine/pmSubtractionMatch.h	(revision 14580)
@@ -0,0 +1,37 @@
+#ifndef PM_SUBTRACTION_MATCH_H
+#define PM_SUBTRACTION_MATCH_H
+
+#include <pslib.h>
+
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmSubtractionKernels.h"
+
+
+bool pmSubtractionMatch(pmReadout *convolved, ///< Output convolved data
+                        const pmReadout *reference, ///< Reference data
+                        const pmReadout *input, ///< Input data
+                        // Stamp parameters
+                        int footprint,  ///< Stamp half-size
+                        float regionSize, ///< Typical size of iso-kernel regions
+                        float stampSpacing, ///< Typical spacing between stamps
+                        float threshold, ///< Threshold for stamps
+                        psString stampsName, ///< Filename for stamps
+                        float targetWidth, ///< Width of PSF for simulated target
+                        // Kernel parameters
+                        pmSubtractionKernelsType type, ///< Kernel type
+                        int size,       ///< Kernel half-size
+                        int order,      ///< Spatial polynomial order
+                        const psVector *widths, ///< ISIS Gaussian widths
+                        const psVector *orders, ///< ISIS Polynomial orders
+                        int inner,      ///< Inner radius for various kernel types
+                        int ringsOrder, ///< RINGS polynomial order
+                        int binning,    ///< SPAM kernel binning
+                        // Operational parameters
+                        int iter,       ///< Rejection iterations
+                        float rej,      ///< Rejection threshold
+                        psMaskType maskBad, ///< Value to mask
+                        psMaskType maskBlank ///< Mask for blank region
+    );
+
+#endif