/** @file psastroMaskUpdates.c
 *
 *  @brief
 *
 *  @ingroup libpsastro
 *
 *  @author IfA
 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-07 02:03:34 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

# include "psastroInternal.h"

# define ESCAPE { \
  psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate"); \
  psFree (view); \
  return false; \
}

/**
 * create a mask or mask regions based on the collection of reference stars that * are in the vicinity of each chip
 */
bool psastroMaskUpdates (pmConfig *config) {

    bool status;
    pmChip *chip = NULL;
    pmCell *cell = NULL;
    pmReadout *readout = NULL;
    float zeropt, exptime, GHOST_MAX_MAG;

    psImageMaskType ghostMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels
    psImageMaskType glintMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for glint pixels (overload ghost)
    psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels
    psImageMaskType starMaskValue  = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels
    psImageMaskType crosstalkMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for crosstalk ghosts

    // psImageMaskType maskBlank  = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels

    // select the current recipe
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE);
    if (!recipe) {
        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe");
        return false;
    }

    // this step is optional
    bool REFSTAR_MASK                      = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK");
    if (!REFSTAR_MASK) return true;

    // convert star positions to crosstalk artifact positions and add to the readout->analysis data
    bool CROSSTALK_MASK                    = psMetadataLookupBool (&status, recipe, "CROSSTALK_MASK");
    if (CROSSTALK_MASK && !psastroLoadCrosstalk (config)) {
      psError(PSASTRO_ERR_CONFIG, false, "Error loading crosstalk data");
      return(false);
    }

    // convert star positions to ghost positions and add to the readout->analysis data
    if (!psastroLoadGhosts (config)) {
        psError(PSASTRO_ERR_CONFIG, false, "Error loading ghosts");
        return false;
    }
    bool COUNT_GHOSTS = psMetadataLookupF32 (&status, recipe, "REFSTAR_COUNT_GHOSTS");
    int nGhosts = 0;

    // convert star positions to glint positions and add to the fpa->analysis data
    if (!psastroLoadGlints (config)) {
        psError(PSASTRO_ERR_CONFIG, false, "Error loading glints");
        return false;
    }

    psLogMsg ("psastro", PS_LOG_INFO, "generating a bright-star mask");

    bool REFSTAR_MASK_BLEED                = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED");

    double REFSTAR_MASK_MAX_MAG            = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_MAX_MAG");
    double REFSTAR_MASK_SATSTAR_MAG_SLOPE  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_SLOPE");
    double REFSTAR_MASK_SATSTAR_MAG_MAX    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_MAX");
    double REFSTAR_MASK_SATSTAR_POS_ZERO   = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_POS_ZERO");
    double REFSTAR_MASK_SATSPIKE_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_MAG_SLOPE");
    double REFSTAR_MASK_SATSPIKE_MAG_MAX   = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_MAG_MAX");
    double REFSTAR_MASK_SATSPIKE_WIDTH     = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_WIDTH");
    double REFSTAR_MASK_SATSPIKE_OFFSET    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_OFFSET");
    double REFSTAR_MASK_BLEED_MAG_MAX      = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX");
    double REFSTAR_MASK_BLEED_MAG_SLOPE    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE");

    double REFSTAR_MASK_CROSSTALK_MAG_MAX  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX");
    double REFSTAR_MASK_CROSSTALK_MAG_SLOPE= psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_SLOPE");

    // select the input data sources
    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    if (!input) {
        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
        return false;
    }
    pmFPA *fpa = input->fpa;

    // really error-out here?  or just skip?
    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &GHOST_MAX_MAG, fpa, recipe)) {
        psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
        return false;
    }

    // recipe values are given in instrumental magnitudes
    // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime)
    float MagOffset = zeropt + 2.5*log10(exptime);
    REFSTAR_MASK_MAX_MAG += MagOffset;
    REFSTAR_MASK_SATSTAR_MAG_MAX += MagOffset;
    REFSTAR_MASK_SATSPIKE_MAG_MAX += MagOffset;
    REFSTAR_MASK_BLEED_MAG_MAX += MagOffset;
    GHOST_MAX_MAG += MagOffset;

    // select the output mask image :: we mosaic to chip mosaic format
    pmFPAfile *outMask = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUTPUT.MASK");
    if (!outMask) {
        psError(PSASTRO_ERR_CONFIG, true, "Can't find output mask");
        return false;
    }
    pmFPA *fpaMask = outMask->fpa;

    // select the reference mask fpa :: we use this to determine cell boundaries
    pmFPAfile *refMask = psMetadataLookupPtr (NULL, config->files, "PSASTRO.REFMASK");
    if (!refMask) {
        psError(PSASTRO_ERR_CONFIG, true, "Can't find mask reference");
        return false;
    }

    psLogMsg ("psastro", PS_LOG_INFO, "loading output and reference mask files");

    // double POSANGLE = PM_RAD_DEG * psMetadataLookupF64 (&status, fpa->concepts, "FPA.POSANGLE");
    // psAssert (status, "POSANGLE missing");
    double ROTANGLE = PM_RAD_DEG * psMetadataLookupF64 (&status, fpa->concepts, "FPA.ROTANGLE");
    psAssert (status, "ROTANGLE missing");

    // de-activate all files except PSASTRO.INPUT.MASK and PSASTRO.OUTPUT.MASK
    pmFPAfileActivate (config->files, false, NULL);
    pmFPAfileActivate (config->files, true, "PSASTRO.INPUT.MASK");
    pmFPAfileActivate (config->files, true, "PSASTRO.OUTPUT.MASK");
    pmFPAfileActivate (config->files, true, "PSASTRO.REFMASK");

    pmFPAview *view = pmFPAviewAlloc (0);
    pmFPAview *viewMask = pmFPAviewAlloc (0);

    // open/load files as needed
    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE;

    // this loop selects the matched stars for all chips
    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
        if (!chip->process || !chip->file_exists) { continue; }
        if (!chip->fromFPA) { continue; }

        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE;

        pmChip *refChip  = pmFPAviewThisChip (view, refMask->fpa);

        // load sequence for mask corresponding to this chip (XXX this is needed if the input mask is not the same format as the astrometry file
        *viewMask = *view;
        while ((cell = pmFPAviewNextCell (viewMask, fpaMask, 1)) != NULL) {
            psTrace ("psastro", 4, "Mask Cell %d: %x %x\n", viewMask->cell, cell->file_exists, cell->process);
            if (!cell->process || !cell->file_exists) { continue; }
            if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_BEFORE)) ESCAPE;

            while ((readout = pmFPAviewNextReadout (viewMask, fpaMask, 1)) != NULL) {
                if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_BEFORE)) ESCAPE;
                if (! readout->data_exists) { continue; }
            }
        }

        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
            if (!cell->process || !cell->file_exists) { continue; }

            // the input mask is a chip-mosaic image
            // the output mask is a chip-mosaic image
            // we mark the masked pixels in the chip space, BUT
            // we need to find the ends of the cells for the bleeds

            // we mask pixels on the input mask image (chip-mosaic)
            // pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa);
            // pmReadout *readoutMask = NULL;
            // if (cellMask->readouts->n) {
            //     readoutMask = cellMask->readouts->data[0];
            // }

            // process each of the readouts
            // XXX there can only be one readout per chip in astrometry, right?
            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
                if (! readout->data_exists) { continue; }

                // XXX why not do this?
                pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
                if (!readoutMask) continue;

                // select the raw objects for this readout
                // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c:
                // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only
                // applied to the SUBSETs
                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
                if (!refstars) continue;

                // we need to generate the following masks regions:
                // 1) circle around the saturated stars (scaled by magnitude)
                // 2) diffraction spikes in direction ROT - ROTo
                // 3) bleed trail in the direction of the readout
                //    update: with 'burntool' applied to the data, the bleed trail mask is not needed.

                for (int i = 0; i < refstars->n; i++) {
                    pmAstromObj *ref = refstars->data[i];

                    if (COUNT_GHOSTS) {
                        if (ref->Mag > GHOST_MAX_MAG) {
                            nGhosts ++;
                        }
                    }

                    if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue;


                    // the reference magnitudes have been converted from the instrumental
                    // values supplied in the recipe to apparent mags (above)

                    // CIRCLE around the stars (scaled by magnitude)
                    float radius = REFSTAR_MASK_SATSTAR_MAG_SLOPE * (REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag);

                    // XXX for now, assume cell binning is 1x1 relative to chip
                    psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius);

                    for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) {

                        float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO;

                        // LINE for boundaries of the saturation spikes (scaled by magnitude)
                        // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;

                        float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
			float spikeLength = pow(10,REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag)) - REFSTAR_MASK_SATSPIKE_OFFSET;
                        float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
			if (spikeLength < 0.0) {
			  spikeLength = 0.0;
			}
                        // XXX we can make the width depend on the spike as well...
                        // The length should also be a function of the image background level

                        psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
                    }

                    // This masking option was needed for persistent charge trails in GPC1; it
                    // has since been replaced with 'burntool', which is applied upon readout
                    // by the camera software, and therefore is aware of the image sequence.
                    if (REFSTAR_MASK_BLEED) {
                        // convert x,y chip coordinates to cells in maskChip
                        pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);

                        // LINE for boundaries of the bleed lines
                        if (refCell) {
                            float xCell = 0.0;
                            float yCell = 0.0;
                            float xEnd = 0.0;
                            float yEnd = 0.0;
                            // find coordinate of star on cell
                            pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
                            // find coordinate of end-point on chip

                            int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
                            pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);

                            float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
                            psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
                        }
                    }
                }

                // select the ghost object for this readout (loaded in psastroExtractGhosts.c).
                // These differ from the reference stars since the star position is not
                // contained by the readout; instead, the ghost position is predicted based on
                // the ghost model, and the ghost positions associated with a given readout are
                // supplied here.
                psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS");
                if (ghosts) {
                    // mask the ghosts on this readout
                    for (int i = 0; i < ghosts->n; i++) {
                        psastroGhost *ghost = ghosts->data[i];
                        // XXX bright vs faint ghost bits? (OR with SUSPECT)
                        psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
                    }
                }

                // Select the glint mask regions for this readout (loaded in
                // psastroExtractGlints.c).  These glint regions are defined as rectangular
                // boxes and are generated for each chip based on the position of the bright
                // stars beyond the edge of the focal plane.  This masking is currently very
                // GPC1-specific
                psArray *glints = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GLINTS");
                if (glints) {
                    // mask the glints on this readout
                    for (int i = 0; i < glints->n; i++) {
                        psRegion *glint = glints->data[i];
                        psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
                    }
                }

                // Select the crosstalk artifact regions for this readout, and mask a circular region
                // corresponding to the source star

                psArray *crosstalks = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.CROSSTALKS");
                if (crosstalks) {
                  for (int i = 0; i < crosstalks->n; i++) {
                    pmAstromObj *ref = crosstalks->data[i];
                    float radius = REFSTAR_MASK_CROSSTALK_MAG_SLOPE * (REFSTAR_MASK_CROSSTALK_MAG_MAX - ref->Mag);
                    psTrace("psastro.crosstalk",2,"Masking star on Chip %s @ (%f,%f) Magnitude: %f Radius %f\n",
                            psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"),
                            ref->chip->x,ref->chip->y,ref->Mag,radius);
                    // XXX for now, assume cell binning is 1x1 relative to chip
                    psastroMaskCircle (readoutMask->mask, crosstalkMaskValue, ref->chip->x, ref->chip->y, radius, radius);
                  }
                }

                // this probably should move into a function of its own:
                {
                    // select the raw objects for this readout, flag is they fall in a mask
                    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
                    if (!detections) continue;

                    psArray *inSources = detections->allSources;
                    psAssert (inSources, "missing sources?");

                    // create a replacement output array:
                    // psArray *outSources = psAllocArrayEmpty(100);

                    // XXX finish this: raise a bit for stars that land on certain types of masks;
                    // others (eg, bright star core) should be ignored.
                    for (int i = 0; i < inSources->n; i++) {
                        pmSource *source = inSources->data[i];

                        int xChip = source->peak->x;
                        int yChip = source->peak->y;

                        bool onChip = true;
                        onChip &= (xChip >= 0);
                        onChip &= (xChip < readoutMask->mask->numCols);
                        onChip &= (yChip >= 0);
                        onChip &= (yChip < readoutMask->mask->numRows);
                        if (!onChip) {
                            // if the source is off the edge of the chip, raise a different bit?
                            source->mode |= PM_SOURCE_MODE_OFF_CHIP;
                            continue;
                        }

                        psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[yChip][xChip];
                        if (value & ghostMaskValue) {
                            source->mode |= PM_SOURCE_MODE_ON_GHOST;
                        }
                        // XXX note that for now, glint and ghost are identical
                        pmSourceMode PM_SOURCE_MODE_ON_GLINT = PM_SOURCE_MODE_ON_GHOST;
                        if (value & glintMaskValue) {
                            source->mode |= PM_SOURCE_MODE_ON_GLINT;
                        }
                        if (value & spikeMaskValue) {
                            source->mode |= PM_SOURCE_MODE_ON_SPIKE;
                        }
                    }
                }
            }
        }

        // output sequence for mask corresponding to this chip (XXX this may not be needed...)
        *viewMask = *view;
        while ((cell = pmFPAviewNextCell (viewMask, outMask->fpa, 1)) != NULL) {
            psTrace ("psastro", 4, "Mask Cell %d: %x %x\n", viewMask->cell, cell->file_exists, cell->process);
            if (!cell->process || !cell->file_exists) { continue; }

            while ((readout = pmFPAviewNextReadout (viewMask, outMask->fpa, 1)) != NULL) {
                if (! readout->data_exists) { continue; }
                if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE;
            }
            if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE;
        }
        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    }
    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;

    if (COUNT_GHOSTS) {
        // save nGhosts to update header.
        psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
        if (!updates) {
            updates = psMetadataAlloc ();
            psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
            psFree (updates);
        }
        psMetadataAddS32 (updates, PS_LIST_TAIL, "NGHOSTS", PS_META_REPLACE, "total expected ghosts", nGhosts);
    }

    // deactivate all files
    pmFPAfileActivate (config->files, false, NULL);

    psFree (view);
    psFree (viewMask);
    return true;
}

