Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 28405)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 28667)
@@ -52,15 +52,9 @@
 
     // Take the square of the normal kernel
-    double sumVariance = 0.0; // Sum of the variance kernels
     for (int v = yMin; v <= yMax; v++) {
         for (int u = xMin; u <= xMax; u++) {
-            sumVariance += out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
-        }
-    }
-
-    // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
-    // This is required to keep the relative scaling between the image and the variance map
-    psBinaryOp(out->image, out->image, "*", psScalarAlloc(1.0 / sumVariance, PS_TYPE_F32));
-
+            out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
+        }
+    }
     return out;
 }
@@ -271,12 +265,12 @@
 // Convolve an image using FFT
 static void convolveVarianceFFT(psImage *target,// Place the result in here
-                              psImage *variance, // Variance map to convolve
-                              psImage *kernelErr, // Kernel error image
-                              psImage *mask, // Mask image
-                              psImageMaskType maskVal, // Value to mask
-                              const psKernel *kernel, // Kernel by which to convolve
-                              psRegion region,// Region of interest
-                              int size        // Size of (square) kernel
-                              )
+                                psImage *variance, // Variance map to convolve
+                                psImage *kernelErr, // Kernel error image
+                                psImage *mask, // Mask image
+                                psImageMaskType maskVal, // Value to mask
+                                const psKernel *kernel, // Kernel by which to convolve
+                                psRegion region,// Region of interest
+                                int size        // Size of (square) kernel
+                                )
 {
     psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
@@ -348,4 +342,5 @@
                                   psImage *image, // Image to convolve
                                   psImage *variance, // Variance map to convolve, or NULL
+                                  const psKernel *covar,               // Covariance, or NULL
                                   psImage *kernelErr, // Kernel error image, or NULL
                                   psImage *subMask, // Subtraction mask
@@ -393,4 +388,14 @@
         if (variance) {
             convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size);
+        }
+    }
+
+    if (variance && covar) {
+        // Apply covariance factor to variance map, to allow for spatial variation
+        float factor = psImageCovarianceCalculateFactor(*kernelImage, covar); // Factor to apply
+        for (int y = region.y0; y < region.y1; y++) {
+            for (int x = region.x0; x < region.x1; x++) {
+                convVariance->data.F32[y][x] *= factor;
+            }
         }
     }
@@ -1085,11 +1090,11 @@
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out1->image, out1->variance, out1->mask, &kernelImage, &kernelVariance,
-                       ro1->image, ro1->variance, kernelErr1, subMask, kernels, polyValues, background,
-                       *region, maskBad, maskPoor, poorFrac, useFFT, false);
+                       ro1->image, ro1->variance, ro1->covariance, kernelErr1, subMask, kernels,
+                       polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT, false);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out2->image, out2->variance, out2->mask, &kernelImage, &kernelVariance,
-                       ro2->image, ro2->variance, kernelErr2, subMask, kernels, polyValues, background,
-                       *region, maskBad, maskPoor, poorFrac, useFFT,
+                       ro2->image, ro2->variance, ro2->covariance, kernelErr2, subMask, kernels,
+                       polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT,
                        kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
@@ -1325,7 +1330,7 @@
 
     // Calculate covariances
-    // 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.
+    // This can be fairly involved, so we only do it for a small number of instances
     float position[NUM_COVAR_POS] = { -1.0, -0.5, 0.0, +0.5, +1.0 }; // Positions for covariance calculations
+    // Enable threads for covariance calculation, since we're not threading on top of it
     oldThreads = psImageCovarianceSetThreads(true);
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -1342,4 +1347,7 @@
         out1->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
+        // Remove covariance factor from covariance, since we've put it in the variance map already
+        float factor = psImageCovarianceFactor(out1->covariance);
+        psBinaryOp(out1->covariance->image, out1->covariance->image, "/", psScalarAlloc(factor, PS_TYPE_F32));
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -1356,4 +1364,7 @@
         out2->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
+        // Remove covariance factor from covariance, since we've put it in the variance map already
+        float factor = psImageCovarianceFactor(out2->covariance);
+        psBinaryOp(out2->covariance->image, out2->covariance->image, "/", psScalarAlloc(factor, PS_TYPE_F32));
     }
     psImageCovarianceSetThreads(oldThreads);