/** @file  psAstrometry.c
*
*  @brief This file defines the basic types for astronomical coordinate
*  transformation
*
*  @ingroup AstroImage
*
*  @author GLG, MHPCC
*
* XXX: We should review the extent of the warning messages on these functions
* when the transformations are not successful.
*
* XXX: Should we implement non-linear cell->chip transforms?
*
*  @version $Revision: 1.12.4.1 $ $Name: not supported by cvs2svn $
*  @date $Date: 2006-02-17 17:13:41 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

/******************************************************************************/
/*  INCLUDE FILES                                                             */
/******************************************************************************/
#include <string.h>
#include <math.h>
#include "pslib.h"

#include "pmAstrometry.h"
#include "pmConcepts.h"

/*****************************************************************************
checkValidImageCoords(): this is a private function which simply determines if
the supplied x,y coordinates are in the range for the supplied psImage.
 
XXX: What about col0 and row0
XXX: This should return a psBool.
XXX: Macro this for speed.
 *****************************************************************************/
static psS32 checkValidImageCoords(
    double x,
    double y,
    psImage* tmpImage)
{
    PS_ASSERT_IMAGE_NON_NULL(tmpImage, 0);

    // The FLT_EPSILON is because -0.0 was failing this.
    if (((x+FLT_EPSILON) < 0.0) || (x > (double)tmpImage->numCols) ||
            ((y+FLT_EPSILON) < 0.0) || (y > (double)tmpImage->numRows)) {
        return (0);
    }

    return (1);
}

/******************************************************************************
 *****************************************************************************/
#define PS_FREE_HIERARCHY 1
static void readoutFree(pmReadout *readout)
{
    if (readout != NULL) {
        psFree(readout->image);
        psFree(readout->mask);
        psFree(readout->weight);
        psFree(readout->analysis);
        #if 0

        psFree(readout->parent);
        #endif

        readout->parent = NULL;
    }
}

static void cellFree(pmCell *cell)
{
    if (cell != NULL) {
        psFree(cell->toChip);
        psFree(cell->toFPA);
        psFree(cell->toSky);
        psFree(cell->concepts);
        psFree(cell->analysis);
        psFree(cell->camera);
        //
        // Set the parent to NULL in all cell->readouts before psFree(cell->readouts)
        // in order to avoid memory reference counter problems.
        //
        #if 0

        for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
            pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
            tmpReadout->parent = NULL;
            if (PS_FREE_HIERARCHY == 1) {
                psFree(tmpReadout);
            }
        }
        psFree(cell->parent);
        #endif

        cell->parent = NULL;

        psFree(cell->readouts);
        psFree(cell->hdu);

    }
}

static void chipFree(pmChip* chip)
{
    if (chip != NULL) {
        psFree(chip->toFPA);
        psFree(chip->fromFPA);
        psFree(chip->concepts);
        psFree(chip->analysis);
        //
        // Set the parent to NULL in all chip->cells before psFree(chip->cells)
        // in order to avoid memory reference counter problems.
        //
        #if 0

        for (psS32 i = 0 ; i < chip->cells->n ; i++) {
            pmCell *tmpCell = (pmCell *) chip->cells->data[i];
            tmpCell->parent = NULL;
            if (PS_FREE_HIERARCHY == 1) {
                psFree(tmpCell);
            }
        }
        psFree(chip->parent);
        #endif

        chip->parent = NULL;
        psFree(chip->cells);
        psFree(chip->hdu);
    }
}


static void FPAFree(pmFPA *fpa)
{
    if (fpa != NULL) {
        psFree(fpa->fromTangentPlane);
        psFree(fpa->toTangentPlane);
        psFree(fpa->projection);
        psFree(fpa->concepts);
        psFree(fpa->analysis);
        psFree(fpa->camera);
        //
        // Set the parent to NULL in all fpa->chips before psFree(fpa->chips)
        // in order to avoid memory reference counter problems.
        //
        #if 0

        for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
            pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
            tmpChip->parent = NULL;
            if (PS_FREE_HIERARCHY == 1) {
                psFree(tmpChip);
            }
        }
        #endif
        psFree(fpa->chips);
        psFree(fpa->hdu);
        psFree(fpa->phu);
    }
}

void p_pmHDUFree(p_pmHDU *hdu)
{
    if (hdu) {
        psFree(hdu->extname);
        psFree(hdu->header);
        psFree(hdu->images);
        psFree(hdu->masks);
        psFree(hdu->weights);
    }
}

// XXX: Verify these default values for row0, col0, rowBins, colBins
// PAP: These values may disappear in the future in favour of values in parent->concepts?
pmReadout *pmReadoutAlloc(pmCell *cell)
{
    pmReadout *tmpReadout = (pmReadout *) psAlloc(sizeof(pmReadout));

    tmpReadout->col0 = 0;
    tmpReadout->row0 = 0;
    tmpReadout->colBins = 0;
    tmpReadout->rowBins = 0;
    tmpReadout->image = NULL;
    tmpReadout->mask = NULL;
    tmpReadout->weight = NULL;
    tmpReadout->analysis = psMetadataAlloc();
    tmpReadout->parent = cell;
    if (cell != NULL) {
        cell->readouts = psArrayAdd(cell->readouts, 1, (psPtr) tmpReadout);
    }
    psMemSetDeallocator(tmpReadout, (psFreeFunc) readoutFree);
    return(tmpReadout);
}

// XXX: Verify these default values for row0, col0.
// PAP: These values may disappear in the future in favour of values in the "concepts"?
pmCell *pmCellAlloc(
    pmChip *chip,
    psMetadata *cameraData,
    const char *name)
{
    pmCell *tmpCell = (pmCell *) psAlloc(sizeof(pmCell));

    tmpCell->col0 = 0;
    tmpCell->row0 = 0;
    tmpCell->toChip = NULL;
    tmpCell->toFPA = NULL;
    tmpCell->toSky = NULL;
    tmpCell->concepts = psMetadataAlloc();
    psBool rc = psMetadataAddStr(tmpCell->concepts, PS_LIST_HEAD, "CELL.NAME", 0, NULL, name);
    if (rc == false) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CELL.NAME to metadata.\n");
    }
    tmpCell->camera = psMemIncrRefCounter(cameraData);
    tmpCell->analysis = psMetadataAlloc();
    tmpCell->readouts = psArrayAlloc(0);
    tmpCell->parent = chip;
    if (chip != NULL) {
        chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell);
    }
    tmpCell->process = true;            // All cells are processed by default
    tmpCell->exists = false;            // Not yet read in
    tmpCell->hdu = NULL;

    pmConceptsBlankCell(tmpCell);

    psMemSetDeallocator(tmpCell, (psFreeFunc) cellFree);
    return(tmpCell);
}

// XXX: Verify these default values for row0, col0.
// PAP: row0, col0 may disappear in the future in favour of storing values in the "concepts".
pmChip *pmChipAlloc(
    pmFPA *fpa,
    const char *name)
{
    pmChip *tmpChip = (pmChip *) psAlloc(sizeof(pmChip));

    tmpChip->col0 = 0;
    tmpChip->row0 = 0;
    tmpChip->toFPA = NULL;
    tmpChip->fromFPA = NULL;
    tmpChip->concepts = psMetadataAlloc();
    psBool rc = psMetadataAddStr(tmpChip->concepts, PS_LIST_HEAD, "CHIP.NAME", 0, NULL, name);
    if (rc == false) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CHIP.NAME to metadata.\n");
    }
    tmpChip->analysis = psMetadataAlloc();
    tmpChip->cells = psArrayAlloc(0);
    tmpChip->parent = fpa;
    if (fpa != NULL) {
        fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip);
    }
    tmpChip->process = true;            // Work on all chips, by default
    tmpChip->exists = false;            // Not read in yet
    tmpChip->hdu = NULL;

    pmConceptsBlankChip(tmpChip);

    psMemSetDeallocator(tmpChip, (psFreeFunc) chipFree);
    return(tmpChip);
}

pmFPA *pmFPAAlloc(const psMetadata *camera)
{
    pmFPA *tmpFPA = (pmFPA *) psAlloc(sizeof(pmFPA));

    tmpFPA->fromTangentPlane = NULL;
    tmpFPA->toTangentPlane = NULL;
    tmpFPA->projection = NULL;
    tmpFPA->concepts = psMetadataAlloc();
    tmpFPA->analysis = NULL;
    tmpFPA->camera = psMemIncrRefCounter((psPtr)camera);
    tmpFPA->chips = psArrayAlloc(0);
    tmpFPA->hdu = NULL;
    tmpFPA->phu = NULL;

    pmConceptsBlankFPA(tmpFPA);

    psMemSetDeallocator(tmpFPA, (psFreeFunc) FPAFree);
    return(tmpFPA);
}

p_pmHDU *p_pmHDUAlloc(const char *extname)
{
    p_pmHDU *hdu = psAlloc(sizeof(p_pmHDU));
    psMemSetDeallocator(hdu, (psFreeFunc)p_pmHDUFree);

    hdu->extname = psStringCopy(extname);
    hdu->header = NULL;
    hdu->images = NULL;
    hdu->masks = NULL;
    hdu->weights = NULL;

    return hdu;
}

static psBool cellCheckParents(pmCell *cell)
{
    if (cell == NULL) {
        return(true);
    }
    psBool flag = true;

    for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
        pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
        PS_ASSERT_PTR_NON_NULL(tmpReadout, false);
        if (tmpReadout->parent != cell) {
            tmpReadout->parent = cell;
            flag = false;
        }
    }
    return(flag);
}

static psBool chipCheckParents(pmChip *chip)
{
    if (chip == NULL) {
        return(true);
    }
    psBool flag = true;

    for (psS32 i = 0 ; i < chip->cells->n ; i++) {
        pmCell *tmpCell = (pmCell *) chip->cells->data[i];
        PS_ASSERT_PTR_NON_NULL(tmpCell, false);
        if (tmpCell->parent != chip) {
            tmpCell->parent = chip;
            flag = false;
        }

        flag &= cellCheckParents(tmpCell);
    }
    return(flag);
}

psBool pmFPACheckParents(pmFPA *fpa)
{
    if (fpa == NULL) {
        return(true);
    }
    psBool flag = true;

    for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
        pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
        PS_ASSERT_PTR_NON_NULL(tmpChip, false);
        if (tmpChip->parent != fpa) {
            tmpChip->parent = fpa;
            flag = false;
        }

        flag &= chipCheckParents(tmpChip);
    }
    return(flag);
}



/*****************************************************************************
 *****************************************************************************/

// Set cells within a chip to be processed or not
static bool setCellsProcess(const pmChip *chip, // Chip of interest
                            bool process  // Process this chip?
                           )
{
    PS_ASSERT_PTR_NON_NULL(chip, false);

    psArray *cells = chip->cells;       // Component cells
    if (! cells) {
        return false;
    }
    for (int i = 0; i < cells->n; i++) {
        pmCell *tmpCell = cells->data[i]; // Cell of interest
        if (tmpCell) {
            tmpCell->process = process;
        }
    }

    return true;
}

/*****************************************************************************
 *****************************************************************************/
bool pmFPASelectChip(
    pmFPA *fpa,
    int chipNum)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);

    psArray *chips = fpa->chips;        // Component chips
    if ((chips == NULL) || (chipNum >= chips->n)) {
        return(false);
    }

    for (int i = 0 ; i < chips->n ; i++) {
        pmChip *tmpChip = (pmChip *) chips->data[i];
        if (tmpChip == NULL) {
            continue;
        }
        if (i == chipNum) {
            tmpChip->process = true;
            setCellsProcess(tmpChip, true);
        } else {
            tmpChip->process = false;
            setCellsProcess(tmpChip, false);
        }

    }

    return true;
}


/*****************************************************************************
XXX: The SDRS is ambiguous on a few things:
    Whether or not the other chips should be set process=true. [PAP: No]
    Should we return the number of chip process=true before or after they're set, [PAP: After]
 *****************************************************************************/
/**
 *
 * pmFPAExcludeChip shall set process to false only for the specified chip
 * number (chipNum). In the event that the specified chip number does not exist
 * within the fpa, the function shall generate a warning, and perform no action.
 * The function shall return the number of chips within the fpa that have process
 * set to true.
 *
 */
int pmFPAExcludeChip(
    pmFPA *fpa,
    int chipNum)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);

    psArray *chips = fpa->chips;        // Component chips
    if (chips == NULL) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n");
        return(0);
    }
    if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum);
        return(0);
    }

    int numChips = 0;                   // Number of chips to be processed
    for (int i = 0 ; i < chips->n ; i++) {
        pmChip *tmpChip = (pmChip *) chips->data[i]; // Chip of interest
        if (tmpChip != NULL) {
            if (i == chipNum) {
                tmpChip->process = false;
                setCellsProcess(tmpChip, false); // Wipe out the cell as well
            } else if (tmpChip->process) {
                numChips++;
            }
        }
    }

    return(numChips);
}


bool pmCellSetWeights(pmCell *cell // Cell for which to set weights
                     )
{
    float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise
    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section

    p_pmHDU *hdu = cell->hdu;           // The data unit, containing the weight and mask originals
    if (!hdu) {
        pmChip *chip = cell->parent;    // The parent chip
        if (chip->hdu) {
            hdu = chip->hdu;
        } else {
            pmFPA *fpa = chip->parent;  // The parent FPA
            hdu = fpa->hdu;
            if (!hdu) {
                psError(PS_ERR_UNKNOWN, true, "Unable to find the HDU in the hierarchy!\n");
                return false;
            }
        }
    }

    psArray *pixels = hdu->images;      // Array of images
    psArray *weights = hdu->weights;    // Array of weight images
    psArray *masks = hdu->masks;        // Array of mask images
    // Generate the weights and masks if required
    if (! weights) {
        weights = psArrayAlloc(pixels->n);
        for (int i = 0; i < pixels->n; i++) {
            psImage *image = pixels->data[i];
            weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
            psImageInit(weights->data[i], 0.0);
        }
        hdu->weights = weights;
    }
    if (! masks) {
        masks = psArrayAlloc(pixels->n);
        for (int i = 0; i < pixels->n; i++) {
            psImage *image = pixels->data[i];
            masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
            psImageInit(masks->data[i], 0);
        }
        hdu->masks = masks;
    }

    // Set the pixels
    psArray *readouts = cell->readouts; // Array of readouts
    for (int i = 0; i < readouts->n; i++) {
        pmReadout *readout = readouts->data[i]; // The readout of interest

        if (! readout->weight) {
            readout->weight = weights->data[i];
        }
        if (! readout->mask) {
            readout->mask = masks->data[i];
        }

        // Mask is already set to 0

        // Set weight image to the variance = g*f + rn^2
        psImage *image = psImageSubset(readout->image, *trimsec); // The pixels
        psImage *weight = psImageSubset(readout->weight, *trimsec); // The weight map
        psBinaryOp(weight, image, "*", psScalarAlloc(gain, PS_TYPE_F32));
        psBinaryOp(weight, weight, "+", psScalarAlloc(readnoise*readnoise, PS_TYPE_F32));
    }

    return true;
}


