/** @file  pmFPA.c
*
*  @brief This file defines the basic types for the FPA hierarchy
*
*  @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.1.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 <assert.h>
#include "pslib.h"

#include "pmFPA.h"
#include "pmConcepts.h"
#include "pmMaskBadPixels.h"

/******************************************************************************
 *****************************************************************************/
#define PS_FREE_HIERARCHY 1
#define PARENT_LINKS 0

static void readoutFree(pmReadout *readout)
{
    if (readout != NULL) {
        psTrace(__func__, 9, "Removing readout %lx from cell %lx...\n", (size_t) readout, (size_t) readout->parent);
        if (readout->parent) {
            psArray *readouts = readout->parent->readouts;
            for (int i = 0; i < readouts->n; i++) {
                if (readouts->data[i] == readout) {
                    //                    pmReadout *tmpReadout = readouts->data[i];
                    readouts->data[i] = NULL;
                    #if PARENT_LINKS

                    psFree(tmpReadout);
                    #endif

                    break;
                }
            }
        }
        psTrace(__func__, 9, "Freeing readout %lx\n", (size_t) readout);

        #if PARENT_LINKS

        psFree(readout->parent);
        #endif

        psFree(readout->image);
        psFree(readout->mask);
        psFree(readout->weight);
        psFree(readout->analysis);
        psFree(readout->bias);
    }
}

static void cellFree(pmCell *cell)
{
    if (cell != NULL) {
        psTrace(__func__, 9, "Removing cell %lx from chip %lx...\n", (size_t)cell, (size_t)cell->parent);
        if (cell->parent) {
            psArray *cells = cell->parent->cells;
            for (int i = 0; i < cells->n; i++) {
                if (cells->data[i] == cell) {
                    //                    pmCell *tmpCell = cells->data[i];
                    cells->data[i] = NULL;
                    #if PARENT_LINKS

                    psFree(tmpCell);
                    #endif

                    break;
                }
            }
        }
        psTrace(__func__, 9, "Freeing cell %lx\n", (size_t)cell);

        pmCellFreeReadouts(cell);
        psFree(cell->readouts);
        #if PARENT_LINKS

        psFree(cell->parent);
        #endif

        psFree(cell->toChip);
        psFree(cell->toFPA);
        psFree(cell->toSky);
        psFree(cell->concepts);
        psFree(cell->analysis);
        psFree(cell->camera);
        psFree(cell->hdu);
    }
}

static void chipFree(pmChip* chip)
{
    if (chip != NULL) {
        psTrace(__func__, 9, "Removing chip %lx from fpa %lx...\n", (size_t)chip, (size_t)chip->parent);
        if (chip->parent) {
            psArray *chips = chip->parent->chips;
            for (int i = 0; i < chips->n; i++) {
                if (chips->data[i] == chip) {
                    //                    pmChip *tmpChip = chips->data[i];
                    chips->data[i] = NULL;
                    #if PARENT_LINKS

                    psFree(tmpChip);
                    #endif

                    break;
                }
            }
        }
        psTrace(__func__, 9, "Freeing chip %lx\n", (size_t)chip);

        pmChipFreeCells(chip);
        psFree(chip->cells);
        #if PARENT_LINKS

        psFree(chip->parent);
        #endif

        psFree(chip->toFPA);
        psFree(chip->fromFPA);
        psFree(chip->concepts);
        psFree(chip->analysis);
        psFree(chip->hdu);
    }
}


static void FPAFree(pmFPA *fpa)
{
    if (fpa != NULL) {
        psTrace(__func__, 9, "Freeing fpa %lx\n", (size_t)fpa);
        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.
        //
        psArray *chips = fpa->chips;
        for (psS32 i = 0 ; i < chips->n ; i++) {
            pmChip *tmpChip = chips->data[i];
            if (! tmpChip) {
                continue;
            }
            chips->data[i] = NULL;
            tmpChip->parent = NULL;
            if (PS_FREE_HIERARCHY == 1) {
                psFree(tmpChip);
            }
        }
        psFree(fpa->chips);
        psFree(fpa->hdu);
        psFree(fpa->phu);
    }
}

void pmCellFreeReadouts(pmCell *cell)
{
    //
    // Set the parent to NULL in all cell->readouts before psFree(cell->readouts)
    // in order to avoid memory reference counter problems.
    //
    psArray *readouts = cell->readouts;
    for (psS32 i = 0 ; i < readouts->n ; i++) {
        pmReadout *tmpReadout = readouts->data[i];
        if (! tmpReadout) {
            continue;
        }
        readouts->data[i] = NULL;
        tmpReadout->parent = NULL;
        psTrace(__func__, 9, "Will now free readout %lx...\n", (size_t)tmpReadout);
        psFree(tmpReadout); // Drop the readout
        #if PARENT_LINKS

        psFree(cell); // Decrement reference counter on cell, since readout->parent isn't held any more
        #endif

    }
    cell->readouts = psArrayRealloc(cell->readouts, 0);
    cell->readouts->n = 0;
}


void pmChipFreeCells(pmChip *chip)
{
    //
    // Set the parent to NULL in all chip->cells before psFree(chip->cells)
    // in order to avoid memory reference counter problems.
    //
    psArray *cells = chip->cells;
    for (int i = 0 ; i < cells->n ; i++) {
        pmCell *tmpCell = cells->data[i];
        if (! tmpCell) {
            continue;
        }
        cells->data[i] = NULL;
        tmpCell->parent = NULL;
        pmCellFreeReadouts(tmpCell);// Drop all readouts the cell holds
        psFree(tmpCell);            // Drop the cell
        #if PARENT_LINKS

        psFree(chip); // Decrement reference counter on chip, since cell->parent isn't held any more
        #endif

    }
    chip->cells = psArrayRealloc(chip->cells, 0);
    chip->cells->n = 0;
}



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->bias = psListAlloc(NULL);
    tmpReadout->analysis = psMetadataAlloc();
    #if PARENT_LINKS

    tmpReadout->parent = psMemIncrRefCounter(cell);
    #else

    tmpReadout->parent = cell;
    #endif

    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);
    #if PARENT_LINKS

    tmpCell->parent = psMemIncrRefCounter(chip);
    #else

    tmpCell->parent = chip;
    #endif

    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);
    #if PARENT_LINKS

    tmpChip->parent = psMemIncrRefCounter(fpa);
    #else

    tmpChip->parent = fpa;
    #endif

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

bool pmChipSelectCell(pmChip *chip,
                      int cellNum
                     )
{
    assert(chip);

    psArray *cells = chip->cells;       // Component cells
    if (!cells || cellNum > cells->n) {
        return false;
    }

    for (int i = 0; i < cells->n; i++) {
        pmCell *cell = cells->data[i];
        if (!cell) {
            continue;
        }
        cell->process = (i == cellNum);
    }

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

int pmChipExcludeCell(pmChip *chip,
                      int cellNum
                     )
{
    assert(chip);

    psArray *cells = chip->cells;       // The component cells
    if (!cells || cellNum > cells->n) {
        return 0;
    }

    int numCells = 0;                   // Number of cells to be processed
    for (int i = 0; i < cells->n; i++) {
        pmCell *cell = cells->data[i];
        if (!cell) {
            continue;
        }
        if (i == cellNum) {
            cell->process = false;
        } else {
            numCells++;
        }
    }

    return numCells;
}


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
    float saturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION"); // Saturation level
    float bad = psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD"); // Bad level

    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 = psMemIncrRefCounter(psImageSubset(weights->data[i], *trimsec));
        }
        if (! readout->mask) {
            readout->mask = psMemIncrRefCounter(psImageSubset(masks->data[i], *trimsec));
        }

        // Set up the mask
        psImage *image = readout->image;// Pixels
        psImage *mask = readout->mask;  // Mask image
        for (int i = 0; i < image->numRows; i++) {
            for (int j = 0; j < image->numCols; j++) {
                if (image->data.F32[i][j] > saturation) {
                    mask->data.U8[i][j] = PM_MASK_SAT;
                }
                if (image->data.F32[i][j] < bad) {
                    mask->data.U8[i][j] = PM_MASK_BAD;
                }
            }
        }

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

    return true;
}


