Index: /branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c	(revision 21302)
+++ /branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c	(revision 21303)
@@ -997,6 +997,5 @@
 
 // XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data
-static bool subtractionConvolvePatch(psKernel **covar1, psKernel **covar2, // Covariance matrices
-                                     int numCols, int numRows, // Size of image
+static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image
                                      int x0, int y0, // Offsets for image
                                      pmReadout *out1, pmReadout *out2, // Output readouts
@@ -1031,5 +1030,4 @@
                        ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region,
                        maskBad, maskPoor, poorFrac, useFFT, false);
-        *covar1 = psImageCovarianceCalculate(kernelImage, ro1->covariance);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -1037,5 +1035,4 @@
                        ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region,
                        maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
-        *covar2 = psImageCovarianceCalculate(kernelImage, ro2->covariance);
     }
 
@@ -1093,16 +1090,6 @@
     bool useFFT = PS_SCALAR_VALUE(args->data[18], PS_TYPE_IMAGE_MASK_DATA); // Use FFT for convolution?
 
-    psKernel *covar1 = NULL, *covar2 = NULL; // Covariance matrices to return
-
-    if (!subtractionConvolvePatch(&covar1, &covar2, numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,
-                                  sys1, sys2, subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG,
-                                  useFFT)) {
-        return false;
-    }
-    job->results = psArrayAlloc(2);
-    job->results->data[0] = covar1;
-    job->results->data[1] = covar2;
-
-    return true;
+    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
+                                    subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
 }
 
@@ -1264,6 +1251,4 @@
         stride = 2 * size + 1;
     }
-
-    psList *covariances1 = psListAlloc(NULL), *covariances2 = psListAlloc(NULL); // List of covariances
 
     for (int j = yMin; j < yMax; j += stride) {
@@ -1321,16 +1306,7 @@
                 psFree(job);
             } else {
-                psKernel *covar1 = NULL, *covar2 = NULL; // Covariance matrices
-                subtractionConvolvePatch(&covar1, &covar2, numCols, numRows, x0, y0, out1, out2, convMask,
-                                         ro1, ro2, sys1, sys2, subMask, maskBad, maskPoor, poorFrac,
-                                         subRegion, kernels, doBG, useFFT);
-                if (covar1) {
-                    psListAdd(covariances1, PS_LIST_TAIL, covar1);
-                }
-                if (covar2) {
-                    psListAdd(covariances2, PS_LIST_TAIL, covar2);
-                }
-                psFree(covar1);
-                psFree(covar2);
+                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
+                                         subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG,
+                                         useFFT);
             }
             psFree(subRegion);
@@ -1352,11 +1328,4 @@
             psAssert(strcmp(job->type, "PSMODULES_SUBTRACTION_CONVOLVE") == 0,
                      "Job has incorrect type: %s", job->type);
-            psAssert(job->results->n == 2, "Job has insufficient results: %ld", job->results->n);
-            if (job->results->data[0]) {
-                psListAdd(covariances1, PS_LIST_TAIL, job->results->data[0]);
-            }
-            if (job->results->data[1]) {
-                psListAdd(covariances2, PS_LIST_TAIL, job->results->data[1]);
-            }
             psFree(job);
         }
@@ -1377,16 +1346,18 @@
     psFree(sys2);
 
-    if (psListLength(covariances1) > 0) {
-        psArray *covar = psListToArray(covariances1);
-        out1->covariance = psImageCovarianceAverage(covar);
-        psFree(covar);
-    }
-    if (psListLength(covariances2) > 0) {
-        psArray *covar = psListToArray(covariances2);
-        out2->covariance = psImageCovarianceAverage(covar);
-        psFree(covar);
-    }
-    psFree(covariances1);
-    psFree(covariances2);
+    // Calculate covariances
+    // This can take a while, so we only do it for a single instance
+    // XXX psImageCovarianceCalculate could be multithreaded
+    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+        out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance);
+        psFree(kernel);
+    }
+    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0,
+                                               kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel
+        out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance);
+        psFree(kernel);
+    }
 
     // Copy anything that wasn't convolved
@@ -1397,4 +1368,5 @@
             out2->variance = psMemIncrRefCounter(ro2->variance);
             out2->mask = psMemIncrRefCounter(ro2->mask);
+            out2->covariance = psMemIncrRefCounter(ro2->covariance);
         }
         break;
@@ -1404,4 +1376,5 @@
             out1->variance = psMemIncrRefCounter(ro1->variance);
             out1->mask = psMemIncrRefCounter(ro1->mask);
+            out1->covariance = psMemIncrRefCounter(ro1->covariance);
         }
         break;