/** @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;

    psImageMaskType ghostMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels
    psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels
    psImageMaskType starMaskValue  = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels

    // 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 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");
    double GHOST_MAX_MAG = psMetadataLookupF32 (&status, recipe, "GHOST_MAX_MAG");
    int nGhosts = 0;

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

    // 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, 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
                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 = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag);
			float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
			// 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);
                    }

		    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 raw objects for this readout (loaded in psastroExtract.c)
                psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS");
                if (ghosts == NULL) { continue; }

                // identify the bright stars of interest
                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 raw objects for this readout, flag is they fall in a mask
		psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
		if (rawstars == NULL) return false;
		
		// XXX finish this:
		for (int i = 0; i < rawstars->n; i++) {
		    pmAstromObj *raw = rawstars->data[i];
		    psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[(int)(raw->chip->x)][(int)(raw->chip->y)];
		    if (value) continue;
		}
            }
        }

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

