# include "psphotInternal.h"

// for now, let's store the detections on the readout->analysis for each readout
bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view)
{
    bool status = true;

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

    // perform full non-linear fits / extended source analysis?
    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) {
	psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
	return true;
    }

    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");

    // loop over the available readouts
    for (int i = 0; i < num; i++) {
	if (!psphotExtendedSourceAnalysisReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
            psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for PSPHOT.INPUT entry %d", i);
	    return false;
	}
    }
    return true;
}

// aperture-like measurements for extended sources
bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {

    bool status;
    int Next = 0;
    int Npetro = 0;
    int Nisophot = 0;
    int Nannuli = 0;
    int Nkron = 0;

    psTimerStart ("psphot.extended");

    // find the currently selected readout
    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // 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?");

    if (!sources->n) {
	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
	return true;
    }

    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    assert (maskVal);

    // XXX require petrosian analysis for non-linear fits? 

    // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");

    // S/N limit to perform full non-linear fits
    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");

    // which extended source analyses should we perform?
    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");

# if (0)
    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
    // XXX use this to set skynoise
    pmReadout *backModel = psphotSelectBackground (config, view);
    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
# endif

    // source analysis is done in S/N order (brightest first)
    sources = psArraySort (sources, pmSourceSortBySN);

    // option to limit analysis to a specific region
    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");

    // choose the sources of interest
    for (int i = 0; i < sources->n; i++) {

	pmSource *source = sources->data[i];

	// skip PSF-like and non-astronomical objects
	if (source->type == PM_SOURCE_TYPE_STAR) continue;
	if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
	if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
	if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
	if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;

	// limit selection to some SN limit
	assert (source->peak); // how can a source not have a peak?
	if (source->peak->SN < SN_LIM) continue;

	// limit selection by analysis region
	if (source->peak->x < AnalysisRegion.x0) continue;
	if (source->peak->y < AnalysisRegion.y0) continue;
	if (source->peak->x > AnalysisRegion.x1) continue;
	if (source->peak->y > AnalysisRegion.y1) continue;

	// replace object in image
	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
	}
	Next ++;

	// force source image to be a bit larger...
	float radius = source->peak->xf - source->pixels->col0;
	radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
	radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
	radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);

	// if we request any of these measurements, we require the radial profile
	if (doPetrosian || doIsophotal || doAnnuli || doKron) {
	    if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
		// all measurements below require the radial profile; skip them all
		// re-subtract the object, leave local sky
		psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
		pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
		continue;
	    }
	    source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
	}

	// Petrosian Mags
	if (doPetrosian) {
	    if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
		psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
	    } else {
		psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
		Npetro ++;
		source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
	    }
	}

# if (0)
	// Isophotal Mags
	if (doIsophotal) {
	    if (!psphotIsophotal (source, recipe, maskVal)) {
		psTrace ("psphot", 5, "failed to measure isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
	    } else {
		psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
		Nisophot ++;
		source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
	    }
	}
	// Kron Mags
	if (doKron) {
	    if (!psphotKron (source, recipe, maskVal)) {
		psTrace ("psphot", 5, "failed to measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
	    } else {
		psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
		Nkron ++;
		source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
	    }
	}
# endif

	// re-subtract the object, leave local sky
	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);

	if (source->extpars) {
	    psFree(source->extpars->radFlux);
	    psFree(source->extpars->ellipticalFlux);
	}
    }

    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
    psLogMsg ("psphot", PS_LOG_INFO, "  %d isophotal\n", Nisophot);
    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
    psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);

    psphotVisualShowResidualImage (readout);

    if (doPetrosian) {
	psphotVisualShowPetrosians (sources);
    }

    return true;
}
