- Timestamp:
- Jan 7, 2010, 12:06:22 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/20091201/psModules/src/objects/pmFootprintCullPeaks.c
r24888 r26533 29 29 */ 30 30 31 # define IN_PEAK 1 31 # define IN_PEAK 1 32 32 psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint 33 const psImage *weight,// corresponding variance image34 35 36 37 33 const psImage *weight, // corresponding variance image 34 pmFootprint *fp, // Footprint containing mortal peaks 35 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive 36 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance 37 const float min_threshold) { // minimum permitted coll height 38 38 assert (img != NULL); assert (img->type.type == PS_TYPE_F32); 39 39 assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); … … 42 42 43 43 if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do 44 45 } 46 47 psRegion subRegion; 44 return PS_ERR_NONE; 45 } 46 47 psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) 48 48 subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; 49 49 subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; … … 55 55 psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); 56 56 57 // We need a psArray of peaks brighter than the current peak. 57 // We need a psArray of peaks brighter than the current peak. 58 58 // We reject peaks which either: 59 59 // 1) are below the local threshold … … 67 67 // The brightest peak is always safe; go through other peaks trying to cull them 68 68 for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 if (idImg->data.S32[y2][x2] == IN_PEAK) 127 128 129 130 131 132 69 const pmPeak *peak = fp->peaks->data[i]; 70 int x = peak->x - subImg->col0; 71 int y = peak->y - subImg->row0; 72 // 73 // Find the level nsigma below the peak that must separate the peak 74 // from any of its friends 75 // 76 assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); 77 78 // const float stdev = sqrt(subWt->data.F32[y][x]); 79 // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; 80 81 const float stdev = sqrt(subWt->data.F32[y][x]); 82 const float flux = subImg->data.F32[y][x]; 83 const float fStdev = fabs(stdev/flux); 84 const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); 85 // if flux is negative, careful with fStdev 86 87 float threshold = flux - nsigma_delta*stdev_pad; 88 89 if (isnan(threshold) || threshold < min_threshold) { 90 // min_threshold is assumed to be below the detection threshold, 91 // so all the peaks are pmFootprint, and this isn't the brightest 92 continue; 93 } 94 95 // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? 96 if (threshold > subImg->data.F32[y][x]) { 97 threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; 98 } 99 100 // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint 101 // perhaps this should alloc a single ID image above and pass it in to be set. 102 103 // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) 104 105 // XXX this is not quite there yet: 106 // pmFootprint *peakFootprint = NULL; 107 // int area = subImg->numCols * subImg->numRows; 108 // if (area > 30000) { 109 110 pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); 111 112 // at this point brightPeaks only has the peaks brighter than the current 113 114 // XXX need to supply the image here 115 // we set the IDs to either 1 (in peak) or 0 (not in peak) 116 pmSetFootprintID (idImg, peakFootprint, IN_PEAK); 117 psFree(peakFootprint); 118 119 // Check if any of the previous (brighter) peaks are within the footprint of this peak 120 // If so, the current peak is bogus; drop it. 121 bool keep = true; 122 for (int j = 0; keep && (j < brightPeaks->n); j++) { 123 const pmPeak *peak2 = fp->peaks->data[j]; 124 int x2 = peak2->x - subImg->col0; 125 int y2 = peak2->y - subImg->row0; 126 if (idImg->data.S32[y2][x2] == IN_PEAK) 127 // There's a brighter peak within the footprint above threshold; so cull our initial peak 128 keep = false; 129 } 130 if (!keep) continue; 131 132 psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); 133 133 } 134 134 … … 151 151 */ 152 152 psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint 153 const psImage *weight,// corresponding variance image154 155 156 157 153 const psImage *weight, // corresponding variance image 154 pmFootprint *fp, // Footprint containing mortal peaks 155 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive 156 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance 157 const float min_threshold) { // minimum permitted coll height 158 158 assert (img != NULL); assert (img->type.type == PS_TYPE_F32); 159 159 assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); … … 162 162 163 163 if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do 164 165 } 166 167 psRegion subRegion; 164 return PS_ERR_NONE; 165 } 166 167 psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) 168 168 subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; 169 169 subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; … … 185 185 // 186 186 for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 #if 1 207 208 209 210 i--;// we moved everything down one211 187 const pmPeak *peak = fp->peaks->data[i]; 188 int x = peak->x - subImg->col0; 189 int y = peak->y - subImg->row0; 190 // 191 // Find the level nsigma below the peak that must separate the peak 192 // from any of its friends 193 // 194 assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); 195 196 // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations) 197 // add additional padding in quadrature: 198 199 const float stdev = sqrt(subWt->data.F32[y][x]); 200 const float flux = subImg->data.F32[y][x]; 201 const float fStdev = stdev/flux; 202 const float stdev_pad = flux * hypot(fStdev, fPadding); 203 float threshold = flux - nsigma_delta*stdev_pad; 204 205 if (isnan(threshold) || threshold < min_threshold) { 206 #if 1 // min_threshold is assumed to be below the detection threshold, 207 // so all the peaks are pmFootprint, and this isn't the brightest 208 // XXX mark peak to be dropped 209 (void)psArrayRemoveIndex(fp->peaks, i); 210 i--; // we moved everything down one 211 continue; 212 212 #else 213 213 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once 214 214 threshold = min_threshold; 215 215 #endif 216 217 218 219 220 221 222 223 224 225 226 const int peak_id = 1;// the ID for the peak of interest227 brightPeaks->n = i;// only stop at a peak brighter than we are228 229 230 231 232 brightPeaks->n = 0;// don't double free233 234 235 236 237 238 239 240 241 242 243 244 245 if (peak2_id == peak_id) {// There's a brighter peak within the footprint above246 ;// threshold; so cull our initial peak247 248 i--;// we moved everything down one249 250 251 252 253 254 255 256 216 } 217 218 // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? 219 if (threshold > subImg->data.F32[y][x]) { 220 threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; 221 } 222 223 // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint 224 // perhaps this should alloc a single ID image above and pass it in to be set. 225 226 const int peak_id = 1; // the ID for the peak of interest 227 brightPeaks->n = i; // only stop at a peak brighter than we are 228 229 // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) 230 231 pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); 232 brightPeaks->n = 0; // don't double free 233 psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id); 234 psFree(peakFootprint); 235 236 // Check if any of the previous (brighter) peaks are within the footprint of this peak 237 // If so, the current peak is bogus; drop it. 238 int j; 239 for (j = 0; j < i; j++) { 240 const pmPeak *peak2 = fp->peaks->data[j]; 241 int x2 = peak2->x - subImg->col0; 242 int y2 = peak2->y - subImg->row0; 243 const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak 244 245 if (peak2_id == peak_id) { // There's a brighter peak within the footprint above 246 ; // threshold; so cull our initial peak 247 (void)psArrayRemoveIndex(fp->peaks, i); 248 i--; // we moved everything down one 249 break; 250 } 251 } 252 if (j == i) { 253 j++; 254 } 255 256 psFree(idImg); 257 257 } 258 258 259 259 brightPeaks->n = 0; psFree(brightPeaks); 260 psFree( (psImage *)subImg);261 psFree( (psImage *)subWt);260 psFree(subImg); 261 psFree(subWt); 262 262 263 263 return PS_ERR_NONE;
Note:
See TracChangeset
for help on using the changeset viewer.
