Index: trunk/ppStack/src/ppStackMatch.c
===================================================================
--- trunk/ppStack/src/ppStackMatch.c	(revision 18918)
+++ trunk/ppStack/src/ppStackMatch.c	(revision 19123)
@@ -43,4 +43,5 @@
     int renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width for renormalisation box
     float penalty = psMetadataLookupF32(NULL, ppsub, "PENALTY"); // Penalty for wideness
+    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
 
     if (!pmReadoutMaskNonfinite(readout, maskIn)) {
@@ -83,11 +84,27 @@
         }
 
+        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, maskIn, rng)) {
+                psWarning("Can't measure background for image.");
+                psErrorClear();
+            } else {
+                minFlux = bg->robustStdev;
+            }
+            psFree(bg);
+        }
+
         // Generate a fake image to match to
-        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;
+        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.1 * powf(10.0, -0.4 * maxMag);
         }
 
@@ -95,5 +112,5 @@
 
         if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL,
-                                      NULL, psf, powf(10.0, -0.4 * maxMag), 0, false)) {
+                                      NULL, psf, minFlux, 0, false)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
             psFree(fake);
@@ -110,4 +127,8 @@
         }
 #endif
+
+        if (threads > 0) {
+            pmSubtractionThreadsInit(output, NULL, readout, fake);
+        }
 
         // Do the image matching
@@ -124,4 +145,8 @@
         psFree(fake);
         psFree(optWidths);
+
+        if (threads > 0) {
+            pmSubtractionThreadsFinalize(output, NULL, readout, fake);
+        }
 
         // Set the variance factor