Index: trunk/psLib/src/image/psImageStats.c
===================================================================
--- trunk/psLib/src/image/psImageStats.c	(revision 767)
+++ trunk/psLib/src/image/psImageStats.c	(revision 773)
@@ -1,2 +1,21 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdarg.h>
+#include <float.h>
+#include <math.h>
+
+#include "psMemory.h"
+#include "psVector.h"
+#include "psTrace.h"
+#include "psError.h"
+#include "psAbort.h"
+#include "psStats.h"
+#include "psSort.h"
+#include "psImage.h"
+#include "psFunctions.h"
+
+#include "float.h"
+#include <math.h>
 /// This routine must determine the various statistics for the image.
 psStats *
@@ -4,28 +23,116 @@
                 psStats *stats)
 {
-    if (input->type.typr == PS_TYPE_UINT8) {
-        return(psArrayStats(input->data.u8, NULL, 0xffffffff, stats));
-    } else if (input->type.typr == PS_TYPE_UINT16) {
-        return(psArrayStats(input->data.u16, NULL, 0xffffffff, stats));
-    } else if (input->type.typr == PS_TYPE_FLOAT) {
-        return(psArrayStats(input->data.F32, NULL, 0xffffffff, stats));
-    } else if (input->type.typr == PS_TYPE_DOUBLE) {
-        return(psArrayStats(input->data.F64, NULL, 0xffffffff, stats));
-    }
+    psVector *junk=NULL;
+
+    junk = psVectorAlloc(input->numRows * input->numCols, input->type.type);
+    return(psArrayStats(junk, NULL, 0xffffffff, stats));
 }
-
 
 psHistogram *
 psImageHistogram(psHistogram *hist,
                  const psImage *input)
-{}
+{
+    return(hist);
+}
 
+
+/*****************************************************************************
+    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] 
+ 
+ *****************************************************************************/
 psPolynomial2D *
 psImageFitPolynomial(const psImage *input,
-                     psPolynomial2D *coeffs);
-{}
+                     psPolynomial2D *coeffs)
+{
+    int x = 0;
+    int y = 0;
+    int i = 0;
+    int j = 0;
+    float **sums = NULL;
+    psPolynomial1D **chebPolys = NULL;
+    int maxChebyPoly = 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;
+        }
+    }
+
+    // Determine how many Chebyshev polynomials are needed, then create them.
+    maxChebyPoly = coeffs->nX;
+    if (coeffs->nX > coeffs->nY) {
+        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++) {
+                    sums[i][j]+= input->data.F32[x][y] *
+                                 psEvalPolynomial1D((float) x, chebPolys[i]) *
+                                 psEvalPolynomial1D((float) y, chebPolys[j]);
+                }
+            }
+        }
+    }
+
+
+    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]);
+    }
+    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]);
+    }
+    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];
+        }
+    }
+    return(coeffs);
+}
 
 int
 psImageEvalPolynomial(const psImage *input,
-                      const psPolynomial2D *coeffs);
-{}
+                      const psPolynomial2D *coeffs)
+{
+    return(0);
+}