/** @file  pmNonLinear.c
 *
 *  Provides polynomial or table lookup non-linearity corrections to readouts.
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-12-05 20:49:40 $
 *
 *  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"
#include "pmSubtractBias.h"

// XXX: Remove, autoconf must be
#define PS_WARN_PTR_NON_NULL(NAME) \
if ((NAME) == NULL) { \
    psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \
} \
/******************************************************************************
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);
    PS_WARN_PTR_NON_NULL(inputReadout->parent);
    if (inputReadout->parent != NULL) {
        PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
    }

    //
    // Determine trimmed image from metadata.
    //
    psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);

    for (psS32 i=0;i<trimmedImg->numRows;i++) {
        for (psS32 j=0;j<trimmedImg->numCols;j++) {
            trimmedImg->data.F32[i][j] = psPolynomial1DEval(input1DPoly,
                                         trimmedImg->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.
 
XXX: Must assert that filename exists.  This should probably happen in
the lookup files.
 *****************************************************************************/
pmReadout *pmNonLinearityLookup(
    pmReadout *inputReadout,
    const char *filename
)
{
    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_WARN_PTR_NON_NULL(inputReadout->parent);
    if (inputReadout->parent != NULL) {
        PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
    }
    //
    // Determine trimmed image from metadata.
    //
    psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);

    psLookupTable *tmpLT = psLookupTableAlloc(filename, "%f %f", 0);
    psS32 numLines = psLookupTableRead(tmpLT);
    psS32 numPixels = 0;
    if (numLines < 2) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: Lookup Table is too small.  Returning original pmReadout.\n");
    } else {
        for (psS32 i=0;i<trimmedImg->numRows;i++) {
            for (psS32 j=0;j<trimmedImg->numCols;j++) {
                psF64 tmpD = psLookupTableInterpolate(tmpLT, trimmedImg->data.F32[i][j], 1);
                if (!isnan(tmpD)) {
                    trimmedImg->data.F32[i][j] = tmpD;
                } else {
                    numPixels++;
                }
            }
        }
        if (numPixels > 0) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels);
        }
    }

    return(inputReadout);
}
