Index: trunk/psLib/src/image/psImageStats.c
===================================================================
--- trunk/psLib/src/image/psImageStats.c	(revision 1073)
+++ trunk/psLib/src/image/psImageStats.c	(revision 1076)
@@ -164,12 +164,13 @@
 float *p_psCalcScaleFactorsEval(int n)
 {
+    int i = 0;
+    float tmp = 0.0;
     return p_psCalcScaleFactorsFit(n);
-    float tmp = 0.0;
-    int i = 0;
-    //    float tmp = 0.0;
+
+    printf("Should not get here\n");
     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 = (float) (n - i);
         tmp = (M_PI * (tmp - 0.5)) / ((float) n);
         scalingFactors[i] = cos(tmp);
@@ -224,6 +225,4 @@
                      psPolynomial2D *coeffs)
 {
-    int k = 0;
-    int l = 0;
     int x = 0;
     int y = 0;
@@ -245,6 +244,6 @@
 
     // We scale the pixel positions to values between -1.0 and 1.0
-    cScalingFactors = p_psCalcScaleFactorsFit(input->numRows);
-    rScalingFactors = p_psCalcScaleFactorsFit(input->numCols);
+    rScalingFactors = p_psCalcScaleFactorsFit(input->numRows);
+    cScalingFactors = p_psCalcScaleFactorsFit(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -264,5 +263,5 @@
 
             }
-            printf("SUM(Cheby(%d) * Cheby(%d)) is %f\n", i, j, tmp);
+            //printf("SUM(Cheby(%d) * Cheby(%d)) is %f\n", i, j, tmp);
         }
     }
@@ -282,22 +281,9 @@
     }
 
-    // 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=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);
+            coeffs->coeff[i][j]/= (float) (input->numRows * input->numCols);
+
             if ((i != 0) && (j != 0)) {
                 coeffs->coeff[i][j]*= 4.0;
@@ -307,33 +293,6 @@
                 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]);
-        }
-    }
-
-    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.
@@ -385,6 +344,6 @@
 
     // We scale the pixel positions to values between -1.0 and 1.0
-    cScalingFactors = p_psCalcScaleFactorsEval(input->numRows);
-    rScalingFactors = p_psCalcScaleFactorsEval(input->numCols);
+    rScalingFactors = p_psCalcScaleFactorsEval(input->numRows);
+    cScalingFactors = p_psCalcScaleFactorsEval(input->numCols);
 
     // Determine how many Chebyshev polynomials are needed, then create them.
@@ -427,2 +386,67 @@
     return(0);
 }
+
+/*****************************************************************************
+p_psImagePixelInterpolation(image, x, y): this routine takes as input an
+image and coordinates (x, y) and produces as output the corresponding pixel
+value at the those coordinates.  For fractional corrdinates (x, y), 2-D
+linear interpolation is performed on the image.
+ *****************************************************************************/
+float psImagePixelInterpolation(psImage *input,
+                                float x,
+                                float y)
+{
+    float floorX = 0.0;
+    float floorY = 0.0;
+    float fracX = 0.0;
+    float fracY = 0.0;
+    int intFloorX = 0;
+    int intFloorY = 0;
+    float x1 = 0.0;
+    float x2 = 0.0;
+    float pixel = 0.0;
+
+    if ((x < 0.0) ||
+            (x > ((float) input->numRows-1)) ||
+            (y < 0.0) ||
+            (y > ((float) input->numCols-1))) {
+        psAbort(__func__,
+                "Fractional coordinates (%f %f) outside image range.", x, y);
+    }
+
+    floorX = floorf(x);
+    intFloorX = (int) floorX;
+    fracX = x - floorX;
+
+    floorY = floorf(y);
+    intFloorY = (int) floorY;
+    fracY = y - floorY;
+
+    if (intFloorX == (input->numRows-1)) {
+        pixel = input->data.F32[intFloorX][intFloorY];
+        if (intFloorY < (input->numCols-1)) {
+            pixel+= fracY * (input->data.F32[intFloorX][intFloorY+1] -
+                             input->data.F32[intFloorX][intFloorY]);
+        }
+        return(pixel);
+    } else if (intFloorX == (input->numCols-1)) {
+        if (intFloorX < (input->numRows-1)) {
+            pixel+= fracX * (input->data.F32[intFloorX+1][intFloorY] -
+                             input->data.F32[intFloorX][intFloorY]);
+        }
+        return(pixel);
+    }
+
+
+    x1 = input->data.F32[intFloorX][intFloorY];
+    x1+= fracY * (input->data.F32[intFloorX][intFloorY+1] -
+                  input->data.F32[intFloorX][intFloorY]);
+
+    x2 = input->data.F32[intFloorX+1][intFloorY];
+    x2+= fracY * (input->data.F32[intFloorX+1][intFloorY+1] -
+                  input->data.F32[intFloorX+1][intFloorY]);
+    pixel = x1;
+    pixel+= fracX * (x2 - x1);
+
+    return(pixel);
+}