Index: trunk/psLib/src/imageops/psImageStats.c
===================================================================
--- trunk/psLib/src/imageops/psImageStats.c	(revision 6959)
+++ trunk/psLib/src/imageops/psImageStats.c	(revision 7535)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.95 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-22 20:52:30 $
+ *  @version $Revision: 1.96 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-13 20:02:05 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -225,5 +225,5 @@
     chebPolys[i][j]
     sums[i][j]: This will contain the sum of
-                input->data.F32[x][y] *
+                input->data.F32[y][x] *
                 psPolynomial1DEval(
 chebPolys[i],
@@ -289,15 +289,15 @@
             //            for (x = 0; x < input->numRows; x++) {
             //                for (y = 0; y < input->numCols; y++) {
-            for (x = input->row0; x < (input->row0 + input->numRows); x++) {
-                for (y = input->col0; y < (input->col0 + input->numCols); y++) {
+            for (y = input->row0; y < (input->row0 + input->numRows); y++) {
+                for (x = input->col0; x < (input->col0 + input->numCols); x++) {
                     double pixel = 0.0;
                     if (input->type.type == PS_TYPE_S8) {
-                        pixel = (double) input->data.S8[x][y];
+                        pixel = (double) input->data.S8[y][x];
                     } else if (input->type.type == PS_TYPE_U16) {
-                        pixel = (double) input->data.U16[x][y];
+                        pixel = (double) input->data.U16[y][x];
                     } else if (input->type.type == PS_TYPE_F32) {
-                        pixel = (double) input->data.F32[x][y];
+                        pixel = (double) input->data.F32[y][x];
                     } else if (input->type.type == PS_TYPE_F64) {
-                        pixel = input->data.F64[x][y];
+                        pixel = input->data.F64[y][x];
                     }
                     sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) *
@@ -341,155 +341,4 @@
 }
 
-
-
-
-/*****************************************************************************
-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.  This is a TEST version of the code.  It is not used by anything.
-  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* psImageFitPolynomialTest(psPolynomial2D* coeffs,
-        const psImage* input)
-{
-    PS_ASSERT_IMAGE_NON_NULL(input, NULL);
-    PS_ASSERT_IMAGE_NON_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_ASSERT_POLY_NON_NULL(coeffs, NULL);
-    PS_ASSERT_POLY_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);
-
-    double min = -1.0;
-    double max = 1.0;
-    double bma = 0.5 * (max-min);  // 1
-    double bpa = 0.5 * (max+min);  // 0
-    // We must calculate the value of the image at the nodes where the
-    // Chebyshev polynomials are 0.
-    for (x = 0; x < input->numRows; x++) {
-        double xTmp = cos(M_PI * (0.5 + ((float) x)) / ((float) input->numRows));
-        double xNode = - ((xTmp + bma + bpa) - 1.0);
-        double xOrig = ((float) input->numRows) * (xNode - min) / (max - min);
-
-        for (y = 0; y < input->numCols; y++) {
-            double yTmp = cos(M_PI * (0.5 + ((float) y)) / ((float) input->numCols));
-            double yNode = - ((yTmp + bma + bpa) - 1.0);
-            double yOrig = ((float) input->numCols) * (yNode - min) / (max - min);
-
-            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yNode, xNode, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
-            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yTmp, xTmp, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
-            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yOrig, xOrig, NULL,
-                                    0, 0.0, PS_INTERPOLATE_BILINEAR);
-        }
-    }
-
-    // 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((1 + coeffs->nX) * sizeof(double *));
-    for (i = 0; i < (1 + coeffs->nX); i++) {
-        sums[i] = (double *)psAlloc((1 + 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 + 1);
-
-    // Compute the sums[][] data structure.
-    for (i = 0; i < (1 + coeffs->nX); i++) {
-        for (j = 0; j < (1 + 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 (0) {
-                        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];
-                        }
-                    }
-                    pixel = nodes->data.F64[x][y];
-                    sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) *
-                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
-                }
-            }
-        }
-    }
-
-    for (i = 0; i < (1 + coeffs->nX); i++) {
-        for (j = 0; j < (1 + 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 + 1; i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-
-    // Free some data
-    for (i = 0; i < (1 + coeffs->nX); i++) {
-        psFree(sums[i]);
-    }
-    psFree(sums);
-    psFree(cScalingFactors);
-    psFree(rScalingFactors);
-    psFree(nodes);
-
-    return (coeffs);
-}
-
 /*****************************************************************************
 XXX: Use static variables for Chebyshev polynomials and scaling factors.
@@ -500,13 +349,8 @@
     PS_ASSERT_POLY_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
 
-    psS32 x = 0;
-    psS32 y = 0;
-    psS32 i = 0;
-    psS32 j = 0;
     psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
+    long maxChebyPoly = 0;
     double *cScalingFactors = NULL;
     double *rScalingFactors = NULL;
-    float polySum = 0.0;
 
     // We scale the pixel positions to values between -1.0 and 1.0
@@ -526,11 +370,9 @@
     chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1);
 
-    //    for (x = 0; x < input->numRows; x++) {
-    //        for (y = 0; y < input->numCols; y++) {
-    for (x = input->row0; x < (input->row0 + input->numRows); x++) {
-        for (y = input->col0; y < (input->col0 + input->numCols); y++) {
-            polySum = 0.0;
-            for (i = 0; i < (1 + coeffs->nX); i++) {
-                for (j = 0; j < (1 + coeffs->nY); j++) {
+    for (long y = input->row0; y < (input->row0 + input->numRows); y++) {
+        for (long x = input->col0; x < (input->col0 + input->numCols); x++) {
+            double polySum = 0.0;
+            for (long i = 0; i < (1 + coeffs->nX); i++) {
+                for (long j = 0; j < (1 + coeffs->nY); j++) {
                     polySum +=
                         psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) *
@@ -541,11 +383,11 @@
 
             if (input->type.type == PS_TYPE_S8) {
-                input->data.S8[x][y] = (char) polySum;
+                input->data.S8[y][x] = (char) polySum;
             } else if (input->type.type == PS_TYPE_U16) {
-                input->data.U16[x][y] = (short int) polySum;
+                input->data.U16[y][x] = (short int) polySum;
             } else if (input->type.type == PS_TYPE_F32) {
-                input->data.F32[x][y] = (float) polySum;
+                input->data.F32[y][x] = (float) polySum;
             } else if (input->type.type == PS_TYPE_F64) {
-                input->data.F64[x][y] = polySum;
+                input->data.F64[y][x] = polySum;
             }
         }
@@ -555,5 +397,5 @@
     // created in this routine.
     // XXX: Use static data structures here.
-    for (i = 0; i < maxChebyPoly + 1; i++) {
+    for (long i = 0; i < maxChebyPoly + 1; i++) {
         psFree(chebPolys[i]);
     }
@@ -750,9 +592,9 @@
         return -1;
     }/* else if (col0 == col1 && row0 == row1) {
-                            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                                    "Invalid psRegion specified.  Region contains only 1 pixel.\n");
-                            return -1;
-                        }
-                    */
+                                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                                        "Invalid psRegion specified.  Region contains only 1 pixel.\n");
+                                return -1;
+                            }
+                        */
     x0 = col0;
     x1 = col1;