# include "psphotInternal.h"

psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, pmDetections *detections)
{
    bool     status  = false;
    psArray *sources = NULL;
    float BIG_RADIUS;

    psTimerStart ("psphot");

    // bit-masks to test for good/bad pixels
    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
    assert (maskVal);

    // bit-mask to mark pixels not used in analysis
    psMaskType markVal = psMetadataLookupU8(&status, recipe, "MARK.PSPHOT");
    assert (markVal);

    // maskVal is used to test for rejected pixels, and must include markVal
    maskVal |= markVal;

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

    psArray *peaks = detections->peaks;
    if (!peaks) {
        psError(PS_ERR_UNEXPECTED_NULL, false, "No peaks found!");
        return NULL;
    }

    sources = psArrayAllocEmpty (peaks->n);

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

        pmPeak *peak = peaks->data[i];
        if (peak->assigned) continue;

        // create a new source, add peak
        pmSource *source = pmSourceAlloc();
        source->peak = psMemIncrRefCounter(peak);

        // 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")) {
            peak->assigned = true;
            psArrayAdd (sources, 100, source);
            psFree (source);
            continue;
        }

        // skip faint sources for moments measurement
        if (source->peak->SN < MIN_SN) {
            peak->assigned = true;
            psArrayAdd (sources, 100, source);
            psFree (source);
            continue;
        }

        // measure a local sky value
        // the local sky is now ignored; kept here for reference only
        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
        if (!status) {
          psFree (source);
          Nfail ++;
          psErrorClear();
          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, maskVal, markVal);
        if (!status) {
          psFree (source);
          Nfail ++;
          psErrorClear();
          continue;
        }

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

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

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

    psphotVisualShowMoments (sources);

    return (sources);
}

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