/* @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 "pmFootprintSpans.h"
#include "pmPeaks.h"

bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp);

 /*
  * 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]);

    // allocate the peakFootprint and peakFPSpans containers -- these are re-used by pmFootprintsFindAtPoint to minimize allocs in this function
    pmFootprint *peakFootprint = pmFootprintAlloc(fp->nspans, subImg);
    pmFootprintSpans *peakFPSpans = pmFootprintSpansAlloc(2*fp->nspans);

    // allocate empty spans for the footprints
    pmFootprintAllocEmptySpans(peakFootprint, fp->nspans);

    psImage *subMask = psImageCopy(NULL, subImg, PS_TYPE_IMAGE_MASK);

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

	// init peakFootprint here?
        pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x);
	if (peakFPSpans->nStartSpans > 2000) {
	    // dumpfootprints(peakFootprint, peakFPSpans);
	    // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans);
	    fprintf (stderr, "big test footprint: %f %f to %f %f (%d pix)\n", peakFootprint->bbox.x0, peakFootprint->bbox.y0, peakFootprint->bbox.x1, peakFootprint->bbox.y1, peakFootprint->npix);
	}

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

        // we set the IDs to either 1 (in peak) or 0 (not in peak)
        pmSetFootprintID (idImg, peakFootprint, IN_PEAK);

        // 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);
    psFree(subMask);
    psFree(peakFootprint);
    psFree(peakFPSpans);

    return PS_ERR_NONE;
}


bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp) {

    FILE *f1 = fopen ("fp.dat", "w");
    if (!f1) return false;

    for (int i = 0; i < fp->nspans; i++) {
	pmSpan *sp = fp->spans->data[i];
	fprintf (f1, "%d - %d : %d\n", sp->x0, sp->x1, sp->y);
    }
    fclose (f1);

    FILE *f2 = fopen ("fpSp.dat", "w");
    if (!f2) return false;

    for (int i = 0; i < fpSp->nStartSpans; i++) {
	pmStartSpan *sp = fpSp->startspans->data[i];
	if (!sp->span) continue;
	fprintf (f2, "%d - %d : %d\n", sp->span->x0, sp->span->x1, sp->span->y);
    }
    fclose (f2);
    return true;
}
