Index: /branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c	(revision 21296)
+++ /branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c	(revision 21297)
@@ -997,5 +997,6 @@
 
 // XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data
-static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image
+static bool subtractionConvolvePatch(psKernel **covar1, psKernel **covar2, // Covariance matrices
+                                     int numCols, int numRows, // Size of image
                                      int x0, int y0, // Offsets for image
                                      pmReadout *out1, pmReadout *out2, // Output readouts
@@ -1030,4 +1031,5 @@
                        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) {
@@ -1035,4 +1037,5 @@
                        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);
     }
 
@@ -1090,6 +1093,17 @@
     bool useFFT = PS_SCALAR_VALUE(args->data[18], PS_TYPE_IMAGE_MASK_DATA); // Use FFT for convolution?
 
-    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
-                                    subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
+    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;
+    }
+    psArrayAdd(job->results, 1, covar1);
+    psArrayAdd(job->results, 1, covar2);
+    psFree(covar1);
+    psFree(covar2);
+
+    return true;
 }
 
@@ -1251,4 +1265,6 @@
         stride = 2 * size + 1;
     }
+
+    psList *covariances1 = psListAlloc(NULL), *covariances2 = psListAlloc(NULL); // List of covariances
 
     for (int j = yMin; j < yMax; j += stride) {
@@ -1306,7 +1322,16 @@
                 psFree(job);
             } else {
-                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,
-                                         sys1, sys2, subMask, maskBad, maskPoor, poorFrac, subRegion,
-                                         kernels, doBG, useFFT);
+                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);
             }
             psFree(subRegion);
@@ -1328,4 +1353,11 @@
             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);
         }
@@ -1341,9 +1373,21 @@
         }
     }
-
     psImageConvolveSetThreads(oldThreads);
 
     psFree(sys1);
     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);
 
     // Copy anything that wasn't convolved