Changeset 26347

Timestamp:

Dec 6, 2009, 8:41:50 AM (16 years ago)

Author:

eugene

Message:

adding SVD analysis of basis function correlations

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (modified) (1 diff)

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c

@@ -973 +973 @@
         }
 
+# define SVD_ANALYSIS 1
         // re-do this section with SVD:
-        // psMatrixSVD ();
+        if (SVD_ANALYSIS) {
+            // We have sumVector and sumMatrix.  we are trying to solve the equation: sumMatrix * x = sumVector.
+            // we can use any standard matrix inversion to solve this.  However, the basis
+            // functions in general have substantial correlation, so that the solution may be
+            // somewhat poorly determiend or unstable.  If not numerically ill-conditioned, the
+            // system of equations may be statistically ill-conditioned.  Noise in the image
+            // will drive insignificant but correlated terms in the solution.  To avoid these
+            // problems, we can use SVD to identify numerically unconstrained values and to
+            // avoid statistically badly determined value.
+
+            // A = sumMatrix, B = sumVector
+            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
+            // x = V (1/w) (U^T B)
+            // \sigma_x = sqrt(diag(A^{-1}))
+            // solve for x and A^{-1} to get x & dx
+            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
+            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
+
+            // I get confused by the index values between the image vs matrix usage:  In terms
+            // of the elements of an image A(x,y) = A->data.F32[y][x] = A_x,y, a matrix
+            // multiplication is: A_k,j * B_i,k = C_i,j
+
+            psImage *U = NULL;
+            psImage *V = NULL;
+            psVector *w = NULL;
+            if (!psMatrixSVD (&U, &w, &V, sumMatrix)) {
+                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
+                return NULL;
+            }
+
+            // calculate A_inverse and x_svd
+            // A'_i,j = \sum_k (V_k,j * w_k * Ut_i,k)
+            // but U^T_i,j = U_j,i, so:
+            // A'_i,j = \sum_k (V_k,j * w_k * U_k,i)
+
+            // v1_j = \sum_k (UT_k,j * B_k)
+            // v1_j = w_j * \sum_k (U_j,k * B_k)
+            // x_i = \sum_j V_j,i * w_j * \sum_k (U_j,k * B_k)
+            
+            // A'_i,j = \sum_k (V_k,j * w_k * U_k,i)
+
+            psImage  *Ainv = psImageAlloc(A->numCols, A->numRows, PS_TYPE_F32);
+            psVector *Xsvd = psVectorAlloc(A->numCols, PS_TYPE_F32);
+            for (int iy = 0; iy < A->numRows; iy++) {
+                double vsum = 0.0;
+                for (int ix = 0; ix < A->numCols; ix++) {
+                    double msum = 0.0;
+                    double vsum1 = 0.0;
+                    for (int k = 0; k < A->numCols; k++) {
+                        if (fabs(w->data.F32[k]) < FLT_EPSILON) continue;
+                        msum += V->data.F32[iy][k] * U->data.F32[ix][k] / w->data.F32[k];
+                        vsum1 += U->data.F32[k][iy] * B->data.F32[k];
+                    }
+                    Ainv->data.F32[iy][ix] = sum;
+                    if (fabs(w->data.F32[ix]) < FLT_EPSILON) continue;
+                    vsum += V->data.F32[iy][ix] * vsum1 / w->data.F32[ix];
+                }
+                Xsvd->data.F32[iy] = vsum;
+            }
+
+            // compare Xsvd[k] to dXsvd and zero out w entries that are not significant
+            # define COEFF_SIG 1.0
+            for (int k = 0; k < A->numRows; k++) {
+                float dXsvd = sqrt(Ainv->data.F32[k][k]);
+                if (fabs(Xsvd->data.F32[k] / dXsvd) < COEFF_SIG) {
+                    w->data.F32[k] = 0.0;
+                }
+            }
+            
+            // 
+
+
+        }
 
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition