Index: /trunk/psphot/src/psphotCullPeaks.c =================================================================== --- /trunk/psphot/src/psphotCullPeaks.c (revision 24876) +++ /trunk/psphot/src/psphotCullPeaks.c (revision 24877) @@ -18,25 +18,14 @@ nsigma_min = 0; } + float fPadding = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_PAD"); + if (!status) { + fPadding = 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 + const float min_threshold = nsigma_min*skyStdev; + 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) { + if (pmFootprintCullPeaks(image, weight, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) { return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); } @@ -45,220 +34,2 @@ 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; -} -