Index: /branches/pap_branch_20090108/ppStack/src/ppStackArguments.c
===================================================================
--- /branches/pap_branch_20090108/ppStack/src/ppStackArguments.c	(revision 21149)
+++ /branches/pap_branch_20090108/ppStack/src/ppStackArguments.c	(revision 21150)
@@ -152,4 +152,5 @@
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-poor-frac", 0, "Fraction of weight for poor pixels", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-deconv-limit", 0, "Maximum deconvolution fraction limit", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-image-rej", 0,
                      "Pixel rejection fraction threshold for rejecting entire image", NAN);
@@ -236,4 +237,5 @@
     VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
     VALUE_ARG_RECIPE_FLOAT("-image-rej",      "IMAGE.REJ",      F32);
+    VALUE_ARG_RECIPE_FLOAT("-deconv-limit",   "DECONV.LIMIT",   F32);
     VALUE_ARG_RECIPE_INT("-rows",             "ROWS",           S32, 0);
     VALUE_ARG_RECIPE_FLOAT("-poor-frac",      "POOR.FRACTION",  F32);
Index: /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c
===================================================================
--- /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c	(revision 21149)
+++ /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c	(revision 21150)
@@ -173,4 +173,6 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
+    float deconvLimit = psMetadataLookupF32(NULL, recipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
+
     // Look up appropriate values from the ppSub recipe
     psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
@@ -204,31 +206,29 @@
         assert(outName);
         // Read convolution kernel
-        {
-            psString filename = NULL;   // Output filename
-            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
-            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
-            psFree(filename);
-            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
-            psFree(resolved);
-            if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
-                psFitsClose(fits);
-                return false;
-            }
+        psString filename = NULL;   // Output filename
+        psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+        psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+        psFree(filename);
+        psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+        psFree(resolved);
+        if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
+            psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
             psFitsClose(fits);
-
-            // Add in variance factor
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
-            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-            if (!isfinite(vf)) {
-                vf = 1.0;
-            }
-            if (isfinite(vfItem->data.F32)) {
-                vfItem->data.F32 *= vf;
-            } else {
-                vfItem->data.F32 = vf;
-            }
+            return false;
+        }
+        psFitsClose(fits);
+
+        // Add in variance factor
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
+                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+        float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
+        psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+        if (!isfinite(vf)) {
+            vf = 1.0;
+        }
+        if (isfinite(vfItem->data.F32)) {
+            vfItem->data.F32 *= vf;
+        } else {
+            vfItem->data.F32 = vf;
         }
 
@@ -253,4 +253,10 @@
         psFree(maskName);
         psFree(weightName);
+
+        psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
+                                               PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
+
+        pmSubtractionAnalysis(readout->analysis, kernels, region,
+                              readout->image->numCols, readout->image->numRows);
     } else {
 #endif
@@ -535,4 +541,14 @@
     }
 
+    // Reject image completely if the maximum deconvolution fraction exceeds the limit
+    float deconv = psMetadataLookupF32(NULL, output->analysis,
+                                       PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
+    if (deconv > deconvLimit) {
+        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
+                  deconv, deconvLimit);
+        psFree(output);
+        return NULL;
+    }
+
     // Renormalise the variances if desired
     if (renorm) {