Changeset 9734 for trunk/psphot/src/psphotBasicDeblend.c
- Timestamp:
- Oct 24, 2006, 2:07:03 PM (20 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/psphotBasicDeblend.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/psphotBasicDeblend.c
r9270 r9734 2 2 3 3 // 2006.02.07 : no leaks 4 bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) { 4 bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) { 5 5 6 6 int N; … … 25 25 // we use an index for this so the spatial sorting is kept 26 26 psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32); 27 SN->n = SN->nalloc;28 27 for (int i = 0; i < SN->n; i++) { 29 30 28 source = sources->data[i]; 29 SN->data.F32[i] = source->peak->counts; 31 30 } 32 31 psVector *index = psVectorSortIndex (NULL, SN); … … 35 34 // examine sources in decreasing SN order 36 35 for (int i = sources->n - 1; i >= 0; i--) { 37 38 36 N = index->data.U32[i]; 37 source = sources->data[N]; 39 38 40 39 if (source->mode & PM_SOURCE_MODE_BLEND) continue; 41 40 42 // temporary array for overlapping objects we find 43 psArray *overlap = psArrayAlloc (100); 44 // DROP overlap->n = 0; 41 // temporary array for overlapping objects we find 42 psArray *overlap = psArrayAllocEmpty (100); 45 43 46 47 48 49 50 51 52 53 54 55 56 57 44 // search backwards for overlapping sources 45 for (int j = N - 1; j >= 0; j--) { 46 testSource = sources->data[j]; 47 if (testSource->peak->x < source->pixels->col0) continue; 48 if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; 49 if (testSource->peak->y < source->pixels->row0) break; 50 if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) { 51 fprintf (stderr, "warning: invalid condition\n"); 52 continue; 53 } 54 psArrayAdd (overlap, 100, testSource); 55 } 58 56 59 60 61 62 63 64 65 66 67 68 69 70 57 // search forwards for overlapping sources 58 for (int j = N + 1; j < sources->n; j++) { 59 testSource = sources->data[j]; 60 if (testSource->peak->x < source->pixels->col0) continue; 61 if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; 62 if (testSource->peak->y < source->pixels->row0) { 63 fprintf (stderr, "warning: invalid condition\n"); 64 continue; 65 } 66 if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break; 67 psArrayAdd (overlap, 100, testSource); 68 } 71 69 72 73 74 75 70 if (overlap->n == 0) { 71 psFree (overlap); 72 continue; 73 } 76 74 77 // this source has overlapping neighbors, check for actual blends 78 // generate source contour (1/4 peak counts) 79 // set the threshold based on user inputs 80 81 // threshold is fraction of the source peak flux 82 // image is background subtracted; source->moments->Sky should always be 0.0 83 threshold = FRACTION * source->moments->Peak; 84 // threshold is no less than NSIGMA above the local median sigma? 85 threshold = PS_MAX (threshold, NSIGMA*sqrt(source->moments->dSky)); 75 // this source has overlapping neighbors, check for actual blends 76 // generate source contour (1/4 peak counts) 77 // set the threshold based on user inputs 86 78 87 // generate a basic contour (set of x,y coordinates at-or-below flux level) 88 psArray *contour = pmSourceContour (source->pixels, source->peak->x, source->peak->y, threshold); 89 if (contour == NULL) { 90 psFree (overlap); 91 continue; 92 } 79 // threshold is fraction of the source peak flux 80 // image is background subtracted; source->moments->Sky should always be 0.0 81 threshold = FRACTION * source->moments->Peak; 82 // threshold is no less than NSIGMA above the local median sigma? 83 threshold = PS_MAX (threshold, NSIGMA*sqrt(source->moments->dSky)); 93 84 94 // the source contour consists of two vectors, xv and yv. the contour is 95 // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]). both 96 // coordinate pairs have the same yv[] value, with xv[i] corresponding to 97 // the left boundary and xv[i+1] corresponding to the right boundary 85 // generate a basic contour (set of x,y coordinates at-or-below flux level) 86 psArray *contour = pmSourceContour (source->pixels, source->peak->x, source->peak->y, threshold); 87 if (contour == NULL) { 88 psFree (overlap); 89 continue; 90 } 98 91 99 psVector *xv = contour->data[0]; 100 psVector *yv = contour->data[1]; 101 for (int k = 0; k < overlap->n; k++) { 102 testSource = overlap->data[k]; 103 if (testSource->peak->counts > source->peak->counts) continue; 104 for (int j = 0; j < xv->n; j+=2) { 105 if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue; 106 if (xv->data.F32[j+0] > testSource->peak->x) break; 107 if (xv->data.F32[j+1] < testSource->peak->x) break; 92 // the source contour consists of two vectors, xv and yv. the contour is 93 // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]). both 94 // coordinate pairs have the same yv[] value, with xv[i] corresponding to 95 // the left boundary and xv[i+1] corresponding to the right boundary 108 96 109 testSource->mode |= PM_SOURCE_MODE_BLEND; 97 psVector *xv = contour->data[0]; 98 psVector *yv = contour->data[1]; 99 for (int k = 0; k < overlap->n; k++) { 100 testSource = overlap->data[k]; 101 if (testSource->peak->counts > source->peak->counts) continue; 102 for (int j = 0; j < xv->n; j+=2) { 103 if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue; 104 if (xv->data.F32[j+0] > testSource->peak->x) break; 105 if (xv->data.F32[j+1] < testSource->peak->x) break; 110 106 111 // add this to the list of source->blends 112 if (source->blends == NULL) { 113 source->blends = psArrayAlloc (16); 114 // DROP source->blends->n = 0; 115 } 116 psArrayAdd (source->blends, 16, testSource); 107 testSource->mode |= PM_SOURCE_MODE_BLEND; 117 108 118 Nblend ++; 119 j = xv->n; 120 } 121 } 122 psFree (overlap); 123 psFree (contour); 109 // add this to the list of source->blends 110 if (source->blends == NULL) { 111 source->blends = psArrayAllocEmpty (16); 112 } 113 psArrayAdd (source->blends, 16, testSource); 114 115 Nblend ++; 116 j = xv->n; 117 } 118 } 119 psFree (overlap); 120 psFree (contour); 124 121 } 125 122 psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f sec)\n", Nblend, psTimerMark ("psphot"));
Note:
See TracChangeset
for help on using the changeset viewer.
