Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 24888) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 26893) @@ -29,11 +29,11 @@ */ -# define IN_PEAK 1 +# 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 + 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); @@ -42,8 +42,8 @@ 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) + 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; @@ -55,5 +55,5 @@ psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); - // We need a psArray of peaks brighter than the current peak. + // We need a psArray of peaks brighter than the current peak. // We reject peaks which either: // 1) are below the local threshold @@ -67,68 +67,73 @@ // 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 + 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; } - // 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); - - // 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]); + psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); } @@ -151,9 +156,9 @@ */ 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 + 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); @@ -162,8 +167,8 @@ 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) + 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; @@ -185,79 +190,79 @@ // 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; + 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; + 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); + } + + // 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); + psFree(subImg); + psFree(subWt); return PS_ERR_NONE;