# include "psphot.h"

bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe) {

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

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

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

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

    // 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 doPSFConvModel = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PSF_CONVOLVED_MODEL");
    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");

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

    // XXX some init functions for the extended source recipe options?

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

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

        // XXX this should use peak?
        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;

        // if model is NULL, we don't have a starting guess
	// XXX this is probably not needed for objects where we only measure aperture-like terms
        if (source->modelEXT == NULL) continue;

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

	if (doPSFConvModel && !psphotPSFConvModel (source, recipe, maskVal)) {
	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in PSF Convolved Model fit");
	    return false;
	} else {
	  // XXX why am I caching the model?
	    pmSourceCacheModel (source, maskVal); // XXX put this in the source model function?
	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
            Nconvolve ++;
	}

	// all of these below require the radial profile
	// XXX push this as a test and call in each of the functions below?
	// XXX this constructs a pmSourceExtendedParameters element
	if (doPetrosian || doIsophotal || doAnnuli || doKron) {
	  if (!psphotRadialProfile (source, recipe, maskVal)) {
	    psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate radial profile");
	    return false;
	  }
	}

	if (doPetrosian && !psphotPetrosian (source, recipe, maskVal)) {
	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Petrosian analysis");
	    return false;
	} else {
	    psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
            Npetro ++;
	}

	if (doIsophotal && !psphotIsophotal (source, recipe, maskVal)) {
	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Isophotal analysis");
	    return false;
	} else {
	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
            Nisophot ++;
	}

	if (doAnnuli && !psphotAnnuli (source, recipe, maskVal)) {
	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis");
	    return false;
	} else {
	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
            Nannuli ++;
	}

	if (doKron && !psphotKron (source, recipe, maskVal)) {
	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Kron analysis");
	    return false;
	} else {
	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
            Nkron ++;
	}

        // re-subtract the object, leave local sky
        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
    }

    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
    psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models\n", Nconvolve);
    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);
    return true;
}
