# include "psphotInternal.h"

// we have 3 possible real filesets:
# define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */

// we have 3 files on which we operate:
// DET (detection image)       : nominally RAW (optionally CNV?)
// SRC (source analysis image) : nominally CNV (optionally RAW)
// OUT (psf-matched images)    : always OUT

bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {

    bool status = false;

    int num = psphotFileruleCount(config, filerule);

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);

    // loop over the available readouts
    for (int i = 0; i < num; i++) {

        // find the currently selected readout
        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
        psAssert (file, "missing file?");

        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
        psAssert (readout, "missing readout?");

        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
        psAssert (detections, "missing detections?");

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

	psphotVisualShowResidualImage (readout, true);
	psphotVisualShowObjectRegions (readout, recipe, sources);
    }
    return true;
}

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

    psArray *objects = NULL; // used below after 'pass1finish' label

    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
    // by the multiple threads, not the total time used by all threads.
    psTimerStart ("psphotReadout");

    pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); // allow models to have ugly fits (eg, central cusp)

    // 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;
    }
    // optional break-point for processing
    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    psAssert (breakPt, "configuration error: set BREAK_POINT");

    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
    char *STACK_DET = STACK_RAW;

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

    // Generate the mask and weight images (if not supplied) and set mask bits
    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!strcasecmp (breakPt, "NOTHING")) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // generate a background model (median, smoothed image)
    if (!psphotModelBackground (config, view, STACK_DET)) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!psphotSubtractBackground (config, view, STACK_DET)) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!psphotSubtractBackground (config, view, STACK_SRC)) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!strcasecmp (breakPt, "BACKMDL")) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // also make the chisq detection image
    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
        return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!strcasecmp (breakPt, "CHISQ")) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // find the detections (by peak and/or footprint) in the image.
    // This finds the detections on Chisq image as well as the individuals
    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
        // this only happens if we had an error in psphotFindDetections
        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
        return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
    // is just a copy of the container pointer; the sources on both DET and SRC are the same
    // memory objects
    if (strcmp(STACK_SRC, STACK_DET)) {
	if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
	    psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
	    return psphotReadoutCleanup (config, view, STACK_SRC);
	}
    }

    // construct sources and measure basic stats (saved on detections->newSources)
    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
        return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!strcasecmp (breakPt, "PEAKS")) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    // psphotDumpTest (config, view, STACK_SRC);
    psMemDump("sourcestats");

    // classify sources based on moments, brightness
    // only run this on detections from the input images, not chisq image
    if (!psphotRoughClass (config, view, STACK_SRC)) {
        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    // only run this on detections from the input images, not chisq image
    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
        return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!strcasecmp (breakPt, "MOMENTS")) {
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // use bright stellar objects to measure PSF
    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
        psLogMsg ("psphot", 3, "failure to construct a psf model");
        return psphotReadoutCleanup (config, view, STACK_SRC);
    }
    if (!strcasecmp (breakPt, "PSFMODEL")) {
        return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // merge the newly selected sources into the existing list
    // NOTE: merge OLD and NEW
    psphotMergeSources (config, view, STACK_SRC);

    // Construct an initial model for each object, set the radius to fitRadius, set circular
    // fit mask.  NOTE: only applied to sources without guess models
    psphotGuessModels (config, view, STACK_SRC);

    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    psphotFitSourcesLinear (config, view, STACK_SRC, false);
    psphotStackVisualFilerule(config, view, STACK_SRC);

    // measure the radial profiles to the sky
    psphotRadialProfileWings (config, view, STACK_SRC);

    // re-measure the kron mags with models subtracted.  this pass starts with a circular
    // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
    // but iterates to an appropriately larger size
    psphotKronIterate(config, view, STACK_SRC);

    // identify CRs and extended sources
    psphotSourceSize (config, view, STACK_SRC, true);

    // non-linear PSF and EXT fit to brighter sources
    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)

    // replace all sources (do NOT ignore subtraction state)
    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)

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

    // linear fit to include all sources (subtract again)
    // NOTE : apply to ALL sources (extended + psf)
    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)

    // NOTE: possibly re-measure background model here with objects subtracted / or masked

    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
    { 
	// if DET and SRC are different images, generate children sources for all sources in
	// the SRC image.  This operation replaces the existing DETECTION container on DET
	// which is currently a view to the one on SRC).  children sources go to
	// det->allSources
	if (strcmp(STACK_SRC, STACK_DET)) {
	    psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 

	    //  subtract all sources from DET (this will subtract using the psf model for SRC, which
	    //  will somewhat oversubtract the sources -- this is OK
	    psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
	}

	// add noise for subtracted objects
	psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)

	// find fainter sources
	// NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
	psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)

	// remove noise for subtracted objects (ie, return to normal noise level)
	// NOTE: this needs to operate only on the OLD sources
	psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)

	// if DET and SRC are different images, copy the detections from DET to SRC 
	// (this operation just ensures the metadata container has a view on SRC as well
	if (strcmp(STACK_SRC, STACK_DET)) {
	    // replace all sources in DET
	    psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources

	    // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
	    if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
		psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
		return psphotReadoutCleanup (config, view, STACK_SRC);
	    }
	}

	// define new sources based on only the new peaks & measure moments
	// NOTE: new sources are saved on detections->newSources
	psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)

	// set source type
	// NOTE: apply only to detections->newSources
	if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
	    psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
	    return psphotReadoutCleanup (config, view, STACK_SRC);
	}

	// replace all sources so fit below applies to all at once
	// NOTE: apply only to OLD sources (which have been subtracted)
	psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2

	// merge the newly selected sources into the existing list
	// NOTE: merge OLD and NEW
	// XXX check on free of sources...
	psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)

	// Construct an initial model for each object, set the radius to fitRadius, set circular
	// fit mask.  NOTE: only applied to sources without guess models
	psphotGuessModels (config, view, STACK_SRC);
    }

pass1finish:

    // generate the objects (objects unify the sources from the different images) NOTE: could
    // this just match the detections for the chisq image, and not bother measuring the source
    // stats in that case...?
    objects = psphotMatchSources (config, view, STACK_SRC);
    psMemDump("matchsources");

    // Construct an initial model for each object, set the radius to fitRadius, set circular
    // fit mask.  NOTE: only applied to sources without guess models
    psphotGuessModels (config, view, STACK_SRC);

    psphotStackObjectsUnifyPosition (objects);

    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);

    // NOTE: apply to ALL sources
    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)

    // measure the radial profiles to the sky (only measures new objects)
    psphotRadialProfileWings (config, view, STACK_SRC);

    // re-measure the kron mags with models subtracted
    // psphotKronMasked(config, view, STACK_SRC);
    psphotKronIterate(config, view, STACK_SRC);

    // measure source size for the remaining sources
    // NOTE: applies only to NEW (unmeasured) sources
    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)

    psMemDump("psfstats");

    // measure elliptical apertures, petrosians (objects sorted by S/N)
    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)

    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)

    // create source children for the OUT filerule (for radial aperture photometry) 
    psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
    if (!objectsRadial) {
	psFree(objects);
	psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // measure circular, radial apertures (objects sorted by S/N)
    // this forces photometry on the undetected sources from other images
    psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
    psMemDump("extmeas");

    int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);

    for (int entry = 1; entry < nRadialEntries; entry++) {
	// NOTE: entry 0 is the unmatched image set

	// re-measure the PSF for the smoothed image (using entries in 'allSources')
	psphotChoosePSF (config, view, STACK_OUT, false);

	// this is necessary to update the models based on the new PSF
	psphotResetModels (config, view, STACK_OUT);

	// this is necessary to get the right normalization for the new models
	psphotFitSourcesLinear (config, view, STACK_OUT, false);

	// measure circular, radial apertures (objects sorted by S/N)
	// entry 0 == unmatched? pass entry + 1?
	psphotRadialApertures (config, view, STACK_OUT, entry); 

	// replace the flux in the image so it is returned to its original state
	psphotReplaceAllSources (config, view, STACK_OUT, false);

	// smooth to the next FWHM, or set 'smoothAgain' to false if no more 
	psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
	psMemDump("matched");
    }

    // measure aperture photometry corrections
    if (!psphotApResid (config, view, STACK_SRC)) {
	psFree (objects);
        psLogMsg ("psphot", 3, "failed on psphotApResid");
	return psphotReadoutCleanup (config, view, STACK_SRC);
    }

    // calculate source magnitudes
    psphotMagnitudes(config, view, STACK_SRC);

    // XXX NOTE: this function wants to have the PSF of the image, but we (so far) only measure the 
    // PSF of the SRC image.  can we fake it by generating the PSF for DET as well (up above)? 
    if (false && !psphotEfficiency(config, view, STACK_DET)) {
        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
        psErrorClear();
    }

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

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

    // drop the references to the image pixels held by each source
    psphotSourceFreePixels (config, view, STACK_OUT);
    psphotSourceFreePixels (config, view, STACK_SRC);

    // remove chisq image from config->file:PSPHOT.INPUT
    psphotStackRemoveChisqFromInputs(config, STACK_DET);
    if (strcmp(STACK_SRC, STACK_DET)) {
	psphotStackRemoveChisqFromInputs(config, STACK_SRC);
    }

    psFree (objects);
    psFree (objectsRadial);

    // create the exported-metadata and free local data
    return psphotReadoutCleanup (config, view, STACK_SRC);
}

/* here is the process:

 * we have three(*) images:
 * RAW : unconvolved image stack 
 * CNV : input convolved image

 * OUT : psf-matched output image (there may be more than one of
 * these.  we will generate the first matched image by selecting the
 * target PSF and doing a full psf-maching process (as used by ppStack
 * and ppSub).  But, additional target output files should use a
 * simple gaussian convolution kernel determind from therms of the
 * current and the target).

 * the output should be / could be one of the matched images, but not
 * all.  should we ensure the first gets written out, and ot save the
 * others (or only optionally).

 * by default, we probably only sve the cmf ffile outputs.

 * load the RAW image (unconvolved stacks)
 * add photcode to the output headers / readout->analysis
 * generate mask and variance image (this is probably never needed in
   practice: we always load an input mask & var.
 * generate & subtract a model background for ?? (RAW? CNV? OUT? all?)
 * load a PSF (probably not yet working)

 * generate the CHISQ image from the RAW input images (why save on OUT?)

 * find detections on RAW

 * copy detections to OUT

 * generate source stats (moments) for OUT

 * match sources across inputs (on OUT?)

 * generate source stats for the new constructions

 * rough class (star, galaxy, cosmic, etc)

 * Image quality

 * generate PSF

 * guess models

 * merge sources (new -> old)

 * linear fit to the psf

 * find ApResid

 * assign common positions

 * radial apertures (** this should be on the PSF-matched images

 * extended analysis (elliptical profile & petrosian)

 * extended fits (sersic, etc)

 * psphot magnitudes


 ******

 the above is all wrong:  first, we should be doing the full
 morphology analysis (ExtendedAnalysis & ExtendedFits) on the CNV or
 RAW image (as desired optionally), etc.

 In the discussion below, 'BST' (best) means optionally RAW or CNV

 * detection : RAW & CHISQ (of RAW)
 * moments : used by psf analysis & classification (BST)
 * rough class : uses moments, not pixels
 * image quality : uses moments as well
 * generate PSF : (BST)
 * guess models (BST)
 * linear fit (BST)
 * find ApResid (BST) -- uses sources not pixels
 * extended analysis (BST)
 * extended fits (BST)
 * detection efficiency (BST)
 
 * somehow need to copy the sources so they point at the pixels on the
 * OUT image 

 * foreach target PSF
   * radial aperture
   * convolve to next target PSF

   * somehow need to organize the output file to have the values from
   * the different PSFs in separate tables (with header info to
   * specify the size of that PSF)

   */
