#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

# include "psphot.h"

// select objects fitted with PSF model
// re-fit all of them with the non-poisson errors
// only allow the normalization to vary 
// XXX eventually, we should be able to do this with linear fitting...
bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 

    int Nfit = 0;
    bool status;
    psF32 *PARp, *PARc;
    psF32 x, y;

    psTimerStart ("psphot");

    // source analysis is done in S/N order (brightest first)
    sources = psArraySort (sources, psphotSortBySN);
    
    // set fitting method to use non-poisson errors (local sky error)
    pmSourceFitModelInit (15, 0.1, false);
    
    // 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 ("psphot", "analysis region mis-defined");

    FILE *f = fopen ("bias.dat", "w");
    if (f == NULL) psAbort ("psphot", "can't open output file");

    for (int i = 0; i < sources->n; i++) {
	pmSource *source = sources->data[i];

	// skip lower-quality objects
	if (source->type != PM_SOURCE_TYPE_STAR) continue; 
	if (source->mode &  PM_SOURCE_MODE_POOR) continue;
	if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
	if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
	if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
	if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;

	// if model is NULL, we don't have a starting guess
	if (source->modelPSF == NULL) continue;

	if (source->moments->x < AnalysisRegion.x0) continue;
	if (source->moments->y < AnalysisRegion.y0) continue;
	if (source->moments->x > AnalysisRegion.x1) continue;
	if (source->moments->y > AnalysisRegion.y1) continue;

	// replace object in image
	pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false);

	// make a temporary model (we don't keep the result of this analysis)
	pmModel *PSF = pmModelCopy (source->modelPSF);

	// extend source radius as needed
	psphotCheckRadiusPSF (readout, source, PSF);

	x = PSF->params->data.F32[2];
	y = PSF->params->data.F32[3];

	// fit PSF model (set/unset the pixel mask)
	psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
	pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM);
	psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
    
	// re-subtract PSF for object, leave local sky
	pmModelSub (source->pixels, source->mask, source->modelPSF, false, false);

	PARp = source->modelPSF->params->data.F32;
	PARc = PSF->params->data.F32;
	fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky);
	Nfit ++;
    }
    fclose (f);
    psLogMsg ("psphot", 3, "measure PSF weighting bias for %d objects: %f sec\n", Nfit, psTimerMark ("psphot"));

    return (true);
}
