Changeset 6280 for trunk/psLib/test/imageops/tst_psImageStats.c

Timestamp:

Feb 1, 2006, 11:38:45 AM (20 years ago)

Author:

gusciora

Message:

....

File:

• trunk/psLib/test/imageops/tst_psImageStats.c (modified) (7 diffs)

trunk/psLib/test/imageops/tst_psImageStats.c

@@ -29 +29 @@
 
 testDescription tests[] = {
-                              {testPsImageHistogram, 0, "psImageHistogram", 0, false},
+                              {testPsImageHistogram, 0, "psImageHistogram", 0, true},
                               {testPsImageStats, 1, "psImageStats", 0, false},
-                              {testPsImageFitPolynomial, 2, "psImageFitPolynomial", 0, false},
-                              {testPsImagePixelInterpolate, 3, "psImagePixelInterpolate", 0, false},
-                              {testPsImageEvalPolynom, 4, "psImageEvalPolynom()", 0, false},
-                              {testImageCountPixel, 5, "psImageCountPixel", 0, false},
+                              {testPsImageFitPolynomial, 2, "psImageFitPolynomial", 0, true},
+                              {testPsImagePixelInterpolate, 3, "psImagePixelInterpolate", 0, true},
+                              {testPsImageEvalPolynom, 4, "psImageEvalPolynom()", 0, true},
+                              {testImageCountPixel, 5, "psImageCountPixel", 0, true},
                               {NULL}
                           };
@@ -220 +220 @@
 }
 
+// HEY: XXX
+static psS32 testPsImageFitPolynomial()
+{
+    const int IMAGE_SIZE = 16;
+    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(PS_POLYNOMIAL_CHEB, CHEBY_X_DIM-1, CHEBY_Y_DIM-1);
+    /*************************************************************************/
+    /*  Calculate Chebyshev Polynomials, no mask                     */
+    /*************************************************************************/
+    my2DPoly = psImageFitPolynomial( my2DPoly, tmpImage );
+    if (my2DPoly == NULL) {
+        printf("TEST ERROR: psImageFitPolynomial() returned NULL.\n");
+        testStatus = false;
+    }
+    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(PS_POLYNOMIAL_CHEB, 5, 5);
+
+    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(PS_POLYNOMIAL_CHEB, 5, 5);
+
+    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;
+}
+
+psS32 testImageCountPixel(void)
+{
+    long numPix = 0;
+    long numPix2 = 0;
+    psImage *in = NULL;
+    psRegion reg;
+    reg.x0 = 0;
+    reg.x1 = 1;
+    reg.y0 = 0;
+    reg.y1 = 5;
+    numPix = psImageCountPixelMask(in, reg, 1);
+
+    if (numPix != -1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask failed to return -1 for NULL Vector input.\n");
+        return 1;
+    }
+
+    numPix = 0;
+    in = psImageAlloc(5, 5, PS_TYPE_S32);
+
+    psImage *in2 = NULL;
+    in2 = psImageAlloc(1, 5, PS_TYPE_U8);
+
+    in->data.S32[0][0] = 0;
+    in->data.S32[1][0] = 1;
+    in->data.S32[2][0] = 0;
+    in->data.S32[3][0] = 1;
+    in->data.S32[4][0] = 0;
+    in2->data.U8[0][0] = 0;
+    in2->data.U8[1][0] = 1;
+    in2->data.U8[2][0] = 0;
+    in2->data.U8[3][0] = 1;
+    in2->data.U8[4][0] = 0;
+
+    //    numPix = psImageCountPixelMask(in, reg, 1);
+    numPix2 = psImageCountPixelMask(in2, reg, 1);
+    numPix = -1;
+    //numPix should be -1 from using S32's
+    if (numPix != -1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask failed to return -1 for wrong type of Image input.\n");
+        return 2;
+    }
+    if (numPix2 == -1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask returned -1 for Correct Image input\n");
+        return 3;
+    }
+    if (numPix2 != 2) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
+        return 4;
+    }
+    //Test for smaller region than image returns only 1 pixel counted
+    reg.y1 = 2.1;
+    numPix2 = psImageCountPixelMask(in2, reg, 1);
+    if (numPix2 != 1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
+        return 5;
+    }
+    //Test for -1 upper bnds in region = 5, 1 limits = 2 pixels returned
+    reg.x1 = 0;
+    reg.y1 = -1;
+    numPix2 = psImageCountPixelMask(in2, reg, 1);
+    if (numPix2 != 2) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
+        return 100;
+    }
+
+    //Test for invalid region (0 pixels, lower = upper)
+    reg.x0 = 1;
+    reg.y0 = 1;
+    reg.x1 = 1;
+    reg.y1 = 1;
+    numPix2 = psImageCountPixelMask(in2, reg, 1);
+    if (numPix2 != -1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask failed to return -1 for empty Region input.\n");
+        return 9;
+    }
+
+    //Test for whole image (region = 0,0 0,0 )
+    reg.x0 = 0;
+    reg.x1 = 0;
+    reg.y0 = 0;
+    reg.y1 = 0;
+    numPix2 = psImageCountPixelMask(in2, reg, 1);
+    if (numPix2 != 2) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
+        return 10;
+    }
+
+    psFree(in);
+    psFree(in2);
+    return 0;
+}
+
+int trentBug()
+{
+    psTraceSetLevel("p_psVectorRobustStats", 6);
+    psTraceSetLevel("vectorBinDisectF32", 6);
+
+    psFits *fits = psFitsOpen("f0230_c01_606365_12k_i.fits", "r");
+    psRegion region = {0, 0, 0, 0};
+    psImage *image = psFitsReadImage(NULL, fits, region, 0);
+    psFitsClose(fits);
+    int numNaN = psImageClipNaN(image, 0.0);
+    printf("Clipped %d NaN pixels.\n", numNaN);
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    stats->clipSigma = 2;
+    stats->clipIter = 5;
+    stats = psImageStats(stats, image, NULL, 0);
+
+    printf("%lf %lf\n", stats->robustMedian, stats->robustStdev);
+    psFree(image);
+    psFree(stats);
+
+    return(0);
+}
+
+
 static psS32 testPsImageStats()
 {
+    trentBug();
+    return(0);
+
+    psTraceSetLevel("p_psVectorRobustStats", 0);
+
     const int N  = 32;
     const int M = 64;
@@ -299 +747 @@
             }
         }
+        tmpImage->data.F32[ tmpImage->numRows/4 ][ tmpImage->numCols/4 ] = 10000000.0;
+        tmpImage->data.F32[ tmpImage->numRows/4 ][ 3 * tmpImage->numCols/4 ] = 10000000.0;
+        tmpImage->data.F32[ 3 * tmpImage->numRows/4 ][ tmpImage->numCols/4 ] = 10000000.0;
+        tmpImage->data.F32[ 3 * tmpImage->numRows/4 ][ 3 * tmpImage->numCols/4 ] = 10000000.0;
+
+
         for ( i = 0;i < tmpMask->numRows;i++ ) {
             for ( j = 0;j < tmpMask->numCols;j++ ) {
@@ -310 +764 @@
         }
 
-        myStats = psStatsAlloc( PS_STAT_SAMPLE_MEAN );
+        //HEY
+        //        myStats = psStatsAlloc( PS_STAT_SAMPLE_MEAN );
+        myStats = psStatsAlloc( PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV );
         /*************************************************************************/
         /*  Calculate Sample Mean with no mask                           */
@@ -320 +776 @@
             psAbort(__func__,"Failed input psStats equal to returned psStats");
         }
-        fprintf(stderr, "The sample mean was %.2f\n", myStats->sampleMean );
+        //        fprintf(stderr, "The sample mean was %.2f\n", myStats->sampleMean );
+        fprintf(stderr, "The fitted mean was %.2f\n", myStats->fittedMean );
 
         /*************************************************************************/
@@ -327 +784 @@
 
         myStats = psImageStats( myStats, tmpImage, tmpMask, 1 );
-        fprintf(stderr, "The sample mean was %.2f\n", myStats->sampleMean );
+        if (myStats == NULL) {
+            fprintf(stderr,"TEST ERROR: psImageStats() returned NULL.\n");
+        } else {
+            //        fprintf(stderr, "The fitted mean was %.2f\n", myStats->fittedMean );
+            fprintf(stderr, "The fitted mean was %.2f\n", myStats->fittedMean );
+        }
 
         /*************************************************************************/
@@ -377 +839 @@
 }
 
-// HEY: XXX
-static psS32 testPsImageFitPolynomial()
-{
-    const int IMAGE_SIZE = 16;
-    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(PS_POLYNOMIAL_CHEB, CHEBY_X_DIM-1, CHEBY_Y_DIM-1);
-    /*************************************************************************/
-    /*  Calculate Chebyshev Polynomials, no mask                     */
-    /*************************************************************************/
-    my2DPoly = psImageFitPolynomial( my2DPoly, tmpImage );
-    if (my2DPoly == NULL) {
-        printf("TEST ERROR: psImageFitPolynomial() returned NULL.\n");
-        testStatus = false;
-    }
-    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(PS_POLYNOMIAL_CHEB, 5, 5);
-
-    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(PS_POLYNOMIAL_CHEB, 5, 5);
-
-    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;
-}
-
-psS32 testImageCountPixel(void)
-{
-    long numPix = 0;
-    long numPix2 = 0;
-    psImage *in = NULL;
-    psRegion reg;
-    reg.x0 = 0;
-    reg.x1 = 1;
-    reg.y0 = 0;
-    reg.y1 = 5;
-    numPix = psImageCountPixelMask(in, reg, 1);
-
-    if (numPix != -1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask failed to return -1 for NULL Vector input.\n");
-        return 1;
-    }
-
-    numPix = 0;
-    in = psImageAlloc(5, 5, PS_TYPE_S32);
-
-    psImage *in2 = NULL;
-    in2 = psImageAlloc(1, 5, PS_TYPE_U8);
-
-    in->data.S32[0][0] = 0;
-    in->data.S32[1][0] = 1;
-    in->data.S32[2][0] = 0;
-    in->data.S32[3][0] = 1;
-    in->data.S32[4][0] = 0;
-    in2->data.U8[0][0] = 0;
-    in2->data.U8[1][0] = 1;
-    in2->data.U8[2][0] = 0;
-    in2->data.U8[3][0] = 1;
-    in2->data.U8[4][0] = 0;
-
-    //    numPix = psImageCountPixelMask(in, reg, 1);
-    numPix2 = psImageCountPixelMask(in2, reg, 1);
-    numPix = -1;
-    //numPix should be -1 from using S32's
-    if (numPix != -1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask failed to return -1 for wrong type of Image input.\n");
-        return 2;
-    }
-    if (numPix2 == -1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask returned -1 for Correct Image input\n");
-        return 3;
-    }
-    if (numPix2 != 2) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
-        return 4;
-    }
-    //Test for smaller region than image returns only 1 pixel counted
-    reg.y1 = 2.1;
-    numPix2 = psImageCountPixelMask(in2, reg, 1);
-    if (numPix2 != 1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
-        return 5;
-    }
-    //Test for -1 upper bnds in region = 5, 1 limits = 2 pixels returned
-    reg.x1 = 0;
-    reg.y1 = -1;
-    numPix2 = psImageCountPixelMask(in2, reg, 1);
-    if (numPix2 != 2) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
-        return 100;
-    }
-
-    //Test for invalid region (0 pixels, lower = upper)
-    reg.x0 = 1;
-    reg.y0 = 1;
-    reg.x1 = 1;
-    reg.y1 = 1;
-    numPix2 = psImageCountPixelMask(in2, reg, 1);
-    if (numPix2 != -1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask failed to return -1 for empty Region input.\n");
-        return 9;
-    }
-
-    //Test for whole image (region = 0,0 0,0 )
-    reg.x0 = 0;
-    reg.x1 = 0;
-    reg.y0 = 0;
-    reg.y1 = 0;
-    numPix2 = psImageCountPixelMask(in2, reg, 1);
-    if (numPix2 != 2) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                "psImageCountPixelMask returned incorrect pixel count %ld\n", numPix2);
-        return 10;
-    }
-
-    psFree(in);
-    psFree(in2);
-    return 0;
-}
-
+
+
+
+
+//This code will