# include "psphotInternal.h"

bool psphotReadout(pmConfig *config, const pmFPAview *view) {

    bool dump = (psTraceGetLevel("psphot") >= 6);

    // select the current recipe
    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    if (!recipe) {
        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
        return false;
    }

    // Interpret the mask values
    const char *maskValStr = psMetadataLookupStr(NULL, recipe, "MASKVAL"); // String with mask names
    if (!maskValStr) {
        psError(PSPHOT_ERR_CONFIG, false, "Missing recipe item: MASKVAL(STR)");
        return false;
    }
    psMaskType maskVal = pmConfigMask(maskValStr, config); // Mask values to mask against
    psMaskType maskMark = pmConfigMask("MARK", config); // Mask value for marking
    psMaskType maskSat = pmConfigMask("SAT", config); // Mask value for saturated pixels
    psMaskType maskBad = pmConfigMask("BAD", config); // Mask value for bad pixels

    // find the currently selected readout
    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    PS_ASSERT_PTR_NON_NULL (readout, false);

    // optional break-point for processing
    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    PS_ASSERT_PTR_NON_NULL (breakPt, false);

    // Use the new pmFootprints approach?
    const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS");

    // generate mask & weight images if they don't already exit
    if (!readout->mask) {
        if (!pmReadoutGenerateMask(readout, maskSat, maskBad)) {
            psError (PSPHOT_ERR_CONFIG, false, "trouble creating mask");
            return false;
        }
    }
    if (!readout->weight) {
        if (!pmReadoutGenerateWeight(readout, true)) {
            psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight");
            return false;
        }
    }

    // set the photcode for this image
    if (!psphotAddPhotcode (recipe, config, view)) {
        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
        return false;
    }

    // I have a valid mask, now mask in the analysis region of interest
    psphotMaskReadout (readout, recipe, maskBad);

    if (psTraceGetLevel("psphot") >= 5) {
        psphotSaveImage (NULL, readout->image,  "image.fits");
        psphotSaveImage (NULL, readout->mask,   "mask.fits");
        psphotSaveImage (NULL, readout->weight, "weight.fits");
    }

    if (!strcasecmp (breakPt, "NOTHING")) {
        return psphotReadoutCleanup(config, readout, recipe, NULL, NULL);
    }

    // generate a background model (median, smoothed image)
    if (!psphotImageMedian (config, view, maskVal)) {
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL);
    }

    if (!strcasecmp (breakPt, "BACKMDL")) {
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL);
    }

    // run a single-model test if desired
    psphotModelTest (config, view, recipe, maskVal, maskMark);

    // find the peaks in the image
    psArray *peaks;
    psArray *footprints = NULL;
    if (useFootprints) {
       footprints = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal);
       int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS");
       if (growRadius > 0) {
           psArray *tmp = pmGrowFootprintArray(footprints, growRadius);
           psFree(footprints);
           footprints = tmp;
       }
       psphotCullPeaks(readout->image, readout->weight, recipe, footprints);

       peaks = pmFootprintArrayToPeaks(footprints);
    } else {
       peaks = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal);
    }

    if (!peaks) {
        psLogMsg ("psphot", 3, "unable to find peaks in this image");
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL);
    }

    // psArray *bounds = psphotGetBounds ();
    // psArray *sources = psphotBoundsToSources ();

    // construct sources and measure basic stats
    psArray *sources = psphotSourceStats (readout, recipe, peaks, maskVal, maskMark);
    if (!sources) return false;
    psFree (peaks);

    if (!strcasecmp (breakPt, "PEAKS")) {
        return psphotReadoutCleanup(config, readout, recipe, NULL, sources);
    }

    psphotDeblendSatstars (sources, recipe);

    // mark blended peaks PS_SOURCE_BLEND
    if (!psphotBasicDeblend (sources, recipe)) {
        psLogMsg ("psphot", 3, "failed on deblend analysis");
        return psphotReadoutCleanup (config, readout, recipe, NULL, sources);
    }

    // classify sources based on moments, brightness
    if (!psphotRoughClass (sources, recipe, true, maskSat)) {
        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
        return psphotReadoutCleanup (config, readout, recipe, NULL, sources);
    }
    if (!strcasecmp (breakPt, "MOMENTS")) {
        return psphotReadoutCleanup(config, readout, recipe, NULL, sources);
    }

    // load the psf model, if suppled
    pmPSF *psf = psphotLoadPSF (config, view, recipe);
    if (!psf) {
        // use bright stellar objects to measure PSF
        psf = psphotChoosePSF (readout, sources, recipe, maskVal, maskMark);
        if (psf == NULL) {
            psLogMsg ("psphot", 3, "failure to construct a psf model");
            return psphotReadoutCleanup (config, readout, recipe, psf, sources);
        }
    }

    // XXX Test dump of PSF parameters across image field.

    if (!strcasecmp (breakPt, "PSFMODEL")) {
        return psphotReadoutCleanup (config, readout, recipe, psf, sources);
    }

    // include externally-supplied sources
    psphotLoadExtSources (config, view, sources);

    // construct an initial model for each object
    psphotGuessModels (readout, sources, recipe, psf, maskVal);

    // XXX test output of models
    // psphotTestSourceOutput (readout, sources, recipe, psf);

    if (dump) psphotSaveImage (NULL, readout->image,  "image.v0.fits");

    // linear PSF fit to peaks
    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE, maskVal);
    if (dump) psphotSaveImage (NULL, readout->image,  "image.v1.fits");

    // XXX test the CR/EXT measurement
    // psphotSourceSize (readout, sources, recipe);

    if (!strcasecmp (breakPt, "ENSEMBLE")) {
        goto finish;
    }

    // non-linear PSF and EXT fit to brighter sources
    psphotBlendFit (readout, sources, recipe, psf, maskVal);
    if (dump) psphotSaveImage (NULL, readout->image,  "image.v2.fits");

    // replace all sources
    psphotReplaceAll (sources, maskVal);
    if (dump) psphotSaveImage (NULL, readout->image,  "image.v3.fits");

    // linear PSF fit to remaining peaks
    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal);
    if (dump) psphotSaveImage (NULL, readout->image,  "image.v4.fits");
    if (!strcasecmp (breakPt, "PASS1")) {
        goto finish;
    }

    // XXX for the moment, drop the re-calc of the background (prove this works)
    // replace background in residual image
    // psphotSkyReplace (config, view);
    // re-measure background model (median, smoothed image)
    // psphotImageMedian (config, view);

    // add noise for subtracted objects
    psphotAddNoise (readout, sources, recipe, true, maskVal);

    // find the peaks in the image
    psArray *newPeaks;
    if (useFootprints) {
       psArray *newFootprints = psphotFindPeaks (readout, recipe, useFootprints, 2, maskVal);

       int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS_2");
       if (growRadius > 0) {
           psArray *tmp = pmGrowFootprintArray(newFootprints, growRadius);
           psFree(newFootprints);
           newFootprints = tmp;
       }

       // merge in old peaks
       const int includePeaks = 0x0 | 0x2; // i.e. just from newFootprints
       psArray *mergedFootprints = pmMergeFootprintArrays(footprints, newFootprints, includePeaks);
       psFree(footprints);
       psFree(newFootprints);

       psphotCullPeaks(readout->image, readout->weight, recipe, mergedFootprints);

       newPeaks = pmFootprintArrayToPeaks(mergedFootprints);

       psFree(mergedFootprints);
    } else {
       newPeaks = psphotFindPeaks(readout, recipe, useFootprints, 2, maskVal);
    }

    // remove noise for subtracted objects
    psphotAddNoise (readout, sources, recipe, false, maskVal);

    // define new sources based on the new peaks
    psArray *newSources = psphotSourceStats (readout, recipe, newPeaks, maskVal, maskMark);
    psFree (newPeaks);

    // set source type
    psphotRoughClass (newSources, recipe, false, maskSat);

    // create full input models
    psphotGuessModels (readout, newSources, recipe, psf, maskVal);

    // replace all sources
    psphotReplaceAll (sources, maskVal);
    if (dump) psphotSaveImage (NULL, readout->image,  "image.v5.fits");

    // merge the newly selected peaks into the existing list
    psphotMergeSources (sources, newSources);
    psFree (newSources);

    // linear PSF fit to remaining peaks
    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal);
    if (dump) psphotSaveImage (NULL, readout->image,  "image.v6.fits");

    psphotExtendedSources (readout, sources, recipe, maskVal);

finish:

    // plot positive sources
    // psphotSourcePlots (readout, sources, recipe);

    // measure aperture photometry corrections
    if (!psphotApResid (readout, sources, recipe, psf, maskVal, maskMark)) {
        psTrace ("psphot", 4, "failure on psphotApResid");
        psError(PSPHOT_ERR_PHOTOM, false, "Measure aperture photometry corrections");
        return false;
    }

    // calculate source magnitudes
    pmReadout *background = psphotSelectBackground (config, view, false);
    psphotMagnitudes(sources, recipe, psf, background, maskVal, maskMark);

    // replace failed sources?
    // psphotReplaceUnfit (sources);

    // replace background in residual image
    psphotSkyReplace (config, view);

    // drop the references to the image pixels held by each source
    psphotSourceFreePixels (sources);

    // create the exported-metadata and free local data
    return psphotReadoutCleanup(config, readout, recipe, psf, sources);
}
