Index: trunk/ppStack/src/ppStackMatch.c
===================================================================
--- trunk/ppStack/src/ppStackMatch.c	(revision 20754)
+++ trunk/ppStack/src/ppStackMatch.c	(revision 20755)
@@ -80,10 +80,10 @@
             continue;
         }
-        coords->data.F32[0] = source->peak->xf;
-        coords->data.F32[1] = source->peak->yf;
+        coords->data.F64[0] = source->peak->xf;
+        coords->data.F64[1] = source->peak->yf;
 
         long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
-        psTrace("ppStack", 5, "Source at %.0f,%.0f has %ld sources in exclusion zone",
-                coords->data.F32[0], coords->data.F32[1], numWithin);
+        psTrace("ppStack", 5, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
+                coords->data.F64[0], coords->data.F64[1], numWithin);
         if (numWithin == 1) {
             // Only itself inside the exclusion zone
@@ -101,4 +101,55 @@
 }
 
+#define BG_SIZE 64                      // Large box half-size in which to measure background
+
+// Generate a background model of the readout we're matching
+psImage *stackBackgroundModel(pmReadout *ro, psMaskType maskVal, const psArray *sources, int size)
+{
+    psImage *image = ro->image, *mask = ro->mask; // Image and mask of readout
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    psImage *bgImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Background image
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for measuring background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+
+        int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source
+
+        int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
+        int yMin = PS_MAX(y - BG_SIZE, 0), yMax = PS_MIN(y + BG_SIZE, numCols);
+
+        psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
+        psImage *subImage = psImageSubset(image, region); // Subimage
+        psImage *subMask = psImageSubset(mask, region); // Subimage mask
+        if (!psImageBackground(stats, NULL, subImage, subMask, maskVal, rng)) {
+            // Can't do anything about it
+            psErrorClear();
+            continue;
+        }
+        psFree(subImage);
+        psFree(subMask);
+
+        float value = stats->robustMedian; // Value to set background to
+
+        int uMin = PS_MAX(x - size, 0), uMax = PS_MIN(x + size, numCols);
+        int vMin = PS_MAX(y - size, 0), vMax = PS_MIN(y + size, numCols);
+
+        for (int v = vMin; v < vMax; v++) {
+            for (int u = uMin; u < uMax; u++) {
+                bgImage->data.F32[v][u] = value;
+            }
+        }
+    }
+
+    psFree(stats);
+    psFree(rng);
+
+    return bgImage;
+}
 
 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
@@ -242,9 +293,10 @@
                     psErrorClear();
                 } else {
-                    minFlux = 0.1 * bg->robustStdev;
+                    minFlux = bg->robustStdev;
                 }
                 psFree(bg);
             }
 
+#if 0
             float maxMag = -INFINITY;   // Maximum magnitude of sources
             for (int i = 0; i < sources->n; i++) {
@@ -256,9 +308,13 @@
             }
             minFlux = PS_MIN(FAINT_SOURCE_FRAC * powf(10.0, -0.4 * maxMag), minFlux);
+#endif
 
             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)) {
+            // For the sake of stamps, remove nearby sources
+            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
+
+            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
+                                          stampSources, NULL, NULL, psf, minFlux, 0, false)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
@@ -267,4 +323,10 @@
                 return false;
             }
+
+            // Add the background into the target image
+            psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
+            psBinaryOp(fake->image, fake->image, "+", bgImage);
+            psFree(bgImage);
+
 
 #ifdef TESTING
@@ -301,7 +363,4 @@
                 }
             }
-
-            // For the sake of stamps, remove nearby sources
-            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
 
             if (threads > 0) {
