Index: trunk/psLib/src/imageops/psImageCovariance.c
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.c	(revision 28405)
+++ 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 28405)
+++ 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);
 }