/** @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.43 $ $Name: not supported by cvs2svn $
*  @date $Date: 2004-10-14 01:06:44 $
*
*  Copyright 2004 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.
/*****************************************************************************
XXX: verify that image/mask have the correct types, and sizes.
 *****************************************************************************/
psStats* psImageStats(psStats* stats,
                      psImage* in,
                      psImage* mask,
                      int maskVal)
{
    psVector* junkData = NULL;
    psVector* junkMask = NULL;

    if (stats == NULL) {
        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageStats",
                   PS_ERR_BAD_PARAMETER_NULL, true,
                   PS_ERRORTEXT_psImage_STAT_NULL);
        return NULL;
    }

    if (in == NULL) {
        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageStats",
                   PS_ERR_BAD_PARAMETER_NULL, true,
                   PS_ERRORTEXT_psImage_IMAGE_NULL);
        return NULL;
    }

    if (stats->options == 0) {
        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageStats",
                   PS_ERR_BAD_PARAMETER_VALUE, true,
                   PS_ERRORTEXT_psImage_NO_STAT_OPTIONS);
        return stats;
    }
    // stuff the image data into a psVector struct.
    junkData = psAlloc(sizeof(psVector));
    junkData->type = in->type;
    junkData->nalloc = in->numRows * in->numCols;
    junkData->n = junkData->nalloc;
    junkData->data.V = in->data.V[0];      // since psImage data is contiguous...

    if (mask != NULL) {
        if (mask->type.type != PS_TYPE_MASK) {
            char* typeStr;
            PS_TYPE_NAME(typeStr,mask->type.type);
            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageStats",
                       PS_ERR_BAD_PARAMETER_TYPE, true,
                       PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE,
                       typeStr, PS_TYPE_MASK_NAME);
            psFree(junkData);
            return NULL;
        }
        // stuff the mask data into a psVector struct.
        junkMask = psAlloc(sizeof(psVector));
        junkMask->type = mask->type;
        junkMask->nalloc = mask->numRows * mask->numCols;
        junkMask->n = junkMask->nalloc;
        junkMask->data.V = mask->data.V[0];
    }

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

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

/*****************************************************************************
    NOTE: We assume that the psHistogram structure out has already been
    allocated and initialized.
    NOTE: verify that image/mask have the, correct types and  sizes.
 *****************************************************************************/
psHistogram* psImageHistogram(psHistogram* out,
                              psImage* in,
                              psImage* mask,
                              unsigned int maskVal)
{
    psVector* junkData = NULL;
    psVector* junkMask = NULL;

    // NOTE: Verify this action.
    PS_CHECK_NULL_PTR_RETURN_NULL(out);
    PS_CHECK_NULL_PTR_RETURN_NULL(in);

    junkData = psAlloc(sizeof(psVector));
    junkData->type = in->type;
    junkData->nalloc = in->numRows * in->numCols;
    junkData->n = junkData->nalloc;
    junkData->data.V = in->data.V[0];  // since psImage data is contiguous...

    if (mask != NULL) {
        if (mask->type.type != PS_TYPE_MASK) {
            char* typeStr;
            PS_TYPE_NAME(typeStr,mask->type.type);
            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageHistogram",
                       PS_ERR_BAD_PARAMETER_TYPE, true,
                       PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE,
                       typeStr, PS_TYPE_MASK_NAME);
            psFree(out);
            psFree(junkData);
            return NULL;
        }
        // stuff the mask data into a psVector struct.
        junkMask = psAlloc(sizeof(psVector));
        junkMask->type = mask->type;
        junkMask->nalloc = mask->numRows * mask->numCols;
        junkMask->n = junkMask->nalloc;
        junkMask->data.V = mask->data.V[0];
    }

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

    psFree(junkMask);
    psFree(junkData);

    return (out);
}

// XXX: Why do we have this function?
float *p_psCalcScaleFactorsFit(int n)
{
    int i = 0;
    float tmp = 0.0;
    float *scalingFactors = (float *)psAlloc(n * sizeof(float));

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

    return (scalingFactors);
}

/*****************************************************************************
p_psCalcScaleFactorsEval(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 scales 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.
 *****************************************************************************/
float *p_psCalcScaleFactorsEval(int n)
{
    int i = 0;
    float tmp = 0.0;

    return p_psCalcScaleFactorsFit(n);

    printf("Should not get here\n");
    float *scalingFactors = (float *)psAlloc(n * sizeof(float));

    for (i = 0; i < n; i++) {
        // scalingFactors[i] = ((2.0 * (float) i)
        // / ((float) (n-1))) - 1.0;
        tmp = (float)(n - i);
        tmp = (M_PI * (tmp - 0.5)) / ((float)n);
        scalingFactors[i] = cos(tmp);
    }
    return (scalingFactors);
}

// XXX: Use a static array of CHebyshev polynomials.
psPolynomial1D **p_psCreateChebyshevPolys(int maxChebyPoly)
{
    if (maxChebyPoly < 1) {
        return(NULL);
    }
    psPolynomial1D **chebPolys = NULL;
    int i = 0;
    int 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((float) x, chebPolys[i]) *
                psPolynomial1DEval((float) y, chebPolys[j]);
        over all pixels (x,y) in the image.
 
XXX: Currently we only support F32.
 *****************************************************************************/
psPolynomial2D* psImageFitPolynomial(psPolynomial2D* coeffs,
                                     const psImage* input)
{
    int x = 0;
    int y = 0;
    int i = 0;
    int j = 0;
    float **sums = NULL;
    psPolynomial1D* *chebPolys = NULL;
    int maxChebyPoly = 0;
    float *cScalingFactors = NULL;
    float *rScalingFactors = NULL;

    // Check for null inputs
    if (input == NULL) {
        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFitPolynomial",
                   PS_ERR_BAD_PARAMETER_NULL, true,
                   PS_ERRORTEXT_psImage_IMAGE_NULL);
        return NULL;
    }

    if (coeffs == NULL) {
        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFitPolynomial",
                   PS_ERR_BAD_PARAMETER_NULL, true,
                   PS_ERRORTEXT_psImage_COEFF_NULL);
        return 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 = (float **)psAlloc(coeffs->nX * sizeof(float *));
    for (i = 0; i < coeffs->nX; i++) {
        sums[i] = (float *)psAlloc(coeffs->nY * sizeof(float));
    }
    // We scale the pixel positions to values
    // between -1.0 and 1.0
    rScalingFactors = p_psCalcScaleFactorsFit(input->numRows);
    cScalingFactors = p_psCalcScaleFactorsFit(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);


    /*
        // Sanity check for the Chebyshevs.
        for (i = 0; i < coeffs->nX; i++) {
            for (j = 0; j < coeffs->nY; j++) {
                tmp = 0.0;
                for (x = 0; x < input->numRows; x++) {
                    tmp +=
                        psPolynomial1DEval
                        (rScalingFactors[x], chebPolys[i]) * psPolynomial1DEval(rScalingFactors[x], chebPolys[j]);
     
                }
                // printf("SUM(Cheby(%d) * Cheby(%d))
                // is %f\n", i, j, tmp);
            }
        }
    */

    // 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++) {
                    sums[i][j] +=
                        input->data.F32[x][y] *
                        psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) *
                        psPolynomial1DEval(cScalingFactors[y], chebPolys[j]);
                }
            }
        }
    }

    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] /= (float)(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);
}

/*****************************************************************************
XXX: Use static variables for chebyshev polynomials and scaling factors. 
 
XXX: Currently we only support F32.
 *****************************************************************************/
psImage* psImageEvalPolynomial(psImage* input, const psPolynomial2D* coeffs)
{
    int x = 0;
    int y = 0;
    int i = 0;
    int j = 0;
    //    float **sums = NULL;
    psPolynomial1D* *chebPolys = NULL;
    int maxChebyPoly = 0;
    float *cScalingFactors = NULL;
    float *rScalingFactors = NULL;
    float polySum = 0.0;

    // 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 = (float **)psAlloc(coeffs->nX * sizeof(float *));
    //    for (i = 0; i < coeffs->nX; i++) {
    //        sums[i] = (float *)psAlloc(coeffs->nY * sizeof(float));
    //    }
    //    for (i = 0; i < coeffs->nX; i++) {
    //        for (j = 0; j < coeffs->nY; j++) {
    //            sums[i][j] = 0.0;
    //        }
    //    }


    // Check for null inputs
    if (input == NULL) {
        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageEvalPolynomial",
                   PS_ERR_BAD_PARAMETER_NULL, true,
                   PS_ERRORTEXT_psImage_IMAGE_NULL);
        return NULL;
    }

    if (coeffs == NULL) {
        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageEvalPolynomial",
                   PS_ERR_BAD_PARAMETER_NULL, true,
                   PS_ERRORTEXT_psImage_COEFF_NULL);
        return NULL;
    }

    // We scale the pixel positions to values between -1.0 and 1.0
    // Use static data structures here.
    rScalingFactors = p_psCalcScaleFactorsEval(input->numRows);
    cScalingFactors = p_psCalcScaleFactorsEval(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(rScalingFactors[x], chebPolys[i]) *
                        psPolynomial1DEval(cScalingFactors[y], chebPolys[j]) *
                        coeffs->coeff[i][j];

                }
            }
            input->data.F32[x][y] = polySum;
        }
    }

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

    psFree(cScalingFactors);
    psFree(rScalingFactors);

    return input;
}
