Index: trunk/psLib/src/imageops/psImageConvolve.c
===================================================================
--- trunk/psLib/src/imageops/psImageConvolve.c	(revision 26641)
+++ trunk/psLib/src/imageops/psImageConvolve.c	(revision 26892)
@@ -246,4 +246,84 @@
     return kernel;
 }
+
+bool psKernelTruncate(psKernel *kernel, float frac)
+{
+    PS_ASSERT_KERNEL_NON_NULL(kernel, false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(frac, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN(frac, 1.0, false);
+
+    if (frac == 0.0) {
+        // Nothing to do
+        return true;
+    }
+
+    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
+    int maxSize = PS_MAX(PS_MAX(PS_MAX(-xMin, xMax), -yMin), yMax); // Maximum size
+
+    // Determine the threshold
+    // Summing absolute values because large negative deviations have power as well
+    double sumKernel = 0.0;             // Sum of the kernel
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sumKernel += fabsf(kernel->kernel[y][x]);
+        }
+    }
+
+    float threshold = sumKernel * (1.0 - frac); // Threshold for truncation
+
+    // Find truncation size
+    int truncate = 0;                     // Truncation radius
+    for (int radius = 1; truncate == 0 && radius < maxSize; radius++) {
+        int uMin = PS_MAX(-radius, xMin);
+        int uMax = PS_MIN(radius, xMax);
+        int vMin = PS_MAX(-radius, yMin);
+        int vMax = PS_MIN(radius, yMax);
+        int r2 = PS_SQR(radius);
+        double sum = 0.0;
+        for (int v = vMin; v <= vMax; v++) {
+            int v2 = PS_SQR(v);
+            for (int u = uMin; u <= uMax; u++) {
+                int u2 = PS_SQR(u);
+                if (u2 + v2 <= r2) {
+                    sum += fabsf(kernel->kernel[v][u]);
+                }
+            }
+        }
+        if (sum > threshold) {
+            // This is the truncation radius
+            truncate = radius;
+        }
+    }
+    if (truncate == maxSize) {
+        // No truncation possible
+        return true;
+    }
+
+    // Truncate the kernel
+    {
+        int uMin = PS_MAX(-truncate, xMin);
+        int uMax = PS_MIN(truncate, xMax);
+        int vMin = PS_MAX(-truncate, yMin);
+        int vMax = PS_MIN(truncate, yMax);
+        int r2 = PS_SQR(truncate);
+        for (int v = vMin; v <= vMax; v++) {
+            int v2 = PS_SQR(v);
+            for (int u = uMin; u <= uMax; u++) {
+                int u2 = PS_SQR(u);
+                if (u2 + v2 > r2) {
+                    kernel->kernel[v][u] = 0.0;
+                }
+            }
+        }
+    }
+    kernel->xMin = PS_MAX(-truncate, kernel->xMin);
+    kernel->xMax = PS_MIN(truncate, kernel->xMax);
+    kernel->yMin = PS_MAX(-truncate, kernel->yMin);
+    kernel->yMax = PS_MIN(truncate, kernel->yMax);
+
+    return true;
+    }
+
+
 
 psImage *psImageConvolveDirect(psImage *out, const psImage *in, const psKernel *kernel)
Index: trunk/psLib/src/imageops/psImageConvolve.h
===================================================================
--- trunk/psLib/src/imageops/psImageConvolve.h	(revision 26641)
+++ trunk/psLib/src/imageops/psImageConvolve.h	(revision 26892)
@@ -138,4 +138,14 @@
 );
 
+/// Truncate a kernel
+///
+/// Truncates the outer parts of the kernel where the contribution is below the nominated fraction of the
+/// total kernel.
+bool psKernelTruncate(
+    psKernel *in,                       ///< Kernel to be truncated
+    float frac                          ///< Fraction for truncation threshold
+    );
+
+
 /// Convolve an image with a kernel, using a direct convolution
 ///
Index: trunk/psLib/src/imageops/psImageCovariance.c
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.c	(revision 26641)
+++ trunk/psLib/src/imageops/psImageCovariance.c	(revision 26892)
@@ -14,6 +14,11 @@
 #include "psTrace.h"
 #include "psBinaryOp.h"
+#include "psScalar.h"
+#include "psThread.h"
 
 #include "psImageCovariance.h"
+
+static bool threaded = false;           // Run threaded?
+
 
 psKernel *psImageCovarianceNone(void)
@@ -24,54 +29,12 @@
 }
 
-
-psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
-{
-    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
-
-    // See http://en.wikipedia.org/wiki/Error_propagation
-    //
-    // If
-    //     f_k = sum_i A_ik x_i
-    // is a set of functions, then the covariance matrix for f is given by:
-    //     M^f_ij = sum_k sum_l A_ik M^x_kl A_lj
-    // 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.
-
-    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 NULL;
-            }
-        }
-    }
-    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 above (double) sum for the covariance matrix means that, for each point in the output covariance
-    // matrix, we need to work out all combinations of getting to the central point via a kernel, input
-    // covariance matrix and another kernel.  This means that the resultant covariance matrix has twice the
-    // size of the kernel plus the size of the input covariance matrix.
-    int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax;
-    int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
-    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
+/// Calculation of covariance matrix element when convolving
+static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix
+                                      const psKernel *kernel, // Convolution kernel
+                                      int x, int y            // Coordinates in output covariance matrix
+                                      )
+{
+    psAssert(covar, "Require covariance matrix");
+    psAssert(kernel, "Require kernel");
 
     // Need to go:
@@ -89,44 +52,234 @@
     // from the source coordinate), we take the smallest possible (because everything else is zero outside).
 
-    double total = 0.0;                 // Total covariance
+    // 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;
+        }
+    }
+
+    return sum;
+}
+
+/// Thread entry point for calculation of covariance matrix element when convolving
+static bool imageCovarianceCalculateThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "No job arguments");
+    psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n);
+
+    psKernel *out = job->args->data[0]; // Output covariance matrix
+    const psKernel *covar = job->args->data[1]; // Input covariance matrix
+    const psKernel *kernel = job->args->data[2]; // Convolution kernel
+    int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix
+    int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix
+
+    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
+
+    return true;
+}
+
+
+
+psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
+{
+    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
+
+    // See http://en.wikipedia.org/wiki/Error_propagation
+    //
+    // If
+    //     f_k = sum_i A_ik x_i
+    // is a set of functions, then the covariance matrix for f is given by:
+    //     M^f_ij = sum_k sum_l A_ik M^x_kl A_lj
+    // 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.
+
+    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 NULL;
+            }
+        }
+    }
+    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 above (double) sum for the covariance matrix means that, for each point in the output covariance
+    // matrix, we need to work out all combinations of getting to the central point via a kernel, input
+    // covariance matrix and another kernel.  This means that the resultant covariance matrix has twice the
+    // size of the kernel plus the size of the input covariance matrix.
+    int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax;
+    int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
+    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
+
     for (int y = yMin; y <= yMax; y++) {
-        // 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);
         for (int x = xMin; x <= xMax; x++) {
-            // 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;
+            if (threaded) {
+                psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE");
+                psArrayAdd(job->args, 1, out);
+                psArrayAdd(job->args, 1, covar);
+                psArrayAdd(job->args, 1, (psKernel*)kernel); // Casting away const
+                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
+                if (!psThreadJobAddPending(job)) {
+                    psFree(covar);
+                    return NULL;
                 }
-            }
-            out->kernel[y][x] = sum;
-            total += sum;
-        }
-    }
-    psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
-
+                psFree(job);
+            } else {
+                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
+            }
+        }
+    }
     psFree(covar);
+
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        return false;
+    }
+
     return out;
 }
+
+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;
+            }
+        }
+    }
+
+    float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
+    psFree(covar);
+    return factor;
+}
+
+// Calculation of covariance matrix element when binning
+static float imageCovarianceBin(const psKernel *covar, // Original covariance matrix
+                                int bin,               // Binning factor
+                                float binVal,          // Convolution kernel value for binning
+                                int x, int y           // Coordinates in output covariance matrix
+    )
+{
+    psAssert(covar, "Require covariance matrix");
+    psAssert(bin > 0 && binVal > 0, "Require binning: %d %f", bin, binVal);
+
+    int binMin = -(bin - 1) / 2, binMax = bin / 2; // Range of "kernel"
+
+    // Range for v
+    int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
+    int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
+    // 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;
+        }
+    }
+
+    return sum;
+}
+
+/// Thread entry point for calculation of covariance matrix element when binning
+static bool imageCovarianceBinThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "No job arguments");
+    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
+
+    psKernel *out = job->args->data[0]; // Output covariance matrix
+    const psKernel *covar = job->args->data[1]; // Input covariance matrix
+    int bin = PS_SCALAR_VALUE(job->args->data[2], S32); // Binning factor
+    float binVal = PS_SCALAR_VALUE(job->args->data[3], F32); // Convolution kernel value for binning
+    int x = PS_SCALAR_VALUE(job->args->data[4], S32); // x coordinate in output covariance matrix
+    int y = PS_SCALAR_VALUE(job->args->data[5], S32); // y coordinate in output covariance matrix
+
+    out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
+
+    return true;
+}
+
 
 psKernel *psImageCovarianceBin(int bin, const psKernel *covariance, bool average)
@@ -165,43 +318,31 @@
     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;
+            if (threaded) {
+                psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_BIN");
+                psArrayAdd(job->args, 1, out);
+                psArrayAdd(job->args, 1, covar);
+                PS_ARRAY_ADD_SCALAR(job->args, bin, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, binVal, PS_TYPE_F32);
+                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
+                if (!psThreadJobAddPending(job)) {
+                    psFree(covar);
+                    return NULL;
                 }
-            }
-            out->kernel[y][x] = sum;
-            total += sum;
-        }
-    }
-    psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
-
+                psFree(job);
+            } else {
+                out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
+            }
+        }
+    }
     psFree(covar);
 
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        return false;
+    }
+
     return out;
 }
@@ -219,15 +360,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];
+        }
     }
 
@@ -428,2 +569,32 @@
     return true;
 }
+
+
+bool psImageCovarianceSetThreads(bool set)
+{
+    bool old = threaded;                // Old value
+    if (set && !threaded) {
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
+            task->function = &imageCovarianceCalculateThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_BIN", 6);
+            task->function = &imageCovarianceBinThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+    } else if (!set && threaded) {
+        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_CALCULATE");
+        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_BIN");
+    }
+    threaded = set;
+    return old;
+}
+
+bool psImageCovarianceGetThreads(void)
+{
+    return threaded;
+}
Index: trunk/psLib/src/imageops/psImageCovariance.h
===================================================================
--- trunk/psLib/src/imageops/psImageCovariance.h	(revision 26641)
+++ trunk/psLib/src/imageops/psImageCovariance.h	(revision 26892)
@@ -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);
@@ -81,4 +90,12 @@
     );
 
+/// Control threading for image covariance functions
+///
+/// Returns old threading status
+bool psImageCovarianceSetThreads(bool threaded ///< Run image covariance functions threaded?
+    );
+
+/// Return whether image covariance functions are threaded
+bool psImageCovarianceGetThreads(void);
 
 /// @}
Index: trunk/psLib/src/imageops/psImageInterpolate.c
===================================================================
--- trunk/psLib/src/imageops/psImageInterpolate.c	(revision 26641)
+++ trunk/psLib/src/imageops/psImageInterpolate.c	(revision 26892)
@@ -197,10 +197,10 @@
 {
     // Casting away const
-    psFree((psImage*)interp->image);
-    psFree((psImage*)interp->mask);
-    psFree((psImage*)interp->variance);
-    psFree((psImage*)interp->kernel);
-    psFree((psImage*)interp->kernel2);
-    psFree((psVector*)interp->sumKernel2);
+    psFree(interp->image);
+    psFree(interp->mask);
+    psFree(interp->variance);
+    psFree(interp->kernel);
+    psFree(interp->kernel2);
+    psFree(interp->sumKernel2);
 }
 
Index: trunk/psLib/src/math/psBinaryOp.c
===================================================================
--- trunk/psLib/src/math/psBinaryOp.c	(revision 26641)
+++ trunk/psLib/src/math/psBinaryOp.c	(revision 26892)
@@ -379,5 +379,5 @@
 }
 
-psMathType* psBinaryOp(psPtr out, const psPtr in1, const char *op, const psPtr in2)
+psMathType* psBinaryOp(psPtr out, psPtr in1, const char *op, psPtr in2)
 {
 
Index: trunk/psLib/src/math/psBinaryOp.h
===================================================================
--- trunk/psLib/src/math/psBinaryOp.h	(revision 26641)
+++ trunk/psLib/src/math/psBinaryOp.h	(revision 26892)
@@ -55,7 +55,7 @@
 psMathType* psBinaryOp(
     psPtr out,                         ///< Output type, either psImage or psVector.
-    const psPtr in1,                   ///< First input, either psImage or psVector.
+    psPtr in1,                   ///< First input, either psImage or psVector.
     const char *op,                    ///< Operator.
-    const psPtr in2                    ///< Second input, either psImage or psVector.
+    psPtr in2                    ///< Second input, either psImage or psVector.
 );
 
Index: trunk/psLib/src/math/psHistogram.c
===================================================================
--- trunk/psLib/src/math/psHistogram.c	(revision 26641)
+++ trunk/psLib/src/math/psHistogram.c	(revision 26892)
@@ -127,5 +127,5 @@
 static void histogramFree(psHistogram* myHist)
 {
-    psFree((void *)myHist->bounds);
+    psFree(myHist->bounds);
     psFree(myHist->nums);
 }
Index: trunk/psLib/src/math/psMatrix.c
===================================================================
--- trunk/psLib/src/math/psMatrix.c	(revision 26641)
+++ trunk/psLib/src/math/psMatrix.c	(revision 26892)
@@ -1033,5 +1033,5 @@
 }
 
-psVector *psMatrixSolveSVD(psVector *out, const psImage *matrix, const psVector *vector)
+psVector *psMatrixSolveSVD(psVector *out, const psImage *matrix, const psVector *vector, float thresh)
 {
     #define psMatrixSolveSVD_EXIT {psFree(out); return NULL; }
@@ -1063,4 +1063,29 @@
     gsl_vector_free(work);
 
+    if (isfinite(thresh) && thresh > 0.0) {
+        // Trim the singular values
+        double total = 0.0;             // Total of singular values
+        for (int i = 0; i < numCols; i++) {
+            total += gsl_vector_get(S, i);
+        }
+        thresh *= total;
+        for (int i = 0; i < numCols; i++) {
+            double value = gsl_vector_get(S, i); // Singular value
+            if (value < thresh) {
+                psTrace("psLib.math", 5, "Trimming singular value %d: %lg", i, value);
+                gsl_vector_set(S, i, 0.0);
+#if 0
+                for (int j = 0; j < numCols; j++) {
+                    // Being thorough; probably unnecessary
+                    gsl_matrix_set(V, j, i, 0.0);
+                    gsl_matrix_set(A, j, i, 0.0);
+                }
+#endif
+            } else {
+                psTrace("psLib.math", 5, "Singular value %d: %lg", i, value);
+            }
+        }
+    }
+
     // Solve system (or minimise least-squares if overconstrained): Ax = b
     gsl_vector *b = gsl_vector_alloc(numCols); // Vector b
@@ -1093,8 +1118,6 @@
 }
 
-
-
 // This code supplied by Andy Becker (becker@astro.washington.edu)
-psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in)
+psImage *psMatrixSVD_old(psImage* evec, psVector* eval, const psImage* in)
 {
     #define psMatrixSVD_EXIT {psFree(evec); psFree(eval); return NULL;}
@@ -1145,2 +1168,56 @@
     return evec;
 }
+
+// this is basically a wrapper for the gsl function: gsl_linalg_SV_decomp() SVD decomposes
+// matrix A based on the following equation: A = U w V^T .  This function (as usual for SVD
+// implementations) returns V not V^T.  U and V are returned to images; w is returned to a
+// vector representing the diagonal of w.  The input image A is not modified.  U, V, and w may
+// be supplied as NULL or may be allocated; their lengths are set here to match the
+// dimensionality of A.  XXX there is no error handling for the gsl functions (anywhere in
+// psMatrix.c)
+bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A)
+{
+    // Error checks  Missing one for eval
+    PS_ASSERT_PTR_NON_NULL(U, false);
+    PS_ASSERT_PTR_NON_NULL(w, false);
+    PS_ASSERT_PTR_NON_NULL(V, false);
+    PS_ASSERT_PTR_NON_NULL(A, false);
+
+    // A is provided with size Nx,Ny = numCols,numRows
+    // U has size Nx,Ny
+    // V has size Nx,Nx
+    // w has size Nx
+
+    // Initialize data
+    int numRows = A->numRows;
+    int numCols = A->numCols;
+
+    *U = psImageRecycle(*U,  numCols, numRows, A->type.type);
+    *V = psImageRecycle(*V,  numCols, numCols, A->type.type);
+    *w = psVectorRecycle(*w, numCols, A->type.type);
+
+    gsl_matrix *Agsl = gsl_matrix_alloc(numRows, numCols);
+    gsl_matrix *Vgsl = gsl_matrix_alloc(numCols, numCols);
+    gsl_vector *Sgsl = gsl_vector_alloc(numCols);
+    gsl_vector *work = gsl_vector_alloc(numCols);
+
+    // Copy psImage data into GSL matrix data
+    matrixPStoGSL(Agsl, A);
+
+    // Calculate SVD decomposition
+    gsl_linalg_SV_decomp(Agsl, Vgsl, Sgsl, work);
+
+    // Copy GSL matrix data to psImage data
+    matrixGSLtoPS(*V, Vgsl);
+    matrixGSLtoPS(*U, Agsl);  // gsl_linalg_SV_decomp replaces A with U
+    vectorGSLtoPS(*w, Sgsl);
+
+    // Free GSL data
+    gsl_matrix_free(Agsl);
+    gsl_matrix_free(Vgsl);
+    gsl_vector_free(Sgsl);
+    gsl_vector_free(work);
+
+    return true;
+}
+
Index: trunk/psLib/src/math/psMatrix.h
===================================================================
--- trunk/psLib/src/math/psMatrix.h	(revision 26641)
+++ trunk/psLib/src/math/psMatrix.h	(revision 26892)
@@ -192,10 +192,10 @@
     psVector *solution,                 ///< Solution to output, or NULL
     const psImage *matrix,              ///< Matrix to be solved
-    const psVector *vector              ///< Vector of values
+    const psVector *vector,             ///< Vector of values
+    float thresh                        ///< Threshold relative to maximum for trimming singular values
     );
 
-
-/// Single value decomposition, provided by Andy Becker
-psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in);
+/// Single value decomposition (original by Andy Becker, updated by EAM)
+bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A);
 
 /// @}
Index: trunk/psLib/src/math/psUnaryOp.c
===================================================================
--- trunk/psLib/src/math/psUnaryOp.c	(revision 26641)
+++ trunk/psLib/src/math/psUnaryOp.c	(revision 26892)
@@ -211,5 +211,5 @@
 }
 
-psMathType* psUnaryOp(psPtr out, const psPtr in, const char *op)
+psMathType* psUnaryOp(psPtr out, psPtr in, const char *op)
 {
     #define psUnaryOp_EXIT { \
Index: trunk/psLib/src/math/psUnaryOp.h
===================================================================
--- trunk/psLib/src/math/psUnaryOp.h	(revision 26641)
+++ trunk/psLib/src/math/psUnaryOp.h	(revision 26892)
@@ -59,5 +59,5 @@
 psMathType* psUnaryOp(
     psPtr out,                         ///< Output type, either psImage or psVector.
-    const psPtr in,                    ///< Input, either psImage or psVector.
+    psPtr in,			       ///< Input, either psImage or psVector.
     const char *op                     ///< Operator.
 );
Index: trunk/psLib/src/mathtypes/psVector.c
===================================================================
--- trunk/psLib/src/mathtypes/psVector.c	(revision 26641)
+++ trunk/psLib/src/mathtypes/psVector.c	(revision 26892)
@@ -729,5 +729,5 @@
     char line[1024];
 
-    sprintf (line, "vector: %s\n", name);
+    sprintf (line, "# vector: %s\n", name);
     if (write(fd, line, strlen(line))) {;} //ignore return value
 
Index: trunk/psLib/src/sys/psErrorCodes.c.in
===================================================================
--- trunk/psLib/src/sys/psErrorCodes.c.in	(revision 26641)
+++ trunk/psLib/src/sys/psErrorCodes.c.in	(revision 26892)
@@ -67,5 +67,5 @@
 static void freeErrorDescription(psErrorDescription* err)
 {
-    psFree((void *)err->description);
+    psFree(err->description);
 }
 
Index: trunk/psLib/src/sys/psMemory.h
===================================================================
--- trunk/psLib/src/sys/psMemory.h	(revision 26641)
+++ trunk/psLib/src/sys/psMemory.h	(revision 26892)
@@ -326,6 +326,7 @@
 
 /** Free memory.  This operates much like free().
- *
+ *  
  *  @see psAlloc, psRealloc
+ *  note: we cast ptr to (void *) in case we are supplied a const pointer.
  */
 #ifdef DOXYGEN
@@ -336,5 +337,5 @@
 #ifndef SWIG
 #define psFree(ptr) \
-        psMemDecrRefCounter(ptr)
+    ptr = psMemDecrRefCounter((void *)ptr);
 #endif // ifndef SWIG
 #endif // ifdef DOXYGEN
Index: trunk/psLib/src/sys/psTrace.c
===================================================================
--- trunk/psLib/src/sys/psTrace.c	(revision 26641)
+++ trunk/psLib/src/sys/psTrace.c	(revision 26892)
@@ -119,5 +119,5 @@
 
     psMemSetPersistent((psPtr)comp->name,false);
-    psFree((void *)comp->name);
+    psFree(comp->name);
 }
 
Index: trunk/psLib/src/types/psArray.c
===================================================================
--- trunk/psLib/src/types/psArray.c	(revision 26641)
+++ trunk/psLib/src/types/psArray.c	(revision 26892)
@@ -166,5 +166,5 @@
 // drop an item from the array and free it
 bool psArrayRemoveData(psArray* array,
-                       const psPtr data)
+                       psPtr data)
 {
     PS_ASSERT_ARRAY_NON_NULL(array, false);
Index: trunk/psLib/src/types/psArray.h
===================================================================
--- trunk/psLib/src/types/psArray.h	(revision 26641)
+++ trunk/psLib/src/types/psArray.h	(revision 26892)
@@ -186,5 +186,5 @@
 bool psArrayRemoveData(
     psArray* array,                    ///< array to operate on
-    const psPtr data                   ///< the data pointer to remove from psArray
+    psPtr data			       ///< the data pointer to remove from psArray
 );
 
Index: trunk/psLib/src/types/psMetadata.c
===================================================================
--- trunk/psLib/src/types/psMetadata.c	(revision 26641)
+++ trunk/psLib/src/types/psMetadata.c	(revision 26892)
@@ -622,5 +622,6 @@
 
 // may need to extend this to change the keyname in the copy
-bool psMetadataItemSupplement(psMetadata *out,
+bool psMetadataItemSupplement(bool *status,
+                              psMetadata *out,
                               const psMetadata *in,
                               const char *key)
@@ -632,5 +633,9 @@
     psMetadataItem *item = psMetadataLookup(in, key);
     if (!item) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key);
+        if (status) {
+            *status = false;
+        } else {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key);
+        }
         return false;
     }
@@ -968,5 +973,12 @@
 }
 
-psMetadataItem* psMetadataLookup(const psMetadata *md,
+// Get entry by index
+psMetadataItem *psMetadataGetIndex(psMetadata *md, long location)
+{
+    PS_ASSERT_METADATA_NON_NULL(md, false);
+    return psListGet(md->list, location);
+}
+
+psMetadataItem *psMetadataLookup(const psMetadata *md,
                                  const char *key)
 {
Index: trunk/psLib/src/types/psMetadata.h
===================================================================
--- trunk/psLib/src/types/psMetadata.h	(revision 26641)
+++ trunk/psLib/src/types/psMetadata.h	(revision 26892)
@@ -545,4 +545,5 @@
  */
 bool psMetadataItemSupplement(
+    bool *status,			///< if supplied, returns true/false if key is found (suppresses the error)
     psMetadata *out,                   ///< output Metadata container for copying.
     const psMetadata *in,              ///< Metadata collection from which to copy.