//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
// one that was being worked on.
//////////////////////////////////////////////////////////////////////////////////////////////////////////////


/** @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.12.1 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-02-17 17:13:41 $
 *
 *  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 "pmFlatField.h"
#include "pmMaskBadPixels.h"
#include "pmFlatFieldErrors.h"


bool pmFlatField(pmReadout *in, const pmReadout *flat)
{
    int i = 0;
    int j = 0;
    int totOffCol = 0;
    int totOffRow = 0;
    psElemType inType;
    psElemType flatType;
    psElemType maskType;

    // 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;
    }

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

    // 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");

    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;
    }

    // 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;
    }
    if (inMask && (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 + y0in - flatImage->col0 - y0flat;
    totOffRow = inImage->row0 + x0in - flatImage->row0 - x0flat;

    // 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(inMask && (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(inMask && inMask->type.type != 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 */           \
                if (inMask) {                                                                                \
                    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 && !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;
}
