Index: trunk/archive/noise_model/simulate.c
===================================================================
--- trunk/archive/noise_model/simulate.c	(revision 28666)
+++ trunk/archive/noise_model/simulate.c	(revision 28667)
@@ -6,9 +6,8 @@
 #define DET_SIZE 2000
 #define SKY_SIZE 1000
-#define SIZE 1024
 #define OUTPUT_ROOT "test"
-#define SCALE 0.654321
+#define SCALE 0.7654321
 #define ROT M_PI_2
-#define INTERPOLATION PS_INTERPOLATE_LANCZOS4
+#define INTERPOLATION PS_INTERPOLATE_LANCZOS3
 #define OFFSET 16
 #define SMOOTH_SIGMA 6.54321
@@ -16,4 +15,5 @@
 #define DUAL_KERNEL "sub.subkernel"
 #define WARP_NUM 10000
+#define CONV_NUM 1000
 
 static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 };
@@ -21,8 +21,13 @@
 static const char *rootNames[] = { "o5298g0209o.fake.warp", "o5298g0210o.fake.warp", "o5298g0211o.fake.warp",
                                    "o5298g0212o.fake.warp", "o5298g0213o.fake.warp", "o5298g0214o.fake.warp",
+
                                    "o5298g0215o.fake.warp", "o5298g0216o.fake.warp" };
+#if 1
 static const char *kernels[] = { "test.5.kernel", "test.11.kernel", "test.17.kernel", "test.23.kernel",
                                  "test.29.kernel", "test.35.kernel", "test.41.kernel", "test.47.kernel" };
-
+#else
+static const char *kernels[] = { "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel",
+                                 "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel" };
+#endif
 
 void writeImage(const psImage *image, const char *suffix)
@@ -76,10 +81,8 @@
 float meanVar(const psImage *var, const psImage *mask, const psKernel *covar)
 {
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    psImageBackground(stats, NULL, var, mask, 0xFF, rng);
+    psImageStats(stats, var, mask, 0xFF);
     float variance = stats->sampleMedian * psImageCovarianceFactor(covar);
     psFree(stats);
-    psFree(rng);
     return variance;
 }
@@ -93,4 +96,7 @@
         for (int x = 0; x < image->numCols; x++) {
             sn->data.F32[y][x] = image->data.F32[y][x] / sqrtf(var->data.F32[y][x] * varFactor);
+            if (!isfinite(sn->data.F32[y][x])) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0xFF;
+            }
         }
     }
@@ -100,4 +106,11 @@
     float offset = stats->sampleMean;
     psFree(stats);
+
+    static int number = 0;
+
+    psString snName = NULL;
+    psStringAppend(&snName, "sn.%d.fits", number);
+    writeImage(sn, snName);
+    fprintf(stderr, "Writing S/N image %d\n", number);
 
     psHistogram *hist = psHistogramAlloc(-5, +5, 101);
@@ -121,8 +134,7 @@
     psFree(sn);
 
-    static int number = 0;
     psString name = NULL;
     fprintf(stderr, "Writing histogram %d\n", number);
-    psStringAppend(&name, "hist_%d.dat", number++);
+    psStringAppend(&name, "hist_%d.dat", number);
     FILE *file = fopen(name, "w");
     psFree(name);
@@ -136,4 +148,6 @@
     psFree(hist);
 
+    number++;
+
     return noise;
 }
@@ -144,14 +158,7 @@
 {
     psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
-    double sum = 0.0, sum2 = 0.0;
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            sum += kernel->kernel[y][x];
-            sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
-        }
-    }
-    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
-        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            kernel2->kernel[y][x] /= sum2;
+            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -176,21 +183,23 @@
 
     psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
+    double sum2 = 0.0;                                               // Sum of kernel squared
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sum2 += PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    float factor = 1.0 / sum2;
     psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar);
     psFree(kernel);
-    psImageCovarianceTransfer(smoothVariance, smoothCovar);
+
+    // Apply square root of significance image scaling factor to image
+    psBinaryOp(smoothImage, smoothImage, "*", psScalarAlloc(sqrtf(factor), PS_TYPE_F32));
+
 #else
     psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
     psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
-    double sum = 0.0, sum2 = 0.0;
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            sum += kernel->kernel[y][x];
-            sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
-        }
-    }
-    fprintf(stderr, "Kernel sum: %f %f\n", sum, sum2);
-    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
-        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
-            kernel2->kernel[y][x] /= sum2;
+            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -332,4 +341,19 @@
 }
 
+float covarianceSum(const psKernel *covar)
+{
+    if (!covar) {
+        return 1.0;
+    }
+    double sum = 0.0;
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            sum += covar->kernel[y][x];
+        }
+    }
+    return sum;
+}
+
+
 int main(int argc, char *argv[])
 {
@@ -383,9 +407,10 @@
         }
 
-        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(inImage, inMask, inVariance, inCovar),
                 psImageCovarianceFactor(inCovar),
-                meanVar(inVariance, inMask, inCovar));
+                meanVar(inVariance, inMask, inCovar),
+                covarianceSum(inCovar));
 
         phot(inImage, inMask, inVariance, inCovar);
@@ -417,9 +442,10 @@
         }
 
-        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(warpImage, warpMask, warpVariance, warpCovar),
                 psImageCovarianceFactor(warpCovar),
-                meanVar(warpVariance, warpMask, warpCovar));
+                meanVar(warpVariance, warpMask, warpCovar),
+                covarianceSum(warpCovar));
 
         phot(warpImage, warpMask, warpVariance, warpCovar);
@@ -434,4 +460,11 @@
         pmReadoutReadSubtractionKernels(ro, fits);
         pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels);
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels);
+#if 1
+//        kernels->solution1->data.F64[normIndex] += 1.0;
+        kernels->solution1->data.F64[bgIndex] = 0.0;
+#endif
+        fprintf(stderr, "Norm: %f BG: %f\n", kernels->solution1->data.F64[normIndex], kernels->solution1->data.F64[bgIndex]);
 #if 0
         for (int i = 0; i < kernels->num; i++) {
@@ -455,4 +488,36 @@
         psFitsClose(fits);
 
+#if 0
+        {
+            psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+            psArray *covariances = psArrayAlloc(CONV_NUM);
+            psVector *factors = psVectorAlloc(CONV_NUM, PS_TYPE_F32);
+            double mean = 0.0;
+            for (int i = 0; i < CONV_NUM; i++) {
+                float x, y;
+                p_pmSubtractionPolynomialNormCoords(&x, &y, psRandomUniform(rng) * SKY_SIZE,
+                                                    psRandomUniform(rng) * SKY_SIZE,
+                                                    kernels->xMin, kernels->xMax,
+                                                    kernels->yMin, kernels->yMax);
+                psKernel *kernel = pmSubtractionKernel(kernels, x, y, false);
+                psKernel *covar = covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
+                psFree(kernel);
+                mean += factors->data.F32[i] = psImageCovarianceFactor(covar);
+            }
+            psFree(rng);
+            psFree(covariances);
+
+            mean /= WARP_NUM;
+
+            double stdev = 0.0;
+            for (int i = 0; i < CONV_NUM; i++) {
+                stdev += PS_SQR(factors->data.F32[i] - mean);
+            }
+            stdev = sqrt(stdev/(CONV_NUM-1));
+            fprintf(stderr, "Conv covariance mean: %f stdev: %f\n", mean, stdev);
+            psFree(factors);
+        }
+#endif
+
         pmReadout *conv = pmReadoutAlloc(NULL);
 #if 1
@@ -460,5 +525,5 @@
         conv->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
         conv->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
-        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.01,
+        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 2 * kernels->size + 1, 0.0, 0.01,
                                        0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
             psErrorStackPrint(stderr, "Error:");
@@ -494,13 +559,15 @@
         }
 
-        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(conv->image, conv->mask, conv->variance, conv->covariance),
                 psImageCovarianceFactor(conv->covariance),
-                meanVar(conv->variance, conv->mask, conv->covariance));
+                meanVar(conv->variance, conv->mask, conv->covariance),
+                covarianceSum(conv->covariance));
 
         phot(conv->image, conv->mask, conv->variance, conv->covariance);
         readouts->data[i] = conv;
 
+        //        exit(1);
     }
 
@@ -624,6 +691,4 @@
                 meanVar(diffVariance, diffMask, diffCovar));
 
-        phot(diffImage, diffMask, diffVariance, diffCovar);
-
         writeImage(diffImage, "ssdiff.image.fits");
         writeImage(diffMask, "ssdiff.mask.fits");
Index: trunk/psLib/src/imageops/psImageCovariance.c
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.c	(revision 28666)
+++ trunk/psLib/src/imageops/psImageCovariance.c	(revision 28667)
@@ -34,5 +34,6 @@
 static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix
                                       const psKernel *kernel, // Convolution kernel
-                                      int x, int y            // Coordinates in output covariance matrix
+                                      int x, int y,           // Coordinates in output covariance matrix
+                                      float scale             // Scale to apply
                                       )
 {
@@ -83,5 +84,5 @@
     }
 
-    return sum;
+    return scale * sum;
 }
 
@@ -91,5 +92,5 @@
     PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     psAssert(job->args, "No job arguments");
-    psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n);
+    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
 
     psKernel *out = job->args->data[0]; // Output covariance matrix
@@ -98,6 +99,7 @@
     int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix
     int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix
-
-    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
+    float scale = PS_SCALAR_VALUE(job->args->data[5], F32); // Scaling to apply
+
+    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
 
     return true;
@@ -127,4 +129,5 @@
 
     // Check for non-finite elements
+    double sumKernel = 0.0, sumKernel2 = 0.0; // Sum of the kernel
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
@@ -135,4 +138,6 @@
                 return NULL;
             }
+            sumKernel += kernel->kernel[y][x];
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -155,4 +160,5 @@
     int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
     psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
+    float scale = 1.0 / sumKernel2;          // Scaling to apply
 
     for (int y = yMin; y <= yMax; y++) {
@@ -165,4 +171,5 @@
                 PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
                 PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, scale, PS_TYPE_F32);
                 if (!psThreadJobAddPending(job)) {
                     psFree(covar);
@@ -170,9 +177,8 @@
                 }
             } else {
-                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
-            }
-        }
-    }
-    psFree(covar);
+                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
+            }
+        }
+    }
 
     if (threaded && !psThreadPoolWait(true)) {
@@ -180,4 +186,6 @@
         return false;
     }
+
+    psFree(covar);
 
     return out;
@@ -194,4 +202,5 @@
 
     // Check for non-finite elements
+    double sumKernel2 = 0.0; // Sum of the squared kernel
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
@@ -202,4 +211,5 @@
                 return NAN;
             }
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -215,5 +225,6 @@
     }
 
-    float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
+    float scale = 1.0 / sumKernel2;     // Scale to apply
+    float factor = imageCovarianceCalculate(covar, kernel, 0, 0, scale); // Covariance factor
     psFree(covar);
     return factor;
@@ -338,10 +349,9 @@
         }
     }
-    psFree(covar);
-
     if (threaded && !psThreadPoolWait(true)) {
         psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
         return false;
     }
+    psFree(covar);
 
     return out;
@@ -616,5 +626,5 @@
     if (set && !threaded) {
         {
-            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 6);
             task->function = &imageCovarianceCalculateThread;
             psThreadTaskAdd(task);
Index: trunk/psLib/src/imageops/psImageInterpolate.c
===================================================================
--- trunk/psLib/src/imageops/psImageInterpolate.c	(revision 28666)
+++ trunk/psLib/src/imageops/psImageInterpolate.c	(revision 28667)
@@ -427,32 +427,42 @@
 
 // Determine the result of the interpolation after all the math has been done
-#define INTERPOLATE_RESULT() \
-    psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; /* Status of interpolation */ \
-    *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage; \
-    if (wantVariance) { \
-        *varianceValue = sumVariance / (sumKernel2 - sumBad); \
-    } \
-    if (sumKernel == 0.0) { \
-        /* No kernel contributions */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->badMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_BAD; \
-    } else if (sumBad == 0) { \
-        /* Completely good pixel */ \
-        status = PS_INTERPOLATE_STATUS_GOOD; \
-    } else if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { \
-        /* Some pixels masked: poor pixel */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->poorMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_POOR; \
-    } else { \
-        /* Many pixels (or a few important pixels) masked: bad pixel */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->badMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_BAD; \
-    }
+static psImageInterpolateStatus interpolateResult(const psImageInterpolation *interp,
+                                                  double *imageValue, double *varianceValue,
+                                                  psImageMaskType *maskValue,
+                                                  double sumImage, double sumVariance, double sumBad,
+                                                  double sumKernel, double sumKernel2,
+                                                  bool wantVariance, bool haveMask)
+{
+    *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage;
+    if (wantVariance) {
+        if (sumBad > 0) {
+            sumVariance *= sumKernel2 / (sumKernel2 - sumBad);
+        }
+        *varianceValue = sumVariance / PS_SQR(sumKernel);
+    }
+    if (sumKernel == 0.0) {
+        // No kernel contributions at all
+        if (haveMask && maskValue) {
+            *maskValue |= interp->badMask;
+        }
+        return PS_INTERPOLATE_STATUS_BAD;
+    }
+    if (sumBad == 0) {
+        // Completely good pixel
+        return PS_INTERPOLATE_STATUS_GOOD;
+    }
+    if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) {
+        // Some pixels masked: poor pixel
+        if (haveMask && maskValue) {
+            *maskValue |= interp->poorMask;
+        }
+        return PS_INTERPOLATE_STATUS_POOR;
+    }
+    // Many pixels (or a few important pixels) masked: bad pixel
+    if (haveMask && maskValue) {
+        *maskValue |= interp->badMask;
+    }
+    return PS_INTERPOLATE_STATUS_BAD;
+}
 
 // Interpolation engine for separable interpolation kernels
@@ -703,6 +713,4 @@
     }
 
-    INTERPOLATE_RESULT();
-
     psFree(xKernelNew);
     psFree(yKernelNew);
@@ -710,5 +718,6 @@
     psFree(yKernel2New);
 
-    return status;
+    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
+                             sumKernel, sumKernel2, wantVariance, haveMask);
 }
 
@@ -861,7 +870,6 @@
     }
 
-    INTERPOLATE_RESULT();
-
-    return status;
+    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
+                             sumKernel, sumKernel2, wantVariance, haveMask);
 }
 
Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 28666)
+++ 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);
Index: trunk/psphot/src/psphotSignificanceImage.c
===================================================================
--- trunk/psphot/src/psphotSignificanceImage.c	(revision 28666)
+++ trunk/psphot/src/psphotSignificanceImage.c	(revision 28667)
@@ -76,5 +76,11 @@
     // Calculate correction factor for the covariance produced by the (potentially multiple) smoothing
     psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing
-    float factor = 1.0 / psImageCovarianceCalculateFactor(kernel, readout->covariance);
+    double sum2 = 0.0;                                               // Sum of kernel squared
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sum2 += PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    float factor = 1.0 / (sum2 * psImageCovarianceCalculateFactor(kernel, readout->covariance));
     psFree(kernel);