/*****************************************************************************
This routine must ensure that the 
    psImageFitPolynomial()
    psImageEvalPolynomial()
functions work correctly.
 
XXX: We do not check image pixels.
 *****************************************************************************/
#include <stdio.h>
#include "pslib.h"
#include "psTest.h"
#include "psMemory.h"
#include "psImage.h"
#include "psImageStats.h"
#include <float.h>
#include <math.h>
#define IMAGE_SIZE 5
#define NUM_COLS IMAGE_SIZE
#define NUM_ROWS IMAGE_SIZE
#define ERROR_TOL 0.1
#define A 1.0
#define B 2.0
#define C 3.0
#define D 4.0
#define E 5.0
#define F 6.0

#define I_POLY(X, Y) \
((A) + (B * X) + (C * Y) + (D * X * Y) + (E * X * X) + (F * Y * Y)) \

psS32 FitChebyF32(int numCols, int numRows)
{
    psImage *tmpImage     = NULL;
    psImage *outImage     = NULL;
    psS32 testStatus      = true;
    psS32 memLeaks        = 0;
    psS32 i               = 0;
    psS32 j               = 0;
    psS32 currentId       = psMemGetId();
    float chi2 = 0.0;


    printf("psImageFitPolynomial(), psImageEvalPolynom(): (%d by %d)\n", numRows, numCols);

    /*************************************************************************/
    /*  Allocate and initialize data structures                      */
    /*************************************************************************/

    printPositiveTestHeader(stdout,
                            "psImageStats functions",
                            "psImageFitPolynomial(), psImageEvalPolynom()");

    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) {
                printf("pixel [%d][%d] is %.2f should be %.2f\n", i, j, actual, expect);
            }
        }
    }
    chi2/= ((float) (numCols * numRows));
    printf("The chi-squared per pixel is %f\n", chi2);

    psFree(tmpImage);
    psFree(outImage);
    psFree(chebys);
    psMemCheckCorruption(1);
    printFooter(stdout,
                "psImageStats functions",
                "psImagePixelInterpolation()",
                testStatus);

    /*************************************************************************/
    /*  Deallocate data structures                                   */
    /*************************************************************************/
    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }

    return(!testStatus);
}

psS32 FitChebyF64(int numCols, int numRows)
{
    psImage *tmpImage     = NULL;
    psImage *outImage     = NULL;
    psS32 testStatus      = true;
    psS32 memLeaks        = 0;
    psS32 i               = 0;
    psS32 j               = 0;
    psS32 currentId       = psMemGetId();
    double chi2 = 0.0;


    printf("psImageFitPolynomial(), psImageEvalPolynom(): (%d by %d)\n", numRows, numCols);

    /*************************************************************************/
    /*  Allocate and initialize data structures                      */
    /*************************************************************************/

    printPositiveTestHeader(stdout,
                            "psImageStats functions",
                            "psImageFitPolynomial(), psImageEvalPolynom()");

    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) {
                printf("pixel [%d][%d] is %.2f should be %.2f\n", i, j, actual, expect);
            }
        }
    }
    chi2/= ((double) (numCols * numRows));
    printf("The chi-squared per pixel is %f\n", chi2);

    psFree(tmpImage);
    psFree(outImage);
    psFree(chebys);
    psMemCheckCorruption(1);
    printFooter(stdout,
                "psImageStats functions",
                "psImagePixelInterpolation()",
                testStatus);

    /*************************************************************************/
    /*  Deallocate data structures                                   */
    /*************************************************************************/
    psMemCheckCorruption(1);
    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    if (0 != memLeaks) {
        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
    }

    return(!testStatus);
}

int main()
{
    FitChebyF32(1, 1);
    FitChebyF32(1, NUM_ROWS);
    FitChebyF32(1, 1);
    FitChebyF32(NUM_COLS, 1);
    FitChebyF32(NUM_COLS, NUM_ROWS);
    FitChebyF64(1, 1);
    FitChebyF64(1, NUM_ROWS);
    FitChebyF64(1, 1);
    FitChebyF64(NUM_COLS, 1);
    FitChebyF64(NUM_COLS, NUM_ROWS);
}
