/** @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.13 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-07-10 20:54:52 $
 *
 *  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 <strings.h>

#include <pslib.h>
#include "pmFPA.h"
#include "pmFlatField.h"
#include "pmFPAMaskWeight.h"
#include "pmFlatFieldErrors.h"


bool pmFlatField(pmReadout *in, const pmReadout *flat)
{
    PS_ASSERT_PTR_NON_NULL(in, false);
    PS_ASSERT_PTR_NON_NULL(in->image, false);
    PS_ASSERT_IMAGE_NON_EMPTY(in->image, false);
    PS_ASSERT_PTR_NON_NULL(flat, false);
    PS_ASSERT_PTR_NON_NULL(flat->image, false);
    PS_ASSERT_IMAGE_NON_EMPTY(flat->image, false);
    if (in->mask) {
        PS_ASSERT_IMAGE_TYPE(in->mask, PS_TYPE_MASK, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(in->mask, in->image, false);
    }
    PS_ASSERT_IMAGE_TYPE(flat->image, in->image->type.type, false);

    psImage *inImage   = in->image;     // Input image
    psImage *inMask    = in->mask;      // Mask for input image
    psImage *flatImage = flat->image;   // Flat-field image

    // Check input image is not larger than flat image; mask is the same size as the input
    if (inImage->numRows > flatImage->numRows || inImage->numCols > flatImage->numCols) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Input image (%dx%d) is larger than flat-field image "
                "(%dx%d).\n", inImage->numCols, inImage->numRows, flatImage->numCols, flatImage->numRows);
        return false;
    }

    // Offsets on the chip
    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
    int x0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.X0");
    int y0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.Y0");

    // Determine offset based on image offset with chip offset: input frame to flat frame
    int yOffset = in->row0 + y0in - flat->col0 - y0flat;
    int xOffset = in->row0 + x0in - flat->col0 - x0flat;

    // Check that offsets are within image limits
    if (inImage->numRows + yOffset > flatImage->numRows ||
            inImage->numCols + xOffset > flatImage->numCols) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Input image (%dx%d) with offsets (%d,%d) is larger than "
                "flat-field image (%dx%d).\n", inImage->numCols, inImage->numRows, xOffset, yOffset,
                flatImage->numCols, flatImage->numRows);
        return false;
    }

    // Macro for all PS types
    #define FLAT_DIVISION_CASE(TYPE, SPECIAL) \
case PS_TYPE_##TYPE: \
    for (int j = 0; j < inImage->numRows; j++) { \
        for (int i = 0; i < inImage->numCols; i++) { \
            if (flatImage->data.TYPE[j][i] <= 0.0) { \
                if (inMask) { \
                    inMask->data.U8[j][i] |= PM_MASK_FLAT; \
                } \
                inImage->data.TYPE[j][i] = SPECIAL; \
            } else { \
                inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \
            } \
        } \
    } \
    break;

    switch (inImage->type.type) {
        FLAT_DIVISION_CASE(U8,  0);
        FLAT_DIVISION_CASE(U16, 0);
        FLAT_DIVISION_CASE(U32, 0);
        FLAT_DIVISION_CASE(U64, 0);
        FLAT_DIVISION_CASE(S8,  0);
        FLAT_DIVISION_CASE(S16, 0);
        FLAT_DIVISION_CASE(S32, 0);
        FLAT_DIVISION_CASE(S64, 0);
        FLAT_DIVISION_CASE(F32, NAN);
        FLAT_DIVISION_CASE(F64, NAN);
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unsupported type for input image: %x\n", inImage->type.type);
        return false;
    }

    return true;
}
