/** @file  pmFlatField.c
 *
 *  @brief Given an input image and a flat field image, pmFlatField shall divide the input image by the flat
 *  field image.
 *
 *  The input image, in, and the flat field image, flat, need not be the same size, since the input image may
 *  already have been trimmed (following overscan subtraction), but the function shall use the offsets in the
 *  image (in->x0 and in->y0) to determine the appropriate offsets to obtain the correct pixel on the flat
 *  field. In the event that the flat image is too small (i.e., pixels on the input image refer to pixels
 *  outside the range of the flat image), the function shall generate an error. Pixels which are negative or
 *  zero in the flat shall be masked in the input image with the value PM_MASK_FLAT. Negative pixels in the
 *  flat may be set to zero so that they are treated identically to zeroes. Any pixels masked in the flat
 *  shall be masked with corresponding values in the output. The function shall not normalize the flat; this
 *  responsibility is left to the caller. This function is basically equivalent to a divide (with psImageOp),
 *  but with care for the region that is divided, checking for negative pixels, and copying of the mask from
 *  the flat to the output.
 *
 *  @author Ross Harman, MHPCC
 *
 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-11-18 19:43:14 $
 *
 *  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 <string.h>

#include "pslib.h"
#include "pmFlatField.h"
#include "pmMaskBadPixels.h"
#include "pmFlatFieldErrors.h"
#include "pmSubtractBias.h"

// XXX: This should be removed when the autoconf stuff handles psConstants.h correctly.
#define PS_WARN_PTR_NON_NULL(NAME) \
if ((NAME) == NULL) { \
    psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \
} \

bool pmFlatField(
    pmReadout *in,
    const pmReadout *flat)
{
    // XXX: Use the proper image and readout asserts.
    PS_ASSERT_PTR_NON_NULL(in, true);
    PS_ASSERT_PTR_NON_NULL(in->image, false);
    PS_ASSERT_PTR_NON_NULL(in->mask, false);
    PS_ASSERT_PTR_NON_NULL(flat, false);
    PS_ASSERT_PTR_NON_NULL(flat->image, false);
    if (in == NULL)
        printf("XXX: NULL\n");

    // XXX: Not sure if this is correct.  Must consult with IfA.
    PS_ASSERT_PTR_NON_NULL(in->mask, false);

    PS_WARN_PTR_NON_NULL(in->parent);
    if (in->parent != NULL) {
        PS_WARN_PTR_NON_NULL(in->parent->concepts);
    }
    int i = 0;
    int j = 0;
    int totOffCol = 0;
    int totOffRow = 0;
    psElemType inType;
    psElemType flatType;
    psElemType maskType;
    psImage *inMask = NULL;
    psImage *flatImage = NULL;

    //
    // Determine trimmed image from metadata.
    //
    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    flatImage = flat->image;
    inMask = in->mask;

    // Check input image and its mask are not larger than flat image

    if (trimmedImg == NULL)
        printf("XXX: 00\n");
    if (flatImage == NULL)
        printf("XXX 01\n");

    if (trimmedImg->numRows>flatImage->numRows || trimmedImg->numCols>flatImage->numCols) {
        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                 PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE,
                 trimmedImg->numRows, trimmedImg->numCols, flatImage->numRows, flatImage->numCols);
        return false;
    }
    if (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols) {
        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                 PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE,
                 inMask->numRows, inMask->numCols, flatImage->numRows, flatImage->numCols);
        return false;
    }

    // Determine total offset based on image offset with chip offset
    totOffCol = trimmedImg->col0 + in->col0;
    totOffRow = trimmedImg->row0 + in->row0;

    // Check that offsets are within image limits
    if(totOffRow>=flatImage->numRows || totOffCol>=flatImage->numCols) {
        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                 PS_ERRORTEXT_pmFlatField_OFFSET_FLAT_IMAGE,
                 totOffRow, totOffCol, flatImage->numRows, flatImage->numCols);
        return false;
    } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) {
        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                 PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE,
                 totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols);
        return false;
    } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {
        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                 PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE,
                 totOffRow, totOffCol, inMask->numRows, inMask->numCols);
        return false;
    }

    // Check for incorrect types
    inType = trimmedImg->type.type;
    flatType = flatImage->type.type;
    maskType = inMask->type.type;
    if(PS_IS_PSELEMTYPE_COMPLEX(inType)) {
        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_pmFlatField_TYPE_INPUT_IMAGE,
                 inType);
        return false;
    } else if(PS_IS_PSELEMTYPE_COMPLEX(flatType)) {
        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_pmFlatField_TYPE_FLAT_IMAGE,
                 flatType);
        return false;
    } else if(maskType != PS_TYPE_MASK) {
        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE,
                 maskType);
        return false;
    } else if(inType != flatType) {
        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_pmFlatField_TYPE_MISMATCH,
                 inType, flatType);
        return false;
    }

    // Macro for all PS types
    #define PM_FLAT_DIVISION(TYPE)                                                                           \
case PS_TYPE_##TYPE:                                                                                         \
    /* Per Eugene's request, use two sets of loops: first to fill mask, second to avoid div with bad pix */  \
    for(j = totOffRow; j < trimmedImg->numRows; j++) {                                                          \
        for(i = totOffCol; i < trimmedImg->numCols; i++) {                                                      \
            if(flatImage->data.TYPE[j][i] <= 0.0) {                                                          \
                /* Negative or zero flat pixels shall be masked in input image as  PM_MASK_FLAT */           \
                inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT;                                        \
                flatImage->data.TYPE[j][i] = 0.0;                                                            \
            }                                                                                                \
        }                                                                                                    \
    }                                                                                                        \
    for(j = totOffRow; j < trimmedImg->numRows; j++) {                                                          \
        for(i = totOffCol; i < trimmedImg->numCols; i++) {                                                      \
            if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                                      \
                /* Module shall divide the input image by the flat-fielded image */                          \
                trimmedImg->data.TYPE[j][i] /= flatImage->data.TYPE[j][i];                                      \
            }                                                                                                \
        }                                                                                                    \
    }                                                                                                        \
    break;

    switch(inType) {
        PM_FLAT_DIVISION(U8);
        PM_FLAT_DIVISION(U16);
        PM_FLAT_DIVISION(U32);
        PM_FLAT_DIVISION(U64);
        PM_FLAT_DIVISION(S8);
        PM_FLAT_DIVISION(S16);
        PM_FLAT_DIVISION(S32);
        PM_FLAT_DIVISION(S64);
        PM_FLAT_DIVISION(F32);
        PM_FLAT_DIVISION(F64);
    default:
        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_pmFlatField_TYPE_UNSUPPORTED,
                 inType);
    }

    return true;
}
