1 /*
2 refine.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: functions to refine extraction of objects.
11 *
12 * Last modify: 27/11/2003
13 *
14 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 */
16
17 #ifdef HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20
21 #include <math.h>
22 #include <stdlib.h>
23 #include <string.h>
24
25 #include "define.h"
26 #include "globals.h"
27 #include "prefs.h"
28 #include "plist.h"
29 #include "extract.h"
30
31 #ifndef RAND_MAX
32 #define RAND_MAX 2147483647
33 #endif
34 #define NSONMAX 1024 /* max. number per level */
35 #define NBRANCH 16 /* starting number per branch */
36
37 static objliststruct *objlist;
38 static short *son, *ok;
39
40 /******************************** parcelout **********************************
41 PROTO parcelout(objliststruct *objlistin, objliststruct *objlistout)
42 PURPOSE Divide a list of isophotal detections in several parts (deblending).
43 INPUT input objlist,
44 output objlist,
45 OUTPUT RETURN_OK if success, RETURN_FATAL_ERROR otherwise (memory overflow).
46 NOTES Even if the object is not deblended, the output objlist threshold is
47 recomputed if a variable threshold is used.
48 AUTHOR E. Bertin (IAP, Leiden & ESO)
49 VERSION 15/03/98
50 ***/
parcelout(objliststruct * objlistin,objliststruct * objlistout)51 int parcelout(objliststruct *objlistin, objliststruct *objlistout)
52
53 {
54 objstruct *obj;
55 static objliststruct debobjlist, debobjlist2;
56 double dthresh, dthresh0, value0;
57 int h,i,j,k,l,m,
58 xn,
59 nbm = NBRANCH,
60 out;
61
62
63 out = RETURN_OK;
64
65 xn = prefs.deblend_nthresh;
66
67 /* ---- initialize lists of objects */
68
69 debobjlist.obj = debobjlist2.obj = NULL;
70 debobjlist.plist = debobjlist2.plist = NULL;
71 debobjlist.nobj = debobjlist2.nobj = 0;
72 debobjlist.npix = debobjlist2.npix = 0;
73 objlistout->thresh = debobjlist2.thresh = objlistin->thresh;
74 memset(objlist, 0, (size_t)xn*sizeof(objliststruct));
75
76 for (l=0; l<objlistin->nobj && out==RETURN_OK; l++)
77 {
78 dthresh0 = objlistin->obj[l].dthresh;
79
80 objlistout->dthresh = debobjlist2.dthresh = dthresh0;
81 if ((out = addobj(l, objlistin, &objlist[0])) == RETURN_FATAL_ERROR)
82 goto exit_parcelout;
83 if ((out = addobj(l, objlistin, &debobjlist2)) == RETURN_FATAL_ERROR)
84 goto exit_parcelout;
85 value0 = objlist[0].obj[0].fdflux*prefs.deblend_mincont;
86 ok[0] = (short)1;
87 for (k=1; k<xn; k++)
88 {
89 /*------ Calculate threshold */
90 dthresh = objlistin->obj[l].fdpeak;
91 if (dthresh>0.0)
92 {
93 if (prefs.detect_type == PHOTO)
94 debobjlist.dthresh= dthresh0 + (dthresh-dthresh0) * (double)k/xn;
95 else
96 debobjlist.dthresh = dthresh0 * pow(dthresh/dthresh0,(double)k/xn);
97 }
98 else
99 debobjlist.dthresh = dthresh0;
100
101 /*--------- Build tree (bottom->up) */
102 if (objlist[k-1].nobj>=NSONMAX)
103 {
104 out = RETURN_FATAL_ERROR;
105 goto exit_parcelout;
106 }
107
108 for (i=0; i<objlist[k-1].nobj; i++)
109 {
110 if ((out=lutz(objlistin,l,&objlist[k-1].obj[i], &debobjlist))
111 ==RETURN_FATAL_ERROR)
112 goto exit_parcelout;
113
114 for (j=h=0; j<debobjlist.nobj; j++)
115 if (belong(j, &debobjlist, i, &objlist[k-1]))
116 {
117 debobjlist.obj[j].dthresh = debobjlist.dthresh;
118 m = addobj(j, &debobjlist, &objlist[k]);
119 if (m==RETURN_FATAL_ERROR || m>=NSONMAX)
120 {
121 out = RETURN_FATAL_ERROR;
122 goto exit_parcelout;
123 }
124 if (h>=nbm-1)
125 if (!(son = (short *)realloc(son,
126 xn*NSONMAX*(nbm+=16)*sizeof(short))))
127 {
128 out = RETURN_FATAL_ERROR;
129 goto exit_parcelout;
130 }
131 son[k-1+xn*(i+NSONMAX*(h++))] = (short)m;
132 ok[k+xn*m] = (short)1;
133 }
134 son[k-1+xn*(i+NSONMAX*h)] = (short)-1;
135 }
136 }
137
138 /*------- cut the right branches (top->down) */
139
140 for (k = xn-2; k>=0; k--)
141 {
142 obj = objlist[k+1].obj;
143 for (i=0; i<objlist[k].nobj; i++)
144 {
145 for (m=h=0; (j=(int)son[k+xn*(i+NSONMAX*h)])!=-1; h++)
146 {
147 if (obj[j].fdflux - obj[j].dthresh*obj[j].fdnpix > value0)
148 m++;
149 ok[k+xn*i] &= ok[k+1+xn*j];
150 }
151 if (m>1)
152 {
153 for (h=0; (j=(int)son[k+xn*(i+NSONMAX*h)])!=-1; h++)
154 if (ok[k+1+xn*j] && obj[j].fdflux-obj[j].dthresh*obj[j].fdnpix
155 > value0)
156 {
157 objlist[k+1].obj[j].flag |= OBJ_MERGED /* Merge flag on */
158 | ((OBJ_ISO_PB|OBJ_APERT_PB|OBJ_OVERFLOW)
159 &debobjlist2.obj[0].flag);
160 if ((out = addobj(j, &objlist[k+1], &debobjlist2))
161 == RETURN_FATAL_ERROR)
162 goto exit_parcelout;
163 }
164 ok[k+xn*i] = (short)0;
165 }
166 }
167 }
168
169 if (ok[0])
170 out = addobj(0, &debobjlist2, objlistout);
171 else
172 out = gatherup(&debobjlist2, objlistout);
173
174 exit_parcelout:
175
176 free(debobjlist2.obj);
177 free(debobjlist2.plist);
178
179 for (k=0; k<xn; k++)
180 {
181 free(objlist[k].obj);
182 free(objlist[k].plist);
183 }
184 }
185
186 free(debobjlist.obj);
187 free(debobjlist.plist);
188
189 return out;
190 }
191
192 /******************************* allocparcelout ******************************/
193 /*
194 Allocate the memory allocated by global pointers in refine.c
195 */
allocparcelout(void)196 void allocparcelout(void)
197 {
198 QMALLOC(son, short, prefs.deblend_nthresh*NSONMAX*NBRANCH);
199 QMALLOC(ok, short, prefs.deblend_nthresh*NSONMAX);
200 QMALLOC(objlist, objliststruct, prefs.deblend_nthresh);
201
202 return;
203 }
204
205 /******************************* freeparcelout *******************************/
206 /*
207 Free the memory allocated by global pointers in refine.c
208 */
freeparcelout(void)209 void freeparcelout(void)
210 {
211 QFREE(son);
212 QFREE(ok);
213 QFREE(objlist);
214 return;
215 }
216
217 /********************************* gatherup **********************************/
218 /*
219 Collect faint remaining pixels and allocate them to their most probable
220 progenitor.
221 */
gatherup(objliststruct * objlistin,objliststruct * objlistout)222 int gatherup(objliststruct *objlistin, objliststruct *objlistout)
223
224 {
225 char *bmp;
226 float *amp, *p, dx,dy, drand, dist, distmin;
227 objstruct *objin = objlistin->obj, *objout, *objt;
228
229 pliststruct *pixelin = objlistin->plist, *pixelout, *pixt,*pixt2;
230
231 int i,k,l, *n, iclst, npix, bmwidth,
232 nobj = objlistin->nobj, xs,ys, x,y, out;
233
234 out = RETURN_OK;
235
236 objlistout->dthresh = objlistin->dthresh;
237 objlistout->thresh = objlistin->thresh;
238
239 QMALLOC(amp, float, nobj);
240 QMALLOC(p, float, nobj);
241 QMALLOC(n, int, nobj);
242
243 for (i=1; i<nobj; i++)
244 preanalyse(i, objlistin, ANALYSE_FULL);
245
246 p[0] = 0.0;
247 bmwidth = objin->xmax - (xs=objin->xmin) + 1;
248 npix = bmwidth * (objin->ymax - (ys=objin->ymin) + 1);
249 if (!(bmp = (char *)calloc(1, npix*sizeof(char))))
250 {
251 bmp = 0;
252 out = RETURN_FATAL_ERROR;
253 goto exit_gatherup;
254 }
255
256 for (objt = objin+(i=1); i<nobj; i++, objt++)
257 {
258 /*-- Now we have passed the deblending section, reset thresholds */
259 objt->dthresh = objlistin->dthresh;
260 objt->thresh = objlistin->thresh;
261
262 /* ------------ flag pixels which are already allocated */
263 for (pixt=pixelin+objin[i].firstpix; pixt>=pixelin;
264 pixt=pixelin+PLIST(pixt,nextpix))
265 bmp[(PLIST(pixt,x)-xs) + (PLIST(pixt,y)-ys)*bmwidth] = '\1';
266
267 if ((n[i] = addobj(i, objlistin, objlistout)) == RETURN_FATAL_ERROR)
268 {
269 out = RETURN_FATAL_ERROR;
270 goto exit_gatherup;
271 }
272 dist = objt->fdnpix/(2*PI*objt->abcor*objt->a*objt->b);
273 amp[i] = dist<70.0? objt->thresh*exp(dist) : 4.0*objt->fdpeak;
274
275 /* ------------ limitate expansion ! */
276 if (amp[i]>4.0*objt->fdpeak)
277 amp[i] = 4.0*objt->fdpeak;
278 }
279
280 objout = objlistout->obj; /* DO NOT MOVE !!! */
281
282 if (!(pixelout=(pliststruct *)realloc(objlistout->plist,
283 (objlistout->npix + npix)*plistsize)))
284 {
285 out = RETURN_FATAL_ERROR;
286 goto exit_gatherup;
287 }
288
289 objlistout->plist = pixelout;
290 k = objlistout->npix;
291 iclst = 0; /* To avoid gcc -Wall warnings */
292 for (pixt=pixelin+objin->firstpix; pixt>=pixelin;
293 pixt=pixelin+PLIST(pixt,nextpix))
294 {
295 x = PLIST(pixt,x);
296 y = PLIST(pixt,y);
297 if (!bmp[(x-xs) + (y-ys)*bmwidth])
298 {
299 pixt2 = pixelout + (l=(k++*plistsize));
300 memcpy(pixt2, pixt, (size_t)plistsize);
301 PLIST(pixt2, nextpix) = -1;
302 distmin = 1e+31;
303 for (objt = objin+(i=1); i<nobj; i++, objt++)
304 {
305 dx = x - objt->mx;
306 dy = y - objt->my;
307 dist=0.5*(objt->cxx*dx*dx+objt->cyy*dy*dy+objt->cxy*dx*dy)/objt->abcor;
308 p[i] = p[i-1] + (dist<70.0?amp[i]*exp(-dist) : 0.0);
309 if (dist<distmin)
310 {
311 distmin = dist;
312 iclst = i;
313 }
314 }
315 if (p[nobj-1] > 1.0e-31)
316 {
317 drand = p[nobj-1]*rand()/RAND_MAX;
318 for (i=1; p[i]<drand; i++);
319 }
320 else
321 i = iclst;
322 objout[n[i]].lastpix=PLIST(pixelout+objout[n[i]].lastpix,nextpix)=l;
323 }
324 }
325
326 objlistout->npix = k;
327 if (!(objlistout->plist = (pliststruct *)realloc(pixelout,
328 objlistout->npix*plistsize)))
329 error (-1, "Not enough memory to update pixel list in ", "gatherup()");
330
331 exit_gatherup:
332
333 free(bmp);
334 free(amp);
335 free(p);
336 free(n);
337
338 return out;
339 }
340
341