# include "psphotInternal.h"

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

    // measure the total elapsed time in psphotReadout.  XXX the current threading plan
    // for psphot envisions threading within psphotReadout, not multiple threads calling
    // the same psphotReadout.  In the current plan, this dtime is the elapsed time used
    // jointly by the multiple threads, not the total time used by all threads.
    psTimerStart ("psphotReadout");

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

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

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

    // Generate the mask and weight images, including the user-defined analysis region of interest
    psphotSetMaskAndWeight (config, readout, recipe);
    if (!strcasecmp (breakPt, "NOTHING")) {
        return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);
    }

    // generate a background model (median, smoothed image)
    if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) {
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
    }
    if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) {
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
    }
    if (!strcasecmp (breakPt, "BACKMDL")) {
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
    }

    // run a single-model test if desired (exits from here if test is run)
    psphotModelTest (config, view, recipe);

    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    pmPSF *psf = psphotLoadPSF (config, view, recipe);

    // several functions below behave differently if we have a PSF model already
    bool havePSF = (psf != NULL);

    // find the detections (by peak and/or footprint) in the image.
    pmDetections *detections = psphotFindDetections (NULL, readout, recipe);
    if (!detections) {
        psLogMsg ("psphot", 3, "unable to find detections in this image");
        return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
    }
 
    // XXX test write out the footprint image
    if (0) {
	psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true);
	psphotSaveImage (NULL, footprintImage, "footprints.1.fits");
	psFree (footprintImage);
    }

    // XXX test write out the footprint image
    if (0) {
        psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true);
        psphotSaveImage (NULL, footprintImage, "footprints.1.fits");
        psFree (footprintImage);
    }

    // construct sources and measure basic stats
    psArray *sources = psphotSourceStats (readout, recipe, detections);
    if (!sources) return false;
    if (!strcasecmp (breakPt, "PEAKS")) {
        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    }

    // find blended neighbors of very saturated stars
    // XXX merge this with Basic Deblend?
    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, detections, psf, sources);
    }

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

    // if we were not supplied a PSF, choose one here
    if (psf == NULL) {
        // use bright stellar objects to measure PSF
        // XXX if we do not have enough stars to generate the PSF, build one
        // from the SEEING guess and model class
        psf = psphotChoosePSF (readout, sources, recipe);
        if (psf == NULL) {
            psLogMsg ("psphot", 3, "failure to construct a psf model");
            return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
        }
        havePSF = true;
    }
    if (!strcasecmp (breakPt, "PSFMODEL")) {
        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    }

    // Define source fitting parameters for everything that follows PSF fits
    // This allows different parameters to be used from the PSF fitting
    {
        bool status = true;             // Status of MD lookup
        int fitIter = psMetadataLookupS32(&status, recipe, "EXT_FIT_ITER"); // Max number of fit iterations
        if (!status || fitIter <= 0) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "EXT_FIT_ITER is not positive");
            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
        }
        float fitTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_TOL"); // Fit tolerance
        if (!status || !isfinite(fitTol) || fitTol <= 0) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "EXT_FIT_TOL is not positive");
            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
        }
        bool poisson = psMetadataLookupBool(&status, recipe, "POISSON.ERRORS.PHOT.LMM"); // Poisson errors?
        if (!status) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "POISSON.ERRORS.PHOT.LMM is not defined");
            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
        }
        float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
        if (!status || !isfinite(skySig) || skySig <= 0) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "SKY_SIG is not positive");
            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
        }
        pmSourceFitModelInit(fitIter, fitTol, PS_SQR(skySig), poisson);
    }

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

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

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

    // linear PSF fit to source peaks
    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);

    if (0) {
        FILE *out = fopen ("out.pass1.dat", "w");

        for (int i = 0; i < sources->n; i++) {
            pmSource *source = sources->data[i];
            if (!source) continue;
            pmPeak *peak = source->peak;
            if (!peak) continue;
            pmModel *model = source->modelPSF;
            if (!model) continue;
            if (!model->params) continue;

            fprintf (out, "%d %f %f  %f %f  %f %f\n",
                     i, peak->xf, peak->yf, peak->flux, peak->SN,
                     model->params->data.F32[PM_PAR_I0],
                     model->dparams->data.F32[PM_PAR_I0]);
        }
        fclose (out);
    }

    // identify CRs and extended sources
    psphotSourceSize (config, readout, sources, recipe, 0);
    if (!strcasecmp (breakPt, "ENSEMBLE")) {
        goto finish;
    }

    // non-linear PSF and EXT fit to brighter sources
    psphotBlendFit (readout, sources, recipe, psf);

    // replace all sources
    psphotReplaceAllSources (sources, recipe);

    // linear fit to include all sources
    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);

    // if we only do one pass, skip to extended source analysis
    if (!strcasecmp (breakPt, "PASS1")) {
        goto pass1finish;
    }

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

    // find fainter sources (pass 2)
    detections = psphotFindDetections (detections, readout, recipe);

    // XXX test write out the footprint image
    if (0) {
        psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true);
        psphotSaveImage (NULL, footprintImage, "footprints.2.fits");
        psFree (footprintImage);
    }
    // remove noise for subtracted objects (ie, return to normal noise level)
    psphotSubNoise (readout, sources, recipe);

    // define new sources based on only the new peaks
    psArray *newSources = psphotSourceStats (readout, recipe, detections);

    // set source type
    if (!psphotRoughClass (newSources, recipe, havePSF)) {
        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    }

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

    // replace all sources so fit below applies to all at once
    psphotReplaceAllSources (sources, recipe);

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

    // linear fit to all sources
    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);

    if (0) {
        FILE *out = fopen ("out.pass2.dat", "w");

        for (int i = 0; i < sources->n; i++) {
            pmSource *source = sources->data[i];
            if (!source) continue;
            pmPeak *peak = source->peak;
            if (!peak) continue;
            pmModel *model = source->modelPSF;
            if (!model) continue;
            if (!model->params) continue;

            fprintf (out, "%d %f %f  %f %f  %f %f\n",
                     i, peak->xf, peak->yf, peak->flux, peak->SN,
                     model->params->data.F32[PM_PAR_I0],
                     model->dparams->data.F32[PM_PAR_I0]);
        }
        fclose (out);
    }

pass1finish:

    // measure source size for the remaining sources
    psphotSourceSize (config, readout, sources, recipe, 0);
    if (0) {
        psphotSaveImage (NULL, readout->mask, "mask.fits");
    }
	
    psphotExtendedSourceAnalysis (readout, sources, recipe);

    psphotExtendedSourceFits (readout, sources, recipe);

finish:

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

    // measure aperture photometry corrections
    if (!psphotApResid (readout, sources, recipe, psf)) {
        psLogMsg ("psphot", 3, "failed on psphotApResid");
        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    }

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

    // replace failed sources?
    // psphotReplaceUnfitSources (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, detections, psf, sources);
}
