Index: /trunk/psLib/test/image/tst_psImageStats.c
===================================================================
--- /trunk/psLib/test/image/tst_psImageStats.c	(revision 3310)
+++ /trunk/psLib/test/image/tst_psImageStats.c	(revision 3310)
@@ -0,0 +1,692 @@
+/*
+*  C Implementation: %{MODULE}
+*
+* Description: 
+*
+*
+* Author: %{AUTHOR} <%{EMAIL}>, (C) %{YEAR}
+*
+* Copyright: See COPYING file that comes with this distribution
+*
+*/
+
+#include <stdio.h>
+#include <float.h>
+#include <math.h>
+
+#include "pslib.h"
+#include "psTest.h"
+
+static psS32 testPsImageHistogram(void);
+static psS32 testPsImageStats(void);
+static psS32 testPsImageFitPolynomial(void);
+static psS32 testPsImagePixelInterpolate(void);
+static psS32 testPsImageEvalPolynom(void);
+
+static bool FitChebyF32(int numCols, int numRows);
+static bool FitChebyF64(int numCols, int numRows);
+
+testDescription tests[] = {
+                              {testPsImageHistogram,0,"psImageHistogram",0,false},
+                              {testPsImageStats,1,"psImageStats",0,false},
+                              {testPsImageFitPolynomial,2,"psImageFitPolynomial",0,false},
+                              {testPsImagePixelInterpolate,3,"psImagePixelInterpolate",0,false},
+                              {testPsImageEvalPolynom,4,"psImageEvalPolynom()",0,false},
+                              {NULL}
+                          };
+
+psS32 main(psS32 argc, char* argv[])
+{
+    psLogSetLevel(PS_LOG_INFO);
+
+    return ! runTestSuite(stderr,"psImage",tests,argc,argv);
+}
+
+
+
+static psS32 testPsImageHistogram(void)
+{
+    const int NUM_BINS = 20;
+    const int N = 32;
+    const int M = 64;
+
+    psHistogram * myHist = NULL;
+    psHistogram *myHist2 = NULL;
+    psImage *tmpImage = NULL;
+    psImage *tmpImage2 = NULL;
+    psImage *tmpMask = NULL;
+    psImage *tmpMask2 = NULL;
+    psS32 testStatus = true;
+    psS32 nb = 0;
+    psS32 i = 0;
+    psS32 j = 0;
+    psS32 IMAGE_X_SIZE = 0;
+    psS32 IMAGE_Y_SIZE = 0;
+    psS32 currentId = 0;
+
+    currentId = psMemGetId();
+    for ( nb = 0;nb < 7;nb++ ) {
+        if ( nb == 0 ) {
+            IMAGE_X_SIZE = 1;
+            IMAGE_Y_SIZE = 1;
+        }
+        if ( nb == 1 ) {
+            IMAGE_X_SIZE = 1;
+            IMAGE_Y_SIZE = N;
+        }
+        if ( nb == 2 ) {
+            IMAGE_X_SIZE = N;
+            IMAGE_Y_SIZE = 1;
+        }
+        if ( nb == 3 ) {
+            IMAGE_X_SIZE = N;
+            IMAGE_Y_SIZE = N;
+        }
+        if ( nb == 4 ) {
+            IMAGE_X_SIZE = N;
+            IMAGE_Y_SIZE = M;
+        }
+        if ( nb == 5 ) {
+            IMAGE_X_SIZE = M;
+            IMAGE_Y_SIZE = N;
+        }
+        if ( nb == 6 ) {
+            IMAGE_X_SIZE = M;
+            IMAGE_Y_SIZE = N;
+        }
+        fprintf(stderr, "*******************************\n" );
+        if (nb != 6) {
+            fprintf(stderr, "* IMAGE SIZE is (%d by %d)\n", IMAGE_X_SIZE, IMAGE_Y_SIZE );
+        } else {
+            fprintf(stderr, "* IMAGE SUBSET SIZE is (%d by %d)\n", IMAGE_X_SIZE, IMAGE_Y_SIZE );
+        }
+        fprintf(stderr, "*******************************\n" );
+        /*********************************************************************/
+        /*  Allocate and initialize data structures                      */
+        /*********************************************************************/
+        if (nb != 6) {
+            tmpImage = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_F32 );
+            tmpMask = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_U8 );
+            tmpImage2 = NULL;
+            tmpMask2 = NULL;
+        } else {
+            tmpImage2 = psImageAlloc( IMAGE_X_SIZE*2, IMAGE_Y_SIZE*2, PS_TYPE_F32 );
+            tmpMask2 = psImageAlloc( IMAGE_X_SIZE*2, IMAGE_Y_SIZE*2, PS_TYPE_U8 );
+            tmpImage = psImageSubset(tmpImage2,
+                                     IMAGE_X_SIZE/4, IMAGE_Y_SIZE/4,
+                                     IMAGE_X_SIZE/4+IMAGE_X_SIZE, IMAGE_Y_SIZE/4+IMAGE_Y_SIZE);
+            tmpMask= psImageSubset(tmpMask2,
+                                   IMAGE_X_SIZE/4, IMAGE_Y_SIZE/4,
+                                   IMAGE_X_SIZE/4+IMAGE_X_SIZE, IMAGE_Y_SIZE/4+IMAGE_Y_SIZE);
+        }
+
+
+        for ( i = 0;i < tmpImage->numRows;i++ ) {
+            for ( j = 0;j < tmpImage->numCols;j++ ) {
+                tmpImage->data.F32[ i ][ j ] = ( float ) ( i + j + 0.1 );
+            }
+        }
+
+
+        for ( i = 0;i < tmpMask->numRows;i++ ) {
+            for ( j = 0;j < tmpMask->numCols;j++ ) {
+                if ( ( i > ( tmpMask->numRows / 2 ) ) &&
+                        ( j > ( tmpMask->numCols / 2 ) ) ) {
+                    tmpMask->data.U8[ i ][ j ] = 1;
+                } else {
+                    tmpMask->data.U8[ i ][ j ] = 0;
+                }
+            }
+        }
+
+        /*************************************************************************/
+        /*  Calculate Histogram with no mask                             */
+        /*************************************************************************/
+
+        myHist = psHistogramAlloc( 0.0, ( float ) ( IMAGE_X_SIZE + IMAGE_Y_SIZE ),
+                                   NUM_BINS );
+        myHist = psImageHistogram( myHist, tmpImage, NULL, 0 );
+        for ( i = 0;i < NUM_BINS;i++ ) {
+            fprintf(stderr, "Bin number %d bounds: (%.1f - %.1f) data (%.2f)\n", i,
+                    myHist->bounds->data.F32[ i ],
+                    myHist->bounds->data.F32[ i + 1 ],
+                    myHist->nums->data.F32[ i ] );
+        }
+        psFree( myHist );
+
+        /*************************************************************************/
+        /*  Calculate Histogram with mask                                */
+        /*************************************************************************/
+
+        myHist = psHistogramAlloc( 0.0, ( float ) ( IMAGE_X_SIZE + IMAGE_Y_SIZE ),
+                                   NUM_BINS );
+        myHist = psImageHistogram( myHist, tmpImage, tmpMask, 1 );
+        for ( i = 0;i < NUM_BINS;i++ ) {
+            fprintf( stderr, "Bin number %d bounds: (%.2f - %.2f) data (%.2f)\n", i,
+                     myHist->bounds->data.F32[ i ],
+                     myHist->bounds->data.F32[ i + 1 ],
+                     myHist->nums->data.F32[ i ] );
+        }
+
+        /*************************************************************************/
+        /*  Deallocate data structures                                   */
+        /*************************************************************************/
+        psFree( myHist );
+        psFree( tmpImage );
+        psFree( tmpImage2 );
+        psFree( tmpMask );
+        psFree( tmpMask2 );
+    }
+    tmpImage = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_F32 );
+    myHist = psHistogramAlloc( 0.0, ( float ) ( IMAGE_X_SIZE + IMAGE_Y_SIZE ),
+                               NUM_BINS );
+    myHist2 = psImageHistogram( NULL, tmpImage, NULL, 0 );
+    if ( myHist2 != NULL ) {
+        fprintf(stderr, "ERROR: myHist2 not equal to NULL\n" );
+    }
+
+    myHist2 = psImageHistogram( myHist, NULL, NULL, 0 );
+    myHist2 = psImageHistogram( NULL, tmpImage, NULL, 0 );
+    if ( myHist2 != NULL ) {
+        fprintf(stderr, "ERROR: myHist2 not equal to NULL\n" );
+    }
+
+
+    tmpMask = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_S8 );
+    for ( i = 0;i < tmpMask->numRows;i++ ) {
+        for ( j = 0;j < tmpMask->numCols;j++ ) {
+            if ( ( i > ( tmpMask->numRows / 2 ) ) &&
+                    ( j > ( tmpMask->numCols / 2 ) ) ) {
+                tmpMask->data.S8[ i ][ j ] = 1;
+            } else {
+                tmpMask->data.S8[ i ][ j ] = 0;
+            }
+        }
+    }
+
+    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for invalid mask type.");
+    myHist2 = psImageHistogram(myHist,tmpImage,tmpMask,1);
+    if (myHist2 != NULL) {
+        fprintf(stderr,"ERROR: myHist2 not equal to NULL\n");
+    }
+
+    psFree( tmpMask );
+    psFree(myHist);
+    psFree(myHist2);
+    psFree( tmpImage );
+
+    return ( !testStatus );
+}
+
+static psS32 testPsImageStats()
+{
+    const int N  = 32;
+    const int M = 64;
+
+    psStats * myStats = NULL;
+    psStats *myStats2 = NULL;
+    psImage *tmpImage = NULL;
+    psImage *tmpImage2 = NULL;
+    psImage *tmpMask = NULL;
+    psImage *tmpMask2 = NULL;
+    psS32 testStatus = true;
+    psS32 nb = 0;
+    psS32 i = 0;
+    psS32 j = 0;
+    psS32 IMAGE_X_SIZE = 0;
+    psS32 IMAGE_Y_SIZE = 0;
+    psS32 currentId = 0;
+
+    currentId = psMemGetId();
+    for ( nb = 0;nb < 7;nb++ ) {
+        if ( nb == 0 ) {
+            IMAGE_X_SIZE = 1;
+            IMAGE_Y_SIZE = 1;
+        }
+        if ( nb == 1 ) {
+            IMAGE_X_SIZE = 1;
+            IMAGE_Y_SIZE = N;
+        }
+        if ( nb == 2 ) {
+            IMAGE_X_SIZE = N;
+            IMAGE_Y_SIZE = 1;
+        }
+        if ( nb == 3 ) {
+            IMAGE_X_SIZE = N;
+            IMAGE_Y_SIZE = N;
+        }
+        if ( nb == 4 ) {
+            IMAGE_X_SIZE = N;
+            IMAGE_Y_SIZE = M;
+        }
+        if ( nb == 5 ) {
+            IMAGE_X_SIZE = M;
+            IMAGE_Y_SIZE = N;
+        }
+        if ( nb == 6 ) {
+            IMAGE_X_SIZE = M;
+            IMAGE_Y_SIZE = N;
+        }
+        fprintf(stderr, "*******************************\n" );
+        if (nb != 6) {
+            fprintf(stderr, "* IMAGE SIZE is (%d by %d)\n", IMAGE_X_SIZE, IMAGE_Y_SIZE );
+        } else {
+            fprintf(stderr, "* IMAGE SUBSET SIZE is (%d by %d)\n", IMAGE_X_SIZE, IMAGE_Y_SIZE );
+        }
+        fprintf(stderr, "*******************************\n" );
+
+        /*************************************************************************/
+        /*  Allocate and initialize data structures                      */
+        /*************************************************************************/
+        if (nb != 6) {
+            tmpImage = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_F32 );
+            tmpMask = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_U8 );
+            tmpImage2 = NULL;
+            tmpMask2 = NULL;
+        } else {
+            tmpImage2 = psImageAlloc( IMAGE_X_SIZE*2, IMAGE_Y_SIZE*2, PS_TYPE_F32 );
+            tmpMask2 = psImageAlloc( IMAGE_X_SIZE*2, IMAGE_Y_SIZE*2, PS_TYPE_U8 );
+            tmpImage = psImageSubset(tmpImage2,
+                                     IMAGE_X_SIZE/4, IMAGE_Y_SIZE/4,
+                                     IMAGE_X_SIZE/4+IMAGE_X_SIZE, IMAGE_Y_SIZE/4+IMAGE_Y_SIZE);
+            tmpMask = psImageSubset(tmpMask2,
+                                    IMAGE_X_SIZE/4, IMAGE_Y_SIZE/4,
+                                    IMAGE_X_SIZE/4+IMAGE_X_SIZE, IMAGE_Y_SIZE/4+IMAGE_Y_SIZE);
+        }
+
+        for ( i = 0;i < tmpImage->numRows;i++ ) {
+            for ( j = 0;j < tmpImage->numCols;j++ ) {
+                tmpImage->data.F32[ i ][ j ] = ( float ) ( i + j );
+            }
+        }
+        for ( i = 0;i < tmpMask->numRows;i++ ) {
+            for ( j = 0;j < tmpMask->numCols;j++ ) {
+                if ( ( i > ( tmpMask->numRows / 2 ) ) &&
+                        ( j > ( tmpMask->numCols / 2 ) ) ) {
+                    tmpMask->data.U8[ i ][ j ] = 1;
+                } else {
+                    tmpMask->data.U8[ i ][ j ] = 0;
+                }
+            }
+        }
+
+        myStats = psStatsAlloc( PS_STAT_SAMPLE_MEAN );
+        /*************************************************************************/
+        /*  Calculate Sample Mean with no mask                           */
+        /*************************************************************************/
+        psStats* tmpStats = myStats;
+        myStats = psImageStats( myStats, tmpImage, NULL, 0 );
+        if ( myStats != tmpStats ) {
+            fprintf(stderr,"ERROR: input psStats not equal to return psStats\n");
+            psAbort(__func__,"Failed input psStats equal to returned psStats");
+        }
+        fprintf(stderr, "The sample mean was %.2f\n", myStats->sampleMean );
+
+        /*************************************************************************/
+        /*  Calculate Sample Mean with mask                              */
+        /*************************************************************************/
+
+        myStats = psImageStats( myStats, tmpImage, tmpMask, 1 );
+        fprintf(stderr, "The sample mean was %.2f\n", myStats->sampleMean );
+
+        /*************************************************************************/
+        /*  Deallocate data structures                                   */
+        /*************************************************************************/
+        psFree( myStats );
+        psFree( tmpImage );
+        psFree( tmpMask );
+        psFree( tmpImage2 );
+        psFree( tmpMask2 );
+    }
+
+    /*************************************************************************/
+    /*  Test With Various Null Inputs                                        */
+    /*************************************************************************/
+
+    tmpImage = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_F32 );
+    tmpMask = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_U8 );
+    myStats = psStatsAlloc( 0 );
+
+    myStats2 = psImageStats( myStats, NULL, NULL, 0 );
+    if ( myStats2 != NULL ) {
+        fprintf(stderr, "ERROR: myStats2 = psImageStats(myStats, NULL, NULL, 0) != NULL\n" );
+    }
+
+    myStats2 = psImageStats( NULL, tmpImage, NULL, 0 );
+    if ( myStats2 != NULL ) {
+        fprintf(stderr, "ERROR: myStats2 = psImageStats(NULL, tmpImage, NULL, 0) != NULL\n" );
+    }
+
+    myStats2 = psImageStats( myStats, tmpImage, NULL, 0 );
+
+    psFree( myStats );
+    psFree( tmpImage );
+    psFree( tmpMask );
+
+    tmpImage = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_F32);
+    tmpMask = psImageAlloc( IMAGE_X_SIZE, IMAGE_Y_SIZE, PS_TYPE_U32);
+    myStats = psStatsAlloc( PS_STAT_SAMPLE_MEAN );
+    myStats2 = psImageStats( myStats, tmpImage, tmpMask, 0 );
+    if ( myStats2 != NULL ) {
+        fprintf(stderr, "ERROR: psImageStats did not return null when mask is invalid type.\n");
+    }
+    psFree(tmpImage);
+    psFree(tmpMask);
+    psFree(myStats);
+
+    return ( !testStatus );
+}
+
+static psS32 testPsImageFitPolynomial()
+{
+    const int IMAGE_SIZE = 64;
+    const int CHEBY_X_DIM = 8;
+    const int CHEBY_Y_DIM = 8;
+    const int THRESHOLD = 10;
+
+    psStats * myStats = NULL;
+    psImage *tmpImage = NULL;
+    psImage *outImage = NULL;
+    psImage *nullImage = NULL;
+    psPolynomial2D *my2DPoly = NULL;
+    psPolynomial2D *null2DPoly = NULL;
+    psS32 testStatus = true;
+    psS32 i = 0;
+    psS32 j = 0;
+    psS32 currentId = 0;
+
+    currentId = psMemGetId();
+    /*************************************************************************/
+    /*  Allocate and initialize data structures                      */
+    /*************************************************************************/
+    tmpImage = psImageAlloc( IMAGE_SIZE, IMAGE_SIZE, PS_TYPE_F32 );
+    outImage = psImageAlloc( IMAGE_SIZE, IMAGE_SIZE, PS_TYPE_F32 );
+    for ( i = 0;i < IMAGE_SIZE;i++ ) {
+        for ( j = 0;j < IMAGE_SIZE;j++ ) {
+            tmpImage->data.F32[ i ][ j ] = 4.0;
+            tmpImage->data.F32[ i ][ j ] = ( float ) ( i + j );
+            tmpImage->data.F32[ i ][ j ] = ( float ) ( i + j ) + ( 4.0 * ( float ) i );
+            outImage->data.F32[ i ][ j ] = 0.0;
+        }
+    }
+    my2DPoly = psPolynomial2DAlloc( CHEBY_X_DIM, CHEBY_Y_DIM, PS_POLYNOMIAL_CHEB );
+    /*************************************************************************/
+    /*  Calculate Chebyshev Polynomials, no mask                     */
+    /*************************************************************************/
+    my2DPoly = psImageFitPolynomial( my2DPoly, tmpImage );
+    for ( i = 0;i < CHEBY_X_DIM;i++ ) {
+        for ( j = 0;j < CHEBY_Y_DIM;j++ ) {
+            fprintf(stderr, "Cheby Polynomial (%d, %d) coefficient is %.2f\n", i, j, 0.0001f+my2DPoly->coeff[ i ][ j ] );
+        }
+    }
+
+    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error message for NULL coeffs argument.");
+    null2DPoly = psImageFitPolynomial( NULL, tmpImage );
+    if ( null2DPoly != NULL ) {
+        fprintf(stderr,"ERROR: psImageFitPolynomial did not return NULL with invalid coeffs argument.");
+    }
+
+    psLogMsg(__func__,PS_LOG_INFO,"Following should be a error message for NULL input argument.");
+    null2DPoly = psImageFitPolynomial( null2DPoly, NULL );
+    if ( null2DPoly != NULL ) {
+        fprintf(stderr,"ERROR: psImageFitPolynomial did not return NULL with null input image.");
+    }
+
+    /*************************************************************************/
+    /*  Evaluate Chebyshev Polynomials, no mask                      */
+    /*************************************************************************/
+    outImage = psImageEvalPolynomial( outImage, my2DPoly );
+    for ( i = 0;i < IMAGE_SIZE;i++ ) {
+        for ( j = 0;j < IMAGE_SIZE;j++ ) {
+
+            //             fprintf(stderr,"pixel[%d][%d] is (%.2f, %.2f)\n", i, j,
+            //                     tmpImage->data.F32[i][j],
+            //                     outImage->data.F32[i][j]);
+            if ( fabs( outImage->data.F32[ i ][ j ] - tmpImage->data.F32[ i ][ j ] ) > THRESHOLD ) {
+                fprintf(stderr, "Pixel (%d, %d) is %.2f, should be %.2f\n", i, j,
+                        outImage->data.F32[ i ][ j ],
+                        tmpImage->data.F32[ i ][ j ] );
+            }
+
+        }
+    }
+
+
+    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error message for NULL coeffs argument.");
+    nullImage = psImageEvalPolynomial( outImage, NULL );
+    if ( nullImage != NULL ) {
+        fprintf(stderr,"ERROR: psImageEvalPolynomial did not return NULL with invalid coeffs argument.");
+    }
+
+    psLogMsg(__func__,PS_LOG_INFO,"Following should be a error message for NULL input argument.");
+    nullImage = psImageEvalPolynomial( NULL, my2DPoly );
+    if ( nullImage != NULL ) {
+        fprintf(stderr,"ERROR: psImageEvalPolynomial did not return NULL with null input image.");
+    }
+
+    /*************************************************************************/
+    /*  Deallocate data structures                                   */
+    /*************************************************************************/
+    psFree( myStats );
+    psFree( tmpImage );
+    psFree( outImage );
+    psFree( my2DPoly );
+
+    return ( !testStatus );
+}
+
+static psS32 testPsImagePixelInterpolate()
+{
+    const int IMAGE_SIZE = 10;
+
+    psImage *tmpImage   = NULL;
+    psS32 testStatus      = true;
+    psS32 i               = 0;
+    psS32 j               = 0;
+    float pixel         = 0.0;
+    float x             = 0.0;
+    float y             = 0.0;
+
+    /*************************************************************************/
+    /*  Allocate and initialize data structures                      */
+    /*************************************************************************/
+    tmpImage = psImageAlloc(IMAGE_SIZE, IMAGE_SIZE, PS_TYPE_F32);
+    for (i=0;i<IMAGE_SIZE;i++) {
+        for (j=0;j<IMAGE_SIZE;j++) {
+            tmpImage->data.F32[i][j] = 4.0;
+            tmpImage->data.F32[i][j] = (float) (i + j) + (4.0 * (float) i);
+            tmpImage->data.F32[i][j] = (float) (i + j);
+            fprintf(stderr,"%.1f ", tmpImage->data.F32[i][j]);
+        }
+        fprintf(stderr,"\n");
+    }
+    for (i=0;i<IMAGE_SIZE-1;i++) {
+        for (j=0;j<IMAGE_SIZE-1;j++) {
+            x = 0.2 + (float) i;
+            y = 0.2 + (float) j;
+            pixel = psImagePixelInterpolate(tmpImage,
+                                            x, y,
+                                            NULL, 0,
+                                            0,
+                                            PS_INTERPOLATE_BILINEAR);
+            fprintf(stderr,"%.1f ", pixel);
+        }
+        fprintf(stderr,"\n");
+    }
+
+    for (i=0;i<IMAGE_SIZE-1;i++) {
+        for (j=0;j<IMAGE_SIZE-1;j++) {
+            x = 0.2 + (float) i;
+            y = 0.2 + (float) j;
+            pixel = psImagePixelInterpolate(tmpImage,
+                                            x, y,
+                                            NULL,0,
+                                            0,
+                                            PS_INTERPOLATE_BILINEAR);
+            fprintf(stderr,"image[%.1f][%.1f] is interpolated at %.1f\n", x, y, pixel);
+        }
+    }
+
+    psFree(tmpImage);
+
+    return (!testStatus);
+}
+
+#define I_POLY(X, Y) \
+((A) + (B * X) + (C * Y) + (D * X * Y) + (E * X * X) + (F * Y * Y)) \
+
+static bool FitChebyF32(int numCols, int numRows)
+{
+    const double A = 1.0;
+    const double B = 2.0;
+    const double C = 3.0;
+    const double D = 4.0;
+    const double E = 5.0;
+    const double F = 6.0;
+    const double ERROR_TOL = 0.1;
+
+    psImage *tmpImage     = NULL;
+    psImage *outImage     = NULL;
+    bool testStatus      = true;
+    psS32 i               = 0;
+    psS32 j               = 0;
+    float chi2 = 0.0;
+
+
+    fprintf(stderr,"psImageFitPolynomial(), psImageEvalPolynom(): (%d by %d)\n", numRows, numCols);
+
+    /*************************************************************************/
+    /*  Allocate and initialize data structures                      */
+    /*************************************************************************/
+
+    tmpImage = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    outImage = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    for (i=0;i<numRows;i++) {
+        for (j=0;j<numCols;j++) {
+            tmpImage->data.F32[i][j] = I_POLY(((float) i), ((float) j));
+        }
+    }
+
+    psPolynomial2D *chebys = psPolynomial2DAlloc(5, 5, PS_POLYNOMIAL_CHEB);
+
+    chebys = psImageFitPolynomial(chebys, tmpImage);
+
+    outImage = psImageEvalPolynomial(outImage, chebys);
+
+    chi2 = 0.0;
+    for (i=0;i<numRows;i++) {
+        for (j=0;j<numCols;j++) {
+            float expect = tmpImage->data.F32[i][j];
+            float actual = outImage->data.F32[i][j];
+            chi2+= (actual-expect) * (actual-expect);
+            if (fabs((actual - expect) / expect) > ERROR_TOL) {
+                fprintf(stderr,"pixel [%d][%d] is %.2f should be %.2f\n", i, j, actual, expect);
+            }
+        }
+    }
+    chi2/= ((float) (numCols * numRows));
+    fprintf(stderr,"The chi-squared per pixel is %.2f\n", chi2);
+
+    psFree(tmpImage);
+    psFree(outImage);
+    psFree(chebys);
+
+    /*************************************************************************/
+    /*  Deallocate data structures                                   */
+    /*************************************************************************/
+
+    return(testStatus);
+}
+
+static bool FitChebyF64(int numCols, int numRows)
+{
+    const double A = 1.0;
+    const double B = 2.0;
+    const double C = 3.0;
+    const double D = 4.0;
+    const double E = 5.0;
+    const double F = 6.0;
+    const double ERROR_TOL = 0.1;
+
+    psImage *tmpImage     = NULL;
+    psImage *outImage     = NULL;
+    bool testStatus      = true;
+    psS32 i               = 0;
+    psS32 j               = 0;
+    double chi2 = 0.0;
+
+
+    fprintf(stderr,"psImageFitPolynomial(), psImageEvalPolynom(): (%d by %d)\n", numRows, numCols);
+
+    /*************************************************************************/
+    /*  Allocate and initialize data structures                      */
+    /*************************************************************************/
+
+    tmpImage = psImageAlloc(numCols, numRows, PS_TYPE_F64);
+    outImage = psImageAlloc(numCols, numRows, PS_TYPE_F64);
+    for (i=0;i<numRows;i++) {
+        for (j=0;j<numCols;j++) {
+            tmpImage->data.F64[i][j] = I_POLY(((double) i), ((double) j));
+        }
+    }
+
+    psPolynomial2D *chebys = psPolynomial2DAlloc(5, 5, PS_POLYNOMIAL_CHEB);
+
+    chebys = psImageFitPolynomial(chebys, tmpImage);
+
+    outImage = psImageEvalPolynomial(outImage, chebys);
+
+    chi2 = 0.0;
+    for (i=0;i<numRows;i++) {
+        for (j=0;j<numCols;j++) {
+            double expect = tmpImage->data.F64[i][j];
+            double actual = outImage->data.F64[i][j];
+            chi2+= (actual-expect) * (actual-expect);
+            if (fabs((actual - expect) / expect) > ERROR_TOL) {
+                fprintf(stderr,"pixel [%d][%d] is %.2f should be %.2f\n", i, j, actual, expect);
+            }
+        }
+    }
+    chi2/= ((double) (numCols * numRows));
+    fprintf(stderr,"The chi-squared per pixel is %.2f\n", chi2);
+
+    psFree(tmpImage);
+    psFree(outImage);
+    psFree(chebys);
+
+    /*************************************************************************/
+    /*  Deallocate data structures                                   */
+    /*************************************************************************/
+
+    return(testStatus);
+}
+
+static psS32 testPsImageEvalPolynom()
+{
+    if (!FitChebyF32(1, 1)) {
+        return 1;
+    }
+    if (!FitChebyF32(1, 5)) {
+        return 2;
+    }
+    if (!FitChebyF32(5, 1)) {
+        return 3;
+    }
+    if (!FitChebyF32(5, 5)) {
+        return 4;
+    }
+    if (!FitChebyF64(1, 1)) {
+        return 5;
+    }
+    if (!FitChebyF64(1, 5)) {
+        return 6;
+    }
+    if (!FitChebyF64(5, 1)) {
+        return 7;
+    }
+    if (!FitChebyF64(5, 5)) {
+        return 8;
+    }
+    return 0;
+}
+