Index: trunk/psLib/src/imageops/psImageStats.c
===================================================================
--- trunk/psLib/src/imageops/psImageStats.c	(revision 896)
+++ trunk/psLib/src/imageops/psImageStats.c	(revision 1022)
@@ -19,5 +19,5 @@
 /// This routine must determine the various statistics for the image.
 /*****************************************************************************
-    GUS: verify that image/mask have the correct types, and sizes.
+    NOTE: verify that image/mask have the correct types, and sizes.
  *****************************************************************************/
 psStats *psImageStats(psStats *stats,
@@ -32,10 +32,10 @@
     int j = 0;
 
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if ((stats == NULL) ||
             (in == NULL)) {
         return(NULL);
     }
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if (stats->options == 0) {
         return(stats);
@@ -52,11 +52,11 @@
 
     if (mask != NULL) {
-        // GUS: verify that mask data is PS_TYPE_U8.
-        // GUS: figure out mask types
+        // NOTE: verify that mask data is PS_TYPE_U8.
+        // NOTE: figure out mask types
 
         junkMask = psVectorAlloc(mask->numRows * mask->numCols,
                                  mask->type.type);
         junkMask->n = junkMask->nalloc;
-        // GUS: Is there a more efficient way to do this?
+        // NOTE: Is there a more efficient way to do this?
         ptr=0;
         for (i=0;i<mask->numRows;i++) {
@@ -78,5 +78,5 @@
     NOTE: We assume that the psHistogram structure out has already been
     allocated and initialized.
-    GUS: verify that image/mask have the, correct types and  sizes.
+    NOTE: verify that image/mask have the, correct types and  sizes.
  *****************************************************************************/
 psHistogram *psImageHistogram(psHistogram *out,
@@ -91,5 +91,5 @@
     int j = 0;
 
-    // GUS: Verify this action.
+    // NOTE: Verify this action.
     if ((out == NULL) ||
             (in == NULL)) {
@@ -100,10 +100,10 @@
     junkData->n = junkData->nalloc;
 
-    // GUS: Is there a more efficient way to do this?  memcopy() won't work.
+    // NOTE: Is there a more efficient way to do this?  memcopy() won't work.
     // Can we trick the junkData structure to use the image buffer, then
     // untrick it before we deallocate it (so we won't deallocate that buffer
     // twice?
 
-    // GUS: Make sure you have the numRows/NumCols in the right place.
+    // NOTE: Make sure you have the numRows/NumCols in the right place.
     ptr=0;
     for (i=0;i<in->numRows;i++) {
@@ -114,6 +114,6 @@
 
     if (mask != NULL) {
-        // GUS: verify that mask data is PS_TYPE_U8.
-        // GUS: figure out mask types
+        // NOTE: verify that mask data is PS_TYPE_U8.
+        // NOTE: figure out mask types
 
         junkMask = psVectorAlloc(mask->numRows * mask->numCols,
@@ -121,5 +121,5 @@
         junkMask->n = junkMask->nalloc;
 
-        // GUS: Is there a more efficient way to do this?
+        // NOTE: Is there a more efficient way to do this?
         ptr=0;
         for (i=0;i<mask->numRows;i++) {
@@ -137,4 +137,65 @@
     return(out);
 }
+
+float *p_psCalcScaleFactorsFit(int n)
+{
+    int i = 0;
+    float tmp = 0.0;
+    float *scalingFactors = (float *) psAlloc(n * sizeof(float));
+
+    for (i=0;i<n;i++) {
+        //     ((2.0 * (float) i) / ((float) (n-1))) - 1.0;
+        //        tmp = (float) (i + 1);
+        tmp = (float) (n - i);
+        tmp = (M_PI * (tmp - 0.5)) / ((float) n);
+        scalingFactors[i] = cos(tmp);
+    }
+
+    return(scalingFactors);
+}
+
+float *p_psCalcScaleFactorsEval(int n)
+{
+    return p_psCalcScaleFactorsFit(n);
+
+    int i = 0;
+    //    float tmp = 0.0;
+    float *scalingFactors = (float *) psAlloc(n * sizeof(float));
+    for (i=0;i<n;i++) {
+        //          scalingFactors[i] = ((2.0 * (float) i) / ((float) (n-1))) - 1.0;
+        tmp = (float) (i + 1);
+        tmp = (M_PI * (tmp - 0.5)) / ((float) n);
+        scalingFactors[i] = cos(tmp);
+    }
+    return(scalingFactors);
+}
+
+psPolynomial1D **p_psCreateChebyshevPolys(int maxChebyPoly)
+{
+    psPolynomial1D **chebPolys = NULL;
+    int i = 0;
+    int j = 0;
+
+    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
+                                            sizeof(psPolynomial1D *));
+    for (i=0;i<maxChebyPoly;i++) {
+        chebPolys[i] = psPolynomial1DAlloc(i+1);
+    }
+
+    // Create the Chebyshev polynomials.  Polynomial i has i-th order.
+    chebPolys[0]->coeff[0] = 1;
+    chebPolys[1]->coeff[1] = 1;
+    for (i=2;i<maxChebyPoly;i++) {
+        for (j=0;j<chebPolys[i-1]->n;j++) {
+            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
+        }
+        for (j=0;j<chebPolys[i-2]->n;j++) {
+            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
+        }
+    }
+
+    return(chebPolys);
+}
+
 
 /*****************************************************************************
@@ -146,6 +207,14 @@
     Output:
     Internal Data Structures:
+ 
  chebPolys[i][j] 
-        sums[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.
+ 
+ 
  *****************************************************************************/
 psPolynomial2D *
@@ -153,4 +222,6 @@
                      psPolynomial2D *coeffs)
 {
+    int k = 0;
+    int l = 0;
     int x = 0;
     int y = 0;
@@ -160,18 +231,18 @@
     psPolynomial1D **chebPolys = NULL;
     int maxChebyPoly = 0;
-
+    float *cScalingFactors = NULL;
+    float *rScalingFactors = NULL;
+    float tmp = 0.0;
 
     // Create the sums[][] data structure.  This will hold the LHS of equation
     // 29 in the ADD: sums[k][l] = SUM { image(x,y) * Tk(x) * Tl(y) }
-
     sums = (float **) psAlloc(coeffs->nX * sizeof(float *));
     for (i=0;i<coeffs->nX;i++) {
         sums[i] = (float *) psAlloc(coeffs->nY * sizeof(float));
     }
-    for (i=0;i<coeffs->nX;i++) {
-        for (j=0;j<coeffs->nY;j++) {
-            sums[i][j] = 0.0;
-        }
-    }
+
+    // We scale the pixel positions to values between -1.0 and 1.0
+    cScalingFactors = p_psCalcScaleFactorsFit(input->numRows);
+    rScalingFactors = p_psCalcScaleFactorsFit(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -180,30 +251,26 @@
         maxChebyPoly = coeffs->nY;
     }
-
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
-                                            sizeof(psPolynomial1D *));
-    for (i=0;i<maxChebyPoly;i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i);
-    }
-
-    // Create the Chebyshev polynomials
-    chebPolys[1]->coeff[0] = 1;
-    for (i=2;i<maxChebyPoly;i++) {
-        for (j=0;j<chebPolys[i-1]->n;j++) {
-            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
-        }
-        for (j=0;j<chebPolys[i-2]->n;j++) {
-            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
-        }
-    }
-
-    for (x=0;x<input->numRows;x++) {
-        for (y=0;y<input->numCols;y++) {
-            for (i=0;i<coeffs->nX;i++) {
-                for (j=0;j<coeffs->nY;j++) {
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
+
+    // Sanity check for the Chebyshevs
+    for (i=0;i<coeffs->nX;i++) {
+        for (j=0;j<coeffs->nY;j++) {
+            tmp = 0.0;
+            for (x=0;x<input->numRows;x++) {
+                tmp+= psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                      psEvalPolynomial1D(rScalingFactors[x], chebPolys[j]);
+
+            }
+        }
+    }
+
+    for (i=0;i<coeffs->nX;i++) {
+        for (j=0;j<coeffs->nY;j++) {
+            sums[i][j] = 0.0;
+            for (x=0;x<input->numRows;x++) {
+                for (y=0;y<input->numCols;y++) {
                     sums[i][j]+= input->data.F32[x][y] *
-                                 psEvalPolynomial1D((float) x, chebPolys[i]) *
-                                 psEvalPolynomial1D((float) y, chebPolys[j]);
-
+                                 psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                                 psEvalPolynomial1D(cScalingFactors[y], chebPolys[j]);
                 }
             }
@@ -211,25 +278,47 @@
     }
 
-    // GUS: Check the math on this at a later date.
-    coeffs->coeff[0][0] = sums[0][0] / (coeffs->nX * coeffs->nY);
-    for (i=0;i<coeffs->nX;i++) {
-        coeffs->coeff[i][0] = 0.5 *
-                              ((sums[i][0] * 4.0) / (coeffs->nX * coeffs->nY)) -
-                              (2.0 * coeffs->coeff[0][0]);
+    // 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] = 0.5 *
-                              ((sums[0][j] * 4.0) / (coeffs->nX * coeffs->nY)) -
-                              (2.0 * coeffs->coeff[0][0]);
+        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) / (coeffs->nX * coeffs->nY) -
-                coeffs->coeff[0][0] -
-                coeffs->coeff[i][0] -
-                coeffs->coeff[0][j];
-        }
-    }
+                ((sums[i][0] * 4.0) / ((float) (coeffs->nX * coeffs->nY))) -
+                (coeffs->coeff[0][0] +
+                 coeffs->coeff[i][0] +
+                 coeffs->coeff[0][j]);
+        }
+    }
+
+    for (k=0;k<coeffs->nX;k++) {
+        for (l=0;l<coeffs->nY;l++) {
+            tmp = 0.0;
+            for (i=0;i<coeffs->nX;i++) {
+                for (j=0;j<coeffs->nY;j++) {
+                    for (x=0;x<input->numRows;x++) {
+                        for (y=0;y<input->numCols;y++) {
+                            tmp+= (coeffs->coeff[i][j] *
+                                   psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                                   psEvalPolynomial1D(rScalingFactors[y], chebPolys[j]) *
+                                   psEvalPolynomial1D(rScalingFactors[x], chebPolys[k]) *
+                                   psEvalPolynomial1D(rScalingFactors[y], chebPolys[l]));
+
+                        }
+                    }
+                }
+            }
+        }
+    }
+
 
     // Free the Chebyshev polynomials that were created in this routine.
@@ -244,4 +333,6 @@
     }
     psFree(sums);
+    psFree(cScalingFactors);
+    psFree(rScalingFactors);
 
     return(coeffs);
@@ -255,11 +346,13 @@
                       const psPolynomial2D *coeffs)
 {
-    int i = 0;
-    int j = 0;
     int x = 0;
     int y = 0;
+    int i = 0;
+    int j = 0;
     float **sums = NULL;
     psPolynomial1D **chebPolys = NULL;
     int maxChebyPoly = 0;
+    float *cScalingFactors = NULL;
+    float *rScalingFactors = NULL;
     float polySum = 0.0;
 
@@ -275,4 +368,8 @@
         }
     }
+
+    // We scale the pixel positions to values between -1.0 and 1.0
+    cScalingFactors = p_psCalcScaleFactorsEval(input->numRows);
+    rScalingFactors = p_psCalcScaleFactorsEval(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -282,20 +379,5 @@
     }
 
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly *
-                                            sizeof(psPolynomial1D *));
-    for (i=0;i<maxChebyPoly;i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i);
-    }
-
-    // Create the Chebyshev polynomials
-    chebPolys[1]->coeff[0] = 1;
-    for (i=2;i<maxChebyPoly;i++) {
-        for (j=0;j<chebPolys[i-1]->n;j++) {
-            chebPolys[i]->coeff[j+1] = 2 * chebPolys[i-1]->coeff[j];
-        }
-        for (j=0;j<chebPolys[i-2]->n;j++) {
-            chebPolys[i]->coeff[j]-= chebPolys[i-2]->coeff[j];
-        }
-    }
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
 
     for (x=0;x<input->numRows;x++) {
@@ -304,7 +386,8 @@
             for (i=0;i<coeffs->nX;i++) {
                 for (j=0;j<coeffs->nY;j++) {
-                    polySum+= psEvalPolynomial1D((float)x, chebPolys[i]) *
-                              psEvalPolynomial1D((float)y, chebPolys[j]) *
+                    polySum+= psEvalPolynomial1D(rScalingFactors[x], chebPolys[i]) *
+                              psEvalPolynomial1D(cScalingFactors[y], chebPolys[j]) *
                               coeffs->coeff[i][j];
+
                 }
             }
@@ -324,4 +407,6 @@
     }
     psFree(sums);
+    psFree(cScalingFactors);
+    psFree(rScalingFactors);
 
     return(0);