Index: trunk/psLib/src/imageops/psImageStats.c
===================================================================
--- trunk/psLib/src/imageops/psImageStats.c	(revision 2261)
+++ trunk/psLib/src/imageops/psImageStats.c	(revision 2262)
@@ -9,6 +9,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.47 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-11-02 01:52:43 $
+*  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-11-02 02:03:56 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -248,10 +248,18 @@
                 psPolynomial1DEval((float) y, chebPolys[j]);
         over all pixels (x,y) in the image.
- 
-XXX: Currently we only support F32.
- *****************************************************************************/
+  *****************************************************************************/
 psPolynomial2D* psImageFitPolynomial(psPolynomial2D* coeffs,
                                      const psImage* input)
 {
+    PS_IMAGE_CHECK_NULL(input, NULL);
+    PS_IMAGE_CHECK_EMPTY(input, NULL);
+    if ((input->type.type != PS_TYPE_S8) &&
+            (input->type.type != PS_TYPE_U16) &&
+            (input->type.type != PS_TYPE_F32) &&
+            (input->type.type != PS_TYPE_F64)) {
+        psError(__func__, "Unalowwable image type.\n");
+    }
+    PS_POLY_CHECK_NULL(coeffs, NULL);
+    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
     psS32 x = 0;
     psS32 y = 0;
@@ -308,8 +316,17 @@
             for (x = 0; x < input->numRows; x++) {
                 for (y = 0; y < input->numCols; y++) {
-                    sums[i][j] +=
-                        input->data.F32[x][y] *
-                        psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) *
-                        psPolynomial1DEval(cScalingFactors[y], chebPolys[j]);
+                    double pixel;
+                    if (input->type.type == PS_TYPE_S8) {
+                        pixel = (double) input->data.S8[x][y];
+                    } else if (input->type.type == PS_TYPE_U16) {
+                        pixel = (double) input->data.U16[x][y];
+                    } else if (input->type.type == PS_TYPE_F32) {
+                        pixel = (double) input->data.F32[x][y];
+                    } else if (input->type.type == PS_TYPE_F64) {
+                        pixel = input->data.F64[x][y];
+                    }
+                    sums[i][j] += pixel *
+                                  psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) *
+                                  psPolynomial1DEval(cScalingFactors[y], chebPolys[j]);
                 }
             }
@@ -320,5 +337,5 @@
         for (j = 0; j < coeffs->nY; j++) {
             coeffs->coeff[i][j] = sums[i][j];
-            coeffs->coeff[i][j] /= (float)(input->numRows * input->numCols);
+            coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols);
 
             if ((i != 0) && (j != 0)) {
@@ -352,7 +369,4 @@
 /*****************************************************************************
 XXX: Use static variables for chebyshev polynomials and scaling factors. 
- 
-XXX: Currently we only support F32.
-     S8, U16, F32, F64
  *****************************************************************************/
 psImage* psImageEvalPolynomial(psImage* input, const psPolynomial2D* coeffs)
@@ -373,5 +387,4 @@
     psS32 i = 0;
     psS32 j = 0;
-    //    float **sums = NULL;
     psPolynomial1D* *chebPolys = NULL;
     psS32 maxChebyPoly = 0;
@@ -407,5 +420,14 @@
                 }
             }
-            input->data.F32[x][y] = polySum;
+
+            if (input->type.type == PS_TYPE_S8) {
+                input->data.S8[x][y] = (char) polySum;
+            } else if (input->type.type == PS_TYPE_U16) {
+                input->data.U16[x][y] = (short int) polySum;
+            } else if (input->type.type == PS_TYPE_F32) {
+                input->data.F32[x][y] = (float) polySum;
+            } else if (input->type.type == PS_TYPE_F64) {
+                input->data.F64[x][y] = polySum;
+            }
         }
     }