Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 26893) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 27672) @@ -20,5 +20,8 @@ #include "pmSpan.h" #include "pmFootprint.h" +#include "pmFootprintSpans.h" #include "pmPeaks.h" + +bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp); /* @@ -65,4 +68,13 @@ 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 @@ -98,22 +110,16 @@ } - // 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); + // 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 - // 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 @@ -144,4 +150,7 @@ psFree(subImg); psFree(subWt); + psFree(subMask); + psFree(peakFootprint); + psFree(peakFPSpans); return PS_ERR_NONE; @@ -149,122 +158,24 @@ - /* - * 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); +bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp) { - if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do - return PS_ERR_NONE; + 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); - 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); + FILE *f2 = fopen ("fpSp.dat", "w"); + if (!f2) return false; - // 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); + 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); } - - brightPeaks->n = 0; psFree(brightPeaks); - psFree(subImg); - psFree(subWt); - - return PS_ERR_NONE; + fclose (f2); + return true; } -