Index: trunk/psLib/src/imageops/psImageStats.c
===================================================================
--- trunk/psLib/src/imageops/psImageStats.c	(revision 1025)
+++ trunk/psLib/src/imageops/psImageStats.c	(revision 1072)
@@ -155,4 +155,11 @@
 }
 
+/*****************************************************************************
+p_psCalcScaleFactorsEval(n): The Chebyshev polynomials are defined over the
+interval [-1.0 : 1.0].  Images typically have sizes of 512x512 or more.  In
+order to use Chebyshev polynomials, we must scales the coordinates from
+0:512 to -1:1.  This routine takes as input an integer N and produces as
+output a vector of evenly spaced floating point values between -1.0:1.0.
+ *****************************************************************************/
 float *p_psCalcScaleFactorsEval(int n)
 {
@@ -200,21 +207,16 @@
 
 /*****************************************************************************
-    psImageFitPolynomial():
- This routine takes as input a 2-D image and produces as output
- the coefficients of the Chebyshev polynomials which match that
- input image.
-    Input:
-    Output:
-    Internal Data Structures:
- 
- chebPolys[i][j] 
- 
- sums[i][j]: This will contain the sum of 
-  input->data.F32[x][y] *
+psImageFitPolynomial(): This routine takes as input a 2-D image and produces
+as output the coefficients of the Chebyshev polynomials which match that
+input image.
+  Input:
+  Output:
+  Internal Data Structures:
+    chebPolys[i][j] 
+    sums[i][j]: This will contain the sum of 
+                input->data.F32[x][y] *
                 psEvalPolynomial1D((float) x, chebPolys[i]) *
                 psEvalPolynomial1D((float) y, chebPolys[j]);
-      over all pixels (x,y) in the image.
- 
- 
+        over all pixels (x,y) in the image.
  *****************************************************************************/
 psPolynomial2D *
@@ -253,5 +255,5 @@
     chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
 
-    // Sanity check for the Chebyshevs
+    // Sanity check for the Chebyshevs.
     for (i=0;i<coeffs->nX;i++) {
         for (j=0;j<coeffs->nY;j++) {
@@ -262,7 +264,9 @@
 
             }
-        }
-    }
-
+            printf("SUM(Cheby(%d) * Cheby(%d)) is %f\n", i, j, tmp);
+        }
+    }
+
+    // Compute the sums[][] data structure.
     for (i=0;i<coeffs->nX;i++) {
         for (j=0;j<coeffs->nY;j++) {
@@ -279,23 +283,34 @@
 
     // NOTE: Check the math on this at a later date.
-    coeffs->coeff[0][0] = sums[0][0] / ((float) (coeffs->nX * coeffs->nY));
-    for (i=0;i<coeffs->nX;i++) {
-        coeffs->coeff[i][0] = ((sums[i][0] * 2.0) /
-                               ((float) (coeffs->nX * coeffs->nY))) -
-                              coeffs->coeff[0][0];
-    }
-    for (j=0;j<coeffs->nY;j++) {
-        coeffs->coeff[0][j] =
-            ((sums[0][j] * 2.0) /
-             ((float) (coeffs->nX * coeffs->nY))) -
-            coeffs->coeff[0][0];
-    }
-    for (i=1;i<coeffs->nX;i++) {
-        for (j=1;j<coeffs->nY;j++) {
-            coeffs->coeff[i][j] =
-                ((sums[i][0] * 4.0) / ((float) (coeffs->nX * coeffs->nY))) -
-                (coeffs->coeff[0][0] +
-                 coeffs->coeff[i][0] +
-                 coeffs->coeff[0][j]);
+    //    coeffs->coeff[0][0] = sums[0][0] / ((float) (coeffs->nX * coeffs->nY));
+    //    for (i=0;i<coeffs->nX;i++) {
+    //        coeffs->coeff[i][0] = ((sums[i][0] * 2.0) /
+    //                               ((float) (coeffs->nX * coeffs->nY))) -
+    //                              coeffs->coeff[0][0];
+    //    }
+    //    for (j=0;j<coeffs->nY;j++) {
+    //        coeffs->coeff[0][j] =
+    //            ((sums[0][j] * 2.0) /
+    //             ((float) (coeffs->nX * coeffs->nY))) -
+    //            coeffs->coeff[0][0];
+    //    }
+
+    for (i=0;i<coeffs->nX;i++) {
+        for (j=0;j<coeffs->nY;j++) {
+            coeffs->coeff[i][j] = sums[i][j];
+            coeffs->coeff[i][j]/= ((float) coeffs->nX) * ((float) coeffs->nY);
+            if ((i != 0) && (j != 0)) {
+                coeffs->coeff[i][j]*= 4.0;
+            } else if ((i == 0) && (j == 0)) {
+                coeffs->coeff[i][j]*= 1.0;
+            } else {
+                coeffs->coeff[i][j]*= 2.0;
+            }
+
+            //            coeffs->coeff[i][j] =
+            //                ((sums[i][0] * 4.0) / ((float) (coeffs->nX * coeffs->nY))) -
+            //                (coeffs->coeff[0][0] +
+            //                 coeffs->coeff[i][0] +
+            //                 coeffs->coeff[0][j]);
         }
     }