Index: trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.c	(revision 13396)
+++ trunk/psLib/src/math/psMinimizePolyFit.c	(revision 15254)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-16 16:16:48 $
+ *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-09 19:25:45 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -340,5 +340,5 @@
     psImage *sums = psImageAlloc(numData, numTerms, PS_TYPE_F64);
     for (int i = 0; i < numTerms; i++) {
-        if (myPoly->mask[i]) {
+        if (myPoly->coeffMask[i] & PS_POLY_MASK_BOTH) {
             continue;
         }
@@ -358,5 +358,5 @@
         dataMask = mask->data.U8;
     }
-    psU8 *termMask = myPoly->mask;      // Mask for polynomial terms
+    psU8 *coeffMask = myPoly->coeffMask;      // Mask for polynomial terms
     psF64 *yData = y->data.F64;         // Coordinate data
     psF64 *yErrData = NULL;             // Errors in the coordinate
@@ -380,5 +380,5 @@
 
         for (int i = 0; i < numTerms; i++) {
-            if (termMask[i]) {
+            if (coeffMask[i] & PS_POLY_MASK_BOTH) {
                 matrix[i][i] = 1.0;
                 continue;
@@ -387,5 +387,5 @@
             matrix[i][i] += sumsData[i][k] * sumsData[i][k] * wt; // The diagonal entry
             for (int j = i + 1; j < numTerms; j++) { // The upper diagonal only: we will use symmetry
-                if (termMask[j]) {
+                if (coeffMask[j] & PS_POLY_MASK_BOTH) {
                     continue;
                 }
@@ -533,5 +533,5 @@
         dataMask = mask->data.U8;
     }
-    psU8 *termMask = myPoly->mask;      // Dereferenced version of mask for polynomial terms
+    psU8 *coeffMask = myPoly->coeffMask;      // Dereferenced version of mask for polynomial terms
     psF64 *ordinates = NULL;            // Dereferenced version of ordinate data
     if (x) {
@@ -567,5 +567,5 @@
 
         for (int i = 0; i < nTerm; i++) {
-            if (termMask[i]) {
+            if (coeffMask[i] & PS_POLY_MASK_SET) {
                 matrix[i][i] = 1.0;
                 continue;
@@ -574,5 +574,5 @@
             matrix[i][i] += sums[2 * i] * wt; // The diagonal entry
             for (int j = i + 1; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
-                if (termMask[j]) {
+                if (coeffMask[j] & PS_POLY_MASK_SET) {
                     continue;
                 }
@@ -584,4 +584,32 @@
     }
     psFree(xSums);
+
+    // elements which are masked for fitting need to be subtracted from the vector
+    for (int i = 0; i < nTerm; i++) {
+	if (coeffMask[i] & PS_POLY_MASK_BOTH) {
+	    continue;
+	}
+	for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
+	    if (coeffMask[j] & PS_POLY_MASK_SET) {
+		continue;
+	    }
+	    if (!(coeffMask[j] & PS_POLY_MASK_FIT)) {
+		continue;
+	    }
+	    vector[i] -= matrix[i][j]*myPoly->coeff[j];
+	}
+    }
+    
+    // set the un-fitted and un-set elements to 0 or 1 for pivots
+    for (int i = 0; i < nTerm; i++) {
+	if (coeffMask[i] & PS_POLY_MASK_BOTH) {
+	    for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
+		matrix[i][j] = 0.0;
+		matrix[j][i] = 0.0;
+	    }
+	    matrix[i][i] = 1.0;
+	    continue;
+	}
+    }
 
     if (psTraceGetLevel("psLib.math") >= 4) {
@@ -610,9 +638,9 @@
             // polynomial coefficients.  this is only true for the 1D case
             for (psS32 k = 0; k < nTerm; k++) {
-                myPoly->coeff[k] = B->data.F64[k];
-                myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
-            }
-        }
-
+		if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
+		myPoly->coeff[k] = B->data.F64[k];
+		myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
+            }
+        }
     } else {
         // LUD version of the fit
@@ -633,4 +661,5 @@
             } else {
                 for (psS32 k = 0; k < nTerm; k++) {
+		    if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
                     myPoly->coeff[k] = coeffs->data.F64[k];
                     // XXX LUD does not give inverse of A
@@ -990,5 +1019,5 @@
         psFree(A);
         psFree(B);
-        psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
+        psTrace("psLib.math", 6, "---- %s() End ----\n", __func__);
         return false;
     }
@@ -997,5 +1026,5 @@
     psF64 **matrix = A->data.F64;       // Dereference the least-squares matrix
     psF64 *vector = B->data.F64;        // Dereference the least-squares vector
-    psU8 **termMask = myPoly->mask;     // Dereference mask for polynomial terms
+    psU8 **coeffMask = myPoly->coeffMask;     // Dereference mask for polynomial terms
     psU8 *dataMask = NULL;              // Dereference mask for data
     if (mask) {
@@ -1031,5 +1060,5 @@
             int l = i / nYterm;         // x index
             int m = i % nYterm;         // y index
-            if (termMask[l][m]) {
+            if (coeffMask[l][m] & PS_POLY_MASK_SET) {
                 matrix[i][i] = 1.0;
                 continue;
@@ -1040,5 +1069,5 @@
                 int p = j / nYterm;     // x index
                 int q = j % nYterm;     // y index
-                if (termMask[p][q]) {
+                if (coeffMask[p][q] & PS_POLY_MASK_SET) {
                     continue;
                 }
@@ -1050,4 +1079,53 @@
     }
     psFree(xySums);
+
+    // elements which are masked for fitting need to be subtracted from the vector
+    for (int i = 0; i < nTerm; i++) {
+	int ix = i / nYterm;         // x index
+	int iy = i % nYterm;         // y index
+	if (coeffMask[ix][iy] & PS_POLY_MASK_BOTH) {
+	    continue;
+	}
+	for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
+	    int jx = j / nYterm;         // x index
+	    int jy = j % nYterm;         // y index
+	    if (coeffMask[jx][jy] & PS_POLY_MASK_SET) {
+		continue;
+	    }
+	    if (!(coeffMask[jx][jy] & PS_POLY_MASK_FIT)) {
+		continue;
+	    }
+	    vector[i] -= matrix[i][j]*myPoly->coeff[jx][jy];
+	}
+    }
+    
+    // set the un-fitted and un-set elements to 0 or 1 for pivots
+    for (int i = 0; i < nTerm; i++) {
+	int ix = i / nYterm;         // x index
+	int iy = i % nYterm;         // y index
+	if (coeffMask[ix][iy] & PS_POLY_MASK_BOTH) {
+	    for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
+		matrix[i][j] = 0.0;
+		matrix[j][i] = 0.0;
+	    }
+	    matrix[i][i] = 1.0;
+	    continue;
+	}
+    }
+
+    if (psTraceGetLevel("psLib.math") >= 4) {
+        printf("Least-squares vector:\n");
+        for (int i = 0; i < nTerm; i++) {
+            printf("%f ", B->data.F64[i]);
+        }
+        printf("\n");
+        printf("Least-squares matrix:\n");
+        for (int i = 0; i < nTerm; i++) {
+            for (int j = 0; j < nTerm; j++) {
+                printf("%f ", A->data.F64[i][j]);
+            }
+            printf("\n");
+        }
+    }
 
     if (!psMatrixGJSolve(A, B)) {
@@ -1062,11 +1140,14 @@
         int l = i / nYterm;         // x index
         int m = i % nYterm;         // y index
-        myPoly->coeff[l][m] = B->data.F64[i];
-        myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
+
+	// retain the incoming values if masked on the fit
+	if (coeffMask[l][m] & PS_POLY_MASK_FIT) continue;
+	myPoly->coeff[l][m] = B->data.F64[i];
+	myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
     }
     psFree(A);
     psFree(B);
 
-    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+    psTrace("psLib.math", 6, "---- %s() end ----\n", __func__);
     return true;
 }
@@ -1387,5 +1468,5 @@
         dataMask = mask->data.U8;
     }
-    psU8 ***termMask = myPoly->mask;    // Mask for polynomial terms
+    psU8 ***coeffMask = myPoly->coeffMask;    // Mask for polynomial terms
     int nYZterm = nYterm * nZterm;      // Multiplication of the numbers, to calculate the index
 
@@ -1411,5 +1492,5 @@
             int iy = (i % nYZterm) / nZterm; // y index
             int iz = (i % nYZterm) % nZterm; // z index
-            if (termMask[ix][iy][iz]) {
+            if (coeffMask[ix][iy][iz] & PS_POLY_MASK_BOTH) {
                 matrix[i][i] = 1.0;
                 continue;
@@ -1422,5 +1503,5 @@
                 int jy = (j % nYZterm) / nZterm; // y index
                 int jz = (j % nYZterm) % nZterm; // z index
-                if (termMask[jx][jy][jz]) {
+                if (coeffMask[jx][jy][jz] & PS_POLY_MASK_BOTH) {
                     continue;
                 }
@@ -1455,4 +1536,5 @@
             int iy = (i % nYZterm) / nZterm; // y index
             int iz = (i % nYZterm) % nZterm; // z index
+	    if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
             myPoly->coeff[ix][iy][iz] = B->data.F64[i];
             myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
@@ -1480,4 +1562,5 @@
                         for (psS32 iz = 0; iz < nZterm; iz++) {
                             psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
+			    if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
                             myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
                             // XXX myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
@@ -1839,5 +1922,5 @@
         dataMask = mask->data.U8;
     }
-    psU8 ****termMask = myPoly->mask;    // Mask for polynomial terms
+    psU8 ****coeffMask = myPoly->coeffMask;    // Mask for polynomial terms
     int nYZTterm = nYterm * nZterm * nTterm; // Multiplication of the numbers, for calculating the index
     int nZTterm = nZterm * nTterm;      // Multiplication of the numbers, for calculating the index
@@ -1865,5 +1948,5 @@
             int iz = ((i % (nYZTterm)) % (nZTterm)) / nTterm; // z index
             int it = ((i % (nYZTterm)) % (nZTterm)) % nTterm; // t index
-            if (termMask[ix][iy][iz][it]) {
+            if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_BOTH) {
                 matrix[i][i] = 1.0;
                 continue;
@@ -1877,5 +1960,5 @@
                 int jz = ((j % nYZTterm) % nZTterm) / nTterm; // z index
                 int jt = ((j % nYZTterm) % nZTterm) % nTterm; // t index
-                if (termMask[jx][jy][jz][jt]) {
+                if (coeffMask[jx][jy][jz][jt] & PS_POLY_MASK_BOTH) {
                     continue;
                 }
@@ -1918,4 +2001,5 @@
             int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
             int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
+	    if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
             myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
             myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
@@ -1944,4 +2028,5 @@
                             for (psS32 it = 0; it < nTterm; it++) {
                                 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+				if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
                                 myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
                                 // myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);