Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 24874) +++ trunk/psModules/src/objects/Makefile.am (revision 24875) @@ -10,4 +10,5 @@ pmFootprintArraysMerge.c \ pmFootprintAssignPeaks.c \ + pmFootprintCullPeaks.c \ pmFootprintFind.c \ pmFootprintFindAtPoint.c \ Index: trunk/psModules/src/objects/pmFootprint.h =================================================================== --- trunk/psModules/src/objects/pmFootprint.h (revision 24874) +++ trunk/psModules/src/objects/pmFootprint.h (revision 24875) @@ -51,12 +51,16 @@ void pmSetFootprintArrayIDsForImage(psImage *idImage, const psArray *footprints, // the footprints to insert - const bool relativeIDs); // show IDs starting at 0, not pmFootprint->id + const bool relativeIDs // show IDs starting at 0, not pmFootprint->id + ); psErrorCode pmFootprintsAssignPeaks(psArray *footprints, const psArray *peaks); -psErrorCode pmFootprintArrayCullPeaks(const psImage *img, const psImage *weight, psArray *footprints, - const float nsigma, const float threshold_min); -psErrorCode pmFootprintCullPeaks(const psImage *img, const psImage *weight, pmFootprint *fp, - const float nsigma, const float threshold_min); +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 + ); psArray *pmFootprintArrayToPeaks(const psArray *footprints); Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 24875) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 24875) @@ -0,0 +1,259 @@ +/* @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 +#endif + +#include +#include +#include +#include +#include +#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) + + // 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; +} + + + /* + * 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((psImage *)subImg); + psFree((psImage *)subWt); + + return PS_ERR_NONE; +} +