Index: trunk/psLib/src/image/psImageStats.c
===================================================================
--- trunk/psLib/src/image/psImageStats.c	(revision 2634)
+++ trunk/psLib/src/image/psImageStats.c	(revision 2735)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-04 19:22:25 $
+ *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-16 21:58:25 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -174,127 +174,4 @@
     return (chebPolys);
 }
-
-/*****************************************************************************
-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] *
-                psPolynomial1DEval(
-chebPolys[i],
-(float) x) *
-                psPolynomial1DEval(
-chebPolys[j],
-(float) y, 
-);
-        over all pixels (x,y) in the image.
-  *****************************************************************************/
-psPolynomial2D* psImageFitPolynomialORIG(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(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n");
-    }
-    PS_POLY_CHECK_NULL(coeffs, NULL);
-    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
-    psS32 x = 0;
-    psS32 y = 0;
-    psS32 i = 0;
-    psS32 j = 0;
-    double **sums = NULL;
-    psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
-    double *cScalingFactors = NULL;
-    double *rScalingFactors = NULL;
-
-    // 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 = (double **)psAlloc(coeffs->nX * sizeof(double *));
-    for (i = 0; i < coeffs->nX; i++) {
-        sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double));
-    }
-    // We scale the pixel positions to values
-    // between -1.0 and 1.0
-    rScalingFactors = calcScaleFactors(input->numRows);
-    cScalingFactors = calcScaleFactors(input->numCols);
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = coeffs->nX;
-    if (coeffs->nY > coeffs->nX) {
-        maxChebyPoly = coeffs->nY;
-    }
-    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
-
-    // Compute the sums[][] data structure.
-    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++) {
-                    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(chebPolys[i],rScalingFactors[x]) *
-                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
-                }
-            }
-        }
-    }
-
-    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] /= (double)(input->numRows * input->numCols);
-
-            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;
-            }
-        }
-    }
-
-    // Free the Chebyshev polynomials that were
-    // created in this routine.
-    for (i = 0; i < maxChebyPoly; i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-
-    // Free some data
-    for (i = 0; i < coeffs->nX; i++) {
-        psFree(sums[i]);
-    }
-    psFree(sums);
-    psFree(cScalingFactors);
-    psFree(rScalingFactors);
-
-    return (coeffs);
-}
-
-
-
 
 /*****************************************************************************
@@ -339,4 +216,127 @@
     double *cScalingFactors = NULL;
     double *rScalingFactors = NULL;
+
+    // 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 = (double **)psAlloc(coeffs->nX * sizeof(double *));
+    for (i = 0; i < coeffs->nX; i++) {
+        sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double));
+    }
+    // We scale the pixel positions to values
+    // between -1.0 and 1.0
+    rScalingFactors = calcScaleFactors(input->numRows);
+    cScalingFactors = calcScaleFactors(input->numCols);
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = coeffs->nX;
+    if (coeffs->nY > coeffs->nX) {
+        maxChebyPoly = coeffs->nY;
+    }
+    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
+
+    // Compute the sums[][] data structure.
+    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++) {
+                    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(chebPolys[i],rScalingFactors[x]) *
+                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
+                }
+            }
+        }
+    }
+
+    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] /= (double)(input->numRows * input->numCols);
+
+            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;
+            }
+        }
+    }
+
+    // Free the Chebyshev polynomials that were
+    // created in this routine.
+    for (i = 0; i < maxChebyPoly; i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+
+    // Free some data
+    for (i = 0; i < coeffs->nX; i++) {
+        psFree(sums[i]);
+    }
+    psFree(sums);
+    psFree(cScalingFactors);
+    psFree(rScalingFactors);
+
+    return (coeffs);
+}
+
+
+
+
+/*****************************************************************************
+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] *
+                psPolynomial1DEval(
+chebPolys[i],
+(float) x) *
+                psPolynomial1DEval(
+chebPolys[j],
+(float) y, 
+);
+        over all pixels (x,y) in the image.
+  *****************************************************************************/
+psPolynomial2D* psImageFitPolynomialHMMM(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(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n");
+    }
+    PS_POLY_CHECK_NULL(coeffs, NULL);
+    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
+    psS32 x = 0;
+    psS32 y = 0;
+    psS32 i = 0;
+    psS32 j = 0;
+    double **sums = NULL;
+    psPolynomial1D* *chebPolys = NULL;
+    psS32 maxChebyPoly = 0;
+    double *cScalingFactors = NULL;
+    double *rScalingFactors = NULL;
     psImage *nodes = psImageAlloc(input->numCols, input->numRows, PS_TYPE_F64);