Index: trunk/ppSub/src/ppSubReadout.c
===================================================================
--- trunk/ppSub/src/ppSubReadout.c	(revision 16916)
+++ trunk/ppSub/src/ppSubReadout.c	(revision 17298)
@@ -11,4 +11,6 @@
 
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
+#define TESTING                         // For test output
+
 
 bool ppSubReadout(pmConfig *config, const pmFPAview *view)
@@ -17,6 +19,8 @@
     pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
     pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
+    pmReadout *inConv = pmReadoutAlloc(NULL); // Convolved version of input
+    pmReadout *refConv = pmReadoutAlloc(NULL); // Convolved version of reference
     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
-    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout: subtraction
     pmFPA *outFPA = outCell->parent->parent; // Output FPA
     pmHDU *outHDU = outFPA->hdu; // Output HDU
@@ -63,8 +67,4 @@
 
     pmSubtractionMode mode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtraction mode
-    pmReadout *conv2 = NULL;            // Convolved image, for dual convolution
-    if (dual) {
-        conv2 = pmReadoutAlloc(NULL);
-    }
 
     // Generate masks if they don't exist
@@ -97,10 +97,11 @@
     }
 
-    if (!pmSubtractionMatch(outRO, conv2, inRO, refRO, footprint, regionSize, spacing, threshold, sources,
+    if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, regionSize, spacing, threshold, sources,
                             stampsName, type, size, order, widths, orders, inner, ringsOrder,
                             binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
                             maskBlank, badFrac, mode)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
-        psFree(conv2);
+        psFree(inConv);
+        psFree(refConv);
         psFree(outRO);
         return false;
@@ -109,40 +110,65 @@
 
     // Add kernel descrption to header
-    pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, outRO->analysis,
+    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis,
                                                         "SUBTRACTION.KERNEL"); // The subtraction kernels
+    if (!kernels) {
+        kernels = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL");
+    }
+    if (!kernels) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL");
+        psFree(inConv);
+        psFree(refConv);
+        psFree(outRO);
+        return false;
+    }
     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0,
                      "Subtraction kernel", kernels->description);
 
-    psImage *kernelImage = psMetadataLookupPtr(NULL, outRO->analysis,
+#ifdef TESTING
+    psImage *kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis,
                                                "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels
+    if (!kernelImage) {
+        kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE");
+    }
     psFits *fits = psFitsOpen("kernel.fits", "w");
     psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
     psFitsClose(fits);
+#endif
 
     // Do the subtraction
     {
         // Subtraction is: minuend - subtrahend
-        psImage *minuendImage = outRO->image;
-        psImage *subtrahendImage = (dual ? conv2->image : inRO->image);
-        psImage *minuendWeight = outRO->image;
-        psImage *subtrahendWeight = (dual ? conv2->weight : inRO->weight);
-        psImage *mask2 = (dual ? conv2->mask : inRO->mask);
+        pmReadout *minuend = inConv;
+        pmReadout *subtrahend = refConv;
 
         if (reverse) {
-            psImage *temp = subtrahendImage;
-            subtrahendImage = minuendImage;
-            minuendImage = temp;
-
-            temp = subtrahendWeight;
-            subtrahendWeight = minuendWeight;
-            minuendWeight = temp;
-        }
-
-        psBinaryOp(outRO->image, minuendImage, "-", subtrahendImage);
-        psBinaryOp(outRO->mask, outRO->mask, "|", mask2);
-        if (minuendWeight && subtrahendWeight) {
-            psBinaryOp(outRO->weight, minuendWeight, "+", subtrahendWeight);
-        }
-    }
+            pmReadout *temp = subtrahend;
+            subtrahend = minuend;
+            minuend = temp;
+        }
+
+        outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
+        outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
+        if (minuend->weight && subtrahend->weight) {
+            outRO->weight = (psImage*)psBinaryOp(outRO->weight, minuend->weight, "+", subtrahend->weight);
+        }
+        outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
+
+#ifdef TESTING
+        {
+            psFits *fits = psFitsOpen("minuend.fits", "w");
+            psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
+            psFitsClose(fits);
+        }
+        {
+            psFits *fits = psFitsOpen("subtrahend.fits", "w");
+            psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
+            psFitsClose(fits);
+        }
+#endif
+    }
+
+    psFree(inConv);
+    psFree(refConv);
 
 #ifdef TESTING