# include "psphotInternal.h"

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");
    if (!status) return NULL;
    float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    if (!status) return NULL;
    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    if (!status) return NULL;
    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    if (!status) return NULL;
    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    if (!status) return NULL;

    sources = psArrayAllocEmpty (peaks->n);

    int Nfail = 0;
    int Nmoments = 0;
    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")) {
            psArrayAdd (sources, 100, source);
            psFree (source);
            continue;
        }

        // skip faint sources
	if (source->peak->SN < MIN_SN) {
            psArrayAdd (sources, 100, source);
            psFree (source);
	    continue;
	}

        // 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);
	  Nfail ++;
          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);
	  Nfail ++;
          continue;
        }

        // measure basic source moments
        status = pmSourceMoments (source, RADIUS);
        if (status) {
            // add to the source array
            psArrayAdd (sources, 100, source);
            psFree (source);
	    Nmoments ++;
            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);
	    Nmoments ++;
            continue;
        }

        psFree (source);
	Nfail ++;
        continue;
    }

    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot"));

    return (sources);
}

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