ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambcomp.c
Revision: 2.21
Committed: Mon Dec 31 18:19:42 2007 UTC (16 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0, rad3R9
Changes since 2.20: +5 -2 lines
Log Message:
Fixed memory leak for certain ambient calculations (thanks to D. Brainard)

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: ambcomp.c,v 2.20 2005/10/28 16:16:33 greg Exp $";
3 #endif
4 /*
5 * Routines to compute "ambient" values using Monte Carlo
6 *
7 * Declarations of external symbols in ambient.h
8 */
9
10 #include "copyright.h"
11
12 #include "ray.h"
13
14 #include "ambient.h"
15
16 #include "random.h"
17
18
19 void
20 inithemi( /* initialize sampling hemisphere */
21 register AMBHEMI *hp,
22 COLOR ac,
23 RAY *r,
24 double wt
25 )
26 {
27 double d;
28 register int i;
29 /* set number of divisions */
30 if (ambacc <= FTINY &&
31 wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
32 wt = d; /* avoid ray termination */
33 hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
34 i = ambacc > FTINY ? 3 : 1; /* minimum number of samples */
35 if (hp->nt < i)
36 hp->nt = i;
37 hp->np = PI * hp->nt + 0.5;
38 /* set number of super-samples */
39 hp->ns = ambssamp * wt + 0.5;
40 /* assign coefficient */
41 copycolor(hp->acoef, ac);
42 d = 1.0/(hp->nt*hp->np);
43 scalecolor(hp->acoef, d);
44 /* make axes */
45 VCOPY(hp->uz, r->ron);
46 hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
47 for (i = 0; i < 3; i++)
48 if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
49 break;
50 if (i >= 3)
51 error(CONSISTENCY, "bad ray direction in inithemi");
52 hp->uy[i] = 1.0;
53 fcross(hp->ux, hp->uy, hp->uz);
54 normalize(hp->ux);
55 fcross(hp->uy, hp->uz, hp->ux);
56 }
57
58
59 int
60 divsample( /* sample a division */
61 register AMBSAMP *dp,
62 AMBHEMI *h,
63 RAY *r
64 )
65 {
66 RAY ar;
67 int hlist[3];
68 double spt[2];
69 double xd, yd, zd;
70 double b2;
71 double phi;
72 register int i;
73 /* ambient coefficient for weight */
74 if (ambacc > FTINY)
75 setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
76 else
77 copycolor(ar.rcoef, h->acoef);
78 if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
79 return(-1);
80 if (ambacc > FTINY) {
81 multcolor(ar.rcoef, h->acoef);
82 scalecolor(ar.rcoef, 1./AVGREFL);
83 }
84 hlist[0] = r->rno;
85 hlist[1] = dp->t;
86 hlist[2] = dp->p;
87 multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
88 zd = sqrt((dp->t + spt[0])/h->nt);
89 phi = 2.0*PI * (dp->p + spt[1])/h->np;
90 xd = tcos(phi) * zd;
91 yd = tsin(phi) * zd;
92 zd = sqrt(1.0 - zd*zd);
93 for (i = 0; i < 3; i++)
94 ar.rdir[i] = xd*h->ux[i] +
95 yd*h->uy[i] +
96 zd*h->uz[i];
97 dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
98 rayvalue(&ar);
99 ndims--;
100 multcolor(ar.rcol, ar.rcoef); /* apply coefficient */
101 addcolor(dp->v, ar.rcol);
102 /* use rt to improve gradient calc */
103 if (ar.rt > FTINY && ar.rt < FHUGE)
104 dp->r += 1.0/ar.rt;
105 /* (re)initialize error */
106 if (dp->n++) {
107 b2 = bright(dp->v)/dp->n - bright(ar.rcol);
108 b2 = b2*b2 + dp->k*((dp->n-1)*(dp->n-1));
109 dp->k = b2/(dp->n*dp->n);
110 } else
111 dp->k = 0.0;
112 return(0);
113 }
114
115
116 static int
117 ambcmp( /* decreasing order */
118 const void *p1,
119 const void *p2
120 )
121 {
122 const AMBSAMP *d1 = (const AMBSAMP *)p1;
123 const AMBSAMP *d2 = (const AMBSAMP *)p2;
124
125 if (d1->k < d2->k)
126 return(1);
127 if (d1->k > d2->k)
128 return(-1);
129 return(0);
130 }
131
132
133 static int
134 ambnorm( /* standard order */
135 const void *p1,
136 const void *p2
137 )
138 {
139 const AMBSAMP *d1 = (const AMBSAMP *)p1;
140 const AMBSAMP *d2 = (const AMBSAMP *)p2;
141 register int c;
142
143 if ( (c = d1->t - d2->t) )
144 return(c);
145 return(d1->p - d2->p);
146 }
147
148
149 double
150 doambient( /* compute ambient component */
151 COLOR acol,
152 RAY *r,
153 double wt,
154 FVECT pg,
155 FVECT dg
156 )
157 {
158 double b, d;
159 AMBHEMI hemi;
160 AMBSAMP *div;
161 AMBSAMP dnew;
162 register AMBSAMP *dp;
163 double arad;
164 int divcnt;
165 register int i, j;
166 /* initialize hemisphere */
167 inithemi(&hemi, acol, r, wt);
168 divcnt = hemi.nt * hemi.np;
169 /* initialize */
170 if (pg != NULL)
171 pg[0] = pg[1] = pg[2] = 0.0;
172 if (dg != NULL)
173 dg[0] = dg[1] = dg[2] = 0.0;
174 setcolor(acol, 0.0, 0.0, 0.0);
175 if (divcnt == 0)
176 return(0.0);
177 /* allocate super-samples */
178 if (hemi.ns > 0 || pg != NULL || dg != NULL) {
179 div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP));
180 if (div == NULL)
181 error(SYSTEM, "out of memory in doambient");
182 } else
183 div = NULL;
184 /* sample the divisions */
185 arad = 0.0;
186 if ((dp = div) == NULL)
187 dp = &dnew;
188 divcnt = 0;
189 for (i = 0; i < hemi.nt; i++)
190 for (j = 0; j < hemi.np; j++) {
191 dp->t = i; dp->p = j;
192 setcolor(dp->v, 0.0, 0.0, 0.0);
193 dp->r = 0.0;
194 dp->n = 0;
195 if (divsample(dp, &hemi, r) < 0) {
196 if (div != NULL)
197 dp++;
198 continue;
199 }
200 arad += dp->r;
201 divcnt++;
202 if (div != NULL)
203 dp++;
204 else
205 addcolor(acol, dp->v);
206 }
207 if (!divcnt) {
208 if (div != NULL)
209 free((void *)div);
210 return(0.0); /* no samples taken */
211 }
212 if (divcnt < hemi.nt*hemi.np) {
213 pg = dg = NULL; /* incomplete sampling */
214 hemi.ns = 0;
215 } else if (arad > FTINY && divcnt/arad < minarad) {
216 hemi.ns = 0; /* close enough */
217 } else if (hemi.ns > 0) { /* else perform super-sampling? */
218 comperrs(div, &hemi); /* compute errors */
219 qsort(div, divcnt, sizeof(AMBSAMP), ambcmp); /* sort divs */
220 /* super-sample */
221 for (i = hemi.ns; i > 0; i--) {
222 dnew = *div;
223 if (divsample(&dnew, &hemi, r) < 0) {
224 dp++;
225 continue;
226 }
227 dp = div; /* reinsert */
228 j = divcnt < i ? divcnt : i;
229 while (--j > 0 && dnew.k < dp[1].k) {
230 *dp = *(dp+1);
231 dp++;
232 }
233 *dp = dnew;
234 }
235 if (pg != NULL || dg != NULL) /* restore order */
236 qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
237 }
238 /* compute returned values */
239 if (div != NULL) {
240 arad = 0.0; /* note: divcnt may be < nt*np */
241 for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
242 arad += dp->r;
243 if (dp->n > 1) {
244 b = 1.0/dp->n;
245 scalecolor(dp->v, b);
246 dp->r *= b;
247 dp->n = 1;
248 }
249 addcolor(acol, dp->v);
250 }
251 b = bright(acol);
252 if (b > FTINY) {
253 b = 1.0/b; /* compute & normalize gradient(s) */
254 if (pg != NULL) {
255 posgradient(pg, div, &hemi);
256 for (i = 0; i < 3; i++)
257 pg[i] *= b;
258 }
259 if (dg != NULL) {
260 dirgradient(dg, div, &hemi);
261 for (i = 0; i < 3; i++)
262 dg[i] *= b;
263 }
264 }
265 free((void *)div);
266 }
267 if (arad <= FTINY)
268 arad = maxarad;
269 else
270 arad = (divcnt+hemi.ns)/arad;
271 if (pg != NULL) { /* reduce radius if gradient large */
272 d = DOT(pg,pg);
273 if (d*arad*arad > 1.0)
274 arad = 1.0/sqrt(d);
275 }
276 if (arad < minarad) {
277 arad = minarad;
278 if (pg != NULL && d*arad*arad > 1.0) { /* cap gradient */
279 d = 1.0/arad/sqrt(d);
280 for (i = 0; i < 3; i++)
281 pg[i] *= d;
282 }
283 }
284 if ((arad /= sqrt(wt)) > maxarad)
285 arad = maxarad;
286 return(arad);
287 }
288
289
290 void
291 comperrs( /* compute initial error estimates */
292 AMBSAMP *da, /* assumes standard ordering */
293 register AMBHEMI *hp
294 )
295 {
296 double b, b2;
297 int i, j;
298 register AMBSAMP *dp;
299 /* sum differences from neighbors */
300 dp = da;
301 for (i = 0; i < hp->nt; i++)
302 for (j = 0; j < hp->np; j++) {
303 #ifdef DEBUG
304 if (dp->t != i || dp->p != j)
305 error(CONSISTENCY,
306 "division order in comperrs");
307 #endif
308 b = bright(dp[0].v);
309 if (i > 0) { /* from above */
310 b2 = bright(dp[-hp->np].v) - b;
311 b2 *= b2 * 0.25;
312 dp[0].k += b2;
313 dp[-hp->np].k += b2;
314 }
315 if (j > 0) { /* from behind */
316 b2 = bright(dp[-1].v) - b;
317 b2 *= b2 * 0.25;
318 dp[0].k += b2;
319 dp[-1].k += b2;
320 } else { /* around */
321 b2 = bright(dp[hp->np-1].v) - b;
322 b2 *= b2 * 0.25;
323 dp[0].k += b2;
324 dp[hp->np-1].k += b2;
325 }
326 dp++;
327 }
328 /* divide by number of neighbors */
329 dp = da;
330 for (j = 0; j < hp->np; j++) /* top row */
331 (dp++)->k *= 1.0/3.0;
332 if (hp->nt < 2)
333 return;
334 for (i = 1; i < hp->nt-1; i++) /* central region */
335 for (j = 0; j < hp->np; j++)
336 (dp++)->k *= 0.25;
337 for (j = 0; j < hp->np; j++) /* bottom row */
338 (dp++)->k *= 1.0/3.0;
339 }
340
341
342 void
343 posgradient( /* compute position gradient */
344 FVECT gv,
345 AMBSAMP *da, /* assumes standard ordering */
346 register AMBHEMI *hp
347 )
348 {
349 register int i, j;
350 double nextsine, lastsine, b, d;
351 double mag0, mag1;
352 double phi, cosp, sinp, xd, yd;
353 register AMBSAMP *dp;
354
355 xd = yd = 0.0;
356 for (j = 0; j < hp->np; j++) {
357 dp = da + j;
358 mag0 = mag1 = 0.0;
359 lastsine = 0.0;
360 for (i = 0; i < hp->nt; i++) {
361 #ifdef DEBUG
362 if (dp->t != i || dp->p != j)
363 error(CONSISTENCY,
364 "division order in posgradient");
365 #endif
366 b = bright(dp->v);
367 if (i > 0) {
368 d = dp[-hp->np].r;
369 if (dp[0].r > d) d = dp[0].r;
370 /* sin(t)*cos(t)^2 */
371 d *= lastsine * (1.0 - (double)i/hp->nt);
372 mag0 += d*(b - bright(dp[-hp->np].v));
373 }
374 nextsine = sqrt((double)(i+1)/hp->nt);
375 if (j > 0) {
376 d = dp[-1].r;
377 if (dp[0].r > d) d = dp[0].r;
378 mag1 += d * (nextsine - lastsine) *
379 (b - bright(dp[-1].v));
380 } else {
381 d = dp[hp->np-1].r;
382 if (dp[0].r > d) d = dp[0].r;
383 mag1 += d * (nextsine - lastsine) *
384 (b - bright(dp[hp->np-1].v));
385 }
386 dp += hp->np;
387 lastsine = nextsine;
388 }
389 mag0 *= 2.0*PI / hp->np;
390 phi = 2.0*PI * (double)j/hp->np;
391 cosp = tcos(phi); sinp = tsin(phi);
392 xd += mag0*cosp - mag1*sinp;
393 yd += mag0*sinp + mag1*cosp;
394 }
395 for (i = 0; i < 3; i++)
396 gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI;
397 }
398
399
400 void
401 dirgradient( /* compute direction gradient */
402 FVECT gv,
403 AMBSAMP *da, /* assumes standard ordering */
404 register AMBHEMI *hp
405 )
406 {
407 register int i, j;
408 double mag;
409 double phi, xd, yd;
410 register AMBSAMP *dp;
411
412 xd = yd = 0.0;
413 for (j = 0; j < hp->np; j++) {
414 dp = da + j;
415 mag = 0.0;
416 for (i = 0; i < hp->nt; i++) {
417 #ifdef DEBUG
418 if (dp->t != i || dp->p != j)
419 error(CONSISTENCY,
420 "division order in dirgradient");
421 #endif
422 /* tan(t) */
423 mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
424 dp += hp->np;
425 }
426 phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
427 xd += mag * tcos(phi);
428 yd += mag * tsin(phi);
429 }
430 for (i = 0; i < 3; i++)
431 gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
432 }