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

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

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

    // 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 (readout, sources, recipe, 0);

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