/** @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.17 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-11-06 02:42:03 $
 *
 *  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"


bool pmFlatField(psReadout *in, const psReadout *flat)
{
    int i = 0;
    int j = 0;
    int totOffCol = 0;
    int totOffRow = 0;
    psElemType inType;
    psElemType flatType;
    psElemType maskType;
    psImage *inImage = NULL;
    psImage *inMask = NULL;
    psImage *flatImage = NULL;


    // Check for nulls
    if (in == NULL) {
        return true;       // Readout may not have data in it
    } else if(flat==NULL) {
        psError( PS_ERR_BAD_PARAMETER_NULL, true,
                 PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT);
        return false;
    }

    inImage = in->image;
    flatImage = flat->image;
    if (inImage == NULL) {
        psError( PS_ERR_BAD_PARAMETER_NULL, true,
                 PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE);
        return false;
    } else if(flatImage == NULL) {
        psError( PS_ERR_BAD_PARAMETER_NULL, true,
                 PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE);
        return false;
    } else if(in->mask == NULL) {
        in->mask = psImageAlloc(inImage->numCols, inImage->numRows, PS_TYPE_MASK);
        memset(in->mask->data.V[0], 0, inImage->numCols*inImage->numRows*PSELEMTYPE_SIZEOF(PS_TYPE_MASK));
    }
    inMask = in->mask;

    // Check input image and its mask are not larger than flat image
    if(inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) {
        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                 PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE,
                 inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols);
        return false;
    } else 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 = inImage->col0 + in->col0;
    totOffRow = inImage->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>=inImage->numRows || totOffCol>=inImage->numCols) {
        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                 PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE,
                 totOffRow, totOffCol, inImage->numRows, inImage->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 = inImage->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 < inImage->numRows; j++) {                                                          \
        for(i = totOffCol; i < inImage->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 < inImage->numRows; j++) {                                                          \
        for(i = totOffCol; i < inImage->numCols; i++) {                                                      \
            if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                                      \
                /* Module shall divide the input image by the flat-fielded image */                          \
                inImage->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;
}
