Index: trunk/psLib/src/imageops/psImageCovariance.c
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.c	(revision 27702)
+++ trunk/psLib/src/imageops/psImageCovariance.c	(revision 28006)
@@ -11,4 +11,6 @@
 #include "psMemory.h"
 #include "psConstants.h"
+#include "psImageStructManip.h"
+#include "psImagePixelManip.h"
 #include "psImageConvolve.h"
 #include "psTrace.h"
@@ -16,9 +18,9 @@
 #include "psScalar.h"
 #include "psThread.h"
+#include "psImageInterpolate.h"
 
 #include "psImageCovariance.h"
 
 static bool threaded = false;           // Run threaded?
-
 
 psKernel *psImageCovarianceNone(void)
@@ -530,4 +532,62 @@
 
 
+psKernel *psImageCovarianceScale(const psKernel *in, float scale)
+{
+    // Trivial cases
+    if (!in) {
+        psKernel *out = psKernelAlloc(0, 0, 0, 0); // Output covariance
+        out->kernel[0][0] = 1.0;
+        return out;
+    }
+    PS_ASSERT_KERNEL_NON_NULL(in, NULL);
+    if (scale == 1.0) {
+        psImage *copy = psImageCopy(NULL, in->image, PS_TYPE_F32); // Copy of input covariance
+        psKernel *out = psKernelAllocFromImage(copy, -in->xMin, -in->yMin); // Output covariance
+        psFree(copy);
+        return out;
+    }
+
+    int xMinIn = in->xMin, xMaxIn = in->xMax, yMinIn = in->yMin, yMaxIn = in->yMax; // Input size
+    int xMinOut = (float)xMinIn / scale - 0.5, xMaxOut = (float)xMaxIn / scale + 0.5;     // Output size in x
+    int yMinOut = (float)yMinIn / scale - 0.5, yMaxOut = (float)yMaxIn / scale + 0.5;     // Output size in y
+
+    // Over-fill the covariance matrix so we're not troubled by edge effects
+    psKernel *overfill = psKernelAlloc(xMinIn - 1, xMaxIn + 1, yMinIn - 1, yMaxIn + 1); // Overfilled covar
+    psImageInit(overfill->image, 0.0);
+    int numOverlay = (xMaxIn - xMinIn + 1) * (yMaxIn - yMinIn + 1); // Number of pixels to overlay
+    if (psImageOverlaySection(overfill->image, in->image, 1, 1, "=") != numOverlay) {
+        psError(psErrorCodeLast(), false, "Unable to overfill covariance matrix.");
+        psFree(overfill);
+        return NULL;
+    }
+
+    psImageInterpolation *interp = psImageInterpolationAlloc(PS_INTERPOLATE_BILINEAR, overfill->image,
+                                                             NULL, NULL, 0, NAN, NAN, 0xFF, 0xFF,
+                                                             0.0, 0); // Interpolation
+    psFree(overfill);
+
+    // In transforming the positions, we get +0.5 to account for the centre of the pixels being at 0.5
+    // and +1 to account for the overfill.
+
+    psKernel *out = psKernelAlloc(xMinOut, xMaxOut, yMinOut, yMaxOut); // Output covariance
+    for (int y = yMinOut; y <= yMaxOut; y++) {
+        float yIn = y * scale + 0.5 - yMinIn + 1; // Position on input image (not the kernel)
+        for (int x = xMinOut; x <= xMaxOut; x++) {
+            float xIn = x * scale + 0.5 - xMinIn + 1; // Position on input (not the kernel)
+            double value;                                     // Value on output
+            if (!psImageInterpolate(&value, NULL, NULL, xIn, yIn, interp)) {
+                psError(psErrorCodeLast(), false, "Unable to interpolate kernel.");
+                return false;
+            }
+            out->kernel[y][x] = value;
+        }
+    }
+
+    psFree(interp);
+
+    return out;
+}
+
+
 bool psImageCovarianceSetThreads(bool set)
 {