Index: trunk/psLib/src/imageops/psImageCovariance.c
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.c	(revision 21207)
+++ trunk/psLib/src/imageops/psImageCovariance.c	(revision 21280)
@@ -11,4 +11,11 @@
 
 #include "psImageCovariance.h"
+
+psKernel *psImageCovarianceNone(void)
+{
+    psKernel *covar = psKernelAlloc(0, 0, 0, 0); // Covariance pseudo-matrix
+    covar->kernel[0][0] = 1.0;
+    return covar;
+}
 
 
@@ -25,12 +32,4 @@
     // where M^x is the covariance matrix for x.
     // Note that the errors in f are correlated (covariance) even if the errors in x are not.
-    //
-    // We don't carry the entire covariance matrix for an image (the size goes as N^2, for N pixels, which
-    // makes storage difficult; and if that's not enough, the time to do the calculation is definitely
-    // impractical).  Since there are (generally) lots of zeros in the covariance matrix, and the same basic
-    // pattern repeats (for background pixels), we can just carry that pattern.  We carry this in a psKernel,
-    // since the values are the covariance between the pixel of consideration (at 0,0 in the kernel) and
-    // neighbouring pixels.  Note that this may not be strictly correct near sources, but this is the best we
-    // can do (and much better than most currently do).
 
     psKernel *covar;                    // Covariance matrix to use
@@ -108,5 +107,44 @@
 float psImageCovarianceFactor(const psKernel *covariance)
 {
-    return covariance ? covariance->kernel[0][0] : 1.0;
+    return covariance ? covariance->kernel[0][0] : NAN;
 }
 
+psKernel *psImageCovarianceAverage(const psArray *array)
+{
+    PS_ASSERT_ARRAY_NON_NULL(array, NULL);
+    PS_ASSERT_ARRAY_NON_EMPTY(array, NULL);
+
+    int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance
+    int num = 0;                        // Number of good matrices to average
+    for (int i = 0; i < array->n; i++) {
+        psKernel *covar = array->data[i]; // Covariance matrix
+        if (!covar) {
+            continue;
+        }
+        xMin = PS_MIN(xMin, covar->xMin);
+        xMax = PS_MAX(xMax, covar->xMax);
+        yMin = PS_MIN(yMin, covar->yMin);
+        yMax = PS_MIN(yMax, covar->yMax);
+        num++;
+    }
+    if (num == 0) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for averaging.");
+        return NULL;
+    }
+
+    psKernel *average = psKernelAlloc(xMin, xMax, yMin, yMax); // Average covariance
+    for (int i = 0; i < array->n; i++) {
+        psKernel *covar = array->data[i]; // Covariance matrix
+        if (!covar) {
+            continue;
+        }
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                average->kernel[y][x] += covar->kernel[y][x];
+            }
+        }
+    }
+    psBinaryOp(average->image, average->image, "/", psScalarAlloc(num, PS_TYPE_F32));
+
+    return average;
+}