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