Index: trunk/ppStack/src/ppStackMatch.c
===================================================================
--- trunk/ppStack/src/ppStackMatch.c	(revision 19231)
+++ trunk/ppStack/src/ppStackMatch.c	(revision 19235)
@@ -15,5 +15,5 @@
                      PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
 
-//#define TESTING                         // Enable debugging output
+#define TESTING                         // Enable debugging output
 
 
@@ -83,7 +83,10 @@
     }
 
+    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
+
 #ifdef TESTING
     // Read previously produced kernel
-    static int numInput = 0;            // Index of input file
+    static int numInput = -1;            // Index of input file
+    numInput++;
     if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
         const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
@@ -97,8 +100,7 @@
             psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
             psFree(resolved);
-            if (!fits || !pmReadoutReadSubtractionKernels(readout, fits)) {
+            if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
                 psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
                 psFitsClose(fits);
-                numInput++;
                 return false;
             }
@@ -115,11 +117,10 @@
         psStringAppend(&weightName, "%s.%d.%s", outName, numInput, tempWeight);
 
-        if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
-            !readImage(&readout->weight, weightName, config)) {
+        if (!readImage(&output->image, imageName, config) || !readImage(&output->mask, maskName, config) ||
+            !readImage(&output->weight, weightName, config)) {
             psError(PS_ERR_IO, false, "Unable to read previously produced image.");
             psFree(imageName);
             psFree(maskName);
             psFree(weightName);
-            numInput++;
             return false;
         }
@@ -127,158 +128,161 @@
         psFree(maskName);
         psFree(weightName);
-
-        numInput++;
-        return true;
-    }
-#endif
-
-    // Normal operations here
-    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
-    if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
-        assert(psf);
-        assert(sources);
-
-        int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
-        float regionSize = psMetadataLookupF32(NULL, ppsub, "REGION.SIZE"); // Size of iso-kernel regs
-        float spacing = psMetadataLookupF32(NULL, ppsub, "STAMP.SPACING"); // Typical stamp spacing
-        int footprint = psMetadataLookupS32(NULL, ppsub, "STAMP.FOOTPRINT"); // Stamp half-size
-        float threshold = psMetadataLookupF32(NULL, ppsub, "STAMP.THRESHOLD"); // Threshold for stmps
-        int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations
-        float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
-        pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(
-            psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type
-        psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
-        psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
-        int inner = psMetadataLookupS32(NULL, ppsub, "INNER"); // Inner radius
-        int ringsOrder = psMetadataLookupS32(NULL, ppsub, "RINGS.ORDER"); // RINGS polynomial order
-        int binning = psMetadataLookupS32(NULL, ppsub, "SPAM.BINNING"); // Binning for SPAM kernel
-        float badFrac = psMetadataLookupF32(NULL, ppsub, "BADFRAC"); // Maximum bad fraction
-        bool optimum = psMetadataLookupBool(&mdok, ppsub, "OPTIMUM"); // Derive optimum parameters?
-        float optMin = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MIN"); // Minimum width for search
-        float optMax = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MAX"); // Maximum width for search
-        float optStep = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.STEP"); // Step for search
-        float optThresh = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.TOL"); // Tolerance for search
-        int optOrder = psMetadataLookupS32(&mdok, ppsub, "OPTIMUM.ORDER"); // Order for search
-        float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
-
-        // These values are specified specifically for stacking
-        const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Stamps filename
-
-        psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
-        if (optimum) {
-            optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
-        }
-
-        float minFlux = INFINITY;       // Minimum flux for fake image
-        {
-            psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for bg
-            if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal, rng)) {
-                psWarning("Can't measure background for image.");
-                psErrorClear();
-            } else {
-                minFlux = 0.1 * bg->robustStdev;
-            }
-            psFree(bg);
-        }
-
-        // Generate a fake image to match to
-        if (!isfinite(minFlux)) {
-            float maxMag = -INFINITY;   // Maximum magnitude of sources
-            for (int i = 0; i < sources->n; i++) {
-                pmSource *source = sources->data[i]; // Source of interest
-                if (source->psfMag > maxMag && source->psfMag <= MAG_IGNORE &&
-                    !(source->mode & SOURCE_MASK)) {
-                    maxMag = source->psfMag;
+    } else {
+#endif
+
+        // Normal operations here
+        if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
+            assert(psf);
+            assert(sources);
+
+            int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
+            float regionSize = psMetadataLookupF32(NULL, ppsub, "REGION.SIZE"); // Size of iso-kernel regs
+            float spacing = psMetadataLookupF32(NULL, ppsub, "STAMP.SPACING"); // Typical stamp spacing
+            int footprint = psMetadataLookupS32(NULL, ppsub, "STAMP.FOOTPRINT"); // Stamp half-size
+            float threshold = psMetadataLookupF32(NULL, ppsub, "STAMP.THRESHOLD"); // Threshold for stmps
+            int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations
+            float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
+            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(
+                psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type
+            psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
+            psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
+            int inner = psMetadataLookupS32(NULL, ppsub, "INNER"); // Inner radius
+            int ringsOrder = psMetadataLookupS32(NULL, ppsub, "RINGS.ORDER"); // RINGS polynomial order
+            int binning = psMetadataLookupS32(NULL, ppsub, "SPAM.BINNING"); // Binning for SPAM kernel
+            float badFrac = psMetadataLookupF32(NULL, ppsub, "BADFRAC"); // Maximum bad fraction
+            bool optimum = psMetadataLookupBool(&mdok, ppsub, "OPTIMUM"); // Derive optimum parameters?
+            float optMin = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MIN"); // Minimum width for search
+            float optMax = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MAX"); // Maximum width for search
+            float optStep = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.STEP"); // Step for search
+            float optThresh = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.TOL"); // Tolerance for search
+            int optOrder = psMetadataLookupS32(&mdok, ppsub, "OPTIMUM.ORDER"); // Order for search
+            float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
+
+            // These values are specified specifically for stacking
+            const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Stamps filename
+
+            psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
+            if (optimum) {
+                optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
+            }
+
+            float minFlux = INFINITY;       // Minimum flux for fake image
+            {
+                psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for bg
+                if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal, rng)) {
+                    psWarning("Can't measure background for image.");
+                    psErrorClear();
+                } else {
+                    minFlux = 0.1 * bg->robustStdev;
                 }
-            }
-            minFlux = 0.01 * powf(10.0, -0.4 * maxMag);
-        }
-
-        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
-
-        if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources,
-                                      NULL, NULL, psf, minFlux, 0, false)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
+                psFree(bg);
+            }
+
+            // Generate a fake image to match to
+            if (!isfinite(minFlux)) {
+                float maxMag = -INFINITY;   // Maximum magnitude of sources
+                for (int i = 0; i < sources->n; i++) {
+                    pmSource *source = sources->data[i]; // Source of interest
+                    if (source->psfMag > maxMag && source->psfMag <= MAG_IGNORE &&
+                        !(source->mode & SOURCE_MASK)) {
+                        maxMag = source->psfMag;
+                    }
+                }
+                minFlux = 0.01 * powf(10.0, -0.4 * maxMag);
+            }
+
+            pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
+
+            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources,
+                                          NULL, NULL, psf, minFlux, 0, false)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
+                psFree(fake);
+                psFree(optWidths);
+                psFree(output);
+                return false;
+            }
+
+#ifdef TESTING
+            {
+                psFits *fits = psFitsOpen("fake.fits", "w");
+                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+                psFitsClose(fits);
+            }
+#endif
+
+            if (threads > 0) {
+                pmSubtractionThreadsInit(output, NULL, readout, fake);
+            }
+
+            // Do the image matching
+            if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold,
+                                    sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
+                                    binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej,
+                                    maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_1)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+                psFree(fake);
+                psFree(optWidths);
+                psFree(output);
+                return false;
+            }
             psFree(fake);
             psFree(optWidths);
-            psFree(output);
-            return false;
+
+            if (threads > 0) {
+                pmSubtractionThreadsFinalize(output, NULL, readout, fake);
+            }
+
+            // Set the variance factor
+            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+            float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR);
+            if (!isfinite(vf)) {
+                vf = 1.0;
+            }
+            if (isfinite(vfItem->data.F32)) {
+                vfItem->data.F32 *= vf;
+            } else {
+                vfItem->data.F32 = vf;
+            }
+
+            // Replace original images with convolved
+            psFree(readout->image);
+            psFree(readout->mask);
+            psFree(readout->weight);
+            readout->image  = psMemIncrRefCounter(output->image);
+            readout->mask   = psMemIncrRefCounter(output->mask);
+            readout->weight = psMemIncrRefCounter(output->weight);
+        } else {
+            // Fake the convolution
+            psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
+            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
+                             PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
+            psFree(region);
+            pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty,
+                                                                     PM_SUBTRACTION_MODE_1);
+            // Set solution to delta function
+            kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
+            psVectorInit(kernels->solution1, 0.0);
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            kernels->solution1->data.F64[normIndex] = 1.0;
+            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
+                             PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels);
+            psFree(kernels);
         }
 
 #ifdef TESTING
+        // Write convolution kernel
         {
-            psFits *fits = psFitsOpen("fake.fits", "w");
-            psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+            const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
+            assert(outName);
+
+            psString filename = NULL;   // Output filename
+            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+            psString resolved = pmConfigConvertFilename(filename, config, true, false); // Resolved filename
+            psFree(filename);
+            psFits *fits = psFitsOpen(resolved, "w"); // FITS file for subtraction kernel
+            psFree(resolved);
+            pmReadoutWriteSubtractionKernels(output, fits);
             psFitsClose(fits);
         }
-#endif
-
-        if (threads > 0) {
-            pmSubtractionThreadsInit(output, NULL, readout, fake);
-        }
-
-        // Do the image matching
-        if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold,
-                                sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
-                                binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej,
-                                maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_1)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
-            psFree(fake);
-            psFree(optWidths);
-            psFree(output);
-            return false;
-        }
-        psFree(fake);
-        psFree(optWidths);
-
-        if (threads > 0) {
-            pmSubtractionThreadsFinalize(output, NULL, readout, fake);
-        }
-
-        // Set the variance factor
-        psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-        float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR);
-        if (!isfinite(vf)) {
-            vf = 1.0;
-        }
-        if (isfinite(vfItem->data.F32)) {
-            vfItem->data.F32 *= vf;
-        } else {
-            vfItem->data.F32 = vf;
-        }
-
-        // Replace original images with convolved
-        psFree(readout->image);
-        psFree(readout->mask);
-        psFree(readout->weight);
-        readout->image  = psMemIncrRefCounter(output->image);
-        readout->mask   = psMemIncrRefCounter(output->mask);
-        readout->weight = psMemIncrRefCounter(output->weight);
-    } else {
-        // Fake the convolution
-        psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
-        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
-                         PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
-        psFree(region);
-        pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty,
-                                                                 PM_SUBTRACTION_MODE_1);
-        // Set solution to delta function
-        kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
-        psVectorInit(kernels->solution1, 0.0);
-        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-        kernels->solution1->data.F64[normIndex] = 1.0;
-        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
-                         PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels);
-        psFree(kernels);
-    }
-
-#ifdef TESTING
-    {
-        psString filename = NULL;   // Output filename
-        psStringAppend(&filename, "stack_kernel_%d.fits", numInput);
-        psFits *fits = psFitsOpen(filename, "w"); // FITS file for subtraction kernel
-        psFree(filename);
-        pmReadoutWriteSubtractionKernels(output, fits);
-        psFitsClose(fits);
     }
 #endif