Index: /branches/eam_branches/20091201/psLib/src/math/psMatrix.c
===================================================================
--- /branches/eam_branches/20091201/psLib/src/math/psMatrix.c	(revision 26550)
+++ /branches/eam_branches/20091201/psLib/src/math/psMatrix.c	(revision 26551)
@@ -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,24 @@
     gsl_vector_free(work);
 
+    if (isfinite(thresh) && thresh > 0.0) {
+        // Trim the singular values
+        thresh *= gsl_vector_get(S, 0);
+        psTrace("psLib.math", 8, "Singular value 0: %lf", gsl_vector_get(S, 0));
+        for (int i = 1; i < numCols; i++) {
+            double value = gsl_vector_get(S, i); // Singular value
+            if (value < thresh) {
+                psTrace("psLib.math", 5, "Trimming singular value %d: %lf", i, value);
+                gsl_vector_set(S, i, 0.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);
+                }
+            } else {
+                psTrace("psLib.math", 5, "Singular value %d: %lf", i, value);
+            }
+        }
+    }
+
     // Solve system (or minimise least-squares if overconstrained): Ax = b
     gsl_vector *b = gsl_vector_alloc(numCols); // Vector b
@@ -1196,2 +1216,3 @@
     return true;
 }
+
Index: /branches/eam_branches/20091201/psLib/src/math/psMatrix.h
===================================================================
--- /branches/eam_branches/20091201/psLib/src/math/psMatrix.h	(revision 26550)
+++ /branches/eam_branches/20091201/psLib/src/math/psMatrix.h	(revision 26551)
@@ -192,7 +192,7 @@
     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 (original by Andy Becker, updated by EAM)