# include "psphotInternal.h"

/*
 * Cull a set of peaks contained in a psArray of pmFootprints
 */
psErrorCode
psphotCullPeaks(const psImage *image,   // the image wherein lives the footprint
                const psImage *weight,  // corresponding variance image
                const psMetadata *recipe,
                psArray *footprints) {  // array of pmFootprints
    bool status = false;
    float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA");
    if (!status) {
        nsigma_delta = 0; // min.
    }
    float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN");
    if (!status) {
        nsigma_min = 0;
    }
    const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV");

    return pmFootprintArrayCullPeaks(image, weight, footprints,
                                     nsigma_delta, nsigma_min*skyStdev);
}


/*
 * Cull an entire psArray of pmFootprints
 * XXX drop this intermediate level function?
 */
psErrorCode
pmFootprintArrayCullPeaks(const psImage *img, // the image wherein lives the footprint
			  const psImage *weight,	// corresponding variance image
			  psArray *footprints, // array of pmFootprints
			  const float nsigma_delta, // how many sigma above local background a peak
    					// needs to be to survive
			  const float min_threshold) { // minimum permitted coll height
    for (int i = 0; i < footprints->n; i++) {
	pmFootprint *fp = footprints->data[i];
	if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) {
	    return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
	}
    }
    
    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 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);
	const float stdev = sqrt(subWt->data.F32[y][x]);
	float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
	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((psImage *)subImg);
    psFree((psImage *)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.
  */

# 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 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;
	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)

	// at this point brightPeaks only has the peaks brighter than the current
	pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);

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

	// 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 && (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) continue;

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

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

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

    return PS_ERR_NONE;
}

