ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/source.c
Revision: 2.52
Committed: Tue May 31 18:01:09 2005 UTC (18 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1
Changes since 2.51: +3 -2 lines
Log Message:
Added Russian roulette ray termination and fixed ambient weights & measures

File Contents

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