ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/source.c
Revision: 2.25
Committed: Tue May 28 13:41:41 1996 UTC (27 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.24: +17 -21 lines
Log Message:
improved scattering calculation accuracy

File Contents

# Content
1 /* Copyright (c) 1995 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * source.c - routines dealing with illumination sources.
9 *
10 * 8/20/85
11 */
12
13 #include "ray.h"
14
15 #include "octree.h"
16
17 #include "otypes.h"
18
19 #include "source.h"
20
21 #include "random.h"
22
23 extern double ssampdist; /* scatter sampling distance */
24
25 #ifndef MAXSSAMP
26 #define MAXSSAMP 16 /* maximum samples per ray */
27 #endif
28
29 /*
30 * Structures used by direct()
31 */
32
33 typedef struct {
34 int sno; /* source number */
35 FVECT dir; /* source direction */
36 COLOR coef; /* material coefficient */
37 COLOR val; /* contribution */
38 } CONTRIB; /* direct contribution */
39
40 typedef struct {
41 int sndx; /* source index (to CONTRIB array) */
42 float brt; /* brightness (for comparison) */
43 } CNTPTR; /* contribution pointer */
44
45 static CONTRIB *srccnt; /* source contributions in direct() */
46 static CNTPTR *cntord; /* source ordering in direct() */
47 static int maxcntr = 0; /* size of contribution arrays */
48
49
50 marksources() /* find and mark source objects */
51 {
52 int foundsource = 0;
53 int i;
54 register OBJREC *o, *m;
55 register int ns;
56 /* initialize dispatch table */
57 initstypes();
58 /* find direct sources */
59 for (i = 0; i < nobjects; i++) {
60
61 o = objptr(i);
62
63 if (!issurface(o->otype) || o->omod == OVOID)
64 continue;
65
66 m = objptr(o->omod);
67
68 if (!islight(m->otype))
69 continue;
70
71 if (m->oargs.nfargs != (m->otype == MAT_GLOW ? 4 :
72 m->otype == MAT_SPOT ? 7 : 3))
73 objerror(m, USER, "bad # arguments");
74
75 if (m->otype == MAT_GLOW &&
76 o->otype != OBJ_SOURCE &&
77 m->oargs.farg[3] <= FTINY)
78 continue; /* don't bother */
79 if (m->oargs.farg[0] <= FTINY && m->oargs.farg[1] <= FTINY &&
80 m->oargs.farg[2] <= FTINY)
81 continue; /* don't bother */
82
83 if (sfun[o->otype].of == NULL ||
84 sfun[o->otype].of->setsrc == NULL)
85 objerror(o, USER, "illegal material");
86
87 if ((ns = newsource()) < 0)
88 goto memerr;
89
90 setsource(&source[ns], o);
91
92 if (m->otype == MAT_GLOW) {
93 source[ns].sflags |= SPROX;
94 source[ns].sl.prox = m->oargs.farg[3];
95 if (source[ns].sflags & SDISTANT)
96 source[ns].sflags |= SSKIP;
97 } else if (m->otype == MAT_SPOT) {
98 source[ns].sflags |= SSPOT;
99 if ((source[ns].sl.s = makespot(m)) == NULL)
100 goto memerr;
101 if (source[ns].sflags & SFLAT &&
102 !checkspot(source[ns].sl.s,source[ns].snorm)) {
103 objerror(o, WARNING,
104 "invalid spotlight direction");
105 source[ns].sflags |= SSKIP;
106 }
107 }
108 if (!(source[ns].sflags & SSKIP))
109 foundsource++;
110 }
111 if (!foundsource) {
112 error(WARNING, "no light sources found");
113 return;
114 }
115 markvirtuals(); /* find and add virtual sources */
116 /* allocate our contribution arrays */
117 maxcntr = nsources + MAXSPART; /* start with this many */
118 srccnt = (CONTRIB *)malloc(maxcntr*sizeof(CONTRIB));
119 cntord = (CNTPTR *)malloc(maxcntr*sizeof(CNTPTR));
120 if (srccnt == NULL | cntord == NULL)
121 goto memerr;
122 return;
123 memerr:
124 error(SYSTEM, "out of memory in marksources");
125 }
126
127
128 srcray(sr, r, si) /* send a ray to a source, return domega */
129 register RAY *sr; /* returned source ray */
130 RAY *r; /* ray which hit object */
131 SRCINDEX *si; /* source sample index */
132 {
133 double d; /* distance to source */
134 register SRCREC *srcp;
135
136 rayorigin(sr, r, SHADOW, 1.0); /* ignore limits */
137
138 while ((d = nextssamp(sr, si)) != 0.0) {
139 sr->rsrc = si->sn; /* remember source */
140 srcp = source + si->sn;
141 if (srcp->sflags & SDISTANT) {
142 if (srcp->sflags & SSPOT && spotout(sr, srcp->sl.s))
143 continue;
144 return(1); /* sample OK */
145 }
146 /* local source */
147 /* check proximity */
148 if (srcp->sflags & SPROX && d > srcp->sl.prox)
149 continue;
150 /* check angle */
151 if (srcp->sflags & SSPOT) {
152 if (spotout(sr, srcp->sl.s))
153 continue;
154 /* adjust solid angle */
155 si->dom *= d*d;
156 d += srcp->sl.s->flen;
157 si->dom /= d*d;
158 }
159 return(1); /* sample OK */
160 }
161 return(0); /* no more samples */
162 }
163
164
165 srcvalue(r) /* punch ray to source and compute value */
166 register RAY *r;
167 {
168 register SRCREC *sp;
169
170 sp = &source[r->rsrc];
171 if (sp->sflags & SVIRTUAL) { /* virtual source */
172 /* check intersection */
173 if (!(*ofun[sp->so->otype].funp)(sp->so, r))
174 return;
175 if (!rayshade(r, r->ro->omod)) /* compute contribution */
176 goto nomat;
177 rayparticipate(r);
178 return;
179 }
180 /* compute intersection */
181 if (sp->sflags & SDISTANT ? sourcehit(r) :
182 (*ofun[sp->so->otype].funp)(sp->so, r)) {
183 if (sp->sa.success >= 0)
184 sp->sa.success++;
185 if (!rayshade(r, r->ro->omod)) /* compute contribution */
186 goto nomat;
187 rayparticipate(r);
188 return;
189 }
190 /* we missed our mark! */
191 if (sp->sa.success < 0)
192 return; /* bitched already */
193 sp->sa.success -= AIMREQT;
194 if (sp->sa.success >= 0)
195 return; /* leniency */
196 sprintf(errmsg, "aiming failure for light source \"%s\"",
197 sp->so->oname);
198 error(WARNING, errmsg); /* issue warning */
199 return;
200 nomat:
201 objerror(r->ro, USER, "material not found");
202 }
203
204
205 sourcehit(r) /* check to see if ray hit distant source */
206 register RAY *r;
207 {
208 int first, last;
209 register int i;
210
211 if (r->rsrc >= 0) { /* check only one if aimed */
212 first = last = r->rsrc;
213 } else { /* otherwise check all */
214 first = 0; last = nsources-1;
215 }
216 for (i = first; i <= last; i++)
217 if ((source[i].sflags & (SDISTANT|SVIRTUAL)) == SDISTANT)
218 /*
219 * Check to see if ray is within
220 * solid angle of source.
221 */
222 if (2.0*PI * (1.0 - DOT(source[i].sloc,r->rdir))
223 <= source[i].ss2) {
224 r->ro = source[i].so;
225 if (!(source[i].sflags & SSKIP))
226 break;
227 }
228
229 if (r->ro != NULL) {
230 for (i = 0; i < 3; i++)
231 r->ron[i] = -r->rdir[i];
232 r->rod = 1.0;
233 r->rox = NULL;
234 return(1);
235 }
236 return(0);
237 }
238
239
240 static int
241 cntcmp(sc1, sc2) /* contribution compare (descending) */
242 register CNTPTR *sc1, *sc2;
243 {
244 if (sc1->brt > sc2->brt)
245 return(-1);
246 if (sc1->brt < sc2->brt)
247 return(1);
248 return(0);
249 }
250
251
252 direct(r, f, p) /* add direct component */
253 RAY *r; /* ray that hit surface */
254 int (*f)(); /* direct component coefficient function */
255 char *p; /* data for f */
256 {
257 extern int (*trace)();
258 register int sn;
259 register CONTRIB *scp;
260 SRCINDEX si;
261 int nshadcheck, ncnts;
262 int nhits;
263 double prob, ourthresh, hwt;
264 RAY sr;
265 /* NOTE: srccnt and cntord global so no recursion */
266 if (nsources <= 0)
267 return; /* no sources?! */
268 /* potential contributions */
269 initsrcindex(&si);
270 for (sn = 0; srcray(&sr, r, &si); sn++) {
271 if (sn >= maxcntr) {
272 maxcntr = sn + MAXSPART;
273 srccnt = (CONTRIB *)realloc((char *)srccnt,
274 maxcntr*sizeof(CONTRIB));
275 cntord = (CNTPTR *)realloc((char *)cntord,
276 maxcntr*sizeof(CNTPTR));
277 if (srccnt == NULL | cntord == NULL)
278 error(SYSTEM, "out of memory in direct");
279 }
280 cntord[sn].sndx = sn;
281 scp = srccnt + sn;
282 scp->sno = sr.rsrc;
283 /* compute coefficient */
284 (*f)(scp->coef, p, sr.rdir, si.dom);
285 cntord[sn].brt = bright(scp->coef);
286 if (cntord[sn].brt <= 0.0)
287 continue;
288 VCOPY(scp->dir, sr.rdir);
289 /* compute potential */
290 sr.revf = srcvalue;
291 rayvalue(&sr);
292 copycolor(scp->val, sr.rcol);
293 multcolor(scp->val, scp->coef);
294 cntord[sn].brt = bright(scp->val);
295 }
296 /* sort contributions */
297 qsort(cntord, sn, sizeof(CNTPTR), cntcmp);
298 { /* find last */
299 register int l, m;
300
301 ncnts = l = sn;
302 sn = 0;
303 while ((m = (sn + ncnts) >> 1) != l) {
304 if (cntord[m].brt > 0.0)
305 sn = m;
306 else
307 ncnts = m;
308 l = m;
309 }
310 }
311 if (ncnts == 0)
312 return; /* no contributions! */
313 /* accumulate tail */
314 for (sn = ncnts-1; sn > 0; sn--)
315 cntord[sn-1].brt += cntord[sn].brt;
316 /* compute number to check */
317 nshadcheck = pow((double)ncnts, shadcert) + .5;
318 /* modify threshold */
319 ourthresh = shadthresh / r->rweight;
320 /* test for shadows */
321 for (nhits = 0, hwt = 0.0, sn = 0; sn < ncnts;
322 hwt += (double)source[scp->sno].nhits /
323 (double)source[scp->sno].ntests,
324 sn++) {
325 /* check threshold */
326 if ((sn+nshadcheck>=ncnts ? cntord[sn].brt :
327 cntord[sn].brt-cntord[sn+nshadcheck].brt)
328 < ourthresh*bright(r->rcol))
329 break;
330 scp = srccnt + cntord[sn].sndx;
331 /* test for hit */
332 rayorigin(&sr, r, SHADOW, 1.0);
333 VCOPY(sr.rdir, scp->dir);
334 sr.rsrc = scp->sno;
335 source[scp->sno].ntests++; /* keep statistics */
336 if (localhit(&sr, &thescene) &&
337 ( sr.ro != source[scp->sno].so ||
338 source[scp->sno].sflags & SFOLLOW )) {
339 /* follow entire path */
340 raycont(&sr);
341 rayparticipate(&sr);
342 if (trace != NULL)
343 (*trace)(&sr); /* trace execution */
344 if (bright(sr.rcol) <= FTINY)
345 continue; /* missed! */
346 copycolor(scp->val, sr.rcol);
347 multcolor(scp->val, scp->coef);
348 }
349 /* add contribution if hit */
350 addcolor(r->rcol, scp->val);
351 nhits++;
352 source[scp->sno].nhits++;
353 }
354 /* source hit rate */
355 if (hwt > FTINY)
356 hwt = (double)nhits / hwt;
357 else
358 hwt = 0.5;
359 #ifdef DEBUG
360 sprintf(errmsg, "%d tested, %d untested, %f conditional hit rate\n",
361 sn, ncnts-sn, hwt);
362 eputs(errmsg);
363 #endif
364 /* add in untested sources */
365 for ( ; sn < ncnts; sn++) {
366 scp = srccnt + cntord[sn].sndx;
367 prob = hwt * (double)source[scp->sno].nhits /
368 (double)source[scp->sno].ntests;
369 if (prob > 1.0)
370 prob = 1.0;
371 scalecolor(scp->val, prob);
372 addcolor(r->rcol, scp->val);
373 }
374 }
375
376
377 srcscatter(r) /* compute source scattering into ray */
378 register RAY *r;
379 {
380 int oldsampndx;
381 int nsamps;
382 RAY sr;
383 SRCINDEX si;
384 double t, d;
385 double re, ge, be;
386 COLOR cvext;
387 int i, j;
388
389 if (r->slights == NULL || r->slights[0] == 0
390 || r->gecc >= 1.-FTINY || r->rot >= FHUGE)
391 return;
392 if (ssampdist <= FTINY || (nsamps = r->rot/ssampdist + .5) < 1)
393 nsamps = 1;
394 #if MAXSSAMP
395 else if (nsamps > MAXSSAMP)
396 nsamps = MAXSSAMP;
397 #endif
398 oldsampndx = samplendx;
399 samplendx = random()&0x7fff; /* randomize */
400 for (i = r->slights[0]; i > 0; i--) { /* for each source */
401 for (j = 0; j < nsamps; j++) { /* for each sample position */
402 samplendx++;
403 t = r->rot * (j+frandom())/nsamps;
404 /* extinction */
405 re = t*colval(r->cext,RED);
406 ge = t*colval(r->cext,GRN);
407 be = t*colval(r->cext,BLU);
408 setcolor(cvext, re > 92. ? 0. : exp(-re),
409 ge > 92. ? 0. : exp(-ge),
410 be > 92. ? 0. : exp(-be));
411 if (intens(cvext) <= FTINY)
412 break; /* too far away */
413 sr.rorg[0] = r->rorg[0] + r->rdir[0]*t;
414 sr.rorg[1] = r->rorg[1] + r->rdir[1]*t;
415 sr.rorg[2] = r->rorg[2] + r->rdir[2]*t;
416 sr.rmax = 0.;
417 initsrcindex(&si); /* sample ray to this source */
418 si.sn = r->slights[i];
419 nopart(&si, &sr);
420 if (!srcray(&sr, NULL, &si) ||
421 sr.rsrc != r->slights[i])
422 continue; /* no path */
423 copycolor(sr.cext, r->cext);
424 copycolor(sr.albedo, r->albedo);
425 sr.gecc = r->gecc;
426 rayvalue(&sr); /* eval. source ray */
427 if (bright(sr.rcol) <= FTINY)
428 continue;
429 if (r->gecc <= FTINY) /* compute P(theta) */
430 d = 1.;
431 else {
432 d = DOT(r->rdir, sr.rdir);
433 d = 1. + r->gecc*r->gecc - 2.*r->gecc*d;
434 d = (1. - r->gecc*r->gecc) / (d*sqrt(d));
435 }
436 /* other factors */
437 d *= si.dom * r->rot / (4.*PI*nsamps);
438 multcolor(sr.rcol, r->cext);
439 multcolor(sr.rcol, r->albedo);
440 scalecolor(sr.rcol, d);
441 multcolor(sr.rcol, cvext);
442 addcolor(r->rcol, sr.rcol); /* add it in */
443 }
444 }
445 samplendx = oldsampndx;
446 }
447
448
449 /****************************************************************
450 * The following macros were separated from the m_light() routine
451 * because they are very nasty and difficult to understand.
452 */
453
454 /* illumblock *
455 *
456 * We cannot allow an illum to pass to another illum, because that
457 * would almost certainly constitute overcounting.
458 * However, we do allow an illum to pass to another illum
459 * that is actually going to relay to a virtual light source.
460 * We also prevent an illum from passing to a glow; this provides a
461 * convenient mechanism for defining detailed light source
462 * geometry behind (or inside) an effective radiator.
463 */
464
465 static int weaksrcmod(obj) int obj; /* efficiency booster function */
466 {register OBJREC *o = objptr(obj);
467 return(o->otype==MAT_ILLUM|o->otype==MAT_GLOW);}
468
469 #define illumblock(m, r) (!(source[r->rsrc].sflags&SVIRTUAL) && \
470 r->rod > 0.0 && \
471 weaksrcmod(source[r->rsrc].so->omod))
472
473 /* wrongsource *
474 *
475 * This source is the wrong source (ie. overcounted) if we are
476 * aimed to a different source than the one we hit and the one
477 * we hit is not an illum that should be passed.
478 */
479
480 #define wrongsource(m, r) (r->rsrc>=0 && source[r->rsrc].so!=r->ro && \
481 (m->otype!=MAT_ILLUM || illumblock(m,r)))
482
483 /* distglow *
484 *
485 * A distant glow is an object that sometimes acts as a light source,
486 * but is too far away from the test point to be one in this case.
487 * (Glows with negative radii should NEVER participate in illumination.)
488 */
489
490 #define distglow(m, r, d) (m->otype==MAT_GLOW && \
491 m->oargs.farg[3] >= -FTINY && \
492 d > m->oargs.farg[3])
493
494 /* badcomponent *
495 *
496 * We must avoid counting light sources in the ambient calculation,
497 * since the direct component is handled separately. Therefore, any
498 * ambient ray which hits an active light source must be discarded.
499 * The same is true for stray specular samples, since the specular
500 * contribution from light sources is calculated separately.
501 */
502
503 #define badcomponent(m, r) (r->crtype&(AMBIENT|SPECULAR) && \
504 !(r->crtype&SHADOW || r->rod < 0.0 || \
505 /* not 100% correct */ distglow(m, r, r->rot)))
506
507 /* passillum *
508 *
509 * An illum passes to another material type when we didn't hit it
510 * on purpose (as part of a direct calculation), or it is relaying
511 * a virtual light source.
512 */
513
514 #define passillum(m, r) (m->otype==MAT_ILLUM && \
515 (r->rsrc<0 || source[r->rsrc].so!=r->ro || \
516 source[r->rsrc].sflags&SVIRTUAL))
517
518 /* srcignore *
519 *
520 * The -dv flag is normally on for sources to be visible.
521 */
522
523 #define srcignore(m, r) !(directvis || r->crtype&SHADOW || \
524 distglow(m, r, raydist(r,PRIMARY)))
525
526
527 m_light(m, r) /* ray hit a light source */
528 register OBJREC *m;
529 register RAY *r;
530 {
531 /* check for over-counting */
532 if (badcomponent(m, r))
533 return(1);
534 if (wrongsource(m, r))
535 return(1);
536 /* check for passed illum */
537 if (passillum(m, r)) {
538 if (m->oargs.nsargs && strcmp(m->oargs.sarg[0], VOIDID))
539 return(rayshade(r, modifier(m->oargs.sarg[0])));
540 raytrans(r);
541 return(1);
542 }
543 /* otherwise treat as source */
544 /* check for behind */
545 if (r->rod < 0.0)
546 return(1);
547 /* check for invisibility */
548 if (srcignore(m, r))
549 return(1);
550 /* check for outside spot */
551 if (m->otype==MAT_SPOT && spotout(r, makespot(m)))
552 return(1);
553 /* get distribution pattern */
554 raytexture(r, m->omod);
555 /* get source color */
556 setcolor(r->rcol, m->oargs.farg[0],
557 m->oargs.farg[1],
558 m->oargs.farg[2]);
559 /* modify value */
560 multcolor(r->rcol, r->pcol);
561 return(1);
562 }