ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/source.c
Revision: 2.41
Committed: Fri Jan 9 05:37:11 2004 UTC (20 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.40: +2 -2 lines
Log Message:
Bug fix in source substructuring test and added debug bounds checking

File Contents

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