/** @file psImageStats.c
 *  \brief Routines for calculating statistics on images.
 *  @ingroup ImageStats
 *
 *  This file will hold the prototypes for procedures which calculate
 *  statistic on images, histograms on images, and fit/evaluate Chebyshev
 *  polynomials to images.
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-05-11 22:02:16 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stdarg.h>
#include <float.h>
#include <math.h>
#include "psMemory.h"
#include "psVector.h"
#include "psTrace.h"
#include "psError.h"
#include "psStats.h"
#include "psImage.h"
#include "psFunctions.h"
#include "psImageStats.h"
#include "psConstants.h"
#include "psImageErrors.h"

/// This routine must determine the various statistics for the image.
/*****************************************************************************
psImageStats(stats, in, mask, maskVal): this routine simply calls the
psVectorStats() routine, which does the actual statistical calculation.  In
order to do so, we create dummy psVectors and set their "data" pointer to that
of the input psImages.
 
XXX: use static psVectors
 
XXX: optimize this.  2k vs 4k, sample mean, takes8 seconds on Gene's machine.
Should take .2.
 *****************************************************************************/
psStats* psImageStats(psStats* stats,
                      const psImage* in,
                      const psImage* mask,
                      psS32 maskVal)
{
    psVector *junkData = NULL;
    psVector *junkMask = NULL;

    PS_PTR_CHECK_NULL(stats, NULL);
    PS_INT_CHECK_ZERO(stats->options, NULL);
    PS_IMAGE_CHECK_NULL(in, NULL)
    if (mask != NULL) {
        PS_IMAGE_CHECK_TYPE(mask, PS_TYPE_U8, NULL);
        PS_IMAGE_CHECK_SIZE_EQUAL(in, mask, NULL);
    }

    if (in->parent == NULL) {
        // stuff the image data into a psVector struct.
        junkData = (psVector *) psAlloc(sizeof(psVector));
        junkData->type = in->type;
        *(int*)&junkData->nalloc = in->numRows * in->numCols;
        junkData->n = junkData->nalloc;
        junkData->data.U8 = in->data.V[0];      // since psImage data is contiguous...
    } else {
        // image not necessarily contiguous
        int numRows = in->numRows;
        int numCols = in->numCols;
        int rowSize = numCols * (PSELEMTYPE_SIZEOF(in->type.type));

        junkData = psVectorAlloc(numRows*numCols, in->type.type);
        junkData->n = junkData->nalloc;

        psU8* data = junkData->data.U8;
        for (int row = 0; row < numRows; row++) {
            memcpy(data, in->data.V[row], rowSize);
            data += rowSize;
        }
    }

    if (mask != NULL) {
        if (mask->parent == NULL) {
            // stuff the mask data into a psVector struct.
            junkMask = psAlloc(sizeof(psVector));
            junkMask->type = mask->type;
            *(int*)&junkMask->nalloc = mask->numRows * mask->numCols;
            junkMask->n = junkMask->nalloc;
            junkMask->data.U8 = mask->data.V[0];
        } else {
            // image not necessarily contiguous
            int numRows = mask->numRows;
            int numCols = mask->numCols;
            int rowSize = numCols * (PSELEMTYPE_SIZEOF(mask->type.type));

            junkMask = psVectorAlloc(numRows*numCols, mask->type.type);
            junkMask->n = junkMask->nalloc;

            psU8* data = junkMask->data.U8;
            for (int row = 0; row < numRows; row++) {
                memcpy(data, mask->data.V[row], rowSize);
                data += rowSize;
            }
        }
    }

    stats = psVectorStats(stats, junkData, NULL, junkMask, maskVal);

    psFree(junkMask);
    psFree(junkData);
    return (stats);
}

/*****************************************************************************
NOTE: We assume that the psHistogram structure out has already been allocated
and initialized.
 *****************************************************************************/
psHistogram* psImageHistogram(psHistogram* out,
                              const psImage* in,
                              const psImage* mask,
                              psU32 maskVal)
{
    PS_PTR_CHECK_NULL(out, NULL);
    PS_PTR_CHECK_NULL(in, NULL);
    if (mask != NULL) {
        PS_IMAGE_CHECK_TYPE(mask, PS_TYPE_U8, NULL);
        PS_IMAGE_CHECK_SIZE_EQUAL(in, mask, NULL);
    }
    psVector* junkData = NULL;
    psVector* junkMask = NULL;

    if (in->parent == NULL) {
        // stuff the image data into a psVector struct.
        junkData = (psVector *) psAlloc(sizeof(psVector));
        junkData->type = in->type;
        *(int*)&junkData->nalloc = in->numRows * in->numCols;
        junkData->n = junkData->nalloc;
        junkData->data.U8 = in->data.V[0];      // since psImage data is contiguous...
    } else {
        // image not necessarily contiguous
        int numRows = in->numRows;
        int numCols = in->numCols;
        int rowSize = numCols * (PSELEMTYPE_SIZEOF(in->type.type));

        junkData = psVectorAlloc(numRows*numCols, in->type.type);
        junkData->n = junkData->nalloc;

        psU8* data = junkData->data.U8;
        for (int row = 0; row < numRows; row++) {
            memcpy(data, in->data.V[row], rowSize);
            data += rowSize;
        }
    }

    if (mask != NULL) {
        if (mask->parent == NULL) {
            // stuff the mask data into a psVector struct.
            junkMask = psAlloc(sizeof(psVector));
            junkMask->type = mask->type;
            *(int*)&junkMask->nalloc = mask->numRows * mask->numCols;
            junkMask->n = junkMask->nalloc;
            junkMask->data.U8 = mask->data.V[0];
        } else {
            // image not necessarily contiguous
            int numRows = mask->numRows;
            int numCols = mask->numCols;
            int rowSize = numCols * (PSELEMTYPE_SIZEOF(mask->type.type));

            junkMask = psVectorAlloc(numRows*numCols, mask->type.type);
            junkMask->n = junkMask->nalloc;

            psU8* data = junkMask->data.U8;
            for (int row = 0; row < numRows; row++) {
                memcpy(data, mask->data.V[row], rowSize);
                data += rowSize;
            }
        }
    }

    out = psVectorHistogram(out, junkData, NULL, junkMask, maskVal);

    psFree(junkMask);
    psFree(junkData);

    return (out);
}

/*****************************************************************************
calcScaleFactorsEval(n): The Chebyshev polynomials are defined over the
interval [-1.0 : 1.0].  Images typically have sizes of 512x512 or more.  In
order to use Chebyshev polynomials, we must scale the coordinates from
0:512 to -1:1.  This routine takes as input an integer N and produces as
output a vector of evenly spaced floating point values between -1.0:1.0.
 
XXX: Use the p_psNormalizeVector here?
 *****************************************************************************/
double* calcScaleFactors(psS32 n)
{
    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
    psS32 i = 0;
    double tmp = 0.0;
    double *scalingFactors = (double *)psAlloc(n * sizeof(double));

    for (i = 0; i < n; i++) {
        tmp = (double)(n - i);
        tmp = (M_PI * (tmp - 0.5)) / ((double)n);
        scalingFactors[i] = cos(tmp);
    }

    return (scalingFactors);
}

// XXX: Use a static array of Chebyshev polynomials.
psPolynomial1D **p_psCreateChebyshevPolys(psS32 maxChebyPoly)
{
    PS_INT_CHECK_POSITIVE(maxChebyPoly, NULL);
    psPolynomial1D **chebPolys = NULL;
    psS32 i = 0;
    psS32 j = 0;

    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
    for (i = 0; i < maxChebyPoly; i++) {
        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
    }

    // Create the Chebyshev polynomials.
    // Polynomial i has i-th order.
    chebPolys[0]->coeff[0] = 1;
    chebPolys[1]->coeff[1] = 1;
    for (i = 2; i < maxChebyPoly; i++) {
        for (j = 0; j < chebPolys[i - 1]->n; j++) {
            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
        }
        for (j = 0; j < chebPolys[i - 2]->n; j++) {
            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
        }
    }

    return (chebPolys);
}

/*****************************************************************************
psImageFitPolynomial(): This routine takes as input a 2-D image and produces
as output the coefficients of the Chebyshev polynomials which match that
input image.
  Input:
  Output:
  Internal Data Structures:
    chebPolys[i][j] 
    sums[i][j]: This will contain the sum of 
                input->data.F32[x][y] *
                psPolynomial1DEval(
chebPolys[i],
(float) x) *
                psPolynomial1DEval(
chebPolys[j],
(float) y, 
);
        over all pixels (x,y) in the image.
  *****************************************************************************/
psPolynomial2D* psImageFitPolynomial(psPolynomial2D* coeffs,
                                     const psImage* input)
{
    PS_IMAGE_CHECK_NULL(input, NULL);
    PS_IMAGE_CHECK_EMPTY(input, NULL);
    if ((input->type.type != PS_TYPE_S8) &&
            (input->type.type != PS_TYPE_U16) &&
            (input->type.type != PS_TYPE_F32) &&
            (input->type.type != PS_TYPE_F64)) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n");
    }
    PS_POLY_CHECK_NULL(coeffs, NULL);
    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
    psS32 x = 0;
    psS32 y = 0;
    psS32 i = 0;
    psS32 j = 0;
    double **sums = NULL;
    psPolynomial1D* *chebPolys = NULL;
    psS32 maxChebyPoly = 0;
    double *cScalingFactors = NULL;
    double *rScalingFactors = NULL;

    // Create the sums[][] data structure.  This
    // will hold the LHS of
    // equation
    // 29 in the ADD: sums[k][l] = SUM {
    // image(x,y) * Tk(x) * Tl(y) }
    sums = (double **)psAlloc(coeffs->nX * sizeof(double *));
    for (i = 0; i < coeffs->nX; i++) {
        sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double));
    }
    // We scale the pixel positions to values
    // between -1.0 and 1.0
    rScalingFactors = calcScaleFactors(input->numRows);
    cScalingFactors = calcScaleFactors(input->numCols);

    // Determine how many Chebyshev polynomials
    // are needed, then create them.
    maxChebyPoly = coeffs->nX;
    if (coeffs->nY > coeffs->nX) {
        maxChebyPoly = coeffs->nY;
    }
    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);

    // Compute the sums[][] data structure.
    for (i = 0; i < coeffs->nX; i++) {
        for (j = 0; j < coeffs->nY; j++) {
            sums[i][j] = 0.0;
            for (x = 0; x < input->numRows; x++) {
                for (y = 0; y < input->numCols; y++) {
                    double pixel = 0.0;
                    if (input->type.type == PS_TYPE_S8) {
                        pixel = (double) input->data.S8[x][y];
                    } else if (input->type.type == PS_TYPE_U16) {
                        pixel = (double) input->data.U16[x][y];
                    } else if (input->type.type == PS_TYPE_F32) {
                        pixel = (double) input->data.F32[x][y];
                    } else if (input->type.type == PS_TYPE_F64) {
                        pixel = input->data.F64[x][y];
                    }
                    sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) *
                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
                }
            }
        }
    }

    for (i = 0; i < coeffs->nX; i++) {
        for (j = 0; j < coeffs->nY; j++) {
            coeffs->coeff[i][j] = sums[i][j];
            coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols);

            if ((i != 0) && (j != 0)) {
                coeffs->coeff[i][j] *= 4.0;
            } else if ((i == 0) && (j == 0)) {
                coeffs->coeff[i][j] *= 1.0;
            } else {
                coeffs->coeff[i][j] *= 2.0;
            }
        }
    }

    // Free the Chebyshev polynomials that were
    // created in this routine.
    for (i = 0; i < maxChebyPoly; i++) {
        psFree(chebPolys[i]);
    }
    psFree(chebPolys);

    // Free some data
    for (i = 0; i < coeffs->nX; i++) {
        psFree(sums[i]);
    }
    psFree(sums);
    psFree(cScalingFactors);
    psFree(rScalingFactors);

    return (coeffs);
}




/*****************************************************************************
psImageFitPolynomial(): This routine takes as input a 2-D image and produces
as output the coefficients of the Chebyshev polynomials which match that input
image.  This is a TEST version of the code.  It is not used by anything.
  Input:
  Output:
  Internal Data Structures:
    chebPolys[i][j] 
    sums[i][j]: This will contain the sum of 
                input->data.F32[x][y] *
                psPolynomial1DEval(
chebPolys[i],
(float) x) *
                psPolynomial1DEval(
chebPolys[j],
(float) y, 
);
        over all pixels (x,y) in the image.
  *****************************************************************************/
psPolynomial2D* psImageFitPolynomialTest(psPolynomial2D* coeffs,
        const psImage* input)
{
    PS_IMAGE_CHECK_NULL(input, NULL);
    PS_IMAGE_CHECK_EMPTY(input, NULL);
    if ((input->type.type != PS_TYPE_S8) &&
            (input->type.type != PS_TYPE_U16) &&
            (input->type.type != PS_TYPE_F32) &&
            (input->type.type != PS_TYPE_F64)) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n");
    }
    PS_POLY_CHECK_NULL(coeffs, NULL);
    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
    psS32 x = 0;
    psS32 y = 0;
    psS32 i = 0;
    psS32 j = 0;
    double **sums = NULL;
    psPolynomial1D* *chebPolys = NULL;
    psS32 maxChebyPoly = 0;
    double *cScalingFactors = NULL;
    double *rScalingFactors = NULL;
    psImage *nodes = psImageAlloc(input->numCols, input->numRows, PS_TYPE_F64);

    double min = -1.0;
    double max = 1.0;
    double bma = 0.5 * (max-min);  // 1
    double bpa = 0.5 * (max+min);  // 0
    // We must calculate the value of the image at the nodes where the
    // Chebyshev polynomials are 0.
    for (x = 0; x < input->numRows; x++) {
        double xTmp = cos(M_PI * (0.5 + ((float) x)) / ((float) input->numRows));
        double xNode = - ((xTmp + bma + bpa) - 1.0);
        double xOrig = ((float) input->numRows) * (xNode - min) / (max - min);

        for (y = 0; y < input->numCols; y++) {
            double yTmp = cos(M_PI * (0.5 + ((float) y)) / ((float) input->numCols));
            double yNode = - ((yTmp + bma + bpa) - 1.0);
            double yOrig = ((float) input->numCols) * (yNode - min) / (max - min);

            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yNode, xNode, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yTmp, xTmp, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yOrig, xOrig, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
        }
    }

    // Create the sums[][] data structure.  This
    // will hold the LHS of
    // equation
    // 29 in the ADD: sums[k][l] = SUM {
    // image(x,y) * Tk(x) * Tl(y) }
    sums = (double **)psAlloc(coeffs->nX * sizeof(double *));
    for (i = 0; i < coeffs->nX; i++) {
        sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double));
    }
    // We scale the pixel positions to values
    // between -1.0 and 1.0
    rScalingFactors = calcScaleFactors(input->numRows);
    cScalingFactors = calcScaleFactors(input->numCols);

    // Determine how many Chebyshev polynomials
    // are needed, then create them.
    maxChebyPoly = coeffs->nX;
    if (coeffs->nY > coeffs->nX) {
        maxChebyPoly = coeffs->nY;
    }
    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);

    // Compute the sums[][] data structure.
    for (i = 0; i < coeffs->nX; i++) {
        for (j = 0; j < coeffs->nY; j++) {
            sums[i][j] = 0.0;
            for (x = 0; x < input->numRows; x++) {
                for (y = 0; y < input->numCols; y++) {
                    double pixel;
                    /*
                                        if (input->type.type == PS_TYPE_S8) {
                                            pixel = (double) input->data.S8[x][y];
                                        } else if (input->type.type == PS_TYPE_U16) {
                                            pixel = (double) input->data.U16[x][y];
                                        } else if (input->type.type == PS_TYPE_F32) {
                                            pixel = (double) input->data.F32[x][y];
                                        } else if (input->type.type == PS_TYPE_F64) {
                                            pixel = input->data.F64[x][y];
                                        }
                    */
                    pixel = nodes->data.F64[x][y];
                    sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) *
                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
                }
            }
        }
    }

    for (i = 0; i < coeffs->nX; i++) {
        for (j = 0; j < coeffs->nY; j++) {
            coeffs->coeff[i][j] = sums[i][j];
            coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols);

            if ((i != 0) && (j != 0)) {
                coeffs->coeff[i][j] *= 4.0;
            } else if ((i == 0) && (j == 0)) {
                coeffs->coeff[i][j] *= 1.0;
            } else {
                coeffs->coeff[i][j] *= 2.0;
            }
        }
    }

    // Free the Chebyshev polynomials that were
    // created in this routine.
    for (i = 0; i < maxChebyPoly; i++) {
        psFree(chebPolys[i]);
    }
    psFree(chebPolys);

    // Free some data
    for (i = 0; i < coeffs->nX; i++) {
        psFree(sums[i]);
    }
    psFree(sums);
    psFree(cScalingFactors);
    psFree(rScalingFactors);
    psFree(nodes);

    return (coeffs);
}

/*****************************************************************************
XXX: Use static variables for Chebyshev polynomials and scaling factors. 
 *****************************************************************************/
psImage* p_psImageEvalPolynomialCheb(psImage* input,
                                     const psPolynomial2D* coeffs)
{
    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);

    psS32 x = 0;
    psS32 y = 0;
    psS32 i = 0;
    psS32 j = 0;
    psPolynomial1D* *chebPolys = NULL;
    psS32 maxChebyPoly = 0;
    double *cScalingFactors = NULL;
    double *rScalingFactors = NULL;
    float polySum = 0.0;

    // We scale the pixel positions to values between -1.0 and 1.0
    // Use static data structures here.
    rScalingFactors = calcScaleFactors(input->numRows);
    cScalingFactors = calcScaleFactors(input->numCols);

    // Determine how many Chebyshev polynomials
    // are needed, then create them.
    maxChebyPoly = coeffs->nX;
    if (coeffs->nY > coeffs->nX) {
        maxChebyPoly = coeffs->nY;
    }

    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);

    for (x = 0; x < input->numRows; x++) {
        for (y = 0; y < input->numCols; y++) {
            polySum = 0.0;
            for (i = 0; i < coeffs->nX; i++) {
                for (j = 0; j < coeffs->nY; j++) {
                    polySum +=
                        psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) *
                        psPolynomial1DEval(chebPolys[j], cScalingFactors[y]) *
                        coeffs->coeff[i][j];
                }
            }

            if (input->type.type == PS_TYPE_S8) {
                input->data.S8[x][y] = (char) polySum;
            } else if (input->type.type == PS_TYPE_U16) {
                input->data.U16[x][y] = (short int) polySum;
            } else if (input->type.type == PS_TYPE_F32) {
                input->data.F32[x][y] = (float) polySum;
            } else if (input->type.type == PS_TYPE_F64) {
                input->data.F64[x][y] = polySum;
            }
        }
    }

    // Free the Chebyshev polynomials that were
    // created in this routine.
    // XXX: Use static data structures here.
    for (i = 0; i < maxChebyPoly; i++) {
        psFree(chebPolys[i]);
    }
    psFree(chebPolys);

    psFree(cScalingFactors);
    psFree(rScalingFactors);

    return input;
}

psImage* p_psImageEvalPolynomialOrd(psImage* input,
                                    const psPolynomial2D* coeffs)
{
    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_ORD, NULL);

    for (int row = 0; row < input->numRows ; row++) {
        for (int col = 0; col < input->numCols ; col++) {
            if (input->type.type == PS_TYPE_S8) {
                input->data.S8[row][col] = (psS8) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col);
            } else if (input->type.type == PS_TYPE_U16) {
                input->data.U16[row][col] = (psS16) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col);
            } else if (input->type.type == PS_TYPE_F32) {
                input->data.F32[row][col] = psPolynomial2DEval(coeffs, (psF32) row, (psF32) col);
            } else if (input->type.type == PS_TYPE_F64) {
                input->data.F64[row][col] = (psF64) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col);
            }
        }
    }

    return(input);
}


/*****************************************************************************
XXX: I added normal polynomials to this routine.  Let IfA know, put it in the
psLib SDR.
 *****************************************************************************/
psImage* psImageEvalPolynomial(psImage* input,
                               const psPolynomial2D* coeffs)
{
    PS_IMAGE_CHECK_NULL(input, NULL);
    PS_IMAGE_CHECK_EMPTY(input, NULL);
    if ((input->type.type != PS_TYPE_S8) &&
            (input->type.type != PS_TYPE_U16) &&
            (input->type.type != PS_TYPE_F32) &&
            (input->type.type != PS_TYPE_F64)) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                "Unallowable image type.\n");
    }
    PS_POLY_CHECK_NULL(coeffs, NULL);

    if (coeffs->type == PS_POLYNOMIAL_ORD) {
        return(p_psImageEvalPolynomialOrd(input, coeffs));
    } else if (coeffs->type == PS_POLYNOMIAL_CHEB) {
        return(p_psImageEvalPolynomialCheb(input, coeffs));
    }
    printf("XXX: Error: wrong polynomial type\n");
    // XXX: psError()
    return(NULL);
}

