Index: trunk/psModules/src/imcombine/pmStackReject.c
===================================================================
--- trunk/psModules/src/imcombine/pmStackReject.c	(revision 27346)
+++ trunk/psModules/src/imcombine/pmStackReject.c	(revision 27365)
@@ -103,5 +103,5 @@
 
     if (!pmSubtractionThreaded()) {
-        pmSubtractionThreadsInit(NULL, NULL);
+        pmSubtractionThreadsInit();
     }
 
Index: trunk/psModules/src/imcombine/pmSubtractionMask.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionMask.c	(revision 27346)
+++ trunk/psModules/src/imcombine/pmSubtractionMask.c	(revision 27365)
@@ -42,11 +42,27 @@
                            pmSubtractionMode mode)
 {
-    PM_ASSERT_READOUT_NON_NULL(ro1, NULL);
-    PM_ASSERT_READOUT_IMAGE(ro1, NULL);
-    PM_ASSERT_READOUT_MASK(ro1, NULL);
-    PM_ASSERT_READOUT_NON_NULL(ro2, NULL);
-    PM_ASSERT_READOUT_IMAGE(ro2, NULL);
-    PM_ASSERT_READOUT_MASK(ro2, NULL);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL);
+    int numCols = 0, numRows = 0;       // Size of the images
+    if (ro1) {
+        PM_ASSERT_READOUT_NON_NULL(ro1, NULL);
+        PM_ASSERT_READOUT_IMAGE(ro1, NULL);
+        PM_ASSERT_READOUT_MASK(ro1, NULL);
+        numCols = ro1->image->numCols;
+        numRows = ro1->image->numRows;
+            }
+    if (ro2) {
+        PM_ASSERT_READOUT_NON_NULL(ro2, NULL);
+        PM_ASSERT_READOUT_IMAGE(ro2, NULL);
+        PM_ASSERT_READOUT_MASK(ro2, NULL);
+        numCols = ro2->image->numCols;
+        numRows = ro2->image->numRows;
+    }
+    if (ro1 && ro2) {
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL);
+    }
+    if (!ro1 && !ro2) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "No image provided.");
+        return false;
+    }
+    psAssert(numCols > 0 && numRows > 0, "There should be an image provided");
     PS_ASSERT_INT_NONNEGATIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
@@ -56,10 +72,9 @@
     }
 
-    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of the images
-
     // Dereference inputs for convenience
-    psF32 **imageData1 = ro1->image->data.F32, **imageData2 = ro2->image->data.F32;
-    psImageMaskType **maskData1 = ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA,
-        **maskData2 = ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA;
+    psF32 **imageData1 = ro1 ? ro1->image->data.F32 : NULL;
+    psF32 **imageData2 = ro2 ? ro2->image->data.F32 : NULL;
+    psImageMaskType **maskData1 = ro1 ? ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+    psImageMaskType **maskData2 = ro2 ? ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
 
     // First, a pass through to determine the fraction of bad pixels
@@ -69,9 +84,9 @@
         for (int y = 0; y < numRows; y++) {
             for (int x = 0; x < numCols; x++) {
-                if ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x])) {
+                if (ro1 && ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x]))) {
                     numBad++;
                     continue;
                 }
-                if ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x])) {
+                if (ro2 && ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x]))) {
                     numBad++;
                     continue;
@@ -128,8 +143,8 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            if (maskData1[y][x] & maskVal) {
+            if (ro1 && maskData1[y][x] & maskVal) {
                 maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1;
             }
-            if (maskData2[y][x] & maskVal) {
+            if (ro2 && maskData2[y][x] & maskVal) {
                 maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2;
             }
@@ -144,12 +159,12 @@
 
     // Pixels that will be bad (or poor) if we convolve with a bad reference pixel
-    if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1, PM_SUBTRACTION_MASK_CONVOLVE_1,
-                             -size, size, -size, size)) {
+    if (ro1 && !psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1, PM_SUBTRACTION_MASK_CONVOLVE_1,
+                                    -size, size, -size, size)) {
         psError(psErrorCodeLast(), false, "Unable to convolve bad pixels from mask 1.");
         psFree(mask);
         return NULL;
     }
-    if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2, PM_SUBTRACTION_MASK_CONVOLVE_2,
-                             -size, size, -size, size)) {
+    if (ro2 && !psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2, PM_SUBTRACTION_MASK_CONVOLVE_2,
+                                    -size, size, -size, size)) {
         psError(psErrorCodeLast(), false, "Unable to convolve bad pixels from mask 2.");
         psFree(mask);
@@ -174,6 +189,6 @@
         psAbort("Unsupported subtraction mode: %x", mode);
     }
-    if (!psImageConvolveMask(mask, mask, maskRej, PM_SUBTRACTION_MASK_REJ,
-                             -footprint, footprint, -footprint, footprint)) {
+    if (ro1 && ro2 && !psImageConvolveMask(mask, mask, maskRej, PM_SUBTRACTION_MASK_REJ,
+                                           -footprint, footprint, -footprint, footprint)) {
         psError(psErrorCodeLast(), false, "Unable to convolve bad pixels.");
         psFree(mask);
Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 27346)
+++ trunk/psModules/src/imcombine/pmSubtractionMatch.c	(revision 27365)
@@ -132,4 +132,6 @@
     if (subMode != PM_SUBTRACTION_MODE_2) {
         PM_ASSERT_READOUT_NON_NULL(conv1, false);
+        PM_ASSERT_READOUT_NON_NULL(ro1, false);
+        PM_ASSERT_READOUT_IMAGE(ro1, false);
         if (conv1->image) {
             psFree(conv1->image);
@@ -147,4 +149,6 @@
     if (subMode != PM_SUBTRACTION_MODE_1) {
         PM_ASSERT_READOUT_NON_NULL(conv2, false);
+        PM_ASSERT_READOUT_NON_NULL(ro2, false);
+        PM_ASSERT_READOUT_IMAGE(ro2, false);
         if (conv2->image) {
             psFree(conv2->image);
@@ -161,9 +165,7 @@
     }
 
-    PM_ASSERT_READOUT_NON_NULL(ro1, false);
-    PM_ASSERT_READOUT_NON_NULL(ro2, false);
-    PM_ASSERT_READOUT_IMAGE(ro1, false);
-    PM_ASSERT_READOUT_IMAGE(ro2, false);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
+    if (ro1 && ro2) {
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
+    }
     PS_ASSERT_INT_NONNEGATIVE(stride, false);
     if (isfinite(normFrac)) {
@@ -379,6 +381,45 @@
     }
 
+    int numCols, numRows;       // Size of image
+    if (ro1) {
+        numCols = ro1->image->numCols;
+        numRows = ro1->image->numRows;
+    } else if (ro2) {
+        numCols = ro2->image->numCols;
+        numRows = ro2->image->numRows;
+    } else {
+        psAbort("No input image provided.");
+    }
+
     pmSubtractionMaskInvalid(ro1, maskVal);
     pmSubtractionMaskInvalid(ro2, maskVal);
+
+    // General background subtraction, since this is done in pmSubtractionMatch
+    {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for background
+        if (ro1) {
+            psStatsInit(bg);
+            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
+                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                return false;
+            }
+            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+        }
+        if (ro2) {
+            psStatsInit(bg);
+            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
+                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                return false;
+            }
+            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+        }
+        psFree(bg);
+        psFree(rng);
+    }
 
     psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
@@ -396,5 +437,4 @@
     psMetadata *outHeader = psMetadataAlloc(); // Output header values
 
-    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of image
     subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, mode, numCols, numRows);
 
@@ -404,6 +444,5 @@
         psRegion *region = regions->data[i]; // Region of interest
 
-        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region,
-                                   ro1->image->numCols, ro1->image->numRows)) {
+        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region, numCols, numRows)) {
             psError(psErrorCodeLast(), false, "Unable to generate QA data");
             psFree(outAnalysis);
@@ -454,4 +493,8 @@
         return false;
     }
+
+    // We need both inputs
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
 
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
Index: trunk/psModules/src/imcombine/pmSubtractionThreads.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionThreads.c	(revision 27346)
+++ trunk/psModules/src/imcombine/pmSubtractionThreads.c	(revision 27365)
@@ -17,38 +17,5 @@
 }
 
-// Initialise a mutex in each of the input
-static void subtractionMutexInit(pmReadout *ro)
-{
-    if (!ro) {
-        return;
-    }
-
-    // XXX if (ro->image) {
-    // XXX     psMutexInit(ro->image);
-    // XXX }
-    // XXX if (ro->variance) {
-    // XXX     psMutexInit(ro->variance);
-    // XXX }
-
-    return;
-}
-
-static void subtractionMutexDestroy(pmReadout *ro)
-{
-    if (!ro) {
-        return;
-    }
-
-    // XXX if (ro->image) {
-    // XXX     psMutexDestroy(ro->image);
-    // XXX }
-    // XXX if (ro->variance) {
-    // XXX     psMutexDestroy(ro->variance);
-    // XXX }
-
-    return;
-}
-
-void pmSubtractionThreadsInit(pmReadout *in1, pmReadout *in2)
+void pmSubtractionThreadsInit(void)
 {
     if (threaded) {
@@ -57,7 +24,4 @@
 
     threaded = true;
-
-    subtractionMutexInit(in1);
-    subtractionMutexInit(in2);
 
     {
@@ -86,5 +50,5 @@
 
 
-void pmSubtractionThreadsFinalize(pmReadout *in1, pmReadout *in2)
+void pmSubtractionThreadsFinalize(void)
 {
     if (!threaded) {
@@ -97,6 +61,4 @@
     psThreadTaskRemove("PSMODULES_SUBTRACTION_CONVOLVE");
 
-    subtractionMutexDestroy(in1);
-    subtractionMutexDestroy(in1);
     return;
 }
Index: trunk/psModules/src/imcombine/pmSubtractionThreads.h
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionThreads.h	(revision 27346)
+++ trunk/psModules/src/imcombine/pmSubtractionThreads.h	(revision 27365)
@@ -7,13 +7,12 @@
 /// Set up threading for image matching
 ///
-/// Sets up thread tasks, and initialises mutexes in readouts
-void pmSubtractionThreadsInit(pmReadout *in1, pmReadout *in2 // Input readouts
-    );
+/// Sets up thread tasks
+void pmSubtractionThreadsInit(void);
 
 
 /// Take down threading for image matching
 ///
-/// Destroys thread tasks, and initialises mutexes in readouts
-void pmSubtractionThreadsFinalize(pmReadout *in1, pmReadout *in2);
+/// Destroys thread tasks
+void pmSubtractionThreadsFinalize(void);
 
 #endif