Changeset 31153 for trunk/psModules/src/objects/pmFootprintCullPeaks.c
- Timestamp:
- Apr 4, 2011, 1:04:41 PM (15 years ago)
- Location:
- trunk/psModules
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/objects/pmFootprintCullPeaks.c (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules
- Property svn:ignore
-
old new 28 28 ChangeLog 29 29 psmodules-*.tar.* 30 a.out.dSYM
-
- Property svn:ignore
-
trunk/psModules/src/objects/pmFootprintCullPeaks.c
r30621 r31153 38 38 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive 39 39 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance 40 const float min_threshold) { // minimum permitted coll height 40 const float min_threshold, // minimum permitted coll height 41 const float max_threshold, 42 const bool useSmoothedImage 43 ) { // maximum permitted coll height 41 44 assert (img != NULL); assert (img->type.type == PS_TYPE_F32); 42 45 assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); … … 44 47 assert (fp != NULL); 45 48 46 if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do49 if (fp->peaks == NULL || fp->peaks->n < 2) { // nothing to do 47 50 return PS_ERR_NONE; 48 51 } … … 53 56 54 57 psImage *subImg = psImageSubset((psImage *)img, subRegion); 55 psImage *subWt = psImageSubset((psImage *)weight, subRegion); 56 assert (subImg != NULL && subWt != NULL); 58 psAssert (subImg, "trouble making local subimage"); 57 59 58 60 psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); … … 68 70 psArrayAdd (brightPeaks, 128, fp->peaks->data[0]); 69 71 72 // XXX test point 73 // pmPeak *myPeak = fp->peaks->data[0]; 74 // if ((fabs(myPeak->x - 205) < 100) && (fabs(myPeak->y - 806) < 100)) { 75 // fprintf (stderr, "test peak\n"); 76 // } 77 70 78 // allocate the peakFootprint and peakFPSpans containers -- these are re-used by pmFootprintsFindAtPoint to minimize allocs in this function 71 79 pmFootprint *peakFootprint = pmFootprintAlloc(fp->nspans, subImg); … … 77 85 psImage *subMask = psImageCopy(NULL, subImg, PS_TYPE_IMAGE_MASK); 78 86 79 // The brightest peak is always safe; go through other peaks trying to cull them 80 for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop 81 const pmPeak *peak = fp->peaks->data[i]; 82 int x = peak->x - subImg->col0; 83 int y = peak->y - subImg->row0; 84 // 85 // Find the level nsigma below the peak that must separate the peak 86 // from any of its friends 87 // 88 assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); 89 90 // const float stdev = sqrt(subWt->data.F32[y][x]); 91 // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; 92 93 const float stdev = sqrt(subWt->data.F32[y][x]); 94 const float flux = subImg->data.F32[y][x]; 95 const float fStdev = fabs(stdev/flux); 96 const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); 97 // if flux is negative, careful with fStdev 98 99 float threshold = flux - nsigma_delta*stdev_pad; 100 101 if (isnan(threshold)) { 102 // min_threshold is assumed to be below the detection threshold, 103 // so all the peaks are pmFootprint, and this isn't the brightest 104 continue; 105 } 106 107 // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? 108 if (threshold > subImg->data.F32[y][x]) { 109 threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; 110 } 111 112 if (threshold < min_threshold) { 113 threshold = min_threshold; 114 } 115 116 // init peakFootprint here? 117 pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x); 118 if (peakFPSpans->nStartSpans > 2000) { 119 // dumpfootprints(peakFootprint, peakFPSpans); 120 // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans); 121 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); 122 } 123 124 // at this point brightPeaks only has the peaks brighter than the current 125 126 // we set the IDs to either 1 (in peak) or 0 (not in peak) 127 pmSetFootprintID (idImg, peakFootprint, IN_PEAK); 128 129 // If this peak has not already been assigned to a source, then we can look for any 130 // brighter peaks within its footprint. Check if any of the previous (brighter) peaks 131 // are within the footprint of this peak If so, the current peak is bogus; drop it. 132 bool keep = true; 133 for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) { 134 const pmPeak *peak2 = fp->peaks->data[j]; 135 int x2 = peak2->x - subImg->col0; 136 int y2 = peak2->y - subImg->row0; 137 if (idImg->data.S32[y2][x2] == IN_PEAK) { 138 // There's a brighter peak within the footprint above threshold; so cull our initial peak 139 keep = false; 140 } 141 } 142 if (!keep) { 143 psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak"); // we should not drop any already assigned peaks. 144 continue; 145 } 146 147 psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); 87 // if we have many peaks in a large footprint, we waste a lot of time generating nearly identical footprints 88 // here we attempt to cull peaks drawing a single footprint for all peaks in some threshold range 89 // fprintf (stderr, "footprint: %d x %d : %d pix, %d peaks\n", subImg->numCols, subImg->numRows, fp->npix, (int) fp->peaks->n); 90 if ((fp->npix > 30000) && (fp->peaks->n > 10)) { 91 92 // max flux is above threshold for brightest peak 93 pmPeak *maxPeak = fp->peaks->data[0]; 94 float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux; 95 96 // we have a relationship between the bin and the threshold of: 97 // threshold = 0.25 beta^2 bin^2 + minThreshold 98 // thus, the max bin is: sqrt(4.0*(maxThreshold - minThreshold)/ALPHA^2) 99 100 # define ALPHA 0.1 101 102 float beta = nsigma_delta * ALPHA; 103 float beta2 = PS_SQR(beta); 104 int nBins = sqrt(4.0*(maxFlux - min_threshold)/beta2) + 10; // let's be extra generous 105 106 // create a vector to store the threshold bins used for each peak 107 psVector *threshbins = psVectorAlloc(fp->peaks->n, PS_TYPE_S32); 108 threshbins->data.S32[0] = -1; // we skip the first peak 109 110 /// create a vector to track if a peak has been tried or not: 111 psVector *peaktried = psVectorAlloc(fp->peaks->n, PS_TYPE_U8); 112 psVectorInit(peaktried, 0); 113 peaktried->data.U8[0] = true; // we skip the first peak 114 115 psVector *threshbounds = psVectorAlloc(nBins, PS_TYPE_F32); 116 for (int i = 0; i < nBins; i++) { 117 threshbounds->data.F32[i] = 0.25*beta2*PS_SQR(i) + min_threshold; 118 } 119 psAssert(threshbounds->data.F32[threshbounds->n-1] > maxFlux, "upper limit does not include max flux"); 120 121 psHistogram *threshist = psHistogramAllocGeneric(threshbounds); 122 123 // assign the peaks to the histogram bins based on their nominal thresholsd 124 for (int i = 1; i < fp->peaks->n; i++) { 125 const pmPeak *peak = fp->peaks->data[i]; 126 float flux = useSmoothedImage ? peak->smoothFlux : peak->rawFlux; 127 float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev; 128 129 // if flux is negative, careful with fStdev 130 const float fStdev = fabs(stdev/flux); 131 const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); 132 133 float threshold = flux - nsigma_delta*stdev_pad; 134 psAssert(!isnan(threshold), "impossible"); 135 136 if (threshold <= min_threshold) { 137 threshist->nums->data.F32[0] += 1.0; 138 threshbins->data.S32[i] = 0; 139 continue; 140 } 141 int bin = sqrt(4.0*(threshold - min_threshold)/beta2); 142 psAssert(bin >= 0, "impossible bin"); 143 144 bin = PS_MIN(bin, threshist->nums->n - 1); 145 threshist->nums->data.F32[bin] += 1.0; 146 147 // record the bin selected for this peak 148 threshbins->data.S32[i] = bin; 149 } 150 151 // XXX TEST: did we assign correctly 152 # if (0) 153 int nPeaks = 1; // we don't put the brightest in the histogram 154 for (int i = 0; i < threshist->nums->n; i++) { 155 if (threshist->nums->data.F32[i] > 0) { 156 fprintf (stderr, "%f : %f : %d\n", threshist->bounds->data.F32[i], threshist->bounds->data.F32[i+1], (int)threshist->nums->data.F32[i]); 157 nPeaks += threshist->nums->data.F32[i]; 158 } 159 } 160 fprintf (stderr, "%d peaks vs %d in histogram\n", (int) fp->peaks->n, nPeaks); 161 # endif 162 163 // XXX for the moment, we will use the alternate cull for all peaks -- it might be 164 // faster to use the standard process for the peaks for which the threshold bin 165 // contains < N sources (N ~ 5-10?) 166 167 // loop over the threshold bins from brightest to faintest 168 for (int i = threshist->nums->n - 1; i >= 0; i--) { 169 if (threshist->nums->data.F32[i] == 0) continue; 170 171 // we are going to examine the footprints at this threshold 172 float threshold = 0.5*(threshist->bounds->data.F32[i] + threshist->bounds->data.F32[i+1]); 173 174 // generate all footprints corresponding to this threshold 175 psArray *myFP = pmFootprintsFind(subImg, threshold, 5); 176 if (!myFP) { 177 psWarning ("missing footprint?"); 178 continue; 179 } 180 if (!myFP->n) { 181 psWarning ("empty footprint?"); 182 continue; 183 } 184 185 // an array to track if the footprint has a brighter peak or not 186 psVector *found = psVectorAlloc(myFP->n + 1, PS_TYPE_U8); 187 psVectorInit(found, 0); 188 int nFound = 0; 189 190 // set IDs to distinguish the footprints 191 psImageInit(idImg, 0); 192 pmSetFootprintArrayIDsForImage(idImg, myFP, true); 193 194 // check which footprints contain already-accepted (brighter) peaks 195 // (we can give up if/when we found a peak for all footprints) 196 for (int j = 0; (j < brightPeaks->n) && (nFound < found->n); j++) { 197 const pmPeak *peak = brightPeaks->data[j]; 198 int x = peak->x - subImg->col0; 199 int y = peak->y - subImg->row0; 200 int myID = idImg->data.S32[y][x]; 201 psAssert(myID >= 0, "impossible"); 202 psAssert(myID < found->n, "impossible"); 203 204 if (myID == 0) continue; // bright peak is not in a footprint 205 if (found->data.U8[myID]) continue; // we already know this footprint contains a peak 206 found->data.U8[myID] = true; 207 nFound ++; 208 } 209 210 // check the peaks from this threshold bin: if they land in a footprint which has 211 // been found, we should drop that peak. otherwise, keep it 212 for (int j = 0; j < fp->peaks->n; j++) { 213 pmPeak *peak = fp->peaks->data[j]; 214 if (peak->assigned) continue; // peak is already claimed by a source -- don't cull 215 216 // skip peaks if we've already tried them 217 if (peaktried->data.U8[j]) continue; 218 219 // is this peak in the threshold bin of interest? 220 if (threshbins->data.S32[j] != i) continue; 221 222 // do not try this peak again 223 peaktried->data.U8[j] = true; 224 225 int x = peak->x - subImg->col0; 226 int y = peak->y - subImg->row0; 227 int myID = idImg->data.S32[y][x]; 228 psAssert(myID < found->n, "impossible"); 229 230 // a peak in this threshold bin without a valid footprint comes from a region 231 // with only a handful of pixels (1 or more from the peak itself). It probably 232 // cannot be joined to a neighbor 233 if (myID == 0) { 234 psArrayAdd (brightPeaks, 128, peak); 235 continue; 236 } 237 238 // keep this peak if found->data.U8[myID] == false (no brighter peak in the footprint) 239 if (!found->data.U8[myID]) { 240 // fprintf (stderr, "keeping %d: %d,%d\n", j, peak->x, peak->y); 241 psArrayAdd (brightPeaks, 128, peak); 242 continue; 243 } 244 // fprintf (stderr, "skipping %d: %d,%d\n", j, peak->x, peak->y); 245 } 246 psFree (myFP); 247 psFree (found); 248 } 249 psFree (threshist); 250 psFree (threshbounds); 251 psFree (threshbins); 252 psFree (peaktried); 253 254 } else { 255 256 // The brightest peak is always safe; go through other peaks trying to cull them 257 for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop 258 const pmPeak *peak = fp->peaks->data[i]; 259 260 float flux = useSmoothedImage ? peak->smoothFlux : peak->rawFlux; 261 float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev; 262 263 // if flux is negative, careful with fStdev 264 const float fStdev = fabs(stdev/flux); 265 const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); 266 267 float threshold = flux - nsigma_delta*stdev_pad; 268 269 if (isnan(threshold)) { 270 // min_threshold is assumed to be below the detection threshold, 271 // so all the peaks are pmFootprint, and this isn't the brightest 272 continue; 273 } 274 275 // just in case, force the threshold below the peak source flux 276 if (threshold > flux) { 277 threshold = flux - 10*FLT_EPSILON; 278 } 279 280 if (threshold < min_threshold) { 281 threshold = min_threshold; 282 } 283 if (threshold > max_threshold) { 284 threshold = max_threshold; 285 } 286 287 pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x); 288 // if (peakFPSpans->nStartSpans > 2000) { 289 // // dumpfootprints(peakFootprint, peakFPSpans); 290 // // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans); 291 // // 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); 292 // } 293 294 // at this point brightPeaks only has the peaks brighter than the current 295 296 // we set the IDs to either 1 (in peak) or 0 (not in peak) 297 pmSetFootprintID (idImg, peakFootprint, IN_PEAK); 298 299 // If this peak has not already been assigned to a source, then we can look for any 300 // brighter peaks within its footprint. Check if any of the previous (brighter) peaks 301 // are within the footprint of this peak If so, the current peak is bogus; drop it. 302 bool keep = true; 303 for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) { 304 // const pmPeak *peak2 = fp->peaks->data[j]; XXX isn't this an error? we only care about the kept brighter peak, right? 305 const pmPeak *peak2 = brightPeaks->data[j]; 306 int x2 = peak2->x - subImg->col0; 307 int y2 = peak2->y - subImg->row0; 308 if (idImg->data.S32[y2][x2] == IN_PEAK) { 309 // There's a brighter peak within the footprint above threshold; so cull our initial peak 310 keep = false; 311 } 312 } 313 if (!keep) { 314 psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak"); // we should not drop any already assigned peaks. 315 continue; 316 } 317 psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); 318 } 148 319 } 149 320 … … 153 324 psFree(idImg); 154 325 psFree(subImg); 155 psFree(subWt);156 326 psFree(subMask); 157 327 psFree(peakFootprint);
Note:
See TracChangeset
for help on using the changeset viewer.
