Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 30621) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 31153) @@ -38,5 +38,8 @@ 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 float min_threshold, // minimum permitted coll height + const float max_threshold, + const bool useSmoothedImage + ) { // maximum permitted coll height assert (img != NULL); assert (img->type.type == PS_TYPE_F32); assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); @@ -44,5 +47,5 @@ assert (fp != NULL); - if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do + if (fp->peaks == NULL || fp->peaks->n < 2) { // nothing to do return PS_ERR_NONE; } @@ -53,6 +56,5 @@ psImage *subImg = psImageSubset((psImage *)img, subRegion); - psImage *subWt = psImageSubset((psImage *)weight, subRegion); - assert (subImg != NULL && subWt != NULL); + psAssert (subImg, "trouble making local subimage"); psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); @@ -68,4 +70,10 @@ psArrayAdd (brightPeaks, 128, fp->peaks->data[0]); + // XXX test point + // pmPeak *myPeak = fp->peaks->data[0]; + // if ((fabs(myPeak->x - 205) < 100) && (fabs(myPeak->y - 806) < 100)) { + // fprintf (stderr, "test peak\n"); + // } + // allocate the peakFootprint and peakFPSpans containers -- these are re-used by pmFootprintsFindAtPoint to minimize allocs in this function pmFootprint *peakFootprint = pmFootprintAlloc(fp->nspans, subImg); @@ -77,73 +85,236 @@ 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 - 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)) { - // 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; - } - - if (threshold < min_threshold) { - threshold = min_threshold; - } - - // 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 - - // we set the IDs to either 1 (in peak) or 0 (not in peak) - pmSetFootprintID (idImg, peakFootprint, IN_PEAK); - - // 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; - } - - psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); + // if we have many peaks in a large footprint, we waste a lot of time generating nearly identical footprints + // here we attempt to cull peaks drawing a single footprint for all peaks in some threshold range + // fprintf (stderr, "footprint: %d x %d : %d pix, %d peaks\n", subImg->numCols, subImg->numRows, fp->npix, (int) fp->peaks->n); + if ((fp->npix > 30000) && (fp->peaks->n > 10)) { + + // max flux is above threshold for brightest peak + pmPeak *maxPeak = fp->peaks->data[0]; + float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux; + + // we have a relationship between the bin and the threshold of: + // threshold = 0.25 beta^2 bin^2 + minThreshold + // thus, the max bin is: sqrt(4.0*(maxThreshold - minThreshold)/ALPHA^2) + +# define ALPHA 0.1 + + float beta = nsigma_delta * ALPHA; + float beta2 = PS_SQR(beta); + int nBins = sqrt(4.0*(maxFlux - min_threshold)/beta2) + 10; // let's be extra generous + + // create a vector to store the threshold bins used for each peak + psVector *threshbins = psVectorAlloc(fp->peaks->n, PS_TYPE_S32); + threshbins->data.S32[0] = -1; // we skip the first peak + + /// create a vector to track if a peak has been tried or not: + psVector *peaktried = psVectorAlloc(fp->peaks->n, PS_TYPE_U8); + psVectorInit(peaktried, 0); + peaktried->data.U8[0] = true; // we skip the first peak + + psVector *threshbounds = psVectorAlloc(nBins, PS_TYPE_F32); + for (int i = 0; i < nBins; i++) { + threshbounds->data.F32[i] = 0.25*beta2*PS_SQR(i) + min_threshold; + } + psAssert(threshbounds->data.F32[threshbounds->n-1] > maxFlux, "upper limit does not include max flux"); + + psHistogram *threshist = psHistogramAllocGeneric(threshbounds); + + // assign the peaks to the histogram bins based on their nominal thresholsd + for (int i = 1; i < fp->peaks->n; i++) { + const pmPeak *peak = fp->peaks->data[i]; + float flux = useSmoothedImage ? peak->smoothFlux : peak->rawFlux; + float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev; + + // if flux is negative, careful with fStdev + const float fStdev = fabs(stdev/flux); + const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); + + float threshold = flux - nsigma_delta*stdev_pad; + psAssert(!isnan(threshold), "impossible"); + + if (threshold <= min_threshold) { + threshist->nums->data.F32[0] += 1.0; + threshbins->data.S32[i] = 0; + continue; + } + int bin = sqrt(4.0*(threshold - min_threshold)/beta2); + psAssert(bin >= 0, "impossible bin"); + + bin = PS_MIN(bin, threshist->nums->n - 1); + threshist->nums->data.F32[bin] += 1.0; + + // record the bin selected for this peak + threshbins->data.S32[i] = bin; + } + + // XXX TEST: did we assign correctly +# if (0) + int nPeaks = 1; // we don't put the brightest in the histogram + for (int i = 0; i < threshist->nums->n; i++) { + if (threshist->nums->data.F32[i] > 0) { + fprintf (stderr, "%f : %f : %d\n", threshist->bounds->data.F32[i], threshist->bounds->data.F32[i+1], (int)threshist->nums->data.F32[i]); + nPeaks += threshist->nums->data.F32[i]; + } + } + fprintf (stderr, "%d peaks vs %d in histogram\n", (int) fp->peaks->n, nPeaks); +# endif + + // XXX for the moment, we will use the alternate cull for all peaks -- it might be + // faster to use the standard process for the peaks for which the threshold bin + // contains < N sources (N ~ 5-10?) + + // loop over the threshold bins from brightest to faintest + for (int i = threshist->nums->n - 1; i >= 0; i--) { + if (threshist->nums->data.F32[i] == 0) continue; + + // we are going to examine the footprints at this threshold + float threshold = 0.5*(threshist->bounds->data.F32[i] + threshist->bounds->data.F32[i+1]); + + // generate all footprints corresponding to this threshold + psArray *myFP = pmFootprintsFind(subImg, threshold, 5); + if (!myFP) { + psWarning ("missing footprint?"); + continue; + } + if (!myFP->n) { + psWarning ("empty footprint?"); + continue; + } + + // an array to track if the footprint has a brighter peak or not + psVector *found = psVectorAlloc(myFP->n + 1, PS_TYPE_U8); + psVectorInit(found, 0); + int nFound = 0; + + // set IDs to distinguish the footprints + psImageInit(idImg, 0); + pmSetFootprintArrayIDsForImage(idImg, myFP, true); + + // check which footprints contain already-accepted (brighter) peaks + // (we can give up if/when we found a peak for all footprints) + for (int j = 0; (j < brightPeaks->n) && (nFound < found->n); j++) { + const pmPeak *peak = brightPeaks->data[j]; + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + int myID = idImg->data.S32[y][x]; + psAssert(myID >= 0, "impossible"); + psAssert(myID < found->n, "impossible"); + + if (myID == 0) continue; // bright peak is not in a footprint + if (found->data.U8[myID]) continue; // we already know this footprint contains a peak + found->data.U8[myID] = true; + nFound ++; + } + + // check the peaks from this threshold bin: if they land in a footprint which has + // been found, we should drop that peak. otherwise, keep it + for (int j = 0; j < fp->peaks->n; j++) { + pmPeak *peak = fp->peaks->data[j]; + if (peak->assigned) continue; // peak is already claimed by a source -- don't cull + + // skip peaks if we've already tried them + if (peaktried->data.U8[j]) continue; + + // is this peak in the threshold bin of interest? + if (threshbins->data.S32[j] != i) continue; + + // do not try this peak again + peaktried->data.U8[j] = true; + + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + int myID = idImg->data.S32[y][x]; + psAssert(myID < found->n, "impossible"); + + // a peak in this threshold bin without a valid footprint comes from a region + // with only a handful of pixels (1 or more from the peak itself). It probably + // cannot be joined to a neighbor + if (myID == 0) { + psArrayAdd (brightPeaks, 128, peak); + continue; + } + + // keep this peak if found->data.U8[myID] == false (no brighter peak in the footprint) + if (!found->data.U8[myID]) { + // fprintf (stderr, "keeping %d: %d,%d\n", j, peak->x, peak->y); + psArrayAdd (brightPeaks, 128, peak); + continue; + } + // fprintf (stderr, "skipping %d: %d,%d\n", j, peak->x, peak->y); + } + psFree (myFP); + psFree (found); + } + psFree (threshist); + psFree (threshbounds); + psFree (threshbins); + psFree (peaktried); + + } else { + + // 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]; + + float flux = useSmoothedImage ? peak->smoothFlux : peak->rawFlux; + float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev; + + // if flux is negative, careful with fStdev + const float fStdev = fabs(stdev/flux); + const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); + + float threshold = flux - nsigma_delta*stdev_pad; + + if (isnan(threshold)) { + // min_threshold is assumed to be below the detection threshold, + // so all the peaks are pmFootprint, and this isn't the brightest + continue; + } + + // just in case, force the threshold below the peak source flux + if (threshold > flux) { + threshold = flux - 10*FLT_EPSILON; + } + + if (threshold < min_threshold) { + threshold = min_threshold; + } + if (threshold > max_threshold) { + threshold = max_threshold; + } + + 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 + + // we set the IDs to either 1 (in peak) or 0 (not in peak) + pmSetFootprintID (idImg, peakFootprint, IN_PEAK); + + // 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]; XXX isn't this an error? we only care about the kept brighter peak, right? + const pmPeak *peak2 = brightPeaks->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; + } + psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); + } } @@ -153,5 +324,4 @@ psFree(idImg); psFree(subImg); - psFree(subWt); psFree(subMask); psFree(peakFootprint);