Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 14513)
+++ trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 14514)
@@ -9,9 +9,4 @@
 #include "pmSubtractionStamps.h"
 
-
-#define PM_SUBTRACTION_MASK_BAD_STAMP (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT | \
-                                       PM_SUBTRACTION_MASK_REJ) // Mask indicates a stamp would be bad
-
-
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Private (file-static) functions
@@ -24,6 +19,33 @@
     psFree(stamp->matrix);
     psFree(stamp->vector);
-}
-
+    psFree(stamp->reference);
+    psFree(stamp->input);
+    psFree(stamp->weight);
+}
+
+// Is this region OK?
+static bool checkStampRegion(int x, int y, // Coordinates of stamp
+                             const psRegion *region // Region of interest
+    )
+{
+    if (!region) {
+        return true;
+    }
+    return (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) ?
+        false : true;
+}
+
+// Is this position unmasked?
+static bool checkStampMask(int x, int y, // Coordinates of stamp
+                           const psImage *mask
+    )
+{
+    if (!mask) {
+        return true;
+    }
+    return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER |
+                                                  PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ?
+        false : true;
+}
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -38,4 +60,5 @@
     stamp->x = NAN;
     stamp->y = NAN;
+    stamp->flux = NAN;
     stamp->xNorm = NAN;
     stamp->yNorm = NAN;
@@ -135,7 +158,5 @@
                 for (int y = yStart; y <= yStop ; y++) {
                     for (int x = xStart; x <= xStop ; x++) {
-                        if ((!subMask || !(subMask->data.PS_TYPE_MASK_DATA[y][x] &
-                                           PM_SUBTRACTION_MASK_BAD_STAMP)) &&
-                            image->data.F32[y][x] > fluxBest) {
+                        if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxBest) {
                             fluxBest = image->data.F32[y][x];
                             xBest = x;
@@ -147,4 +168,5 @@
                 stamp->x = xBest;
                 stamp->y = yBest;
+                stamp->flux = fluxBest;
 
                 // Reset the postage stamps since we're making a new stamp
@@ -173,5 +195,5 @@
 
 
-psArray *pmSubtractionSetStamps(const psVector *x, const psVector *y,
+psArray *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux,
                                 const psImage *subMask, const psRegion *region)
 {
@@ -181,4 +203,9 @@
     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
     PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
+    if (flux) {
+        PS_ASSERT_VECTOR_NON_NULL(flux, false);
+        PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, false);
+    }
     if (subMask) {
         PS_ASSERT_IMAGE_NON_NULL(subMask, false);
@@ -191,12 +218,8 @@
     for (int i = 0; i < num; i++) {
         int xStamp = x->data.F32[i] + 0.5, yStamp = y->data.F32[i] + 0.5; // Pixel coords of stamp
-        if (region && (xStamp < region->x0 || xStamp > region->x1 ||
-                       yStamp < region->y0 || yStamp > region->y1)) {
-            // Not interested at this time
-            continue;
-        }
-
-        if (subMask && subMask->data.PS_TYPE_MASK_DATA[yStamp][xStamp] & PM_SUBTRACTION_MASK_BAD_STAMP) {
-            // It's bad
+        if (!checkStampRegion(xStamp, yStamp, region)) {
+            continue;
+        }
+        if (!checkStampMask(xStamp, yStamp, subMask)) {
             continue;
         }
@@ -205,4 +228,7 @@
         stamp->x = x->data.F32[i];
         stamp->y = y->data.F32[i];
+        if (flux) {
+            stamp->flux = flux->data.F32[i];
+        }
         stamps->data[stamps->n++] = stamp;
     }
@@ -278,40 +304,39 @@
 
 
-psArray *pmSubtractionGenerateStamps(const psVector *x, const psVector *y, const psVector *flux,
-                                     float sigma, int footprint, const pmSubtractionKernels *kernels)
-{
-    PS_ASSERT_VECTOR_NON_NULL(x, false);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
-    PS_ASSERT_VECTOR_NON_NULL(y, false);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
-    PS_ASSERT_VECTOR_NON_NULL(flux, false);
-    PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(flux, y, false);
+bool pmSubtractionGenerateStamps(psArray *stamps, float sigma, int footprint,
+                                 const pmSubtractionKernels *kernels)
+{
+    PS_ASSERT_ARRAY_NON_NULL(stamps, false);
     PS_ASSERT_FLOAT_LARGER_THAN(sigma, 0.0, false);
+    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
 
     int size = kernels->size + footprint; // Size of postage stamps
-    int num = x->n;                     // Number of stamps
-    psArray *stamps = psArrayAlloc(num); // Stamps, to return
+    int num = stamps->n;                // Number of stamps
 
     for (int i = 0; i < num; i++) {
-        pmSubtractionStamp *stamp = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_CALCULATE); // Stamp
-        stamps->data[i] = stamp;
-
-        stamp->x = x->data.F32[i];
-        stamp->y = y->data.F32[i];
-
-        float xStamp = x->data.F32[i] - (int)(stamp->x + 0.5); // x coordinate of star in stamp frame
-        float yStamp = y->data.F32[i] - (int)(stamp->y + 0.5); // y coordinate of star in stamp frame
-        float fluxStamp = flux->data.F32[i]; // Flux of star
-
+        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
+        if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) {
+            continue;
+        }
+
+        float x = stamp->x, y = stamp->y; // Coordinates of stamp
+        float flux = stamp->flux; // Flux of star
+        if (!isfinite(flux)) {
+            psWarning("Unable to generate PSF for stamp %d --- bad flux.", i);
+            continue;
+        }
+
+        float xStamp = x - (int)(x + 0.5); // x coordinate of star in stamp frame
+        float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame
+
+        psFree(stamp->input);
         stamp->input = psKernelAlloc(-size, size, -size, size);
         psKernel *input = stamp->input; // Target stamp
 
         // Put in a Gaussian, just for fun!
-        for (int y = -size; y <= size; y++) {
-            for (int x = -size; x <= size; x++) {
-                input->kernel[y][x] = fluxStamp / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 *
-                    expf(0.5 * (PS_SQR(x + xStamp) + PS_SQR(y + yStamp)) / PS_SQR(sigma));
+        for (int v = -size; v <= size; v++) {
+            for (int u = -size; u <= size; u++) {
+                input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 *
+                    expf(0.5 * (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / PS_SQR(sigma));
             }
         }
Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionStamps.h	(revision 14513)
+++ trunk/psModules/src/imcombine/pmSubtractionStamps.h	(revision 14514)
@@ -18,4 +18,5 @@
 typedef struct {
     float x, y;                         ///< Position
+    float flux;                         ///< Flux
     float xNorm, yNorm;                 ///< Normalised position
     psKernel *reference;                ///< Reference image postage stamp
@@ -39,4 +40,5 @@
 psArray *pmSubtractionSetStamps(const psVector *x, ///< x coordinates for each stamp
                                 const psVector *y, ///< y coordinates for each stamp
+                                const psVector *flux, ///< Flux for each stamp, or NULL
                                 const psImage *mask, ///< Mask, or NULL
                                 const psRegion *region ///< Region in which to search for stamps, or NULL
@@ -53,10 +55,8 @@
 
 /// Generate target for stamps based on list
-psArray *pmSubtractionGenerateStamps(const psVector *x, ///< x coordinates for each stamp
-                                     const psVector *y, ///< y coordinates for each stamp
-                                     const psVector *flux, ///< Flux for each stamp
-                                     float sigma, ///< Gaussian width for each stamp
-                                     int footprint, ///< Stamp footprint size
-                                     const pmSubtractionKernels *kernels ///< Kernel basis functions
+bool pmSubtractionGenerateStamps(psArray *stamps, ///< Stamps
+                                 float sigma, ///< Gaussian width for each stamp
+                                 int footprint, ///< Stamp footprint size
+                                 const pmSubtractionKernels *kernels ///< Kernel basis functions
     );