ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/source.c
Revision: 1.31
Committed: Tue May 21 17:41:26 1991 UTC (32 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.30: +7 -3 lines
Log Message:
added stratified random sampling with urand

File Contents

# Content
1 /* Copyright (c) 1990 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 "source.h"
18
19 #include "otypes.h"
20
21 #include "cone.h"
22
23 #include "face.h"
24
25 #include "random.h"
26
27
28 extern double dstrsrc; /* source distribution amount */
29 extern double shadthresh; /* relative shadow threshold */
30 extern double shadcert; /* shadow testing certainty */
31
32 SRCREC *source = NULL; /* our list of sources */
33 int nsources = 0; /* the number of sources */
34
35 static CONTRIB *srccnt; /* source contributions in direct() */
36 static CNTPTR *cntord; /* source ordering in direct() */
37
38
39 marksources() /* find and mark source objects */
40 {
41 register OBJREC *o, *m;
42 register int i;
43
44 for (i = 0; i < nobjects; i++) {
45
46 o = objptr(i);
47
48 if (o->omod == OVOID)
49 continue;
50
51 m = objptr(o->omod);
52
53 if (!islight(m->otype))
54 continue;
55
56 if (m->oargs.nfargs != (m->otype == MAT_GLOW ? 4 :
57 m->otype == MAT_SPOT ? 7 : 3))
58 objerror(m, USER, "bad # arguments");
59
60 if (m->otype == MAT_GLOW &&
61 o->otype != OBJ_SOURCE &&
62 m->oargs.farg[3] <= FTINY)
63 continue; /* don't bother */
64
65 if (source == NULL)
66 source = (SRCREC *)malloc(sizeof(SRCREC));
67 else
68 source = (SRCREC *)realloc((char *)source,
69 (unsigned)(nsources+1)*sizeof(SRCREC));
70 if (source == NULL)
71 goto memerr;
72
73 newsource(&source[nsources], o);
74
75 if (m->otype == MAT_GLOW) {
76 source[nsources].sflags |= SPROX;
77 source[nsources].sl.prox = m->oargs.farg[3];
78 if (o->otype == OBJ_SOURCE)
79 source[nsources].sflags |= SSKIP;
80 } else if (m->otype == MAT_SPOT) {
81 source[nsources].sflags |= SSPOT;
82 source[nsources].sl.s = makespot(m);
83 }
84 nsources++;
85 }
86 if (nsources <= 0) {
87 error(WARNING, "no light sources found");
88 return;
89 }
90 srccnt = (CONTRIB *)malloc(nsources*sizeof(CONTRIB));
91 cntord = (CNTPTR *)malloc(nsources*sizeof(CNTPTR));
92 if (srccnt != NULL && cntord != NULL)
93 return;
94 /* fall through */
95 memerr:
96 error(SYSTEM, "out of memory in marksources");
97 }
98
99
100 newsource(src, so) /* add a source to the array */
101 register SRCREC *src;
102 register OBJREC *so;
103 {
104 double cos(), tan(), sqrt();
105 double theta;
106 FACE *f;
107 CONE *co;
108 int j;
109 register int i;
110
111 src->sflags = 0;
112 src->aimsuccess = 2*AIMREQT-1; /* bitch on second failure */
113 src->nhits = 1; src->ntests = 2; /* start probability = 1/2 */
114 src->so = so;
115
116 switch (so->otype) {
117 case OBJ_SOURCE:
118 if (so->oargs.nfargs != 4)
119 objerror(so, USER, "bad arguments");
120 src->sflags |= SDISTANT;
121 VCOPY(src->sloc, so->oargs.farg);
122 if (normalize(src->sloc) == 0.0)
123 objerror(so, USER, "zero direction");
124 theta = PI/180.0/2.0 * so->oargs.farg[3];
125 if (theta <= FTINY)
126 objerror(so, USER, "zero size");
127 src->ss = theta >= PI/4 ? 1.0 : tan(theta);
128 src->ss2 = 2.0*PI * (1.0 - cos(theta));
129 break;
130 case OBJ_SPHERE:
131 VCOPY(src->sloc, so->oargs.farg);
132 src->ss = so->oargs.farg[3];
133 src->ss2 = PI * src->ss * src->ss;
134 break;
135 case OBJ_FACE:
136 /* get the face */
137 f = getface(so);
138 /* find the center */
139 for (j = 0; j < 3; j++) {
140 src->sloc[j] = 0.0;
141 for (i = 0; i < f->nv; i++)
142 src->sloc[j] += VERTEX(f,i)[j];
143 src->sloc[j] /= f->nv;
144 }
145 if (!inface(src->sloc, f))
146 objerror(so, USER, "cannot hit center");
147 src->ss = sqrt(f->area / PI);
148 src->ss2 = f->area;
149 break;
150 case OBJ_RING:
151 /* get the ring */
152 co = getcone(so, 0);
153 VCOPY(src->sloc, CO_P0(co));
154 if (CO_R0(co) > 0.0)
155 objerror(so, USER, "cannot hit center");
156 src->ss = CO_R1(co);
157 src->ss2 = PI * src->ss * src->ss;
158 break;
159 default:
160 objerror(so, USER, "illegal material");
161 }
162 }
163
164
165 SPOT *
166 makespot(m) /* make a spotlight */
167 register OBJREC *m;
168 {
169 extern double cos();
170 register SPOT *ns;
171
172 if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
173 error(SYSTEM, "out of memory in makespot");
174 ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
175 VCOPY(ns->aim, m->oargs.farg+4);
176 if ((ns->flen = normalize(ns->aim)) == 0.0)
177 objerror(m, USER, "zero focus vector");
178 return(ns);
179 }
180
181
182 double
183 srcray(sr, r, sn) /* send a ray to a source, return domega */
184 register RAY *sr; /* returned source ray */
185 RAY *r; /* ray which hit object */
186 register int sn; /* source number */
187 {
188 register double *norm = NULL; /* plane normal */
189 double ddot; /* (distance times) cosine */
190 FVECT vd;
191 double d;
192 register int i;
193
194 if (source[sn].sflags & SSKIP)
195 return(0.0); /* skip this source */
196
197 rayorigin(sr, r, SHADOW, 1.0); /* ignore limits */
198
199 sr->rsrc = sn; /* remember source */
200 /* get source direction */
201 if (source[sn].sflags & SDISTANT)
202 /* constant direction */
203 VCOPY(sr->rdir, source[sn].sloc);
204 else { /* compute direction */
205 for (i = 0; i < 3; i++)
206 sr->rdir[i] = source[sn].sloc[i] - sr->rorg[i];
207
208 if (source[sn].so->otype == OBJ_FACE)
209 norm = getface(source[sn].so)->norm;
210 else if (source[sn].so->otype == OBJ_RING)
211 norm = getcone(source[sn].so,0)->ad;
212
213 if (norm != NULL && (ddot = -DOT(sr->rdir, norm)) <= FTINY)
214 return(0.0); /* behind surface! */
215 }
216 if (dstrsrc > FTINY) {
217 /* distribute source direction */
218 dimlist[ndims++] = sn;
219 for (i = 0; i < 3; i++) {
220 dimlist[ndims] = i + 8831;
221 vd[i] = dstrsrc * source[sn].ss *
222 (1.0 - 2.0*urand(ilhash(dimlist,ndims+1)+samplendx));
223 }
224 ndims--;
225 if (norm != NULL) { /* project offset */
226 d = DOT(vd, norm);
227 for (i = 0; i < 3; i++)
228 vd[i] -= d * norm[i];
229 }
230 for (i = 0; i < 3; i++) /* offset source direction */
231 sr->rdir[i] += vd[i];
232
233 } else if (source[sn].sflags & SDISTANT)
234 /* already normalized */
235 return(source[sn].ss2);
236
237 if ((d = normalize(sr->rdir)) == 0.0)
238 /* at source! */
239 return(0.0);
240
241 if (source[sn].sflags & SDISTANT)
242 /* domega constant */
243 return(source[sn].ss2);
244
245 /* check proximity */
246 if (source[sn].sflags & SPROX &&
247 d > source[sn].sl.prox)
248 return(0.0);
249 /* compute dot product */
250 if (norm != NULL)
251 ddot /= d;
252 else
253 ddot = 1.0;
254 /* check angle */
255 if (source[sn].sflags & SSPOT) {
256 if (source[sn].sl.s->siz < 2.0*PI *
257 (1.0 + DOT(source[sn].sl.s->aim,sr->rdir)))
258 return(0.0);
259 d += source[sn].sl.s->flen; /* adjust length */
260 }
261 /* compute domega */
262 return(ddot*source[sn].ss2/(d*d));
263 }
264
265
266 sourcehit(r) /* check to see if ray hit distant source */
267 register RAY *r;
268 {
269 int first, last;
270 register int i;
271
272 if (r->rsrc >= 0) { /* check only one if aimed */
273 first = last = r->rsrc;
274 } else { /* otherwise check all */
275 first = 0; last = nsources-1;
276 }
277 for (i = first; i <= last; i++)
278 if (source[i].sflags & SDISTANT)
279 /*
280 * Check to see if ray is within
281 * solid angle of source.
282 */
283 if (2.0*PI * (1.0 - DOT(source[i].sloc,r->rdir))
284 <= source[i].ss2) {
285 r->ro = source[i].so;
286 if (!(source[i].sflags & SSKIP))
287 break;
288 }
289
290 if (r->ro != NULL) {
291 for (i = 0; i < 3; i++)
292 r->ron[i] = -r->rdir[i];
293 r->rod = 1.0;
294 r->rox = NULL;
295 return(1);
296 }
297 return(0);
298 }
299
300
301 static int
302 cntcmp(sc1, sc2) /* contribution compare (descending) */
303 register CNTPTR *sc1, *sc2;
304 {
305 if (sc1->brt > sc2->brt)
306 return(-1);
307 if (sc1->brt < sc2->brt)
308 return(1);
309 return(0);
310 }
311
312
313 direct(r, f, p) /* add direct component */
314 RAY *r; /* ray that hit surface */
315 int (*f)(); /* direct component coefficient function */
316 char *p; /* data for f */
317 {
318 extern double pow();
319 register int sn;
320 int nshadcheck, ncnts;
321 int nhits;
322 double prob, ourthresh, hwt;
323 RAY sr;
324 /* NOTE: srccnt and cntord global so no recursion */
325 if (nsources <= 0)
326 return; /* no sources?! */
327 /* compute number to check */
328 nshadcheck = pow((double)nsources, shadcert) + .5;
329 /* modify threshold */
330 ourthresh = shadthresh / r->rweight;
331 /* potential contributions */
332 for (sn = 0; sn < nsources; sn++) {
333 cntord[sn].sno = sn;
334 cntord[sn].brt = 0.0;
335 /* get source ray */
336 if ((srccnt[sn].dom = srcray(&sr, r, sn)) == 0.0)
337 continue;
338 VCOPY(srccnt[sn].dir, sr.rdir);
339 /* compute coefficient */
340 (*f)(srccnt[sn].val, p, srccnt[sn].dir, srccnt[sn].dom);
341 cntord[sn].brt = bright(srccnt[sn].val);
342 if (cntord[sn].brt <= 0.0)
343 continue;
344 /* compute intersection */
345 if (source[sn].sflags & SDISTANT ?
346 sourcehit(&sr) :
347 (*ofun[source[sn].so->otype].funp)
348 (source[sn].so, &sr)) {
349 if (source[sn].aimsuccess >= 0)
350 source[sn].aimsuccess++;
351 } else {
352 cntord[sn].brt = 0.0;
353 if (source[sn].aimsuccess < 0)
354 continue; /* bitched already */
355 source[sn].aimsuccess -= AIMREQT;
356 if (source[sn].aimsuccess >= 0)
357 continue; /* leniency */
358 sprintf(errmsg,
359 "aiming failure for light source \"%s\"",
360 source[sn].so->oname);
361 error(WARNING, errmsg);
362 continue;
363 }
364 /* compute contribution */
365 raycont(&sr);
366 multcolor(srccnt[sn].val, sr.rcol);
367 cntord[sn].brt = bright(srccnt[sn].val);
368 }
369 /* sort contributions */
370 qsort(cntord, nsources, sizeof(CNTPTR), cntcmp);
371 { /* find last */
372 register int l, m;
373
374 sn = 0; ncnts = l = nsources;
375 while ((m = (sn + ncnts) >> 1) != l) {
376 if (cntord[m].brt > 0.0)
377 sn = m;
378 else
379 ncnts = m;
380 l = m;
381 }
382 }
383 /* accumulate tail */
384 for (sn = ncnts-1; sn > 0; sn--)
385 cntord[sn-1].brt += cntord[sn].brt;
386 /* test for shadows */
387 nhits = 0;
388 for (sn = 0; sn < ncnts; sn++) {
389 /* check threshold */
390 if ((sn+nshadcheck>=ncnts ? cntord[sn].brt :
391 cntord[sn].brt-cntord[sn+nshadcheck].brt)
392 < ourthresh*bright(r->rcol))
393 break;
394 /* get statistics */
395 source[cntord[sn].sno].ntests++;
396 /* test for hit */
397 rayorigin(&sr, r, SHADOW, 1.0);
398 VCOPY(sr.rdir, srccnt[cntord[sn].sno].dir);
399 sr.rsrc = cntord[sn].sno;
400 if (localhit(&sr, &thescene) &&
401 sr.ro != source[cntord[sn].sno].so) {
402 /* check for transmission */
403 raycont(&sr);
404 if (bright(sr.rcol) <= FTINY)
405 continue; /* missed! */
406 (*f)(srccnt[cntord[sn].sno].val, p,
407 srccnt[cntord[sn].sno].dir,
408 srccnt[cntord[sn].sno].dom);
409 multcolor(srccnt[cntord[sn].sno].val, sr.rcol);
410 }
411 /* add contribution if hit */
412 addcolor(r->rcol, srccnt[cntord[sn].sno].val);
413 nhits++;
414 source[cntord[sn].sno].nhits++;
415 }
416 /* surface hit rate */
417 if (sn > 0)
418 hwt = (double)nhits / (double)sn;
419 else
420 hwt = 0.5;
421 #ifdef DEBUG
422 sprintf(errmsg, "%d tested, %d untested, %f hit rate\n",
423 sn, ncnts-sn, hwt);
424 eputs(errmsg);
425 #endif
426 /* add in untested sources */
427 for ( ; sn < ncnts; sn++) {
428 prob = hwt * (double)source[cntord[sn].sno].nhits /
429 (double)source[cntord[sn].sno].ntests;
430 scalecolor(srccnt[cntord[sn].sno].val, prob);
431 addcolor(r->rcol, srccnt[cntord[sn].sno].val);
432 }
433 }
434
435
436 #define wrongsource(m, r) (m->otype!=MAT_ILLUM && \
437 r->rsrc>=0 && \
438 source[r->rsrc].so!=r->ro)
439
440 #define badambient(m, r) ((r->crtype&(AMBIENT|SHADOW))==AMBIENT && \
441 !(r->rtype&REFLECTED) && /* hack! */\
442 !(m->otype==MAT_GLOW&&r->rot>m->oargs.farg[3]))
443
444 #define passillum(m, r) (m->otype==MAT_ILLUM && \
445 !(r->rsrc>=0&&source[r->rsrc].so==r->ro))
446
447
448 m_light(m, r) /* ray hit a light source */
449 register OBJREC *m;
450 register RAY *r;
451 {
452 /* check for over-counting */
453 if (wrongsource(m, r) || badambient(m, r))
454 return;
455 /* check for passed illum */
456 if (passillum(m, r)) {
457
458 if (m->oargs.nsargs < 1 || !strcmp(m->oargs.sarg[0], VOIDID))
459 raytrans(r);
460 else
461 rayshade(r, modifier(m->oargs.sarg[0]));
462
463 /* otherwise treat as source */
464 } else {
465 /* check for behind */
466 if (r->rod < 0.0)
467 return;
468 /* get distribution pattern */
469 raytexture(r, m->omod);
470 /* get source color */
471 setcolor(r->rcol, m->oargs.farg[0],
472 m->oargs.farg[1],
473 m->oargs.farg[2]);
474 /* modify value */
475 multcolor(r->rcol, r->pcol);
476 /* assign distance */
477 r->rt = r->rot;
478 }
479 }