Index: trunk/psLib/src/image/psImageStats.c
===================================================================
--- trunk/psLib/src/image/psImageStats.c	(revision 775)
+++ trunk/psLib/src/image/psImageStats.c	(revision 776)
@@ -40,5 +40,4 @@
     return(hist);
 }
-
 
 /*****************************************************************************
@@ -140,4 +139,48 @@
                       const 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];
+        }
+    }
+
+
     return(0);
 }