/** @file  tst_psAstrometry01.c
*
*  @brief This code will test the code in psAstrometry.[ch]
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-07-20 01:21:13 $
*
* XXX: Must test
* pmSourceAddModel
* pmSourceSubModel
* pmSourceFitModel
* pmSourceContour
* pmSourceModelGuess
* pmSourceRoughClass
* pmSourceMoments
* pmSourceLocalSky
*
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

#include <stdlib.h>

#include "psTest.h"
#include "pslib_strict.h"

static psS32 test1( void );
static psS32 test2( void );
static psS32 test3( void );
static psS32 test4( void );
static psS32 test5( void );
static psS32 test6( void );

testDescription tests[] = {
                              {test1, -1, "Tests psFixedPatternAlloc()", 0, false},
                              {test2, -2, "Tests psGrommitAlloc()", 0, false},
                              {test3, -3, "MISC", 0, false},
                              {test4, -1, "psPlaneTransformCombine()", 0, false},
                              {test5, -1, "psPlaneTransformFit()", 0, false},
                              {test6, -1, "psPlaneTransformInvert()", 0, false},
                              {NULL}
                          };

#define PS_PERCENT_COMPARE(X, Y, PERCENT_FRACTION) \
(fabs((Y)-(X))/fabs(X) < (PERCENT_FRACTION))

#define RA  1.0
#define DEC  2.0
#define HA  3.0
#define ZD  4.0
#define AZ  5.0
#define ROT_ANGLE 6.0
#define TEMPERATURE 7.0
#define PRESSURE 8.0
#define HUMIDITY 20.0
#define EXP_TIME 22.0
#define WAVELENGTH 11.0
#define NAME  "The Kaiser Royal Observatory"
#define LATITUDE 20.7
#define LONGITUDE 156.3
#define HEIGHT  3055.0
#define LAPSERATE 20.0

#define NUM_READOUTS 2
#define READOUT_NUM_ROWS 8
#define READOUT_NUM_COLS 8

#define CELL_WIDTH READOUT_NUM_COLS
#define CELL_HEIGHT READOUT_NUM_ROWS

#define CELL_GAP 2

#define NUM_CELLS 4
#define CELL_MIN_X 0
#define CELL_MAX_X CELL_WIDTH
#define CELL_MIN_Y 0
#define CELL_MAX_Y CELL_HEIGHT

#define NUM_CHIPS 2
#define CHIP_MIN_X 0
#define CHIP_MAX_X ((NUM_CELLS * CELL_WIDTH) + ((NUM_CELLS) * CELL_GAP))
#define CHIP_MIN_Y 0
#define CHIP_MAX_Y CELL_HEIGHT

#define CHIP_WIDTH CHIP_MAX_X
#define CHIP_HEIGHT CHIP_MAX_Y
#define CHIP_GAP 2

#define FPA_MIN_X 0
#define FPA_MAX_X ((NUM_CHIPS * CHIP_WIDTH) + ((NUM_CHIPS) * CHIP_GAP))
#define FPA_MIN_Y 0
#define FPA_MAX_Y CHIP_HEIGHT

#define ROW0  0
#define COL0  0

static void psPlaneDistortFree(psPlaneDistort *pd)
{
    psFree(pd->x);
    psFree(pd->y);
}

static void psPlaneTransformFree(psPlaneTransform *pd)
{
    psFree(pd->x);
    psFree(pd->y);
}

#define PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(NAME) \
{\
    NAME = (psPlaneTransform *) psAlloc(sizeof(psPlaneTransform)); \
    psMemSetDeallocator(NAME, (psFreeFunc) psPlaneTransformFree); \
    NAME->x = psPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD); \
    NAME->y = psPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD); \
    NAME->x->coeff[1][0] = 1.0; \
    NAME->y->coeff[0][1] = 1.0; \
}

#define PS_CREATE_4D_IDENTITY_PLANE_DISTORT(NAME) \
{\
    NAME = (psPlaneDistort *) psAlloc(sizeof(psPlaneDistort)); \
    psMemSetDeallocator(NAME, (psFreeFunc) psPlaneDistortFree); \
    NAME->x = psPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD); \
    NAME->y = psPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD); \
    NAME->x->coeff[1][0][0][0] = 1.0; \
    NAME->y->coeff[0][1][0][0] = 1.0; \
}

/*
    This routine will create a system of FPAs/Chips/Cells/Readouts.  For
    simplicity, an FPA is defined as a linear array of chips, and a chip is
    defined as a linear array of cells, both in the x direction.  The transforms
    between the various layers take into account the cell/chip and the boundaries
    between each.
*/
psFPA *genSystem()
{
    psS32 i;
    psS32 j;
    psS32 k;
    psImage *tmpImage;
    //    psTime* now = psTimeGetNow(PS_TIME_UTC);
    //    psObservatory* obs = psObservatoryAlloc(NAME, LATITUDE, LONGITUDE, HEIGHT, LAPSERATE);
    //    psExposure* exp = psExposureAlloc(RA, DEC, HA, ZD, AZ,
    //                                      now, ROT_ANGLE, TEMPERATURE, PRESSURE, HUMIDITY,
    //                                      EXP_TIME, WAVELENGTH, obs);
    //    psGrommit *grom = psGrommitAlloc(exp);
    psFPA *myFPA = psFPAAlloc(NUM_CHIPS, NULL);
    psChip **chips;
    psCell **cells;
    psReadout **readouts;

    // We create the transforms between the FPA and tangent plane.  The
    // transform is a simple identity transform.
    PS_CREATE_4D_IDENTITY_PLANE_DISTORT(myFPA->fromTangentPlane);
    PS_CREATE_4D_IDENTITY_PLANE_DISTORT(myFPA->toTangentPlane);

    chips = (psChip **) myFPA->chips->data;
    myFPA->chips->n = NUM_CHIPS;

    for (i=0;i<NUM_CHIPS;i++) {
        chips[i] = psChipAlloc(NUM_CELLS, myFPA);
        // We create the transforms between the chip and FPA.  The
        // transform is a simple identity transform.
        PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(chips[i]->toFPA);
        chips[i]->toFPA->x->coeff[0][0] = (double) (i * (CHIP_WIDTH + CHIP_GAP));
        chips[i]->toFPA->x->coeff[1][1] = 0.0;
        PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(chips[i]->fromFPA);
        chips[i]->fromFPA->x->coeff[0][0] = (double) (-i * (CHIP_WIDTH + CHIP_GAP));
        chips[i]->fromFPA->x->coeff[1][1] = 0.0;

        cells = (psCell **) chips[i]->cells->data;
        chips[i]->cells->n = NUM_CELLS;

        for (j=0;j<NUM_CELLS;j++) {
            cells[j] = psCellAlloc(NUM_READOUTS, chips[i]);
            // We create the transforms between the cell and the chip.  The
            // transform is a simple identity transform.
            PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(cells[j]->toChip);
            cells[j]->toChip->x->coeff[0][0] = (double) (j * (CELL_WIDTH + CELL_GAP));
            cells[j]->toChip->x->coeff[1][1] = 0.0;
            PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(cells[j]->fromChip);
            cells[j]->fromChip->x->coeff[0][0] = (double) (-j * (CELL_WIDTH + CELL_GAP));
            cells[j]->fromChip->x->coeff[1][1] = 0.0;
            PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(cells[j]->toFPA);
            cells[j]->toFPA->x->coeff[0][0] = (double) ((i * (CHIP_WIDTH + CHIP_GAP)) +
                                              (j * (CELL_WIDTH + CELL_GAP)));
            cells[j]->toFPA->x->coeff[1][1] = 0.0;
            PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(cells[j]->toTP);
            cells[j]->toTP->x->coeff[0][0] = cells[j]->toFPA->x->coeff[0][0];
            cells[j]->toTP->x->coeff[1][1] = 0.0;

            readouts = (psReadout **) cells[j]->readouts->data;
            cells[j]->readouts->n = NUM_READOUTS;

            for (k=0;k<NUM_READOUTS;k++) {
                tmpImage = psImageAlloc(READOUT_NUM_COLS, READOUT_NUM_ROWS, PS_TYPE_F32);
                for (psS32 row=0;row<READOUT_NUM_ROWS;row++) {
                    for(psS32 col=0;col<READOUT_NUM_COLS;col++) {
                        tmpImage->data.F32[row][col] = (float) ((i * (CHIP_WIDTH + CHIP_GAP)) +
                                                                (j * (CELL_WIDTH + CELL_GAP)));
                    }
                }
                readouts[k] = psReadoutAlloc();
                readouts[k]->image = tmpImage;
                *(psS32*)&readouts[k]->col0 = (j * (CELL_WIDTH + CELL_GAP));
                *(psS32*)&readouts[k]->row0 = 0;


            }
        }
    }

    return(myFPA);
}

psS32 main( psS32 argc, char* argv[] )
{
    psLogSetLevel( PS_LOG_INFO );

    return ( ! runTestSuite( stderr, "psAstrometry", tests, argc, argv ) );
}

#define X_NUM_ROWS 4
#define X_NUM_COLS 5
#define Y_NUM_ROWS 6
#define Y_NUM_COLS 7
#define X_SCALE  1.0
#define Y_SCALE  2.0
#define X0  50.0
#define Y0  60.0

psS32 test1( void )
{
    psS32 i;
    psS32 j;
    psFixedPattern *tmp = NULL;
    float x0 = X0;
    float y0 = Y0;
    float xScale = X_SCALE;
    float yScale = Y_SCALE;
    psImage *x = psImageAlloc(X_NUM_COLS, X_NUM_ROWS, PS_TYPE_F64);
    psImage *y = psImageAlloc(Y_NUM_COLS, Y_NUM_ROWS, PS_TYPE_F64);
    psS32 testStatus = 0;
    psS32 numPixels = 0;

    for (i=0;i<X_NUM_ROWS;i++) {
        for (j=0;j<X_NUM_COLS;j++) {
            x->data.F64[i][j] = (double) (i + j);
        }
    }
    for (i=0;i<Y_NUM_ROWS;i++) {
        for (j=0;j<Y_NUM_COLS;j++) {
            y->data.F64[i][j] = (double) (i - j);
        }
    }

    tmp = psFixedPatternAlloc(x0, y0, xScale, yScale, x, y);

    // Test the x pattern
    for (i=0;i<X_NUM_ROWS;i++) {
        for (j=0;j<X_NUM_COLS;j++) {
            if ((tmp->x)[i][j] != x->data.F64[i][j]) {
                printf("ERROR: (tmp->x)[%d][%d] is %.2f, should be %.2f\n", i, j,
                       (tmp->x)[i][j], x->data.F64[i][j]);
                testStatus = 1;
                numPixels++;
            }
        }
    }
    if (numPixels != 0) {
        printf("ERROR: x: %d incorrect pixels\n", numPixels);
    }

    // Test the y pattern
    numPixels = 0;
    for (i=0;i<Y_NUM_ROWS;i++) {
        for (j=0;j<Y_NUM_COLS;j++) {
            if ((tmp->y)[i][j] != y->data.F64[i][j]) {
                printf("ERROR: (tmp->y)[%d][%d] is %.2f, should be %.2f\n", i, j,
                       (tmp->y)[i][j], y->data.F64[i][j]);

                testStatus = 1;
                numPixels++;
            }
        }
    }
    if (numPixels != 0) {
        printf("ERROR: y: %d incorrect pixels\n", numPixels);
    }

    // Test the x variables
    if (tmp->nX != (X_NUM_ROWS * X_NUM_COLS)) {
        printf("ERROR: tmp->nX is %d, should be %d\n", tmp->nX, (X_NUM_ROWS * X_NUM_COLS));
        testStatus = 1;
    }
    if (tmp->x0 != X0) {
        printf("ERROR: tmp->x0 is %.2f, should be %.2f\n", tmp->x0, X0);
        testStatus = 1;
    }
    if (tmp->xScale != X_SCALE) {
        printf("ERROR: tmp->xScale is %.2f, should be %.2f\n", tmp->xScale, X_SCALE);
        testStatus = 1;
    }

    // Test the y variables
    if (tmp->nY != (Y_NUM_ROWS * Y_NUM_COLS)) {
        printf("ERROR: tmp->nY is %d, should be %d\n", tmp->nY, (Y_NUM_ROWS * Y_NUM_COLS));
        testStatus = 1;
    }
    if (tmp->y0 != Y0) {
        printf("ERROR: tmp->y0 is %.2f, should be %.2f\n", tmp->y0, Y0);
        testStatus = 1;
    }
    if (tmp->yScale != Y_SCALE) {
        printf("ERROR: tmp->yScale is %.2f, should be %.2f\n", tmp->yScale, Y_SCALE);
        testStatus = 1;
    }

    psFree(tmp);
    psFree(x);
    psFree(y);

    return 0;
}

psS32 test2()
{
    psS32 testStatus = 0;

    psTime* now = psTimeGetNow(PS_TIME_TAI);

    psObservatory* obs = psObservatoryAlloc(NAME, LATITUDE, LONGITUDE, HEIGHT, LAPSERATE);

    psExposure* exp = psExposureAlloc(RA, DEC, HA, ZD, AZ,
                                      now, ROT_ANGLE, TEMPERATURE, PRESSURE, HUMIDITY,
                                      EXP_TIME, WAVELENGTH, obs);

    psGrommit *grom = psGrommitAlloc(exp);

    if (grom->height != HEIGHT) {
        printf("ERROR: grom->height is %.2f, should be %.2f\n", grom->height, HEIGHT);
        testStatus = 1;
    }
    if (grom->temperature != TEMPERATURE) {
        printf("ERROR: grom->temperature is %.2f, should be %.2f\n", grom->temperature, TEMPERATURE);
        testStatus = 2;
    }

    if (grom->pressure != PRESSURE) {
        printf("ERROR: grom->pressure is %.2f, should be %.2f\n", grom->pressure, PRESSURE);
        testStatus = 3;
    }

    if (grom->humidity != HUMIDITY) {
        printf("ERROR: grom->humidity is %.2f, should be %.2f\n", grom->humidity, HUMIDITY);
        testStatus = 4;
    }

    if (grom->wavelength != WAVELENGTH) {
        printf("ERROR: grom->wavelength is %.2f, should be %.2f\n", grom->wavelength, WAVELENGTH);
        testStatus = 5;
    }

    if (grom->lapseRate != LAPSERATE) {
        printf("ERROR: grom->lapseRate is %.2f, should be %.2f\n", grom->lapseRate, LAPSERATE);
        testStatus = 6;
    }

    /*
        XXX: Don't know how to verify: latitude, sinLat, cosLat, abberationMag,
             refractA, refractB, longitudeOffset, siderealTime

       if (grom-> != ) {
            printf("ERROR: grom-> is %.2f, should be %.2f\n", grom->, );
            testStatus = 1;
       }
    */

    printf("grom->latitude is %.2f\n", grom->latitude);
    printf("grom->longitude is %.2f\n", grom->longitude);
    printf("grom->height is %.2f\n", grom->height);
    printf("grom->abberationMag is %.2f\n", grom->abberationMag);
    printf("grom->temperature is %.2f\n", grom->temperature);
    printf("grom->pressure is %.2f\n", grom->pressure);
    printf("grom->humidity is %.2f\n", grom->humidity);
    printf("grom->wavelength is %.2f\n", grom->wavelength);
    printf("grom->lapseRate is %.2f\n", grom->lapseRate);
    printf("grom->refractA is %.2f\n", grom->refractA);
    printf("grom->refractB is %.2f\n", grom->refractB);

    // printf("grom->siderealTime is %.2f\n", grom->siderealTime);

    psFree(now);
    psFree(obs);
    psFree(exp);
    psFree(grom);
    psLibFinalize();

    return(testStatus);
}

psS32 printCell(psCell *cell)
{
    psReadout **readouts = (psReadout **) cell->readouts->data;
    psReadout *readout = readouts[0];

    printf("-------------------------------\n");
    for (psS32 i = 0; i < readout->image->numRows ; i++) {
        for (psS32 j = 0; j < readout->image->numCols ; j++) {
            printf("(%.1f) ", readout->image->data.F32[i][j]);
        }
        printf("\n");
    }
    return(0);
}


/******************************************************************************
This routine tests many Astrometry functions:
    psChipInFPA()
    psCellInFPA()
    psCoordFPAToChip()
    psCoordChipToCell()
This routine loops through all valid cell coordinates and maintains a running
set of (x,y) coordinates in the FPA plane.  It calls the first two functions
with each FPA coordinate and determines the chip/cell that that coordinate
corresponds to.  It then calls the last two functions to determine the actual
chip and cell coordinates.  The cell coordinate is verified for correctness.
 
It then calls
    psCoordCellToChip()
to determine if the chip coordinates that were generated agree with the chip
coords returned above.
 
Then, it calls
    psCellInChip()
and flags an ERROR if the cell it returns is different from that of
psCellInFPA().
 
Then, it calls
    psCoordChipToFPA
to determine if the FPA coordinates are correct.
 
 
 *****************************************************************************/
psS32 test3( void )
{
    psS32 x;
    psS32 y;
    psPlane fpaCoord;
    psFPA *myFPA = genSystem();
    psCell *myCell = NULL;
    psPlane chipCoord;
    psPlane cellCoord;
    psPlane testCoord;

    psS32 xReadout = 0;
    psS32 xFPA = 0;
    for (psS32 chip=0;chip<NUM_CHIPS;chip++) {
        for (psS32 cell=0;cell<NUM_CELLS;cell++) {
            for(x=0;x<CELL_WIDTH;x++) {
                for (y=0;y<CELL_HEIGHT;y++) {
                    fpaCoord.x = (double) xFPA;
                    fpaCoord.y = (double) y;
                    //printf("------------------ (%.2f, %.2f) ------------------\n", fpaCoord.x, fpaCoord.y);
                    psChip* tmpChip = psChipInFPA(&fpaCoord, myFPA);
                    myCell = psCellInFPA(&fpaCoord, myFPA);

                    if ((myCell == NULL) || (tmpChip == NULL)) {
                        if (myCell == NULL) {
                            printf("ERROR: myCell is NULL\n");
                        }
                        if (tmpChip == NULL) {
                            printf("ERROR: tmpChip is NULL\n");
                        }
                    } else {
                        psCoordFPAToChip(&chipCoord, &fpaCoord, tmpChip);
                        psCoordChipToCell(&cellCoord, &chipCoord, myCell);

                        if (x != (psS32) cellCoord.x) {
                            printf("ERROR: x coord was %.2f, should be %d\n", cellCoord.x, x);
                        }
                        if (y != (psS32) cellCoord.y) {
                            printf("ERROR: y coord was %.2f, should be %d\n", cellCoord.y, y);
                        }

                        psCoordCellToChip(&testCoord, &cellCoord, myCell);
                        if (testCoord.x != chipCoord.x) {
                            printf("ERROR: psCoordCellToChip() x coord was %.2f, should be %d\n", cellCoord.x, x);
                        }
                        if (testCoord.y != chipCoord.y) {
                            printf("ERROR: psCoordCellToChip() y coord was %.2f, should be %d\n", cellCoord.y, y);
                        }

                        psCell *myCell2 = psCellInChip(&chipCoord, tmpChip);
                        if (myCell2 != myCell) {
                            printf("ERROR: psCellInFPA() != psCellInChip(psChipInFPA()) (%p %p)\n", myCell2, myCell);
                            printCell(myCell2);
                            printCell(myCell);

                        }

                        psCoordChipToFPA(&testCoord, &chipCoord, tmpChip);
                        if (testCoord.x != xFPA) {
                            printf("ERROR: psCoordChipToFPA() x coord was %.2f, should be %d\n", cellCoord.x, x);
                        }
                        if (testCoord.y != y) {
                            printf("ERROR: psCoordChipToFPA() y coord was %.2f, should be %d\n", cellCoord.y, y);
                        }

                        psCoordFPAToTP(&testCoord, &fpaCoord, 0.0, 0.0, myFPA);
                        if (testCoord.x != fpaCoord.x) {
                            printf("ERROR: psCoordFPAToTP() x coord was %.2f, should be %d\n", cellCoord.x, x);
                        }
                        if (testCoord.y != fpaCoord.y) {
                            printf("ERROR: psCoordFPAToTP() y coord was %.2f, should be %d\n", cellCoord.y, y);
                        }

                        psCoordCellToFPA(&testCoord, &cellCoord, myCell);
                        if (testCoord.x != fpaCoord.x) {
                            printf("ERROR: psCoordFPAToTP() x coord was %.2f, should be %d\n", cellCoord.x, x);
                        }
                        if (testCoord.y != fpaCoord.y) {
                            printf("ERROR: psCoordFPAToTP() y coord was %.2f, should be %d\n", cellCoord.y, y);
                        }

                        psCoordTPToFPA(&testCoord, &fpaCoord, 0.0, 0.0, myFPA);
                        if (testCoord.x != fpaCoord.x) {
                            printf("ERROR: psCoordFPAToTP() x coord was %.2f, should be %d\n", cellCoord.x, x);
                        }
                        if (testCoord.y != fpaCoord.y) {
                            printf("ERROR: psCoordFPAToTP() y coord was %.2f, should be %d\n", cellCoord.y, y);
                        }
                    }
                }
                xFPA++;
                xReadout++;
            }
            // We skip the gaps between cells.
            for(x=0;x<CELL_GAP;x++) {
                xFPA++;
                xReadout=0;
            }
        }
        // We skip the gaps between chips.
        for (x=0;x<CHIP_GAP;x++) {
            xFPA++;
            xReadout=0;
        }
    }
    psFree(myFPA);

    return(0);
}

int printTrans(psPlaneTransform *trans)
{
    printf("-------------------- trans->x coeffs --------------------\n");
    for (psS32 i = 0 ; i < trans->x->nX ; i++) {
        for (psS32 j = 0 ; j < trans->x->nY ; j++) {
            printf("(%.1f) ", trans->x->coeff[i][j]);
        }
        printf("\n");
    }
    printf("-------------------- trans->y coeffs --------------------\n");
    for (psS32 i = 0 ; i < trans->y->nX ; i++) {
        for (psS32 j = 0 ; j < trans->y->nY ; j++) {
            printf("(%.1f) ", trans->y->coeff[i][j]);
        }
        printf("\n");
    }
    return(0);
}


/******************************************************************************
TEST04: The strategy behind test04() is to generate a pair of transforms t1
and t2 with a wide variety of sizes for the x/y components of the x/y
transforms, then set the coefficients of t1/t2 to arbitray values, then
generate a new transform t3 via psPlaneTransformCombine() function.  Then, for
several arbitrary input plane coordinates, we tarnsform them via the new
transform t3, as well as individually via t1 then t2, and verify that they
produce identical output coordinates.
 *****************************************************************************/
// These macros determine the number of x/y test points which will be x-formed.
#define TST04_X_MIN 0.0
#define TST04_X_MAX 10.0
#define TST04_Y_MIN 0.0
#define TST04_Y_MAX 10.0
#define TST04_NUM 5.0
// These macros define the max sizes of the various input transforms.
#define TST04_T1_X_X 6
#define TST04_T1_X_Y 3
#define TST04_T1_Y_X 3
#define TST04_T1_Y_Y 3
#define TST04_T2_X_X 3
#define TST04_T2_X_Y 3
#define TST04_T2_Y_X 3
#define TST04_T2_Y_Y 3
/******************************************************************************
tstTransforms(t1, t2, t3): this function generates a set of arbitrary x/y
coordinates, then transforms them into output coordinates via the t3
transform, as well as the t1 followed by the t2 transform.  It verifies the
output coordinates are identical and returns TRUE if so.  Otherwise, it prints
an error message and returns FALSE.
 *****************************************************************************/
bool tstTransforms(psPlaneTransform *t1, psPlaneTransform *t2, psPlaneTransform *t3)
{
    bool testStatus = true;
    psPlane *in = psPlaneAlloc();
    in->xErr = 0.0;
    in->yErr = 0.0;

    for (in->x = TST04_X_MIN ;
            in->x <= TST04_X_MAX ;
            in->x+= (TST04_X_MAX-TST04_X_MIN) / TST04_NUM) {
        for (in->y = TST04_Y_MIN ; in->y <= TST04_Y_MAX ; in->y+= (TST04_Y_MAX-TST04_Y_MIN) / TST04_NUM) {
            // Apply the t1/t2 transforms individually to the in coords.
            psPlane *mid = psPlaneTransformApply(NULL, t1, in);
            if (mid == NULL) {
                printf("TEST ERROR: intermediate psPlane coords are NULL.\n");
                psFree(in);
                return(false);
            }

            psPlane *outT1T2 = psPlaneTransformApply(NULL, t2, mid);
            if (outT1T2 == NULL) {
                printf("TEST ERROR: intermediate psPlane coords are NULL.\n");
                psFree(mid);
                return(false);
            }

            // Apply the t3 transforms individually to the in coords.
            psPlane *outT3 = psPlaneTransformApply(NULL, t3, in);
            if (outT3 == NULL) {
                printf("TEST ERROR: intermediate psPlane coords are NULL.\n");
                psFree(mid);
                psFree(outT1T2);
                return(false);
            }

            // Verify that the results are identical.
            if (fabs(outT3->x - outT1T2->x) > FLT_EPSILON) {
                printf("TEST ERROR: x is %f, should be %f\n", outT3->x, outT1T2->x);
                testStatus = false;
            }
            if (fabs(outT3->y - outT1T2->y) > FLT_EPSILON) {
                printf("TEST ERROR: y is %f, should be %f\n", outT3->y, outT1T2->y);
                testStatus = false;
            }
            psFree(mid);
            psFree(outT1T2);
            psFree(outT3);
        }
    }

    psFree(in);
    return(testStatus);
}

/******************************************************************************
setCoeffs(t1, t2): this function sets the coefficients of the t1 and t2
transforms to somewhat arbitrary values.
 *****************************************************************************/
int setCoeffs(psPlaneTransform *t1, psPlaneTransform *t2)
{
    for (psS32 t1xx = 0 ; t1xx < t1->x->nX ; t1xx++) {
        for (psS32 t1xy = 0 ; t1xy < t1->x->nY ; t1xy++) {
            t1->x->coeff[t1xx][t1xy] = (psF32) (t1xx + t1xy);
        }
    }

    for (psS32 t1yx = 0 ; t1yx < t1->y->nX ; t1yx++) {
        for (psS32 t1yy = 0 ; t1yy < t1->y->nY ; t1yy++) {
            t1->y->coeff[t1yx][t1yy] = (psF32) (t1yx + t1yy - 10);
        }
    }
    for (psS32 t2xx = 0 ; t2xx < t2->x->nX ; t2xx++) {
        for (psS32 t2xy = 0 ; t2xy < t2->x->nY ; t2xy++) {
            t2->x->coeff[t2xx][t2xy] = (psF32) (2 + (t2xx * t2xy));
        }
    }
    for (psS32 t2yx = 0 ; t2yx < t2->y->nX ; t2yx++) {
        for (psS32 t2yy = 0 ; t2yy < t2->y->nY ; t2yy++) {
            t2->y->coeff[t2yx][t2yy] = (psF32) (3 - t2yx + t2yy);
        }
    }

    return(0);
}

/******************************************************************************
test04(): We test psPlaneTransformCombine() with a variety of input
transforms, as well as input coordinates, as well as erroneous input
parameters.
 *****************************************************************************/
psS32 test4( void )
{
    bool testStatus = true;


    //
    // Create a variety of input input transforms, then test them.
    //
    for (psS32 t1xx = 1 ; t1xx < TST04_T1_X_X ; t1xx++) {
        for (psS32 t1xy = 1 ; t1xy < TST04_T1_X_Y ; t1xy++) {
            for (psS32 t1yx = 1 ; t1yx < TST04_T1_Y_X ; t1yx++) {
                for (psS32 t1yy = 1 ; t1yy < TST04_T1_Y_Y ; t1yy++) {
                    for (psS32 t2xx = 1 ; t2xx < TST04_T2_X_X ; t2xx++) {
                        for (psS32 t2xy = 1 ; t2xy < TST04_T2_X_Y ; t2xy++) {
                            for (psS32 t2yx = 1 ; t2yx < TST04_T2_Y_X ; t2yx++) {
                                for (psS32 t2yy = 1 ; t2yy < TST04_T2_Y_Y ; t2yy++) {
                                    psPlaneTransform *trans1 = psPlaneTransformAlloc(PS_MAX(t1xx, t1yx), PS_MAX(t1xy, t1yy));
                                    psPlaneTransform *trans2 = psPlaneTransformAlloc(PS_MAX(t2xx, t2yx), PS_MAX(t2xy, t2yy));
                                    psPlaneTransform *trans3 = NULL;
                                    setCoeffs(trans1, trans2);
                                    trans3 = psPlaneTransformCombine(NULL, trans1, trans2,psRegionSet(NAN,NAN,NAN,NAN),0);
                                    //XXX: the last two parameters are bogus.  Needs to change.  -rdd

                                    if (trans3 == NULL) {
                                        printf("TEST ERROR: psPlaneTransformCombine() returned NULL/\n");
                                        testStatus = false;
                                    } else {
                                        testStatus = tstTransforms(trans1, trans2, trans3);
                                    }
                                    psFree(trans1);
                                    psFree(trans2);
                                    psFree(trans3);
                                }
                            }
                        }
                    }
                }
            }
        }
    }


    //
    // Test erroneous input parameter combinations
    //
    psPlaneTransform *trans1 = psPlaneTransformAlloc(2, 2);
    psPlaneTransform *trans2 = psPlaneTransformAlloc(2, 2);
    psPlaneTransform *trans3 = NULL;

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformCombine with NULL trans1.  Should generate error and return NULL.\n");
    trans3 = psPlaneTransformCombine(NULL, NULL, trans2, psRegionSet(0,0,0,0), 10);
    if (trans3 != NULL) {
        printf("TEST ERROR: psPlaneTransformCombine() returned a non-NULL psPlaneTransform.\n");
        testStatus = false;
        psFree(trans3);
    }

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformCombine with NULL trans2.  Should generate error and return NULL.\n");
    trans3 = psPlaneTransformCombine(NULL, trans1, NULL, psRegionSet(0,0,0,0), 10);
    if (trans3 != NULL) {
        printf("TEST ERROR: psPlaneTransformCombine() returned a non-NULL psPlaneTransform.\n");
        testStatus = false;
        psFree(trans3);
    }

    psFree(trans1);
    psFree(trans2);
    return(!testStatus);
}


#define TST05_X_MIN 0.0
#define TST05_X_MAX 10.0
#define TST05_Y_MIN 0.0
#define TST05_Y_MAX 10.0
#define TST05_NUM_X 5.0
#define TST05_NUM_Y 5.0
#define TST05_X_POLY_ORDER 2
#define TST05_Y_POLY_ORDER 2
#define TST05_NUM_DATA (TST05_NUM_X * TST05_NUM_Y)
/******************************************************************************
XXX: This is only a rudimentary test of the psPlaneTransformFit() function.
It tests a few NULL input parameter conditions, and some simple linear
transformations.
 *****************************************************************************/
psS32 test5( void )
{
    bool testStatus = true;
    psPlaneTransform *trans = psPlaneTransformAlloc(TST05_X_POLY_ORDER, TST05_Y_POLY_ORDER);
    psPlaneTransform *transInit = psPlaneTransformAlloc(TST05_X_POLY_ORDER, TST05_Y_POLY_ORDER);
    psArray *src = psArrayAlloc((int) TST05_NUM_DATA);
    psArray *dst = psArrayAlloc((int) TST05_NUM_DATA);
    bool rc;

    //
    // We set an arbitrary non-linear transformation.
    //
    for (psS32 x = 0 ; x < TST05_X_POLY_ORDER ; x++) {
        for (psS32 y = 0 ; y < TST05_Y_POLY_ORDER ; y++) {
            transInit->x->coeff[x][y] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
            transInit->y->coeff[x][y] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
        }
    }

    // okay, for now, it's linear.
    transInit->x->coeff[1][1] = 0.1;
    transInit->y->coeff[1][1] = 0.1;

    /* HEY
        //
        // We set an arbitrary linear transformation.
        //
        transInit->x->coeff[0][0] = 1.0;
        transInit->x->coeff[0][1] = 3.0;
        transInit->x->coeff[1][0] = 2.0;
        transInit->x->coeff[1][1] = 7.0;
        transInit->y->coeff[0][0] = 4.0;
        transInit->y->coeff[0][1] = 6.0;
        transInit->y->coeff[1][0] = 5.0;
        transInit->y->coeff[1][1] = 3.0;
    */
    //
    // We generate a grid of input data points in the x1,y1 plane, calculate the
    // corresponding values in the transformed plane, and set these to
    // the src and dst psVectors.
    //
    psS32 i = 0;
    for (psF32 x = TST05_X_MIN ; x < TST05_X_MAX ; x+= (TST05_X_MAX-TST05_X_MIN)/TST05_NUM_X) {
        for (psF32 y = TST05_Y_MIN ; y < TST05_Y_MAX ; y+= (TST05_Y_MAX-TST05_Y_MIN)/TST05_NUM_Y) {
            src->data[i] = (psPtr *) psPlaneAlloc();
            ((psPlane *) src->data[i])->x = x;
            ((psPlane *) src->data[i])->y = y;

            dst->data[i] = psPlaneTransformApply(NULL, transInit, ((psPlane *) src->data[i]));
            i++;
        }
    }

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformFit with NULL trans.  Should generate error and return NULL.\n");
    rc = psPlaneTransformFit(NULL, src, dst, 100, 100.0);
    if (rc == true) {
        printf("TEST ERROR: psPlaneTransformFit() returned TRUE.\n");
        testStatus = false;
    }

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformFit with NULL src psArray.  Should generate error and return NULL.\n");
    rc = psPlaneTransformFit(trans, NULL, dst, 100, 100.0);
    if (rc == true) {
        printf("TEST ERROR: psPlaneTransformFit() returned TRUE.\n");
        testStatus = false;
    }

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformFit with NULL dst psArray.  Should generate error and return NULL.\n");
    rc = psPlaneTransformFit(trans, src, NULL, 100, 100.0);
    if (rc == true) {
        printf("TEST ERROR: psPlaneTransformFit() returned TRUE.\n");
        testStatus = false;
    }

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformFit with acceptable data.\n");
    rc = psPlaneTransformFit(trans, src, dst, 100, 100.0);
    if (rc != true) {
        printf("TEST ERROR: psPlaneTransformFit() returned FALSE.\n");
        testStatus = false;
    }

    //
    // For the initial grid of input points, we transform them to output points with
    // the derived transformation, and verify that they are within 10%.
    //

    // XXX: We ignore endpoints.
    for (psS32 i = 1 ; i < src->n-1 ; i++) {
        psPlane *inData = (psPlane *) src->data[i];
        psPlane *outData = (psPlane *) dst->data[i];
        psPlane *outDataDeriv = psPlaneTransformApply(NULL, trans, inData);
        // HERE
        if (!PS_PERCENT_COMPARE(outDataDeriv->x, outData->x, 0.20) ||
                !PS_PERCENT_COMPARE(outDataDeriv->y, outData->y, 0.20)) {
            printf("TEST ERROR: the derived output coords (%d) were (%f, %f) should have been (%f, %f).\n",
                   i, outDataDeriv->x, outDataDeriv->y, outData->x, outData->y);
            testStatus = false;
        }
        psFree(outDataDeriv);
    }

    psFree(transInit);
    psFree(trans);
    psFree(src);
    psFree(dst);

    return(!testStatus);
}

/******************************************************************************
XXX: This is only a rudimentary test of the psPlaneTransformInvert()
function.  It tests a few NULL input parameter conditions.
 
XXX: We call it with acceptable data, but we do not verify the output.
 *****************************************************************************/
psS32 test6( void )
{
    bool testStatus = true;
    psPlaneTransform *trans = psPlaneTransformAlloc(2, 2);
    psPlaneTransform *transOut = NULL;
    psRegion myRegion = psRegionSet(1.0, 5.0, 1.0, 5.0);

    //
    // Set non-linear transformation
    //
    trans->x->coeff[0][0] = 1.0;
    trans->x->coeff[0][1] = 3.0;
    trans->x->coeff[1][0] = 2.0;
    trans->x->coeff[1][1] = 1.0;
    trans->y->coeff[0][0] = 4.0;
    trans->y->coeff[0][1] = 6.0;
    trans->y->coeff[1][0] = 5.0;
    trans->y->coeff[1][1] = 1.0;

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformInvert with NULL trans.  Should generate error and return NULL.\n");
    transOut = psPlaneTransformInvert(NULL, NULL, myRegion, 10);
    if (transOut != NULL) {
        printf("TEST ERROR: psPlaneTransformInvert() returned a non-NULL psPlaneTransform.\n");
        testStatus = false;
        psFree(transOut);
    }

    printf("----------------------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformInvert with zero nSamples.  Should generate error and return NULL.\n");
    transOut = psPlaneTransformInvert(NULL, trans, myRegion, 0);
    if (transOut != NULL) {
        printf("TEST ERROR: psPlaneTransformInvert() returned a non-NULL psPlaneTransform.\n");
        testStatus = false;
        psFree(transOut);
    }

    printf("Calling psPlaneTransformInvert with an acceptable linear transformation.\n");
    for (psS32 n = 0 ; n < 100 ; n++) {
        //
        // We create a rand() linear transformation and hope it is invertible.
        //
        trans->x->coeff[0][0] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
        trans->x->coeff[0][1] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
        trans->x->coeff[1][0] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
        trans->x->coeff[1][1] = 0.0;
        trans->y->coeff[0][0] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
        trans->y->coeff[0][1] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
        trans->y->coeff[1][0] = (psF32)rand() / (psF32)RAND_MAX * 10.0f;
        trans->y->coeff[1][1] = 0.0;

        transOut = psPlaneTransformInvert(NULL, trans, myRegion, 10);
        if (transOut == NULL) {
            printf("TEST ERROR: psPlaneTransformInvert() returned a NULL psPlaneTransform.\n");
            testStatus = false;
        }
        psPlane *in = psPlaneAlloc();
        in->x = 2.0;
        in->y = 3.0;

        psPlane *mid = psPlaneTransformApply(NULL, trans, in);
        psPlane *out = psPlaneTransformApply(NULL, transOut, mid);

        if (((fabs(in->x - out->x) > FLT_EPSILON)) ||
                ((fabs(in->y - out->y) > FLT_EPSILON))) {
            printf("TEST ERROR: in coords were (%f, %f), out coords were (%f, %f).\n",
                   in->x, in->y, out->x, out->y);
            testStatus = false;
        }

        psFree(in);
        psFree(mid);
        psFree(out);
        psFree(transOut);
    }

    psFree(trans);

    fflush(stdout);

    return(!testStatus);
}
