Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 14480)
+++ trunk/psModules/src/imcombine/pmSubtractionStamps.c	(revision 14513)
@@ -9,4 +9,9 @@
 #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
@@ -31,8 +36,8 @@
     psMemSetDeallocator(stamp, (psFreeFunc)subtractionStampFree);
 
-    stamp->x = 0;
-    stamp->y = 0;
-    stamp->xNorm = 0.0;
-    stamp->yNorm = 0.0;
+    stamp->x = NAN;
+    stamp->y = NAN;
+    stamp->xNorm = NAN;
+    stamp->yNorm = NAN;
     stamp->matrix = NULL;
     stamp->vector = NULL;
@@ -131,6 +136,5 @@
                     for (int x = xStart; x <= xStop ; x++) {
                         if ((!subMask || !(subMask->data.PS_TYPE_MASK_DATA[y][x] &
-                                           (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT |
-                                            PM_SUBTRACTION_MASK_REJ))) &&
+                                           PM_SUBTRACTION_MASK_BAD_STAMP)) &&
                             image->data.F32[y][x] > fluxBest) {
                             fluxBest = image->data.F32[y][x];
@@ -143,6 +147,4 @@
                 stamp->x = xBest;
                 stamp->y = yBest;
-                stamp->xNorm = 2.0 * (float)(xBest - numCols/2.0) / (float)numCols;
-                stamp->yNorm = 2.0 * (float)(yBest - numRows/2.0) / (float)numRows;
 
                 // Reset the postage stamps since we're making a new stamp
@@ -156,5 +158,5 @@
                     numFound++;
                     psTrace("psModules.imcombine", 5, "Found stamp in region (%d,%d): %d,%d\n",
-                            i, j, stamp->x, stamp->y);
+                            i, j, (int)stamp->x, (int)stamp->y);
                 } else {
                     stamp->status = PM_SUBTRACTION_STAMP_NONE;
@@ -166,4 +168,43 @@
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps in %dx%d regions",
              numFound, xNumStamps, yNumStamps);
+
+    return stamps;
+}
+
+
+psArray *pmSubtractionSetStamps(const psVector *x, const psVector *y,
+                                const psImage *subMask, const psRegion *region)
+{
+    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);
+    if (subMask) {
+        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
+        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
+    }
+
+    int num = x->n;                     // Number of stamps
+    psArray *stamps = psArrayAllocEmpty(num); // Stamps, to return
+
+    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
+            continue;
+        }
+
+        pmSubtractionStamp *stamp = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_CALCULATE); // Stamp
+        stamp->x = x->data.F32[i];
+        stamp->y = y->data.F32[i];
+        stamps->data[stamps->n++] = stamp;
+    }
 
     return stamps;
@@ -198,7 +239,16 @@
     for (int i = 0; i < stamps->n; i++) {
         pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
-
-
-        int x = stamp->x, y = stamp->y; // Stamp coordinates
+        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
+            continue;
+        }
+
+        if (isnan(stamp->xNorm)) {
+            stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
+        }
+        if (isnan(stamp->yNorm)) {
+            stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
+        }
+
+        int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
         if (x < size || x > numCols - size || y < size || y > numRows - size) {
             psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
@@ -228,2 +278,45 @@
 
 
+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);
+    PS_ASSERT_FLOAT_LARGER_THAN(sigma, 0.0, false);
+
+    int size = kernels->size + footprint; // Size of postage stamps
+    int num = x->n;                     // Number of stamps
+    psArray *stamps = psArrayAlloc(num); // Stamps, to return
+
+    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
+
+        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));
+            }
+        }
+
+    }
+
+    return stamps;
+}