Index: /trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- /trunk/psLib/src/math/psMinimizePolyFit.c	(revision 6526)
+++ /trunk/psLib/src/math/psMinimizePolyFit.c	(revision 6527)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-24 23:43:15 $
+ *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-03-06 20:57:42 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -41,13 +41,12 @@
 #define PS_VECTOR_GEN_CHEBY_INDEX(VEC, SIZE, TYPE) \
 VEC = psVectorAlloc(SIZE, TYPE); \
+VEC->n = VEC->nalloc; \
 if (TYPE == PS_TYPE_F64) { \
     for (psS32 i = 0 ; i < SIZE ; i++) { \
         VEC->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
-        VEC->n++; \
     }\
 } else if (TYPE == PS_TYPE_F32){ \
     for (psS32 i = 0 ; i < SIZE ; i++) { \
         VEC->data.F32[i] = ((2.0 / ((psF32) (SIZE - 1))) * ((psF32) i)) - 1.0; \
-        VEC->n++; \
     }\
 }\
@@ -309,5 +308,5 @@
     linear equations which can be easily solved.  The resulting vector is the
     coefficients of the Chebyshev polys.
- 
+    
     This method is significantly slower than the standard NR algorithm.  It
     was explicitly requested that we not use the NR algorithm.
@@ -370,4 +369,5 @@
     psImage *A = psImageAlloc(numPoly, numPoly, PS_TYPE_F64);
     psVector *B = psVectorAlloc(numPoly, PS_TYPE_F64);
+    B->n = B->nalloc;
     for (psS32 i = 0 ; i < numPoly ; i++) {
         for (psS32 j = 0 ; j < numPoly ; j++) {
@@ -375,5 +375,4 @@
         }
         B->data.F64[i] = ordPoly->coeff[i];
-        B->n++;
     }
 
@@ -817,4 +816,6 @@
         }
         for (psS32 i = 0; i < nTerm; i++) {
+            if (myPoly->mask[i])
+                continue;
             B->data.F64[i] += f->data.F64[k] * xSums->data.F64[i] * wt;
         }
@@ -823,5 +824,9 @@
         // we must handle masked orders
         for (psS32 i = 0; i < nTerm; i++) {
+            if (myPoly->mask[i])
+                continue;
             for (psS32 j = 0; j < nTerm; j++) {
+                if (myPoly->mask[j])
+                    continue;
                 A->data.F64[i][j] += xSums->data.F64[i + j] * wt;
             }
@@ -829,4 +834,16 @@
     }
 
+    // set masked elements in A,B appropriately
+    for (int i = 0; i < nTerm; i++) {
+        if (!myPoly->mask[i])
+            continue;
+        B->data.F64[i] = 0;
+        for (int j = 0; j < nTerm; j++) {
+            A->data.F64[i][j] = (i == j) ? 1 : 0;
+        }
+    }
+
+
+    // XXX: rel10_ifa used psGaussJordan().  However, this failed tests.  So, I'm using psMatrixLUD().
     if (0) {
         // GaussJordan version
@@ -840,4 +857,5 @@
             for (psS32 k = 0; k < nTerm; k++) {
                 myPoly->coeff[k] = B->data.F64[k];
+                myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
             }
         }
@@ -864,4 +882,5 @@
                 for (psS32 k = 0; k < nTerm; k++) {
                     myPoly->coeff[k] = coeffs->data.F64[k];
+                    // XXX LUD does not give inverse of A
                 }
             }
@@ -1245,4 +1264,6 @@
         for (psS32 n = 0; n < nXterm; n++) {
             for (psS32 m = 0; m < nYterm; m++) {
+                if (myPoly->mask[n][m])
+                    continue;
                 B->data.F64[n+m*nXterm] += f->data.F64[k] * Sums->data.F64[n][m] * wt;
             }
@@ -1251,8 +1272,28 @@
         for (psS32 i = 0; i < nXterm; i++) {
             for (psS32 j = 0; j < nYterm; j++) {
+                if (myPoly->mask[i][j])
+                    continue;
                 for (psS32 n = 0; n < nXterm; n++) {
                     for (psS32 m = 0; m < nYterm; m++) {
+                        if (myPoly->mask[n][m])
+                            continue;
                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
                     }
+                }
+            }
+        }
+    }
+
+    // set masked elements appropriately
+    for (int i = 0; i < nXterm; i++) {
+        for (int j = 0; j < nYterm; j++) {
+            if (!myPoly->mask[i][j])
+                continue;
+            int nx = i+j*nXterm;
+            B->data.F64[nx] = 0;
+            for (int n = 0; n < nXterm; n++) {
+                for (int m = 0; m < nYterm; m++) {
+                    int ny = n+m*nXterm;
+                    A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
                 }
             }
@@ -1269,4 +1310,5 @@
             for (psS32 m = 0; m < nYterm; m++) {
                 myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+                myPoly->coeffErr[n][m] = sqrt(A->data.F64[n+m*nXterm][n+m*nXterm]);
             }
         }
@@ -1723,4 +1765,5 @@
     }
 
+    // XXX: rel10_ifa used psGaussJordan().  However, this failed tests.  So, I'm using psMatrixLUD().
     if (0) {
         // does the solution in place
@@ -1775,5 +1818,6 @@
                             psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
                             myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
-                            myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
+                            // XXX myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
+                            // XXX this is wrong: A is not replaced with inverse(A), as for GaussJordan
                         }
                     }
@@ -2278,4 +2322,5 @@
 
 
+    // XXX: rel10_ifa used psGaussJordan().  However, this failed tests.  So, I'm using psMatrixLUD().
     if (0) {
         // does the solution in place
@@ -2336,5 +2381,6 @@
                                 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                                 myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
-                                myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
+                                // myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
+                                // XXX this is wrong: LUD does not supply inverse(A)
                             }
                         }
Index: /trunk/psLib/test/math/tst_psPolyFit4D.c
===================================================================
--- /trunk/psLib/test/math/tst_psPolyFit4D.c	(revision 6526)
+++ /trunk/psLib/test/math/tst_psPolyFit4D.c	(revision 6527)
@@ -521,4 +521,5 @@
     // All Vectors non-NULL
     testStatus &= genericTest(TS00_MASK_U8 | TS00_FERR_F32 | TS00_X_F32 | TS00_Y_F32 | TS00_Z_F32 | TS00_T_F32 | TS00_F_F32 | TS00_POLY_ORD, POLY_ORDER_X, POLY_ORDER_Y, POLY_ORDER_Z, POLY_ORDER_T, NUM_DATA, true);
+
     // Some Vectors NULL
     testStatus &= genericTest(TS00_MASK_U8 | TS00_FERR_NULL | TS00_X_F32 | TS00_Y_F32 | TS00_Z_F32 | TS00_T_F32 | TS00_F_F32 | TS00_POLY_ORD, POLY_ORDER_X, POLY_ORDER_Y, POLY_ORDER_Z, POLY_ORDER_T, NUM_DATA, true);