# include "psphotInternal.h"

// aperture-like measurements for extended sources
// flux in simple, circular apertures

// **** it looks like this function will re-point the source pixels at the specified FILERULE
// **** I need to distinguish PSF-matched images from raw
// **** save (somewhere) the PSF-matched PSF values

// this function measures the radial aperture fluxes for the set of readouts.  this function
// may be called multiple times (presumably with different matched PSF sizes).  we must have
// already added an entry to the readout->analysis identifying the FWHM of this version.

bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF) {

    bool status;
    int Nradial = 0;

    psTimerStart ("psphot.radial");

    // 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, "RADIAL_APERTURES")) {
	psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n");
	return true;
    }

    // number of images used to define sources
    int nImages = psphotFileruleCount(config, filerule);

    // radMax stores the upper bounds of the annuli
    // XXX keep the same name here as for the petrosian / elliptical apertures?
    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
    float outerRadius = radMax->data.F32[radMax->n - 1];

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

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

    // how many target PSFs do we want?
    int nPSFsizes = 0;
    { 
	psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
	if (status) {
	    nPSFsizes = 1;
	} else {
	    psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
	    psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
	    nPSFsizes = fwhmValues->n;
	}
    }
    
    // source analysis is done in S/N order (brightest first)
    objects = psArraySort (objects, pmPhotObjSortByFlux);

    // generate look-up arrays for readouts
    psArray *readouts = psArrayAlloc(nImages);
    for (int i = 0; i < nImages; 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?");

	psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
	if (!fwhmValues) {
	    psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
	    return false;
	}
	if (fwhmValues->n != nMatchedPSF + 1) {
	    psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM sequence is inconsistent this readout");
	    return false;
	}
	psLogMsg ("psphot", PS_LOG_DETAIL, "PSF FWHM of %s : %f pixels\n", file->name, fwhmValues->data.F32[nMatchedPSF]);

	readouts->data[i] = psMemIncrRefCounter(readout);
    }

    // process the objects in order.  
    for (int i = 0; i < objects->n; i++) {
        pmPhotObj *object = objects->data[i];
	if (!object) continue;
	if (!object->sources) continue;

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

	    pmSource *source = object->sources->data[j];

	    // skip PSF-like and non-astronomical objects
	    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;

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

	    int index = source->imageID;
	    if (index >= readouts->n) continue; // skip the sources generated by the chisq image
	    pmReadout *readout = readouts->data[index];

	    // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
	    // psphotVisualShowImage(readout);

	    // allocate pmSourceExtendedParameters, if not already defined
	    if (source->parent) {
		if (!source->parent->radialAper) {
		    source->parent->radialAper = psArrayAlloc(nPSFsizes);
		}
	    } else {
		if (!source->radialAper) {
		    source->radialAper = psArrayAlloc(nPSFsizes);
		}
	    }

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

	    // we need to change the view for the radial aperture analysis, but we want to recover exactly
	    // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
	    psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
	    psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
	    psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
	    psRegion oldRegion    = source->region;

	    Nradial ++;

	    // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
	    // psphotVisualShowImage(readout);

	    // force source image to be a bit larger...
	    pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);

	    if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, nMatchedPSF)) {
		psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
	    } else {
		source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
		if (source->parent) {
		    source->parent->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
		}
	    }

	    pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
	    psFree(source->maskObj);   source->maskObj   = oldMaskObj;
	    psFree(source->modelFlux); source->modelFlux = oldModelFlux;
	    psFree(source->psfImage);  source->psfImage  = oldPSFimage;

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

    psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);

    psFree(readouts);
    return true;
}
