Index: trunk/ppStack/src/ppStackMatch.c
===================================================================
--- trunk/ppStack/src/ppStackMatch.c	(revision 15850)
+++ trunk/ppStack/src/ppStackMatch.c	(revision 16605)
@@ -9,9 +9,19 @@
 #include "ppStack.h"
 
+#define ARRAY_BUFFER 16                 // Number to add to array at a time
+
+
 //#define TESTING
 
-bool ppStackMatch(pmReadout *output, const pmReadout *input, const pmReadout *sourcesRO,
-                  const pmPSF *psf, const pmConfig *config)
+bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
+                  const pmReadout *sourcesRO, const pmPSF *psf, const pmConfig *config)
 {
+    assert(readout);
+    assert(regions && !*regions);
+    assert(kernels && !*kernels);
+    assert(sourcesRO);
+    assert(psf);
+    assert(config);
+
     // Look up appropriate values from the ppSub recipe
     bool mdok;                          // Status of MD lookup
@@ -68,5 +78,5 @@
     pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
 
-    if (!pmReadoutFakeFromSources(fake, input->image->numCols, input->image->numRows, sources, NULL, NULL,
+    if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL, NULL,
                                   psf, powf(10.0, -0.4 * maxMag), 0, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
@@ -85,5 +95,6 @@
 
     // Do the image matching
-    if (!pmSubtractionMatch(output, NULL, input, fake, footprint, regionSize, spacing, threshold,
+    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
+    if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold,
                             sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
                             binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
@@ -92,4 +103,5 @@
         psFree(fake);
         psFree(optWidths);
+        psFree(output);
         return false;
     }
@@ -97,4 +109,51 @@
     psFree(optWidths);
 
+    // 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);
+
+    // Extract the regions and solutions used in the image matching
+    // This stops them from being freed when we iterate back up the FPA
+    *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
+    {
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_REGION);
+            *regions = psArrayAdd(*regions, ARRAY_BUFFER, item->data.V);
+        }
+        psFree(iter);
+    }
+    *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
+    {
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_UNKNOWN);
+            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
+            // images by scans.
+            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
+            kernel->numCols = readout->image->numCols;
+            kernel->numRows = readout->image->numRows;
+
+            *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel);
+        }
+        psFree(iter);
+    }
+    assert((*regions)->n == (*kernels)->n);
+
+
+    psFree(output);
+
     return true;
 }