1 /*
2 growth.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: Make growth curves.
11 *
12 * Last modify: 15/02/2005
13 *
14 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 */
16
17 #ifdef HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25
26 #include "define.h"
27 #include "globals.h"
28 #include "prefs.h"
29 #include "growth.h"
30
31 /*------------------------------- variables ---------------------------------*/
32
33 static double *growth;
34 static int ngrowth;
35 static obj2struct *obj2 = &outobj2;
36
37 /******************************** initgrowth *********************************/
38 /*
39 Allocate memory for growth curve stuff.
40 */
initgrowth()41 void initgrowth()
42 {
43
44 QMALLOC(growth, double, GROWTH_NSTEP);
45 /* Quick (and dirty) fix to allow FLUX_RADIUS support */
46 if (FLAG(obj2.flux_radius) && !prefs.flux_radiussize)
47 {
48 QCALLOC(obj2->flux_radius, float, 1);
49 }
50
51 return;
52 }
53
54
55 /******************************** endgrowth **********************************/
56 /*
57 Free memory occupied by growth curve stuff.
58 */
endgrowth()59 void endgrowth()
60 {
61 free(growth);
62 if (FLAG(obj2.flux_radius) && !prefs.flux_radiussize)
63 free(obj2->flux_radius);
64
65 return;
66 }
67
68
69 /****************************** makeavergrowth *******************************/
70 /*
71 Build growth curve based on averages.
72 */
makeavergrowth(picstruct * field,picstruct * wfield,objstruct * obj)73 void makeavergrowth(picstruct *field, picstruct *wfield, objstruct *obj)
74
75 {
76 float *fgrowth;
77 double *growtht,
78 dx,dx1,dy,dy2,mx,my, r2,r,rlim, d, rextlim2, raper,
79 offsetx,offsety,scalex,scaley,scale2, ngamma, locarea,
80 tv, sigtv, area, pix, var, gain, dpos,step,step2, dg,
81 stepdens, backnoise2, prevbinmargin, nextbinmargin;
82 int i,j,n, x,y, x2,y2, xmin,xmax,ymin,ymax, sx,sy, w,h,
83 fymin,fymax, pflag,corrflag, ipos;
84 LONG pos;
85 PIXTYPE *strip,*stript, *wstrip,*wstript,
86 pdbkg, wthresh;
87
88
89 if (wfield)
90 wthresh = wfield->weight_thresh;
91 else
92 wthresh = 0.0; /* To avoid gcc -Wall warnings */
93
94 /* Clear the growth-curve buffer */
95 memset(growth, 0, (size_t)(GROWTH_NSTEP*sizeof(double)));
96
97 mx = obj->mx;
98 my = obj->my;
99 w = field->width;
100 h = field->stripheight;
101 fymin = field->ymin;
102 fymax = field->ymax;
103 pflag = (prefs.detect_type==PHOTO)? 1:0;
104 corrflag = (prefs.mask_type==MASK_CORRECT);
105 var = backnoise2 = field->backsig*field->backsig;
106 gain = prefs.gain;
107
108 /* Integration radius */
109 rlim = GROWTH_NSIG*obj->a;
110 if (rlim<prefs.autoaper[0])
111 rlim = prefs.autoaper[0];
112 raper = rlim+1.5; /* margin for interpolation */
113 /* External radius */
114 rextlim2 = raper*raper;
115
116 stepdens = GROWTH_NSTEP/rlim;
117 /* It is useless to oversample too much! */
118 if (0/*stepdens>GROWTH_OVERSAMP*/)
119 {
120 ngrowth = (int)(rlim*GROWTH_OVERSAMP);
121 stepdens = ngrowth/rlim;
122 }
123 else
124 ngrowth = GROWTH_NSTEP;
125 step = 1/stepdens;
126
127 /* Critical distances (in pixels) from bin boundaries */
128 prevbinmargin = 0.75; /* > 1/sqrt(2) */
129 nextbinmargin = step - 0.75; /* <step - sqrt(2) */
130 /* For photographic data */
131 if (pflag)
132 {
133 ngamma = field->ngamma;
134 pdbkg = exp(obj->dbkg/ngamma);
135 }
136 else
137 {
138 ngamma = 0.0;
139 pdbkg = 0.0;
140 }
141 tv = sigtv = area = 0.0;
142 scaley = scalex = 1.0/GROWTH_OVERSAMP;
143 scale2 = scalex*scaley;
144 offsetx = 0.5*(scalex-1.0);
145 offsety = 0.5*(scaley-1.0);
146 xmin = (int)(mx-raper+0.499999);
147 xmax = (int)(mx+raper+1.499999);
148 ymin = (int)(my-raper+0.499999);
149 ymax = (int)(my+raper+1.499999);
150
151 if (xmin < 0)
152 {
153 xmin = 0;
154 obj->flag |= OBJ_APERT_PB;
155 }
156 if (xmax > w)
157 {
158 xmax = w;
159 obj->flag |= OBJ_APERT_PB;
160 }
161 if (ymin < fymin)
162 {
163 ymin = fymin;
164 obj->flag |= OBJ_APERT_PB;
165 }
166 if (ymax > fymax)
167 {
168 ymax = fymax;
169 obj->flag |= OBJ_APERT_PB;
170 }
171
172 strip = field->strip;
173 wstrip = wstript = NULL; /* To avoid gcc -Wall warnings */
174 if (wfield)
175 wstrip = wfield->strip;
176 for (y=ymin; y<ymax; y++)
177 {
178 stript = strip + (pos = (y%h)*w + xmin);
179 if (wfield)
180 wstript = wstrip + pos;
181 for (x=xmin; x<xmax; x++, stript++, wstript++)
182 {
183 dx = x - mx;
184 dy = y - my;
185 if ((r2=dx*dx+dy*dy)<rextlim2)
186 {
187 /*------ Here begin tests for pixel and/or weight overflows. Things are a */
188 /*------ bit intricated to have it running as fast as possible in the most */
189 /*------ common cases */
190 if ((pix=*stript)<=-BIG || (wfield && (var=*wstript)>=wthresh))
191 {
192 if (corrflag
193 && (x2=(int)(2*mx+0.49999-x))>=0 && x2<w
194 && (y2=(int)(2*my+0.49999-y))>=fymin && y2<fymax
195 && (pix=*(strip + (pos = (y2%h)*w + x2)))>-BIG)
196 {
197 if (wfield)
198 {
199 var = *(wstrip + pos);
200 if (var>=wthresh)
201 pix = var = 0.0;
202 }
203 }
204 else
205 {
206 pix = 0.0;
207 if (wfield)
208 var = 0.0;
209 }
210 }
211 if (pflag)
212 pix = exp(pix/ngamma) - pdbkg;
213
214 /*------ Check if oversampling is needed (close enough to a bin boundary) */
215 d = fmod(r=sqrt(r2),step);
216 if (d<prevbinmargin || d>nextbinmargin)
217 {
218 dx += offsetx;
219 dy += offsety;
220 locarea = 0.0;
221 for (sy=GROWTH_OVERSAMP; sy--; dy+=scaley)
222 {
223 dx1 = dx;
224 dy2 = dy*dy;
225 for (sx=GROWTH_OVERSAMP; sx--; dx1+=scalex)
226 {
227 j = (int)(sqrt(dx1*dx1+dy2)*stepdens);
228 if (j<ngrowth)
229 {
230 growth[j] += scale2*pix;
231 locarea += scale2;
232 }
233 }
234 }
235 }
236 else
237 {
238 j = (int)(r*stepdens);
239 if (j<ngrowth)
240 {
241 growth[j] += pix;
242 locarea = 1.0;
243 }
244 }
245 area += locarea;
246 /*
247 if (pflag)
248 sigtv += var*locarea*pix*pix;
249 else
250 sigtv += var*locarea;
251 tv += locarea*pix;
252 if (wfield && pix>0.0 && gain>0.0)
253 sigtv += pix/gain*var/backnoise2;
254 */
255 }
256 }
257 }
258
259 /*
260 if (pflag)
261 {
262 tv = ngamma*(tv-area*exp(obj->dbkg/ngamma));
263 sigtv /= ngamma*ngamma;
264 }
265 else
266 {
267 tv -= area*obj->dbkg;
268 if (!wfield && gain > 0.0 && tv>0.0)
269 sigtv += tv/gain;
270 }
271 */
272
273 /* Integrate the growth curve */
274 pix = 0.0;
275 growtht = growth;
276 for (i=ngrowth; i--;)
277 {
278 *growtht += pix;
279 pix = *(growtht++);
280 }
281
282 /* Now let's remap the growth-curve to match user's choice */
283 if (FLAG(obj2.flux_growth))
284 {
285 n = prefs.flux_growthsize;
286 if (FLAG(obj2.flux_growthstep))
287 obj2->flux_growthstep = rlim/n;
288 fgrowth = obj2->flux_growth;
289 step2 = (double)GROWTH_NSTEP/n;
290 j = 1;
291 for (i=n; i--; j++)
292 {
293 ipos = (int)(dpos=step2*j-0.99999);
294 if (dpos<0.0)
295 *(fgrowth++) = (float)(*growth*(0.99999+dpos));
296 else
297 {
298 growtht = growth + ipos;
299 *(fgrowth++) = (float)(*growtht+(*(growtht+1)-*growtht)*(dpos-ipos));
300 }
301 }
302 }
303
304 if (FLAG(obj2.mag_growth))
305 {
306 n = prefs.mag_growthsize;
307 if (FLAG(obj2.mag_growthstep))
308 obj2->mag_growthstep = rlim/n;
309 fgrowth = obj2->mag_growth;
310 step2 = (double)GROWTH_NSTEP/n;
311 j = 1;
312 for (i=n; i--; j++)
313 {
314 ipos = (int)(dpos=step2*j-0.99999);
315 if (dpos<0.0)
316 pix = *growth*(0.99999+dpos);
317 else
318 {
319 growtht = growth + ipos;
320 pix = *growtht+(*(growtht+1)-*growtht)*(dpos-ipos);
321 }
322 *(fgrowth++) = pix>0.0?(prefs.mag_zeropoint-2.5*log10(pix)):99.0;
323 }
324 }
325
326 if (FLAG(obj2.flux_radius))
327 {
328 n = ngrowth-1;
329 for (j=0; j<prefs.nflux_frac; j++)
330 {
331 tv = prefs.flux_frac[j]*obj2->flux_auto;
332 growtht = growth-1;
333 for (i=0; i<n && *(++growtht)<tv; i++);
334 obj2->flux_radius[j] = step
335 *(i? ((dg=*growtht - *(growtht-1)) != 0.0 ?
336 i + (tv - *(growtht-1))/dg
337 : i)
338 : (*growth !=0.0 ?tv/(*growth) : 0.0));
339 }
340 }
341
342 /* Specific to Half-light radius used by other parameters */
343 if (FLAG(obj2.hl_radius))
344 {
345 n = ngrowth-1;
346 tv = 0.5*obj2->flux_auto;
347 growtht = growth-1;
348 for (i=0; i<n && *(++growtht)<tv; i++);
349 obj2->hl_radius = step*(i? ((dg=*growtht - *(growtht-1)) != 0.0 ?
350 i + (tv - *(growtht-1))/dg
351 : i)
352 : (*growth !=0.0 ?tv/(*growth) : 0.0));
353 }
354
355 return;
356 }
357
358
359
360