1 /*
2 winpos.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: Compute windowed barycenter
11 *
12 * Last modify: 22/09/2005
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
24 #include "define.h"
25 #include "globals.h"
26 #include "prefs.h"
27 #include "winpos.h"
28
29 static obj2struct *obj2 = &outobj2;
30
31 /****** compute_winpos ********************************************************
32 PROTO void compute_winpos(picstruct *field, picstruct *wfield,
33 objstruct *obj)
34 PURPOSE Compute windowed source barycenter.
35 INPUT Picture structure pointer,
36 Weight-map structure pointer,
37 object structure.
38 OUTPUT -.
39 NOTES obj->posx and obj->posy are taken as initial centroid guesses.
40 AUTHOR E. Bertin (IAP)
41 VERSION 22/09/2005
42 ***/
compute_winpos(picstruct * field,picstruct * wfield,objstruct * obj)43 void compute_winpos(picstruct *field, picstruct *wfield, objstruct *obj)
44
45 {
46 float r2, raper,raper2, rintlim,rintlim2,rextlim2,
47 dx,dx1,dy,dy2, sig, pdbkg,
48 offsetx,offsety,scalex,scaley,scale2, ngamma, locarea;
49 double tv, tv2, pix, var, backnoise2, gain, locpix,
50 dxpos,dypos, twosig2, err,err2, emx2,emy2,emxy,
51 esum, temp,temp2, mx2, my2,mxy,pmx2, theta, mx,my,
52 mx2ph, my2ph;
53 int i,x,y, x2,y2, xmin,xmax,ymin,ymax, sx,sy, w,h,
54 fymin,fymax, pflag,corrflag, gainflag, errflag,
55 momentflag;
56 long pos;
57 PIXTYPE *strip,*stript, *wstrip,*wstript,
58 wthresh = 0.0;
59
60 if (wfield)
61 wthresh = wfield->weight_thresh;
62 wstrip = wstript = NULL;
63 w = field->width;
64 h = field->stripheight;
65 fymin = field->ymin;
66 fymax = field->ymax;
67 pflag = (prefs.detect_type==PHOTO)? 1:0;
68 corrflag = (prefs.mask_type==MASK_CORRECT);
69 gainflag = wfield && prefs.weightgain_flag;
70 errflag = FLAG(obj2.winposerr_mx2);
71 momentflag = FLAG(obj2.win_mx2) | FLAG(obj2.winposerr_mx2);
72 var = backnoise2 = field->backsig*field->backsig;
73 gain = prefs.gain;
74 sig = obj2->hl_radius*2.0/2.35; /* From half-FWHM to sigma */
75 twosig2 = 2.0*sig*sig;
76
77 /* Integration radius */
78 raper = WINPOS_NSIG*sig;
79
80 /* For photographic data */
81 if (pflag)
82 {
83 ngamma = field->ngamma;
84 pdbkg = exp(obj->dbkg/ngamma);
85 }
86 else
87 {
88 ngamma = 0.0;
89 pdbkg = 0.0;
90 }
91 raper2 = raper*raper;
92 /* Internal radius of the oversampled annulus (<r-sqrt(2)/2) */
93 rintlim = raper - 0.75;
94 rintlim2 = (rintlim>0.0)? rintlim*rintlim: 0.0;
95 /* External radius of the oversampled annulus (>r+sqrt(2)/2) */
96 rextlim2 = (raper + 0.75)*(raper + 0.75);
97 scaley = scalex = 1.0/WINPOS_OVERSAMP;
98 scale2 = scalex*scaley;
99 offsetx = 0.5*(scalex-1.0);
100 offsety = 0.5*(scaley-1.0);
101 /* Use isophotal centroid as a first guess */
102 mx = obj2->posx - 1.0;
103 my = obj2->posy - 1.0;
104
105 for (i=0; i<WINPOS_NITERMAX; i++)
106 {
107 xmin = (int)(mx-raper+0.499999);
108 xmax = (int)(mx+raper+1.499999);
109 ymin = (int)(my-raper+0.499999);
110 ymax = (int)(my+raper+1.499999);
111 mx2ph = mx*2.0 + 0.49999;
112 my2ph = my*2.0 + 0.49999;
113
114 if (xmin < 0)
115 {
116 xmin = 0;
117 obj->flag |= OBJ_APERT_PB;
118 }
119 if (xmax > w)
120 {
121 xmax = w;
122 obj->flag |= OBJ_APERT_PB;
123 }
124 if (ymin < fymin)
125 {
126 ymin = fymin;
127 obj->flag |= OBJ_APERT_PB;
128 }
129 if (ymax > fymax)
130 {
131 ymax = fymax;
132 obj->flag |= OBJ_APERT_PB;
133 }
134
135 tv = esum = emxy = emx2 = emy2 = mx2 = my2 = mxy = 0.0;
136 dxpos = dypos = 0.0;
137 strip = field->strip;
138 wstrip = wstript = NULL; /* To avoid gcc -Wall warnings */
139 if (wfield)
140 wstrip = wfield->strip;
141 for (y=ymin; y<ymax; y++)
142 {
143 stript = strip + (pos = (y%h)*w + xmin);
144 if (wfield)
145 wstript = wstrip + pos;
146 for (x=xmin; x<xmax; x++, stript++, wstript++)
147 {
148 dx = x - mx;
149 dy = y - my;
150 if ((r2=dx*dx+dy*dy)<rextlim2)
151 {
152 if (WINPOS_OVERSAMP>1 && r2> rintlim2)
153 {
154 dx += offsetx;
155 dy += offsety;
156 locarea = 0.0;
157 for (sy=WINPOS_OVERSAMP; sy--; dy+=scaley)
158 {
159 dx1 = dx;
160 dy2 = dy*dy;
161 for (sx=WINPOS_OVERSAMP; sx--; dx1+=scalex)
162 if (dx1*dx1+dy2<raper2)
163 locarea += scale2;
164 }
165 }
166 else
167 locarea = 1.0;
168 locarea *= exp(-r2/twosig2);
169 /*-------- Here begin tests for pixel and/or weight overflows. Things are a */
170 /*-------- bit intricated to have it running as fast as possible in the most */
171 /*-------- common cases */
172 if ((pix=*stript)<=-BIG || (wfield && (var=*wstript)>=wthresh))
173 {
174 if (corrflag
175 && (x2=(int)(mx2ph-x))>=0 && x2<w
176 && (y2=(int)(my2ph-y))>=fymin && y2<fymax
177 && (pix=*(strip + (pos = (y2%h)*w + x2)))>-BIG)
178 {
179 if (wfield)
180 {
181 var = *(wstrip + pos);
182 if (var>=wthresh)
183 pix = var = 0.0;
184 }
185 }
186 else
187 {
188 pix = 0.0;
189 if (wfield)
190 var = 0.0;
191 }
192 }
193 if (pflag)
194 pix=exp(pix/ngamma);
195 dx = x - mx;
196 dy = y - my;
197 locpix = locarea*pix;
198 tv += locpix;
199 dxpos += locpix*dx;
200 dypos += locpix*dy;
201 if (errflag)
202 {
203 err = var;
204 if (pflag)
205 err *= locpix*pix/(ngamma*ngamma);
206 else if (gain>0.0 && pix>0.0)
207 {
208 if (gainflag)
209 err += pix/gain*var/backnoise2;
210 else
211 err += pix/gain;
212 }
213 err2 = locarea*locarea*err;
214 esum += err2;
215 emx2 += err2*(dx*dx+0.0833); /* Finite pixel size */
216 emy2 += err2*(dy*dy+0.0833); /* Finite pixel size */
217 emxy += err2*dx*dy;
218 }
219 if (momentflag)
220 {
221 mx2 += locpix*dx*dx;
222 my2 += locpix*dy*dy;
223 mxy += locpix*dx*dy;
224 }
225 }
226 }
227 }
228
229 if (tv>0.0)
230 {
231 mx += (dxpos /= tv)*WINPOS_GRADFAC;
232 my += (dypos /= tv)*WINPOS_GRADFAC;
233 }
234 else
235 break;
236
237 /*-- Stop here if position does not change */
238 if (dxpos*dxpos+dypos*dypos < WINPOS_STEPMIN*WINPOS_STEPMIN)
239 break;
240 }
241 mx2 = mx2/tv - dxpos*dxpos;
242 my2 = my2/tv - dypos*dypos;
243 mxy = mxy/tv - dxpos*dypos;
244 obj2->winpos_x = mx + 1.0; /* The dreaded 1.0 FITS offset */
245 obj2->winpos_y = my + 1.0; /* The dreaded 1.0 FITS offset */
246 obj2->winpos_niter = i+1;
247
248 /* WINdowed flux */
249 if (FLAG(obj2.flux_win))
250 {
251 obj2->flux_win = tv;
252 obj2->fluxerr_win = sqrt(esum);
253 }
254 temp2=mx2*my2-mxy*mxy;
255 obj2->win_flag = (tv <= 0.0)*4 + (mx2 < 0.0 || my2 < 0.0)*2
256 + (temp2<0.0);
257 if (obj2->win_flag)
258 {
259 /*--- Negative values: revert to isophotal estimates */
260 if (errflag)
261 {
262 obj2->winposerr_mx2 = obj->poserr_mx2;
263 obj2->winposerr_my2 = obj->poserr_my2;
264 obj2->winposerr_mxy = obj->poserr_mxy;
265 if (FLAG(obj2.winposerr_a))
266 {
267 obj2->winposerr_a = obj2->poserr_a;
268 obj2->winposerr_b = obj2->poserr_b;
269 obj2->winposerr_theta = obj2->poserr_theta;
270 }
271 if (FLAG(obj2.winposerr_cxx))
272 {
273 obj2->winposerr_cxx = obj2->poserr_cxx;
274 obj2->winposerr_cyy = obj2->poserr_cyy;
275 obj2->winposerr_cxy = obj2->poserr_cxy;
276 }
277 }
278 if (momentflag)
279 {
280 obj2->win_mx2 = obj->mx2;
281 obj2->win_my2 = obj->my2;
282 obj2->win_mxy = obj->mxy;
283 if (FLAG(obj2.win_cxx))
284 {
285 obj2->win_cxx = obj->cxx;
286 obj2->win_cyy = obj->cyy;
287 obj2->win_cxy = obj->cxy;
288 }
289 if (FLAG(obj2.win_a))
290 {
291 obj2->win_a = obj->a;
292 obj2->win_b = obj->b;
293 obj2->win_polar = obj2->polar;
294 obj2->win_theta = obj->theta;
295 }
296 }
297 }
298 else
299 {
300 if (errflag)
301 {
302 tv2 = tv*tv;
303 emx2 /= tv2;
304 emy2 /= tv2;
305 emxy /= tv2;
306 /*-- Handle fully correlated profiles (which cause a singularity...) */
307 esum *= 0.08333/tv2;
308 if (obj->singuflag && (emx2*emy2-emxy*emxy) < esum*esum)
309 {
310 emx2 += esum;
311 emy2 += esum;
312 }
313
314 obj2->winposerr_mx2 = emx2;
315 obj2->winposerr_my2 = emy2;
316 obj2->winposerr_mxy = emxy;
317 /*---- Error ellipse parameters */
318 if (FLAG(obj2.winposerr_a))
319 {
320 double pmx2,pmy2,temp,theta;
321
322 if (fabs(temp=emx2-emy2) > 0.0)
323 theta = atan2(2.0 * emxy,temp) / 2.0;
324 else
325 theta = PI/4.0;
326
327 temp = sqrt(0.25*temp*temp+ emxy*emxy);
328 pmy2 = pmx2 = 0.5*(emx2+emy2);
329 pmx2+=temp;
330 pmy2-=temp;
331
332 obj2->winposerr_a = (float)sqrt(pmx2);
333 obj2->winposerr_b = (float)sqrt(pmy2);
334 obj2->winposerr_theta = theta*180.0/PI;
335 }
336
337 if (FLAG(obj2.winposerr_cxx))
338 {
339 double temp;
340
341 obj2->winposerr_cxx = (float)(emy2/(temp=emx2*emy2-emxy*emxy));
342 obj2->winposerr_cyy = (float)(emx2/temp);
343 obj2->winposerr_cxy = (float)(-2*emxy/temp);
344 }
345 }
346
347 if (momentflag)
348 {
349 /*-- Handle fully correlated profiles (which cause a singularity...) */
350 if ((temp2=mx2*my2-mxy*mxy)<0.00694)
351 {
352 mx2 += 0.0833333;
353 my2 += 0.0833333;
354 temp2 = mx2*my2-mxy*mxy;
355 }
356 obj2->win_mx2 = mx2;
357 obj2->win_my2 = my2;
358 obj2->win_mxy = mxy;
359
360 if (FLAG(obj2.win_cxx))
361 {
362 obj2->win_cxx = (float)(my2/temp2);
363 obj2->win_cyy = (float)(mx2/temp2);
364 obj2->win_cxy = (float)(-2*mxy/temp2);
365 }
366
367 if (FLAG(obj2.win_a))
368 {
369 if ((fabs(temp=mx2-my2)) > 0.0)
370 theta = atan2(2.0 * mxy,temp) / 2.0;
371 else
372 theta = PI/4.0;
373
374 temp = sqrt(0.25*temp*temp+mxy*mxy);
375 pmx2 = 0.5*(mx2+my2);
376 obj2->win_a = (float)sqrt(pmx2 + temp);
377 obj2->win_b = (float)sqrt(pmx2 - temp);
378 if (FLAG(obj2.win_polar))
379 obj2->win_polar = temp / pmx2;
380 obj2->win_theta = theta*180.0/PI;
381 }
382 }
383 }
384
385 return;
386 }
387
388
389