Index: trunk/psLib/src/imageops/psImageCovariance.c
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.c	(revision 23859)
+++ trunk/psLib/src/imageops/psImageCovariance.c	(revision 23882)
@@ -130,4 +130,81 @@
 }
 
+psKernel *psImageCovarianceBin(int bin, const psKernel *covariance, bool average)
+{
+    psKernel *covar;                    // Covariance matrix to use
+    if (covariance) {
+        covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
+    } else {
+        covar = psImageCovarianceNone();
+    }
+
+    // Check for non-finite elements
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            if (!isfinite(covar->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
+                psFree(covar);
+                return NULL;
+            }
+        }
+    }
+
+    // The following calculation is the same as for psImageCovarianceCalculate(), except that the "convolution
+    // kernel" has a constant value, and the output covariance kernel is binned down, so we unbin it in the
+    // ranges for the inputs.
+
+    // Values for "convolution kernel"
+    int binMin = -(bin - 1) / 2, binMax = bin / 2; // Range of "kernel"
+    int binDiff = binMax - binMin;      // Difference between the max and min
+    float binVal = average ? 1.0 / PS_SQR(bin) : 1.0; // Value of "kernel" pixels
+
+    // Size of the output covariance matrix
+    int xMin = (covar->xMin - binDiff) / bin + 0.5, xMax = (covar->xMax + binDiff) / bin + 0.5;
+    int yMin = (covar->yMin - binDiff) / bin + 0.5, yMax = (covar->yMax + binDiff) / bin + 0.5;
+    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
+
+    double total = 0.0;                 // Total covariance
+    for (int y = yMin; y <= yMax; y++) {
+        // Range for v
+        int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
+        int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
+        for (int x = xMin; x <= xMax; x++) {
+            // Range for u
+            int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin);
+            int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax);
+
+            double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
+            for (int v = vMin; v <= vMax; v++) {
+                // Range for q
+                int qMin = PS_MAX(v + covar->yMin, binMin);
+                int qMax = PS_MIN(v + covar->yMax, binMax);
+                for (int u = uMin; u <= uMax; u++) {
+                    // Range for p
+                    int pMin = PS_MAX(u + covar->xMin, binMin);
+                    int pMax = PS_MIN(u + covar->xMax, binMax);
+
+                    double xyuvValue = binVal; // Value for (x,y) --> (u,v)
+
+                    double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
+                    for (int q = qMin; q <= qMax; q++) {
+                        for (int p = pMin; p <= pMax; p++) {
+                            uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal;
+                        }
+                    }
+                    sum += xyuvValue * uvpqValue;
+                }
+            }
+            out->kernel[y][x] = sum;
+            total += sum;
+        }
+    }
+    psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
+
+    psFree(covar);
+
+    return out;
+}
+
 float psImageCovarianceFactor(const psKernel *covariance)
 {
Index: trunk/psLib/src/imageops/psImageCovariance.h
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.h	(revision 23859)
+++ trunk/psLib/src/imageops/psImageCovariance.h	(revision 23882)
@@ -35,4 +35,11 @@
     );
 
+/// Calculate the covariance pseudo-matrix for binning
+psKernel *psImageCovarianceBin(
+    int bin,                            ///< Binning factor
+    const psKernel *covariance,         ///< Current covariance pseudo-matrix
+    bool average                        ///< Averaging pixels when binning?
+    );
+
 /// Return the pixel-to-pixel covariance factor
 float psImageCovarianceFactor(