Index: /trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- /trunk/psLib/src/math/psMinimizePolyFit.c	(revision 7134)
+++ /trunk/psLib/src/math/psMinimizePolyFit.c	(revision 7135)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-05-18 01:21:16 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-05-18 03:05:47 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -379,4 +379,5 @@
         for (int i = 0; i < numTerms; i++) {
             if (termMask[i]) {
+                matrix[i][i] = 1.0;
                 continue;
             }
@@ -560,4 +561,5 @@
         for (int i = 0; i < nTerm; i++) {
             if (termMask[i]) {
+                matrix[i][i] = 1.0;
                 continue;
             }
@@ -575,19 +577,4 @@
     }
     psFree(xSums);
-
-    #if 0
-    // set masked elements in A,B appropriately
-    // PAP: Is this necessary, given we initialise to zero, and we have already masked on the terms?
-    for (int i = 0; i < nTerm; i++) {
-        if (!termMask[i])
-            continue;
-        vector[i] = 0;
-        matrix[i][i] = 1.0;
-        for (int j = i + 1; j < nTerm; j++) {
-            matrix[i][j] = 0.0;
-            matrix[j][i] = 0.0;
-        }
-    }
-    #endif
 
     if (psTraceGetLevel(__func__) >= 4) {
@@ -1027,7 +1014,8 @@
         // Iterating over the matrix
         for (int i = 0; i < nTerm; i++) {
-            int l = i % nXterm;         // x index
-            int m = i / nXterm;         // y index
+            int l = i / nYterm;         // x index
+            int m = i % nYterm;         // y index
             if (termMask[l][m]) {
+                matrix[i][i] = 1.0;
                 continue;
             }
@@ -1035,6 +1023,6 @@
             matrix[i][i] += sums[2*l][2*m] * wt; // The diagonal entry
             for (int j = i + 1; j < nTerm; j++) { // Doing the upper diagonal only: we will use symmetry
-                int p = j % nXterm;     // x index
-                int q = j / nXterm;     // y index
+                int p = j / nYterm;     // x index
+                int q = j % nYterm;     // y index
                 if (termMask[p][q]) {
                     continue;
@@ -1054,9 +1042,9 @@
     } else {
         // select the appropriate solution entries
-        for (int n = 0; n < nXterm; n++) {
-            for (int 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]);
-            }
+        for (int i = 0; i < nTerm; i++) {
+            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]);
         }
     }
@@ -1451,5 +1439,5 @@
     }
     psU8 ***termMask = myPoly->mask;    // Mask for polynomial terms
-    int nXYterm = nXterm * nYterm;      // Multiplication of the numbers, to calculate the index
+    int nYZterm = nYterm * nZterm;      // Multiplication of the numbers, to calculate the index
 
     // Build the B and A data structs.
@@ -1471,8 +1459,9 @@
 
         for (int i = 0; i < nTerm; i++) {
-            int iz = i / nXYterm; // z index
-            int iy = (i % nXYterm) / nXterm; // y index
-            int ix = (i % nXYterm) % nXterm; // x index
+            int ix = i / nYZterm; // x index
+            int iy = (i % nYZterm) / nZterm; // y index
+            int iz = (i % nYZterm) % nZterm; // z index
             if (termMask[ix][iy][iz]) {
+                matrix[i][i] = 1.0;
                 continue;
             }
@@ -1481,7 +1470,7 @@
             matrix[i][i] += Sums[2*ix][2*iy][2*iz] * wt;
             for (int j = i + 1; j < nTerm; j++) {
-                int jz = j / (nXterm * nYterm); // z index
-                int jy = (j % (nXterm * nYterm)) / nXterm; // y index
-                int jx = (j % (nXterm * nYterm)) % nXterm; // x index
+                int jx = j / (nYZterm); // x index
+                int jy = (j % nYZterm) / nZterm; // y index
+                int jz = (j % nYZterm) % nZterm; // z index
                 if (termMask[jx][jy][jz]) {
                     continue;
@@ -1515,12 +1504,10 @@
         }
         // select the appropriate solution entries
-        for (psS32 ix = 0; ix < nXterm; ix++) {
-            for (psS32 iy = 0; iy < nYterm; iy++) {
-                for (psS32 iz = 0; iz < nZterm; iz++) {
-                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                    myPoly->coeff[ix][iy][iz] = B->data.F64[nx];
-                    myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
-                }
-            }
+        for (int i = 0; i < nTerm; i++) {
+            int ix = i / nYZterm; // x index
+            int iy = (i % nYZterm) / nZterm; // y index
+            int iz = (i % nYZterm) % nZterm; // z index
+            myPoly->coeff[ix][iy][iz] = B->data.F64[i];
+            myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
         }
     } else {
@@ -1933,5 +1920,5 @@
     PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
-    if (mask != NULL) {
+    if (mask) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
@@ -1976,6 +1963,6 @@
     }
     psU8 ****termMask = myPoly->mask;    // Mask for polynomial terms
-    int nXYZterm = nXterm * nYterm * nZterm; // Multiplication of the numbers, for calculating the index
-    int nXYterm = nXterm * nYterm;      // Multiplication of the numbers, for calculating the index
+    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
 
     // Build the B and A data structs.
@@ -1997,9 +1984,10 @@
 
         for (int i = 0; i < nTerm; i++) {
-            int it = i / (nXYZterm); // t index
-            int iz = (i % (nXYZterm)) / (nXYterm); // z index
-            int iy = ((i % (nXYZterm)) % (nXYterm)) / nXterm; // y index
-            int ix = ((i % (nXYZterm)) % (nXYterm)) % nXterm; // x index
+            int ix = i / (nYZTterm); // x index
+            int iy = (i % (nYZTterm)) / (nZTterm); // y index
+            int iz = ((i % (nYZTterm)) % (nZTterm)) / nTterm; // z index
+            int it = ((i % (nYZTterm)) % (nZTterm)) % nTterm; // t index
             if (termMask[ix][iy][iz][it]) {
+                matrix[i][i] = 1.0;
                 continue;
             }
@@ -2008,8 +1996,8 @@
             matrix[i][i] += Sums[2*ix][2*iy][2*iz][2*it] * wt;
             for (int j = i + 1; j < nTerm; j++) {
-                int jt = j / nXYZterm; // t index
-                int jz = (j % nXYZterm) / nXYterm; // z index
-                int jy = ((j % nXYZterm) % nXYterm) / nXterm; // y index
-                int jx = ((j % nXYZterm) % nXYterm) % nXterm; // x index
+                int jx = j / nYZTterm; // x index
+                int jy = (j % nYZTterm) / nZTterm; // y index
+                int jz = ((j % nYZTterm) % nZTterm) / nTterm; // z index
+                int jt = ((j % nYZTterm) % nZTterm) % nTterm; // t index
                 if (termMask[jx][jy][jz][jt]) {
                     continue;
@@ -2047,14 +2035,11 @@
 
         // select the appropriate solution entries
-        for (psS32 ix = 0; ix < nXterm; ix++) {
-            for (psS32 iy = 0; iy < nYterm; iy++) {
-                for (psS32 iz = 0; iz < nZterm; iz++) {
-                    for (psS32 it = 0; it < nTterm; it++) {
-                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
-                        myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx];
-                        myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
-                    }
-                }
-            }
+        for (int i = 0; i < nTerm; i++) {
+            int ix = i / nYZTterm; // x index
+            int iy = (i % nYZTterm) / nZTterm; // y index
+            int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
+            int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
+            myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
+            myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
         }
     } else {
@@ -2136,5 +2121,5 @@
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
+    if (mask) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
@@ -2307,5 +2292,5 @@
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
+    if (mask) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);