# include "psphotInternal.h"
# define RADIUS_TYPE int

static float PSF_FIT_NSIGMA;
static float PSF_FIT_PADDING;
static float PSF_APERTURE = 0;	// radius to use in PSF aperture mags
static float PSF_FIT_RADIUS = 0;	// radius to use in fitting (ignored if <= 0,
					// and a per-object radius is calculated)

bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) {

    bool status = true;

    PSF_FIT_NSIGMA = psMetadataLookupF32(&status, recipe, "PSF_FIT_NSIGMA");
    PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING");

    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS");
    if (!status) {
	PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
    }

    PSF_APERTURE =  psMetadataLookupF32(&status, analysis, "PSF_APERTURE");
    if (!status) {
	PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
    }

    return true;
}

// call this function whenever you (re)-define the PSF model
bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal)
{
    psF32 *PAR = model->params->data.F32;

    // XXX do we have a better value for the sky noise level?  not really...
    pmMoments *moments = source->moments;

    // set the fit radius based on the object flux limit and the model
    float radiusFit = PSF_FIT_RADIUS;
    if (radiusFit <= 0) {		// use fixed radius
	if (moments == NULL) {
	    radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
	} else {
	    radiusFit = model->modelRadius(model->params, 1.0);
	}
	model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
    } else {
	model->fitRadius = radiusFit;
    }
    if (isnan(model->fitRadius)) psAbort("error in radius");
	
    if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	model->fitRadius *= 2;
    }

    // radius used to measure aperture photometry
    source->apRadius = PSF_APERTURE;

    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);

    // set the mask to flag the excluded pixels
    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
    return status;
}

bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal, float dR) {

    psF32 *PAR = model->params->data.F32;

    pmMoments *moments = source->moments;
    if (moments == NULL) return false;

    // set the fit radius based on the object flux limit and the model
    float radiusFit = PSF_FIT_RADIUS;
    if (radiusFit <= 0) {		// use fixed radius
	if (moments == NULL) {
	    radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
	} else {
	    radiusFit = model->modelRadius(model->params, 1.0);
	}
	model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
    } else {
	model->fitRadius = radiusFit;
    }
    if (isnan(model->fitRadius)) psAbort("error in radius");

    // above sets a radius for a single star, bump by blend separation
    model->fitRadius += dR;

    if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
	model->fitRadius *= 2;
    }

    // radius used to measure aperture photometry
    source->apRadius = PSF_APERTURE;

    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);

    // set the mask to flag the excluded pixels
    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
    return status;
}

static float EXT_FIT_NSIGMA;
static float EXT_FIT_PADDING;
static float EXT_FIT_MAX_RADIUS;

bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) {

    bool status;

    EXT_FIT_NSIGMA     = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA");
    EXT_FIT_PADDING    = psMetadataLookupF32 (&status, recipe, "EXT_FIT_PADDING");
    EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");

    return true;
}

// call this function whenever you (re)-define the EXT model
float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) {

    psAssert (source, "source not defined??");
    psAssert (source->peak, "peak not defined??");

    pmPeak *peak = source->peak;

    // set the radius based on the footprint:
    if (!peak->footprint) goto escape;
    pmFootprint *footprint = peak->footprint;
    if (!footprint->spans) goto escape;
    if (footprint->spans->n < 1) goto escape;

    // find the max radius
    float radius = 0.0;
    for (int j = 0; j < footprint->spans->n; j++) {
	pmSpan *span = footprint->spans->data[j];

	float dY  = span->y  - peak->yf;
	float dX0 = span->x0 - peak->xf;
	float dX1 = span->x1 - peak->xf;

	radius = PS_MAX (radius, hypot(dY, dX0));
	radius = PS_MAX (radius, hypot(dY, dX1));
    }

    radius += EXT_FIT_PADDING;
    if (isnan(radius)) psAbort("error in radius");

    radius = PS_MIN (radius, EXT_FIT_MAX_RADIUS);

    // redefine the pixels if needed
    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius);

    // set the mask to flag the excluded pixels
    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal);
    return radius;

escape:
    return NAN;
    // bool result = psphotCheckRadiusEXT (readout, source, model, markVal);
    // return result;
}

// alternative EXT radius based on model guess (for use without footprints)
bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) {

    psAbort ("do not use this function");

    psF32 *PAR = model->params->data.F32;

    pmMoments *moments = source->moments;
    if (moments == NULL) return false;

    // set the fit radius based on the object flux limit and the model
    float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky);

    model->fitRadius = rawRadius + EXT_FIT_PADDING;
    if (isnan(model->fitRadius)) psAbort("error in radius");

    // redefine the pixels if needed
    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);

    // set the mask to flag the excluded pixels
    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
    return status;
}
