Index: trunk/ppSub/src/ppSubReadout.c
===================================================================
--- trunk/ppSub/src/ppSubReadout.c	(revision 13363)
+++ trunk/ppSub/src/ppSubReadout.c	(revision 13371)
@@ -16,6 +16,8 @@
     pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
     pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
+#if 0
     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
     pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+#endif
 
     psImage *input = inRO->image;       // Input image
@@ -87,7 +89,8 @@
     psVector *solution = NULL;          // Solution to match PSF
 
-    for (int i = 0, numRejected = 0; i < iter && numRejected != 0; i++) {
+    int numRejected = -1;               // Number of rejected stamps in each iteration
+    for (int i = 0; i < iter && numRejected != 0; i++) {
         stamps = pmSubtractionFindStamps(stamps, refRO->image, stampMask, MASK_BAD, MASK_STAMP,
-                                         threshold, spacing, footprint);
+                                         threshold, spacing, size + footprint);
         if (!stamps) {
             psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image.");
@@ -116,5 +119,6 @@
     }
 
-    if (!pmSubtractionConvolve(&outRO->image, &outRO->weight, &outRO->mask,
+    psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images
+    if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask,
                                refRO->image, refRO->weight, refRO->mask,
                                MASK_BAD, maskBlank, solution, kernels)) {
@@ -125,15 +129,20 @@
     // Do the subtraction
     if (reverse) {
-        (void)psBinaryOp(outRO->image, outRO->image, "-", inRO->image);
+        (void)psBinaryOp(inRO->image, convImage, "-", inRO->image);
     } else {
-        (void)psBinaryOp(outRO->image, inRO->image, "-", outRO->image);
+        (void)psBinaryOp(inRO->image, inRO->image, "-", convImage);
     }
-    (void)psBinaryOp(outRO->mask, outRO->mask, "|", inRO->mask);
-    (void)psBinaryOp(outRO->weight, outRO->weight, "+", inRO->weight);
+    (void)psBinaryOp(inRO->mask, convMask, "|", inRO->mask);
+    (void)psBinaryOp(inRO->weight, convWeight, "+", inRO->weight);
+
+    psFree(convImage);
+    psFree(convMask);
+    psFree(convWeight);
 
 #if 0
     // Generate image with convolution kernels
-    int fullSize = 2 * size + 1;        // Full size of kernel
+    int fullSize = 2 * size + 1 + 1;    // Full size of kernel
     psImage *convKernels = psImageAlloc(5 * fullSize, 5 * fullSize, PS_TYPE_F32);
+    psImageInit(fullSize, NAN);
     for (int j = -2; j <= 2; j++) {
         for (int i = -2; i <= 2; i++) {
@@ -156,4 +165,9 @@
 
     // XXX What do we do with this image?
+
+    psFits *kernelFile = psFitsOpen("kernel.fits", "w");
+    (void)psFitsWriteImage(kernelFile, NULL, convKernels, 0, NULL);
+    psFitsClose(kernelFile);
+
     psFree(convKernels);
 #endif
@@ -163,5 +177,4 @@
         psFree(stamps);
         psFree(solution);
-        psFree(outRO);
         return true;
 
@@ -170,5 +183,4 @@
         psFree(stamps);
         psFree(solution);
-        psFree(outRO);
         return false;
 }