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

# include "psphot.h"

// 2006.02.02 : no leaks
psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks) 
{
    bool     status  = false;
    psArray *sources = NULL;
    float BIG_RADIUS;

    psTimerStart ("psphot");

    // determine properties (sky, moments) of initial sources
    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");

    sources = psArrayAlloc (peaks->n);

    for (int i = 0; i < peaks->n; i++) {

	// create a new source, add peak
	pmSource *source = pmSourceAlloc();
	source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]);

	// allocate image, weight, mask arrays for each peak (square of radius OUTER)
	pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
	if (!strcasecmp (breakPt, "PEAKS")) continue;

	// XXX skip faint sources?

	// measure a local sky value
	// XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
	status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
	if (!status) {
	  psFree (source);
	  continue;
	}

	// measure the local sky variance (needed if noise is not sqrt(signal))
	// XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
	status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER);
	if (!status) {
	  psFree (source);
	  continue;
	}

	// measure basic source moments
	status = pmSourceMoments (source, RADIUS);
	if (status) {
	    // add to the source array
	    psArrayAdd (sources, 100, source);
	    psFree (source);
	    continue;
	}

	// if no valid pixels, or massive swing, likely saturated source,
	// try a much larger box
	BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
	psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
	status = pmSourceMoments (source, BIG_RADIUS);
	if (status) {
	    // add to the source array
	    psArrayAdd (sources, 100, source);
	    psFree (source);
	    continue;
	}

	psFree (source);
	continue;
    }
     
    psLogMsg ("psphot", 3, "%ld moments: %f sec\n", sources->n, psTimerMark ("psphot"));

    return (sources);
}

// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky)
