/** @file  pmNonLinear.c
 *
 *  Provides polynomial or table lookup non-linearity corrections to readouts.
 *
 *  @author George Gusciora, MHPCC
 *
 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-10-14 21:44:08 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#if HAVE_CONFIG_H
#include <config.h>
#endif

#include<stdio.h>
#include<math.h>

#include "pmNonLinear.h"

/******************************************************************************
pmNonLinearityLookup(): This routine will take an psReadout image as input
and a 1-D polynomial.  For each pixel in the input image, the polynomial will
be evaluated at that pixels value, and the image pixel will then be set to
that value.
 *****************************************************************************/

psReadout *pmNonLinearityPolynomial(psReadout *in,
                                    const psPolynomial1D *coeff)
{
    int i;
    int j;

    PS_CHECK_NULL_PTR_RETURN_NULL(in);
    PS_CHECK_NULL_PTR_RETURN_NULL(in->image);
    PS_CHECK_NULL_PTR_RETURN_NULL(coeff);

    for (i=0;i<(in->image)->numRows;i++) {
        for (j=0;j<(in->image)->numCols;j++) {
            (in->image)->data.F32[i][j] = psPolynomial1DEval((in->image)->data.F32[i][j], coeff);
        }
    }
    return(in);
}


/******************************************************************************
pmNonLinearityLookup(): This routine will take an psReadout image as input
and two input vectors, which constitute a lookup table.  For each pixel in
the input image, that pixels value will be determined in the input vector
inFluxe, and the corresponding value in outFlux.  The image pixel will then
be set to the value from outFlux.
 *****************************************************************************/
psReadout *pmNonLinearityLookup(psReadout *in,
                                const psVector *inFlux,
                                const psVector *outFlux)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(in);
    PS_CHECK_NULL_PTR_RETURN_NULL(in->image);
    PS_CHECK_NULL_PTR_RETURN_NULL(inFlux);
    PS_CHECK_NULL_PTR_RETURN_NULL(outFlux);

    int i;
    int j;
    int tableSize = inFlux->n;
    int binNum;
    psScalar x;
    int numPixels = 0;
    float slope;

    x.type.type = PS_TYPE_F32;
    if (inFlux->n != outFlux->n) {
        // XXX: psWarning("pmNonLinearityLookup(): input vectors have different size");
        tableSize = PS_MIN(inFlux->n, outFlux->n);
    }
    if (inFlux->n < 2) {
        // XXX: psWarning("pmNonLinearityLookup(): input vector less than 2 elements");
        return(in);
    }
    for (i=0;i<(in->image)->numRows;i++) {
        for (j=0;j<(in->image)->numCols;j++) {
            x.data.F32 = (in->image)->data.F32[i][j];

            binNum = p_psVectorBinDisect((psVector *)inFlux, &x);

            if (binNum == -2) {
                // We get here if x is below the table lookup range.
                (in->image)->data.F32[i][j] = outFlux->data.F32[0];
                numPixels++;

            } else if (binNum == -1) {
                // We get here if x is above the table lookup range.
                (in->image)->data.F32[i][j] = outFlux->data.F32[tableSize-1];

                numPixels++;
            } else {
                // Perform linear interpolation.
                slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) /
                        (inFlux->data.F32[binNum+1]  - inFlux->data.F32[binNum]);
                (in->image)->data.F32[i][j] = outFlux->data.F32[binNum] +
                                              ((x.data.F32 - inFlux->data.F32[binNum]) * slope);
            }
        }
    }
    if (numPixels > 0) {
        // XXX: psWarning("Warning: pmNonLinearityLookup(): %d pixels outside table.", numPixels);
    }
    return(in);
}
