/** @file  pmNonLinear.c
 *
 *  Provides polynomial or table lookup non-linearity corrections to readouts.
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-09-28 20:43:52 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 *  XXX: The SDR is silent about image types.  Only F32 was implemented.
 *
 */

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

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

#include "pmNonLinear.h"

/******************************************************************************
pmNonLinearityLookup(): This routine will take an pmReadout 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.
 *****************************************************************************/

pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout,
                                    const psPolynomial1D *input1DPoly)
{
    PS_ASSERT_PTR_NON_NULL(inputReadout, NULL);
    PS_ASSERT_PTR_NON_NULL(inputReadout->image, NULL);
    PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
    PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL);

    psS32 i;
    psS32 j;

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


/******************************************************************************
pmNonLinearityLookup(): This routine will take an pmReadout 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.
 *****************************************************************************/
pmReadout *pmNonLinearityLookup(pmReadout *inputReadout,
                                const psVector *inFlux,
                                const psVector *outFlux)
{
    PS_ASSERT_PTR_NON_NULL(inputReadout,NULL);
    PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL);
    PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
    PS_ASSERT_PTR_NON_NULL(inFlux,NULL);
    if (inFlux->n < 2) {
        psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): input vector less than 2 elements.  Returning inputReadout image.");
        return(inputReadout);
    }
    PS_ASSERT_PTR_NON_NULL(outFlux,NULL);
    psS32 tableSize = inFlux->n;
    if (inFlux->n != outFlux->n) {
        tableSize = PS_MIN(inFlux->n, outFlux->n);
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): input vectors have different sizes (%d, %d)\n", inFlux->n, outFlux->n);
    }
    PS_ASSERT_VECTOR_TYPE(inFlux, PS_TYPE_F32, NULL);
    PS_ASSERT_VECTOR_TYPE(outFlux, PS_TYPE_F32, NULL);

    psS32 i;
    psS32 j;
    psS32 binNum;
    psScalar x;
    psS32 numPixels = 0;
    psF32 slope;

    x.type.type = PS_TYPE_F32;
    for (i=0;i<inputReadout->image->numRows;i++) {
        for (j=0;j<inputReadout->image->numCols;j++) {
            x.data.F32 = inputReadout->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.
                inputReadout->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.
                inputReadout->image->data.F32[i][j] = outFlux->data.F32[tableSize-1];
                numPixels++;

            } else if (binNum < -2) {
                // We get here if there was some other problem.
                psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): Could not perform p_psVectorBinDisect().  Returning inputReadout image.");
                return(inputReadout);
                numPixels++;
            } else {
                // Perform linear interpolation.
                slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) /
                        (inFlux->data.F32[binNum+1]  - inFlux->data.F32[binNum]);
                inputReadout->image->data.F32[i][j] = outFlux->data.F32[binNum] +
                                                      ((x.data.F32 - inFlux->data.F32[binNum]) * slope);
            }
        }
    }
    if (numPixels > 0) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels);
    }
    return(inputReadout);
}
