/* @file  pmFootprint.c
 * low-level pmFootprint functions
 *
 * @author RHL, Princeton & IfA; EAM, IfA
 *
 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
 * @date $Date: 2008-12-08 02:51:14 $
 * Copyright 2006 Institute for Astronomy, University of Hawaii
 */

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

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <strings.h>
#include <pslib.h>
#include "pmSpan.h"
#include "pmFootprint.h"
#include "pmPeaks.h"

 /*
  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
  * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
  * to reach a still higher peak --- and if that coll's more than nsigma DN below your
  * starting point, discard the peak.
  */

# define IN_PEAK 1
psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
                                 const psImage *weight, // corresponding variance image
                                 pmFootprint *fp, // Footprint containing mortal peaks
                                 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
                                 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
                                 const float min_threshold) { // minimum permitted coll height
    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    assert (fp != NULL);

    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
        return PS_ERR_NONE;
    }

    psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;

    psImage *subImg = psImageSubset((psImage *)img, subRegion);
    psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    assert (subImg != NULL && subWt != NULL);

    psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);

    // We need a psArray of peaks brighter than the current peak.
    // We reject peaks which either:
    // 1) are below the local threshold
    // 2) have a brighter peak within their threshold

    // allocate the full-sized array.  if the final array is much smaller, we can realloc
    // at that point.
    psArray *brightPeaks = psArrayAllocEmpty(fp->peaks->n);
    psArrayAdd (brightPeaks, 128, fp->peaks->data[0]);

    // The brightest peak is always safe; go through other peaks trying to cull them
    for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
        const pmPeak *peak = fp->peaks->data[i];
        int x = peak->x - subImg->col0;
        int y = peak->y - subImg->row0;
        //
        // Find the level nsigma below the peak that must separate the peak
        // from any of its friends
        //
        assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);

        // const float stdev = sqrt(subWt->data.F32[y][x]);
        // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;

        const float stdev = sqrt(subWt->data.F32[y][x]);
        const float flux = subImg->data.F32[y][x];
        const float fStdev = fabs(stdev/flux);
        const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
        // if flux is negative, careful with fStdev

        float threshold = flux - nsigma_delta*stdev_pad;

        if (isnan(threshold) || threshold < min_threshold) {
            // min_threshold is assumed to be below the detection threshold,
            // so all the peaks are pmFootprint, and this isn't the brightest
            continue;
        }

        // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
        if (threshold > subImg->data.F32[y][x]) {
            threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
        }

        // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
        // perhaps this should alloc a single ID image above and pass it in to be set.

        // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)

        // XXX this is not quite there yet:
        // pmFootprint *peakFootprint = NULL;
        // int area = subImg->numCols * subImg->numRows;
        // if (area > 30000) {

        pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);

        // at this point brightPeaks only has the peaks brighter than the current

        // XXX need to supply the image here
        // we set the IDs to either 1 (in peak) or 0 (not in peak)
        pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
        psFree(peakFootprint);

        // If this peak has not already been assigned to a source, then we can look for any
        // brighter peaks within its footprint. Check if any of the previous (brighter) peaks
        // are within the footprint of this peak If so, the current peak is bogus; drop it.
        bool keep = true;
        for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) {
            const pmPeak *peak2 = fp->peaks->data[j];
            int x2 = peak2->x - subImg->col0;
            int y2 = peak2->y - subImg->row0;
            if (idImg->data.S32[y2][x2] == IN_PEAK) {
                // There's a brighter peak within the footprint above threshold; so cull our initial peak
                keep = false;
	    }
        }
        if (!keep) {
	    psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak");  // we should not drop any already assigned peaks.
	    continue;
	}

        psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
    }

    psFree (fp->peaks);
    fp->peaks = brightPeaks;

    psFree(idImg);
    psFree(subImg);
    psFree(subWt);

    return PS_ERR_NONE;
}


 /*
  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
  * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
  * to reach a still higher peak --- and if that coll's more than nsigma DN below your
  * starting point, discard the peak.
  */
psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
                                 const psImage *weight, // corresponding variance image
                                 pmFootprint *fp, // Footprint containing mortal peaks
                                 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
                                 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
                                 const float min_threshold) { // minimum permitted coll height
    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    assert (fp != NULL);

    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
        return PS_ERR_NONE;
    }

    psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    const psImage *subImg = psImageSubset((psImage *)img, subRegion);
    const psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    assert (subImg != NULL && subWt != NULL);
    //
    // We need a psArray of peaks brighter than the current peak.  We'll fake this
    // by reusing the fp->peaks but lying about n.
    //
    // We do this for efficiency (otherwise I'd need two peaks lists), and we are
    // rather too chummy with psArray in consequence.  But it works.
    //
    psArray *brightPeaks = psArrayAlloc(0);
    psFree(brightPeaks->data);
    brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
    //
    // The brightest peak is always safe; go through other peaks trying to cull them
    //
    for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
        const pmPeak *peak = fp->peaks->data[i];
        int x = peak->x - subImg->col0;
        int y = peak->y - subImg->row0;
        //
        // Find the level nsigma below the peak that must separate the peak
        // from any of its friends
        //
        assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);

        // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations)
        // add additional padding in quadrature:

        const float stdev = sqrt(subWt->data.F32[y][x]);
        const float flux = subImg->data.F32[y][x];
        const float fStdev = stdev/flux;
        const float stdev_pad = flux * hypot(fStdev, fPadding);
        float threshold = flux - nsigma_delta*stdev_pad;

        if (isnan(threshold) || threshold < min_threshold) {
#if 1     // min_threshold is assumed to be below the detection threshold,
          // so all the peaks are pmFootprint, and this isn't the brightest
            // XXX mark peak to be dropped
            (void)psArrayRemoveIndex(fp->peaks, i);
            i--;                        // we moved everything down one
            continue;
#else
#error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
            threshold = min_threshold;
#endif
        }

        // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
        if (threshold > subImg->data.F32[y][x]) {
            threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
        }

        // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
        // perhaps this should alloc a single ID image above and pass it in to be set.

        const int peak_id = 1;          // the ID for the peak of interest
        brightPeaks->n = i;             // only stop at a peak brighter than we are

        // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)

        pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
        brightPeaks->n = 0;             // don't double free
        psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
        psFree(peakFootprint);

        // Check if any of the previous (brighter) peaks are within the footprint of this peak
        // If so, the current peak is bogus; drop it.
        int j;
        for (j = 0; j < i; j++) {
            const pmPeak *peak2 = fp->peaks->data[j];
            int x2 = peak2->x - subImg->col0;
            int y2 = peak2->y - subImg->row0;
            const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak

            if (peak2_id == peak_id) {  // There's a brighter peak within the footprint above
                ;                       // threshold; so cull our initial peak
                (void)psArrayRemoveIndex(fp->peaks, i);
                i--;                    // we moved everything down one
                break;
            }
        }
        if (j == i) {
            j++;
        }

        psFree(idImg);
    }

    brightPeaks->n = 0; psFree(brightPeaks);
    psFree(subImg);
    psFree(subWt);

    return PS_ERR_NONE;
}

