Index: trunk/psLib/src/imageops/psImageCovariance.c
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.c	(revision 21298)
+++ trunk/psLib/src/imageops/psImageCovariance.c	(revision 21316)
@@ -5,4 +5,5 @@
 #include <stdio.h>
 #include <limits.h>
+#include <string.h>
 
 #include "psMemory.h"
@@ -172,2 +173,42 @@
     return sum;
 }
+
+
+psKernel *psImageCovarianceTruncate(const psKernel *covar, float frac)
+{
+    PS_ASSERT_KERNEL_NON_NULL(covar, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(frac, 0, 1, NULL);
+
+    int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax; // Range
+
+    double sum = 0.0;                   // Sum of covariance
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sum += covar->kernel[y][x];
+        }
+    }
+
+    // Determine radius for truncation
+    double threshold = (1.0 - frac) * sum; // Threshold for truncation
+    int maxRadius = PS_MAX(PS_MAX(PS_MAX(xMax, yMax), -xMin), yMin); // Maximum radius of covariance matrix
+    double enclosed = NAN;              // Enclosed value
+    int radius;                         // Radius at which to truncate
+    for (radius = 0; radius <= maxRadius && enclosed < threshold; radius++) {
+        // Not so efficient to execute, but easy to code; it should be fast anyway
+        enclosed = 0.0;
+        for (int y = PS_MAX(yMin, -radius); y <= PS_MIN(yMax, radius); y++) {
+            for (int x = PS_MAX(xMin, -radius); x <= PS_MIN(xMax, radius); x++) {
+                enclosed += covar->kernel[y][x];
+            }
+        }
+    }
+
+    // Generate truncated version
+    psKernel *trunc = psKernelAlloc(-radius, radius, -radius, radius); // Truncated covariance matrix
+    int numBytes = (2 * radius + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
+    for (int y = -radius; y <= radius; y++) {
+        memcpy(trunc->kernel[y], covar->kernel[y], numBytes);
+    }
+
+    return trunc;
+}
Index: trunk/psLib/src/imageops/psImageCovariance.h
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.h	(revision 21298)
+++ trunk/psLib/src/imageops/psImageCovariance.h	(revision 21316)
@@ -5,6 +5,6 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-02-05 00:42:48 $
+ * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-05 20:15:32 $
  * Copyright 2009 Institute for Astronomy, University of Hawaii
  */
@@ -21,8 +21,8 @@
 // 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).
+// background pixels), we can just carry that pattern ("covariance pseudo-matrix").  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).
 
 /// Allocate a covariance pseudo-matrix with no covariance
@@ -50,4 +50,12 @@
     );
 
+/// Truncate covariance pseudo-matrix
+///
+/// The covariance pseudo-matrix is truncated by removing the outer regions that contribute less than the
+/// nominated fraction of the total.
+psKernel *psImageCovarianceTruncate(
+    const psKernel *covar,              ///< Covariance pseudo-matrix
+    float frac                          ///< Fraction of covariance to truncate
+    );
 
 /// @}