Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26346)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26347)
@@ -973,6 +973,79 @@
 	}
 
+# 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