ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambcomp.c
Revision: 2.12
Committed: Sun Jul 27 22:12:03 2003 UTC (20 years, 8 months ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R6, rad3R6P1
Changes since 2.11: +2 -2 lines
Log Message:
Added grouping parens to reduce ambiguity warnings.

File Contents

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