Changeset 26551

Timestamp:

Jan 8, 2010, 4:45:02 PM (16 years ago)

Author:

Paul Price

Message:

Allow trimming of singular values when solving using SVD.

Location:

branches/eam_branches/20091201/psLib/src/math

Files:

• psMatrix.c (modified) (3 diffs)
• psMatrix.h (modified) (1 diff)

branches/eam_branches/20091201/psLib/src/math/psMatrix.c

@@ -1033 +1033 @@
 }
 
-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 +1063 @@
     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 +1216 @@
     return true;
 }
+

branches/eam_branches/20091201/psLib/src/math/psMatrix.h

@@ -192 +192 @@
     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)