Index: /trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 28149)
+++ /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 28150)
@@ -32,4 +32,5 @@
 #define MIN_SAMPLE_STATS    7           // Minimum number to use sample statistics; otherwise use quartiles
 #define USE_KERNEL_ERR                  // Use kernel error image?
+#define NUM_COVAR_POS 5                 // Number of positions for covariance calculation
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1328,17 +1329,33 @@
     // This can be fairly involved, so we only do it for a single instance
     // Enable threads for covariance calculation, since we're not threading on top of it.
+    float position[NUM_COVAR_POS] = { -1.0, -0.5, 0.0, +0.5, +1.0 }; // Positions for covariance calculations
     oldThreads = psImageCovarianceSetThreads(true);
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
-        psKernelTruncate(kernel, covarFrac);
-        out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance);
-        psFree(kernel);
+        psArray *covars = psArrayAlloc(PS_SQR(NUM_COVAR_POS)); // Covariances
+        for (int y = 0, i = 0; y < NUM_COVAR_POS; y++) {
+            for (int x = 0; x < NUM_COVAR_POS; x++, i++) {
+                psKernel *kernel = pmSubtractionKernel(kernels, position[x], position[y],
+                                                       false); // Convolution kernel
+                psKernelTruncate(kernel, covarFrac);
+                covars->data[i] = psImageCovarianceCalculate(kernel, ro1->covariance);
+                psFree(kernel);
+            }
+        }
+        out1->covariance = psImageCovarianceAverage(covars);
+        psFree(covars);
     }
     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
-        psKernelTruncate(kernel, covarFrac);
-        out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance);
-        psFree(kernel);
+        psArray *covars = psArrayAlloc(PS_SQR(NUM_COVAR_POS)); // Covariances
+        for (int y = 0, i = 0; y < NUM_COVAR_POS; y++) {
+            for (int x = 0; x < NUM_COVAR_POS; x++, i++) {
+                psKernel *kernel = pmSubtractionKernel(kernels, position[x], position[y],
+                                                       kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Convolution kernel
+                psKernelTruncate(kernel, covarFrac);
+                covars->data[i] = psImageCovarianceCalculate(kernel, ro2->covariance);
+                psFree(kernel);
+            }
+        }
+        out2->covariance = psImageCovarianceAverage(covars);
+        psFree(covars);
     }
     psImageCovarianceSetThreads(oldThreads);