Changeset 26893 for trunk/psModules/src/objects/pmFootprintCullPeaks.c
- Timestamp:
- Feb 10, 2010, 7:34:39 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmFootprintCullPeaks.c
r24888 r26893 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 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 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 // If this peak has not already been assigned to a source, then we can look for any 120 // brighter peaks within its footprint. Check if any of the previous (brighter) peaks 121 // are within the footprint of this peak If so, the current peak is bogus; drop it. 122 bool keep = true; 123 for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) { 124 const pmPeak *peak2 = fp->peaks->data[j]; 125 int x2 = peak2->x - subImg->col0; 126 int y2 = peak2->y - subImg->row0; 127 if (idImg->data.S32[y2][x2] == IN_PEAK) { 128 // There's a brighter peak within the footprint above threshold; so cull our initial peak 129 keep = false; 130 } 131 } 132 if (!keep) { 133 psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak"); // we should not drop any already assigned peaks. 92 134 continue; 93 135 } 94 136 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]); 137 psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); 133 138 } 134 139 … … 151 156 */ 152 157 psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint 153 const psImage *weight,// corresponding variance image154 155 156 157 158 const psImage *weight, // corresponding variance image 159 pmFootprint *fp, // Footprint containing mortal peaks 160 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive 161 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance 162 const float min_threshold) { // minimum permitted coll height 158 163 assert (img != NULL); assert (img->type.type == PS_TYPE_F32); 159 164 assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); … … 162 167 163 168 if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do 164 165 } 166 167 psRegion subRegion; 169 return PS_ERR_NONE; 170 } 171 172 psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) 168 173 subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; 169 174 subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; … … 185 190 // 186 191 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 192 const pmPeak *peak = fp->peaks->data[i]; 193 int x = peak->x - subImg->col0; 194 int y = peak->y - subImg->row0; 195 // 196 // Find the level nsigma below the peak that must separate the peak 197 // from any of its friends 198 // 199 assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); 200 201 // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations) 202 // add additional padding in quadrature: 203 204 const float stdev = sqrt(subWt->data.F32[y][x]); 205 const float flux = subImg->data.F32[y][x]; 206 const float fStdev = stdev/flux; 207 const float stdev_pad = flux * hypot(fStdev, fPadding); 208 float threshold = flux - nsigma_delta*stdev_pad; 209 210 if (isnan(threshold) || threshold < min_threshold) { 211 #if 1 // min_threshold is assumed to be below the detection threshold, 212 // so all the peaks are pmFootprint, and this isn't the brightest 213 // XXX mark peak to be dropped 214 (void)psArrayRemoveIndex(fp->peaks, i); 215 i--; // we moved everything down one 216 continue; 212 217 #else 213 218 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once 214 219 threshold = min_threshold; 215 220 #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 221 } 222 223 // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? 224 if (threshold > subImg->data.F32[y][x]) { 225 threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; 226 } 227 228 // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint 229 // perhaps this should alloc a single ID image above and pass it in to be set. 230 231 const int peak_id = 1; // the ID for the peak of interest 232 brightPeaks->n = i; // only stop at a peak brighter than we are 233 234 // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) 235 236 pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); 237 brightPeaks->n = 0; // don't double free 238 psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id); 239 psFree(peakFootprint); 240 241 // Check if any of the previous (brighter) peaks are within the footprint of this peak 242 // If so, the current peak is bogus; drop it. 243 int j; 244 for (j = 0; j < i; j++) { 245 const pmPeak *peak2 = fp->peaks->data[j]; 246 int x2 = peak2->x - subImg->col0; 247 int y2 = peak2->y - subImg->row0; 248 const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak 249 250 if (peak2_id == peak_id) { // There's a brighter peak within the footprint above 251 ; // threshold; so cull our initial peak 252 (void)psArrayRemoveIndex(fp->peaks, i); 253 i--; // we moved everything down one 254 break; 255 } 256 } 257 if (j == i) { 258 j++; 259 } 260 261 psFree(idImg); 257 262 } 258 263 259 264 brightPeaks->n = 0; psFree(brightPeaks); 260 psFree( (psImage *)subImg);261 psFree( (psImage *)subWt);265 psFree(subImg); 266 psFree(subWt); 262 267 263 268 return PS_ERR_NONE;
Note:
See TracChangeset
for help on using the changeset viewer.
