Index: branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.c
===================================================================
--- branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.c	(revision 26703)
+++ branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.c	(revision 26704)
@@ -24,4 +24,69 @@
 }
 
+float psImageCovarianceCalculateFactor(const psKernel *kernel, const psKernel *covariance)
+{
+    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 = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            if (!isfinite(kernel->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+                psFree(covar);
+                return NAN;
+            }
+        }
+    }
+    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 NAN;
+            }
+        }
+    }
+
+    int x = 0, y = 0;                   // Coordinates on the output covariance matrix
+
+    // Range for v
+    int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin);
+    int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax);
+    // Range for u
+    int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin);
+    int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax);
+
+    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, kernel->yMin);
+        int qMax = PS_MIN(v + covar->yMax, kernel->yMax);
+        for (int u = uMin; u <= uMax; u++) {
+            // Range for p
+            int pMin = PS_MAX(u + covar->xMin, kernel->xMin);
+            int pMax = PS_MIN(u + covar->xMax, kernel->xMax);
+
+            double xyuvValue = kernel->kernel[v-y][u-x]; // 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)kernel->kernel[q][p];
+                }
+            }
+            sum += xyuvValue * uvpqValue;
+        }
+    }
+
+    psFree(covar);
+    return sum;
+}
 
 psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
@@ -219,15 +284,15 @@
 
     for (int y = covar->yMin; y <= covar->yMax; y++) {
-	if (y < -radius) continue;
-	if (y > +radius) continue;
-	for (int x = covar->xMin; x <= covar->xMax; x++) {
-	    if (x < -radius) continue;
-	    if (x > +radius) continue;
-
-	    if (hypot(x, y) > radius) continue;
-
-	    psAssert (isfinite(covar->kernel[y][x]), "invalid NAN in covariance matrix");
-	    Sum += covar->kernel[y][x];
-	}
+        if (y < -radius) continue;
+        if (y > +radius) continue;
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            if (x < -radius) continue;
+            if (x > +radius) continue;
+
+            if (hypot(x, y) > radius) continue;
+
+            psAssert (isfinite(covar->kernel[y][x]), "invalid NAN in covariance matrix");
+            Sum += covar->kernel[y][x];
+        }
     }
 
Index: branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.h
===================================================================
--- branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.h	(revision 26703)
+++ branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.h	(revision 26704)
@@ -47,4 +47,13 @@
     );
 
+/// Return the pixel-to-pixel covariance factor following calculation
+///
+/// This doesn't require calculation of the entire covariance matrix, so is much faster.
+float psImageCovarianceCalculateFactor(
+    const psKernel *kernel,             ///< Convolution kernel
+    const psKernel *covariance          ///< Current covariance pseudo-matrix
+    );
+
+
 /// Return the covariance factor for an aperture of a given radius
 float psImageCovarianceFactorForAperture(const psKernel *covar, float radius);