ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 1.11
Committed: Fri Sep 7 08:17:42 1990 UTC (33 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.10: +9 -2 lines
Log Message:
Added check for ambset overflow

File Contents

# User Rev Content
1 greg 1.1 /* Copyright (c) 1986 Regents of the University of California */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * ambient.c - routines dealing with ambient (inter-reflected) component.
9     *
10     * The macro AMBFLUSH (if defined) is the number of ambient values
11     * to wait before flushing to the ambient file.
12     *
13     * 5/9/86
14     */
15    
16     #include "ray.h"
17    
18     #include "octree.h"
19    
20 greg 1.8 #include "otypes.h"
21    
22 greg 1.1 #include "random.h"
23    
24     #define OCTSCALE 0.5 /* ceil((valid rad.)/(cube size)) */
25    
26     extern CUBE thescene; /* contains space boundaries */
27    
28     extern COLOR ambval; /* global ambient component */
29     extern double ambacc; /* ambient accuracy */
30     extern int ambres; /* ambient resolution */
31     extern int ambdiv; /* number of divisions for calculation */
32     extern int ambssamp; /* number of super-samples */
33     extern int ambounce; /* number of ambient bounces */
34     extern char *amblist[]; /* ambient include/exclude list */
35     extern int ambincl; /* include == 1, exclude == 0 */
36    
37 greg 1.11 #define MAXASET 511 /* maximum number of elements in ambient set */
38     OBJECT ambset[MAXASET+1]={0}; /* ambient include/exclude set */
39 greg 1.1
40     double maxarad; /* maximum ambient radius */
41     double minarad; /* minimum ambient radius */
42    
43     typedef struct ambval {
44     FVECT pos; /* position in space */
45     FVECT dir; /* normal direction */
46     int lvl; /* recursion level of parent ray */
47     float weight; /* weight of parent ray */
48     COLOR val; /* computed ambient value */
49     float rad; /* validity radius */
50     struct ambval *next; /* next in list */
51     } AMBVAL; /* ambient value */
52    
53     typedef struct ambtree {
54     AMBVAL *alist; /* ambient value list */
55     struct ambtree *kid; /* 8 child nodes */
56     } AMBTREE; /* ambient octree */
57    
58     typedef struct {
59     float k; /* error contribution per sample */
60     COLOR v; /* ray sum */
61     int n; /* number of samples */
62     short t, p; /* theta, phi indices */
63     } AMBSAMP; /* ambient sample */
64    
65     static AMBTREE atrunk; /* our ambient trunk node */
66    
67     static FILE *ambfp = NULL; /* ambient file pointer */
68    
69     #define newambval() (AMBVAL *)bmalloc(sizeof(AMBVAL))
70    
71     #define newambtree() (AMBTREE *)calloc(8, sizeof(AMBTREE))
72    
73     double sumambient(), doambient(), makeambient();
74    
75    
76     setambient(afile) /* initialize calculation */
77     char *afile;
78     {
79     long ftell();
80     AMBVAL amb;
81 greg 1.8
82 greg 1.1 maxarad = thescene.cusize / 2.0; /* maximum radius */
83     /* minimum radius */
84     minarad = ambres > 0 ? thescene.cusize/ambres : 0.0;
85    
86     /* open ambient file */
87     if (afile != NULL)
88     if ((ambfp = fopen(afile, "r+")) != NULL) {
89 greg 1.9 while (fread((char *)&amb,sizeof(AMBVAL),1,ambfp) == 1)
90 greg 1.1 avinsert(&amb, &atrunk, thescene.cuorg,
91     thescene.cusize);
92     /* align */
93     fseek(ambfp, -(ftell(ambfp)%sizeof(AMBVAL)), 1);
94     } else if ((ambfp = fopen(afile, "w")) == NULL) {
95     sprintf(errmsg, "cannot open ambient file \"%s\"",
96     afile);
97     error(SYSTEM, errmsg);
98 greg 1.8 }
99     }
100    
101    
102     ambnotify(obj) /* record new modifier */
103     OBJECT obj;
104     {
105 greg 1.11 static int hitlimit = 0;
106 greg 1.8 register OBJREC *o = objptr(obj);
107     register char **amblp;
108    
109 greg 1.11 if (hitlimit || !ismodifier(o->otype))
110 greg 1.8 return;
111     for (amblp = amblist; *amblp != NULL; amblp++)
112     if (!strcmp(o->oname, *amblp)) {
113 greg 1.11 if (ambset[0] >= MAXASET) {
114     error(WARNING, "too many modifiers in ambient list");
115     hitlimit++;
116     return; /* should this be fatal? */
117     }
118 greg 1.8 insertelem(ambset, obj);
119     return;
120 greg 1.1 }
121     }
122    
123    
124     ambient(acol, r) /* compute ambient component for ray */
125     COLOR acol;
126     register RAY *r;
127     {
128     static int rdepth = 0; /* ambient recursion */
129     double wsum;
130    
131     rdepth++; /* increment level */
132    
133     if (ambdiv <= 0) /* no ambient calculation */
134     goto dumbamb;
135     /* check number of bounces */
136     if (rdepth > ambounce)
137     goto dumbamb;
138     /* check ambient list */
139     if (ambincl != -1 && r->ro != NULL &&
140     ambincl != inset(ambset, r->ro->omod))
141     goto dumbamb;
142    
143     if (ambacc <= FTINY) { /* no ambient storage */
144     if (doambient(acol, r) == 0.0)
145     goto dumbamb;
146     goto done;
147     }
148     /* get ambient value */
149     setcolor(acol, 0.0, 0.0, 0.0);
150     wsum = sumambient(acol, r, &atrunk, thescene.cuorg, thescene.cusize);
151     if (wsum > FTINY)
152     scalecolor(acol, 1.0/wsum);
153     else if (makeambient(acol, r) == 0.0)
154     goto dumbamb;
155     goto done;
156    
157     dumbamb: /* return global value */
158     copycolor(acol, ambval);
159     done: /* must finish here! */
160     rdepth--;
161     }
162    
163    
164     double
165     sumambient(acol, r, at, c0, s) /* get interpolated ambient value */
166     COLOR acol;
167     register RAY *r;
168     AMBTREE *at;
169     FVECT c0;
170     double s;
171     {
172     extern double sqrt();
173     double d, e1, e2, wt, wsum;
174     COLOR ct;
175     FVECT ck0;
176     int i;
177     register int j;
178     register AMBVAL *av;
179 greg 1.7 /* do this node */
180 greg 1.1 wsum = 0.0;
181     for (av = at->alist; av != NULL; av = av->next) {
182     /*
183     * Ray strength test.
184     */
185     if (av->lvl > r->rlvl || av->weight < r->rweight-FTINY)
186     continue;
187     /*
188     * Ambient radius test.
189     */
190     e1 = 0.0;
191     for (j = 0; j < 3; j++) {
192     d = av->pos[j] - r->rop[j];
193     e1 += d * d;
194     }
195     e1 /= av->rad * av->rad;
196     if (e1 > ambacc*ambacc*1.21)
197     continue;
198     /*
199     * Normal direction test.
200     */
201     e2 = (1.0 - DOT(av->dir, r->ron)) * r->rweight;
202     if (e2 < 0.0) e2 = 0.0;
203     if (e1 + e2 > ambacc*ambacc*1.21)
204     continue;
205     /*
206     * Ray behind test.
207     */
208     d = 0.0;
209     for (j = 0; j < 3; j++)
210     d += (r->rop[j] - av->pos[j]) *
211     (av->dir[j] + r->ron[j]);
212     if (d < -minarad)
213     continue;
214     /*
215     * Jittering final test reduces image artifacts.
216     */
217     wt = sqrt(e1) + sqrt(e2);
218 greg 1.7 wt *= .9 + .2*frandom();
219 greg 1.6 if (wt > ambacc)
220 greg 1.1 continue;
221     if (wt <= 1e-3)
222     wt = 1e3;
223     else
224     wt = 1.0 / wt;
225     wsum += wt;
226     copycolor(ct, av->val);
227     scalecolor(ct, wt);
228     addcolor(acol, ct);
229 greg 1.7 }
230     if (at->kid == NULL)
231     return(wsum);
232     /* do children */
233     s *= 0.5;
234     for (i = 0; i < 8; i++) {
235     for (j = 0; j < 3; j++) {
236     ck0[j] = c0[j];
237     if (1<<j & i)
238     ck0[j] += s;
239     if (r->rop[j] < ck0[j] - OCTSCALE*s)
240     break;
241     if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
242     break;
243     }
244     if (j == 3)
245     wsum += sumambient(acol, r, at->kid+i, ck0, s);
246 greg 1.1 }
247     return(wsum);
248     }
249    
250    
251     double
252     makeambient(acol, r) /* make a new ambient value */
253     COLOR acol;
254     register RAY *r;
255     {
256     AMBVAL amb;
257    
258     amb.rad = doambient(acol, r); /* compute ambient */
259     if (amb.rad == 0.0)
260     return(0.0);
261     /* store it */
262     VCOPY(amb.pos, r->rop);
263     VCOPY(amb.dir, r->ron);
264     amb.lvl = r->rlvl;
265     amb.weight = r->rweight;
266     copycolor(amb.val, acol);
267     /* insert into tree */
268     avinsert(&amb, &atrunk, thescene.cuorg, thescene.cusize);
269     avsave(&amb); /* write to file */
270     return(amb.rad);
271     }
272    
273    
274     double
275     doambient(acol, r) /* compute ambient component */
276     COLOR acol;
277     register RAY *r;
278     {
279     extern int ambcmp();
280     extern double sin(), cos(), sqrt();
281     double phi, xd, yd, zd;
282 greg 1.3 double b, b2;
283 greg 1.1 register AMBSAMP *div;
284     AMBSAMP dnew;
285     RAY ar;
286     FVECT ux, uy;
287     double arad;
288     int ndivs, nt, np, ns, ne, i, j;
289     register int k;
290    
291     setcolor(acol, 0.0, 0.0, 0.0);
292     /* set number of divisions */
293     nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5;
294     np = 2 * nt;
295     ndivs = nt * np;
296     /* check first */
297     if (ndivs == 0 || rayorigin(&ar, r, AMBIENT, 0.5) < 0)
298     return(0.0);
299     /* set number of super-samples */
300     ns = ambssamp * r->rweight + 0.5;
301     if (ns > 0) {
302     div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
303     if (div == NULL)
304     error(SYSTEM, "out of memory in doambient");
305     }
306     /* make axes */
307     uy[0] = uy[1] = uy[2] = 0.0;
308     for (k = 0; k < 3; k++)
309     if (r->ron[k] < 0.6 && r->ron[k] > -0.6)
310     break;
311     uy[k] = 1.0;
312     fcross(ux, r->ron, uy);
313     normalize(ux);
314     fcross(uy, ux, r->ron);
315     /* sample divisions */
316     arad = 0.0;
317     ne = 0;
318     for (i = 0; i < nt; i++)
319     for (j = 0; j < np; j++) {
320     rayorigin(&ar, r, AMBIENT, 0.5); /* pretested */
321     zd = sqrt((i+frandom())/nt);
322     phi = 2.0*PI * (j+frandom())/np;
323     xd = cos(phi) * zd;
324     yd = sin(phi) * zd;
325     zd = sqrt(1.0 - zd*zd);
326     for (k = 0; k < 3; k++)
327     ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
328     rayvalue(&ar);
329     if (ar.rot < FHUGE)
330     arad += 1.0 / ar.rot;
331     if (ns > 0) { /* save division */
332     div[ne].k = 0.0;
333     copycolor(div[ne].v, ar.rcol);
334     div[ne].n = 0;
335     div[ne].t = i; div[ne].p = j;
336     /* sum errors */
337 greg 1.3 b = bright(ar.rcol);
338 greg 1.1 if (i > 0) { /* from above */
339 greg 1.3 b2 = bright(div[ne-np].v) - b;
340     b2 *= b2 * 0.25;
341     div[ne].k += b2;
342 greg 1.1 div[ne].n++;
343 greg 1.3 div[ne-np].k += b2;
344 greg 1.1 div[ne-np].n++;
345     }
346     if (j > 0) { /* from behind */
347 greg 1.3 b2 = bright(div[ne-1].v) - b;
348     b2 *= b2 * 0.25;
349     div[ne].k += b2;
350 greg 1.1 div[ne].n++;
351 greg 1.3 div[ne-1].k += b2;
352 greg 1.1 div[ne-1].n++;
353     }
354     if (j == np-1) { /* around */
355 greg 1.3 b2 = bright(div[ne-(np-1)].v) - b;
356     b2 *= b2 * 0.25;
357     div[ne].k += b2;
358 greg 1.1 div[ne].n++;
359 greg 1.3 div[ne-(np-1)].k += b2;
360 greg 1.1 div[ne-(np-1)].n++;
361     }
362     ne++;
363     } else
364     addcolor(acol, ar.rcol);
365     }
366     for (k = 0; k < ne; k++) { /* compute errors */
367     if (div[k].n > 1)
368     div[k].k /= div[k].n;
369     div[k].n = 1;
370     }
371     /* sort the divisions */
372     qsort(div, ne, sizeof(AMBSAMP), ambcmp);
373     /* skim excess */
374     while (ne > ns) {
375     ne--;
376     addcolor(acol, div[ne].v);
377     }
378     /* super-sample */
379     for (i = ns; i > 0; i--) {
380     rayorigin(&ar, r, AMBIENT, 0.5); /* pretested */
381     zd = sqrt((div[0].t+frandom())/nt);
382     phi = 2.0*PI * (div[0].p+frandom())/np;
383     xd = cos(phi) * zd;
384     yd = sin(phi) * zd;
385     zd = sqrt(1.0 - zd*zd);
386     for (k = 0; k < 3; k++)
387     ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
388     rayvalue(&ar);
389     if (ar.rot < FHUGE)
390     arad += 1.0 / ar.rot;
391     /* recompute error */
392     copycolor(dnew.v, div[0].v);
393     addcolor(dnew.v, ar.rcol);
394     dnew.n = div[0].n + 1;
395     dnew.t = div[0].t; dnew.p = div[0].p;
396 greg 1.3 b2 = bright(dnew.v)/dnew.n - bright(ar.rcol);
397 greg 1.5 b2 = b2*b2 + div[0].k*(div[0].n*div[0].n);
398     dnew.k = b2/(dnew.n*dnew.n);
399 greg 1.1 /* reinsert */
400     for (k = 0; k < ne-1 && dnew.k < div[k+1].k; k++)
401 greg 1.9 copystruct(&div[k], &div[k+1]);
402     copystruct(&div[k], &dnew);
403 greg 1.1
404     if (ne >= i) { /* extract darkest division */
405     ne--;
406     if (div[ne].n > 1)
407     scalecolor(div[ne].v, 1.0/div[ne].n);
408     addcolor(acol, div[ne].v);
409     }
410     }
411     scalecolor(acol, 1.0/ndivs);
412     if (arad <= FTINY)
413     arad = FHUGE;
414     else
415     arad = (ndivs+ns) / arad / sqrt(r->rweight);
416     if (arad > maxarad)
417     arad = maxarad;
418     else if (arad < minarad)
419     arad = minarad;
420     if (ns > 0)
421     free((char *)div);
422     return(arad);
423     }
424    
425    
426     static int
427     ambcmp(d1, d2) /* decreasing order */
428     AMBSAMP *d1, *d2;
429     {
430     if (d1->k < d2->k)
431     return(1);
432     if (d1->k > d2->k)
433     return(-1);
434     return(0);
435     }
436    
437    
438     static
439     avsave(av) /* save an ambient value */
440     AMBVAL *av;
441     {
442     #ifdef AMBFLUSH
443     static int nunflshed = 0;
444     #endif
445     if (ambfp == NULL)
446     return;
447 greg 1.9 if (fwrite((char *)av, sizeof(AMBVAL), 1, ambfp) != 1)
448 greg 1.1 goto writerr;
449     #ifdef AMBFLUSH
450     if (++nunflshed >= AMBFLUSH) {
451     if (fflush(ambfp) == EOF)
452     goto writerr;
453     nunflshed = 0;
454     }
455     #endif
456     return;
457     writerr:
458     error(SYSTEM, "error writing ambient file");
459     }
460    
461    
462     static
463     avinsert(aval, at, c0, s) /* insert ambient value in a tree */
464     AMBVAL *aval;
465     register AMBTREE *at;
466     FVECT c0;
467     double s;
468     {
469     FVECT ck0;
470     int branch;
471     register AMBVAL *av;
472     register int i;
473    
474     if ((av = newambval()) == NULL)
475     goto memerr;
476 greg 1.9 copystruct(av, aval);
477 greg 1.1 VCOPY(ck0, c0);
478     while (s*(OCTSCALE/2) > av->rad*ambacc) {
479     if (at->kid == NULL)
480     if ((at->kid = newambtree()) == NULL)
481     goto memerr;
482     s *= 0.5;
483     branch = 0;
484     for (i = 0; i < 3; i++)
485     if (av->pos[i] > ck0[i] + s) {
486     ck0[i] += s;
487     branch |= 1 << i;
488     }
489     at = at->kid + branch;
490     }
491     av->next = at->alist;
492     at->alist = av;
493     return;
494     memerr:
495     error(SYSTEM, "out of memory in avinsert");
496     }