ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcsupp.c
Revision: 2.1
Committed: Tue Nov 12 17:10:22 1991 UTC (32 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.16: +0 -0 lines
Log Message:
updated revision number for release 2.0

File Contents

# Content
1 /* Copyright (c) 1991 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * Support routines for source objects and materials
9 */
10
11 #include "ray.h"
12
13 #include "otypes.h"
14
15 #include "source.h"
16
17 #include "cone.h"
18
19 #include "face.h"
20
21 #define SRCINC 4 /* realloc increment for array */
22
23 SRCREC *source = NULL; /* our list of sources */
24 int nsources = 0; /* the number of sources */
25
26 SRCFUNC sfun[NUMOTYPE]; /* source dispatch table */
27
28
29 initstypes() /* initialize source dispatch table */
30 {
31 extern VSMATERIAL mirror_vs, direct1_vs, direct2_vs;
32 extern int fsetsrc(), ssetsrc(), sphsetsrc(), cylsetsrc(), rsetsrc();
33 extern int nopart(), flatpart(), cylpart();
34 extern double fgetplaneq(), rgetplaneq();
35 extern double fgetmaxdisk(), rgetmaxdisk();
36 static SOBJECT fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
37 static SOBJECT ssobj = {ssetsrc, nopart};
38 static SOBJECT sphsobj = {sphsetsrc, nopart};
39 static SOBJECT cylsobj = {cylsetsrc, cylpart};
40 static SOBJECT rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};
41
42 sfun[MAT_MIRROR].mf = &mirror_vs;
43 sfun[MAT_DIRECT1].mf = &direct1_vs;
44 sfun[MAT_DIRECT2].mf = &direct2_vs;
45 sfun[OBJ_FACE].of = &fsobj;
46 sfun[OBJ_SOURCE].of = &ssobj;
47 sfun[OBJ_SPHERE].of = &sphsobj;
48 sfun[OBJ_CYLINDER].of = &cylsobj;
49 sfun[OBJ_RING].of = &rsobj;
50 }
51
52
53 int
54 newsource() /* allocate new source in our array */
55 {
56 if (nsources == 0)
57 source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
58 else if (nsources%SRCINC == 0)
59 source = (SRCREC *)realloc((char *)source,
60 (unsigned)(nsources+SRCINC)*sizeof(SRCREC));
61 if (source == NULL)
62 return(-1);
63 source[nsources].sflags = 0;
64 source[nsources].nhits = 1;
65 source[nsources].ntests = 2; /* initial hit probability = 1/2 */
66 return(nsources++);
67 }
68
69
70 setflatss(src) /* set sampling for a flat source */
71 register SRCREC *src;
72 {
73 double mult;
74 register int i;
75
76 src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
77 for (i = 0; i < 3; i++)
78 if (src->snorm[i] < 0.6 && src->snorm[i] > -0.6)
79 break;
80 src->ss[SV][i] = 1.0;
81 fcross(src->ss[SU], src->ss[SV], src->snorm);
82 mult = .5 * sqrt( src->ss2 / DOT(src->ss[SU],src->ss[SU]) );
83 for (i = 0; i < 3; i++)
84 src->ss[SU][i] *= mult;
85 fcross(src->ss[SV], src->snorm, src->ss[SU]);
86 }
87
88
89 fsetsrc(src, so) /* set a face as a source */
90 register SRCREC *src;
91 OBJREC *so;
92 {
93 register FACE *f;
94 register int i, j;
95 double d;
96
97 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
98 src->so = so;
99 /* get the face */
100 f = getface(so);
101 /* find the center */
102 for (j = 0; j < 3; j++) {
103 src->sloc[j] = 0.0;
104 for (i = 0; i < f->nv; i++)
105 src->sloc[j] += VERTEX(f,i)[j];
106 src->sloc[j] /= (double)f->nv;
107 }
108 if (!inface(src->sloc, f))
109 objerror(so, USER, "cannot hit center");
110 src->sflags |= SFLAT;
111 VCOPY(src->snorm, f->norm);
112 src->ss2 = f->area;
113 /* find maximum radius */
114 src->srad = 0.;
115 for (i = 0; i < f->nv; i++) {
116 d = dist2(VERTEX(f,i), src->sloc);
117 if (d > src->srad)
118 src->srad = d;
119 }
120 src->srad = sqrt(src->srad);
121 /* compute size vectors */
122 if (f->nv == 4 || (f->nv == 5 && /* parallelogram case */
123 dist2(VERTEX(f,0),VERTEX(f,4)) <= FTINY*FTINY))
124 for (j = 0; j < 3; j++) {
125 src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
126 src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
127 }
128 else
129 setflatss(src);
130 }
131
132
133 ssetsrc(src, so) /* set a source as a source */
134 register SRCREC *src;
135 register OBJREC *so;
136 {
137 double theta;
138
139 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
140 src->so = so;
141 if (so->oargs.nfargs != 4)
142 objerror(so, USER, "bad arguments");
143 src->sflags |= SDISTANT;
144 VCOPY(src->sloc, so->oargs.farg);
145 if (normalize(src->sloc) == 0.0)
146 objerror(so, USER, "zero direction");
147 theta = PI/180.0/2.0 * so->oargs.farg[3];
148 if (theta <= FTINY)
149 objerror(so, USER, "zero size");
150 src->ss2 = 2.0*PI * (1.0 - cos(theta));
151 /* the following is approximate */
152 src->srad = sqrt(src->ss2/PI);
153 VCOPY(src->snorm, src->sloc);
154 setflatss(src); /* hey, whatever works */
155 src->ss[SW][0] = src->ss[SW][1] = src->ss[SW][2] = 0.0;
156 }
157
158
159 sphsetsrc(src, so) /* set a sphere as a source */
160 register SRCREC *src;
161 register OBJREC *so;
162 {
163 register int i;
164
165 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
166 src->so = so;
167 if (so->oargs.nfargs != 4)
168 objerror(so, USER, "bad # arguments");
169 if (so->oargs.farg[3] <= FTINY)
170 objerror(so, USER, "illegal radius");
171 VCOPY(src->sloc, so->oargs.farg);
172 src->srad = so->oargs.farg[3];
173 src->ss2 = PI * src->srad * src->srad;
174 for (i = 0; i < 3; i++)
175 src->ss[SU][i] = src->ss[SV][i] = src->ss[SW][i] = 0.0;
176 for (i = 0; i < 3; i++)
177 src->ss[i][i] = .7236 * so->oargs.farg[3];
178 }
179
180
181 rsetsrc(src, so) /* set a ring (disk) as a source */
182 register SRCREC *src;
183 OBJREC *so;
184 {
185 register CONE *co;
186
187 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
188 src->so = so;
189 /* get the ring */
190 co = getcone(so, 0);
191 VCOPY(src->sloc, CO_P0(co));
192 if (CO_R0(co) > 0.0)
193 objerror(so, USER, "cannot hit center");
194 src->sflags |= SFLAT;
195 VCOPY(src->snorm, co->ad);
196 src->srad = CO_R1(co);
197 src->ss2 = PI * src->srad * src->srad;
198 setflatss(src);
199 }
200
201
202 cylsetsrc(src, so) /* set a cylinder as a source */
203 register SRCREC *src;
204 OBJREC *so;
205 {
206 register CONE *co;
207 register int i;
208
209 src->sa.success = 4*AIMREQT-1; /* bitch on fourth failure */
210 src->so = so;
211 /* get the cylinder */
212 co = getcone(so, 0);
213 if (CO_R0(co) > .2*co->al) /* heuristic constraint */
214 objerror(so, WARNING, "source aspect too small");
215 src->sflags |= SCYL;
216 for (i = 0; i < 3; i++)
217 src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
218 src->srad = .5*co->al;
219 src->ss2 = 2.*CO_R0(co)*co->al;
220 /* set sampling vectors */
221 for (i = 0; i < 3; i++)
222 src->ss[SU][i] = .5 * co->al * co->ad[i];
223 src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
224 for (i = 0; i < 3; i++)
225 if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
226 break;
227 src->ss[SV][i] = 1.0;
228 fcross(src->ss[SW], src->ss[SV], co->ad);
229 normalize(src->ss[SW]);
230 for (i = 0; i < 3; i++)
231 src->ss[SW][i] *= .8559 * CO_R0(co);
232 fcross(src->ss[SV], src->ss[SW], co->ad);
233 }
234
235
236 SPOT *
237 makespot(m) /* make a spotlight */
238 register OBJREC *m;
239 {
240 register SPOT *ns;
241
242 if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
243 return(NULL);
244 ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
245 VCOPY(ns->aim, m->oargs.farg+4);
246 if ((ns->flen = normalize(ns->aim)) == 0.0)
247 objerror(m, USER, "zero focus vector");
248 return(ns);
249 }
250
251
252 double
253 fgetmaxdisk(ocent, op) /* get center and squared radius of face */
254 FVECT ocent;
255 OBJREC *op;
256 {
257 double maxrad2;
258 double d;
259 register int i, j;
260 register FACE *f;
261
262 f = getface(op);
263 if (f->area == 0.)
264 return(0.);
265 for (i = 0; i < 3; i++) {
266 ocent[i] = 0.;
267 for (j = 0; j < f->nv; j++)
268 ocent[i] += VERTEX(f,j)[i];
269 ocent[i] /= (double)f->nv;
270 }
271 d = DOT(ocent,f->norm);
272 for (i = 0; i < 3; i++)
273 ocent[i] += (f->offset - d)*f->norm[i];
274 maxrad2 = 0.;
275 for (j = 0; j < f->nv; j++) {
276 d = dist2(VERTEX(f,j), ocent);
277 if (d > maxrad2)
278 maxrad2 = d;
279 }
280 return(maxrad2);
281 }
282
283
284 double
285 rgetmaxdisk(ocent, op) /* get center and squared radius of ring */
286 FVECT ocent;
287 OBJREC *op;
288 {
289 register CONE *co;
290
291 co = getcone(op, 0);
292 VCOPY(ocent, CO_P0(co));
293 return(CO_R1(co)*CO_R1(co));
294 }
295
296
297 double
298 fgetplaneq(nvec, op) /* get plane equation for face */
299 FVECT nvec;
300 OBJREC *op;
301 {
302 register FACE *fo;
303
304 fo = getface(op);
305 VCOPY(nvec, fo->norm);
306 return(fo->offset);
307 }
308
309
310 double
311 rgetplaneq(nvec, op) /* get plane equation for ring */
312 FVECT nvec;
313 OBJREC *op;
314 {
315 register CONE *co;
316
317 co = getcone(op, 0);
318 VCOPY(nvec, co->ad);
319 return(DOT(nvec, CO_P0(co)));
320 }
321
322
323 commonspot(sp1, sp2, org) /* set sp1 to intersection of sp1 and sp2 */
324 register SPOT *sp1, *sp2;
325 FVECT org;
326 {
327 FVECT cent;
328 double rad2, cos1, cos2;
329
330 cos1 = 1. - sp1->siz/(2.*PI);
331 cos2 = 1. - sp2->siz/(2.*PI);
332 if (sp2->siz >= 2.*PI-FTINY) /* BIG, just check overlap */
333 return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
334 sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
335 /* compute and check disks */
336 rad2 = intercircle(cent, sp1->aim, sp2->aim,
337 1./(cos1*cos1) - 1., 1./(cos2*cos2) - 1.);
338 if (rad2 <= FTINY || normalize(cent) == 0.)
339 return(0);
340 VCOPY(sp1->aim, cent);
341 sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
342 return(1);
343 }
344
345
346 commonbeam(sp1, sp2, dir) /* set sp1 to intersection of sp1 and sp2 */
347 register SPOT *sp1, *sp2;
348 FVECT dir;
349 {
350 FVECT cent, c1, c2;
351 double rad2, d;
352 register int i;
353 /* move centers to common plane */
354 d = DOT(sp1->aim, dir);
355 for (i = 0; i < 3; i++)
356 c1[i] = sp1->aim[i] - d*dir[i];
357 d = DOT(sp2->aim, dir);
358 for (i = 0; i < 3; i++)
359 c2[i] = sp2->aim[i] - d*dir[i];
360 /* compute overlap */
361 rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
362 if (rad2 <= FTINY)
363 return(0);
364 VCOPY(sp1->aim, cent);
365 sp1->siz = PI*rad2;
366 return(1);
367 }
368
369
370 checkspot(sp, nrm) /* check spotlight for behind source */
371 register SPOT *sp; /* spotlight */
372 FVECT nrm; /* source surface normal */
373 {
374 double d, d1;
375
376 d = DOT(sp->aim, nrm);
377 if (d > FTINY) /* center in front? */
378 return(1);
379 /* else check horizon */
380 d1 = 1. - sp->siz/(2.*PI);
381 return(1.-FTINY-d*d < d1*d1);
382 }
383
384
385 double
386 spotdisk(oc, op, sp, pos) /* intersect spot with object op */
387 FVECT oc;
388 OBJREC *op;
389 register SPOT *sp;
390 FVECT pos;
391 {
392 FVECT onorm;
393 double offs, d, dist;
394 register int i;
395
396 offs = getplaneq(onorm, op);
397 d = -DOT(onorm, sp->aim);
398 if (d >= -FTINY && d <= FTINY)
399 return(0.);
400 dist = (DOT(pos, onorm) - offs)/d;
401 if (dist < 0.)
402 return(0.);
403 for (i = 0; i < 3; i++)
404 oc[i] = pos[i] + dist*sp->aim[i];
405 return(sp->siz*dist*dist/PI/(d*d));
406 }
407
408
409 double
410 beamdisk(oc, op, sp, dir) /* intersect beam with object op */
411 FVECT oc;
412 OBJREC *op;
413 register SPOT *sp;
414 FVECT dir;
415 {
416 FVECT onorm;
417 double offs, d, dist;
418 register int i;
419
420 offs = getplaneq(onorm, op);
421 d = -DOT(onorm, dir);
422 if (d >= -FTINY && d <= FTINY)
423 return(0.);
424 dist = (DOT(sp->aim, onorm) - offs)/d;
425 for (i = 0; i < 3; i++)
426 oc[i] = sp->aim[i] + dist*dir[i];
427 return(sp->siz/PI/(d*d));
428 }
429
430
431 double
432 intercircle(cc, c1, c2, r1s, r2s) /* intersect two circles */
433 FVECT cc; /* midpoint (return value) */
434 FVECT c1, c2; /* circle centers */
435 double r1s, r2s; /* radii squared */
436 {
437 double a2, d2, l;
438 FVECT disp;
439 register int i;
440
441 for (i = 0; i < 3; i++)
442 disp[i] = c2[i] - c1[i];
443 d2 = DOT(disp,disp);
444 /* circle within overlap? */
445 if (r1s < r2s) {
446 if (r2s >= r1s + d2) {
447 VCOPY(cc, c1);
448 return(r1s);
449 }
450 } else {
451 if (r1s >= r2s + d2) {
452 VCOPY(cc, c2);
453 return(r2s);
454 }
455 }
456 a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
457 /* no overlap? */
458 if (a2 <= 0.)
459 return(0.);
460 /* overlap, compute center */
461 l = sqrt((r1s - a2)/d2);
462 for (i = 0; i < 3; i++)
463 cc[i] = c1[i] + l*disp[i];
464 return(a2);
465 }
466
467
468 sourcehit(r) /* check to see if ray hit distant source */
469 register RAY *r;
470 {
471 int first, last;
472 register int i;
473
474 if (r->rsrc >= 0) { /* check only one if aimed */
475 first = last = r->rsrc;
476 } else { /* otherwise check all */
477 first = 0; last = nsources-1;
478 }
479 for (i = first; i <= last; i++)
480 if ((source[i].sflags & (SDISTANT|SVIRTUAL)) == SDISTANT)
481 /*
482 * Check to see if ray is within
483 * solid angle of source.
484 */
485 if (2.0*PI * (1.0 - DOT(source[i].sloc,r->rdir))
486 <= source[i].ss2) {
487 r->ro = source[i].so;
488 if (!(source[i].sflags & SSKIP))
489 break;
490 }
491
492 if (r->ro != NULL) {
493 for (i = 0; i < 3; i++)
494 r->ron[i] = -r->rdir[i];
495 r->rod = 1.0;
496 r->rox = NULL;
497 return(1);
498 }
499 return(0);
500 }
501
502
503 /****************************************************************
504 * The following macros were separated from the m_light() routine
505 * because they are very nasty and difficult to understand.
506 */
507
508 /* wrongillum *
509 *
510 * We cannot allow an illum to pass to another illum, because that
511 * would almost certainly constitute overcounting.
512 * However, we do allow an illum to pass to another illum
513 * that is actually going to relay to a virtual light source.
514 */
515
516 #define wrongillum(m, r) (!(source[r->rsrc].sflags&SVIRTUAL) && \
517 objptr(source[r->rsrc].so->omod)->otype==MAT_ILLUM)
518
519 /* wrongsource *
520 *
521 * This source is the wrong source (ie. overcounted) if we are
522 * aimed to a different source than the one we hit and the one
523 * we hit is not an illum which should be passed.
524 */
525
526 #define wrongsource(m, r) (r->rsrc>=0 && source[r->rsrc].so!=r->ro && \
527 (m->otype!=MAT_ILLUM || wrongillum(m,r)))
528
529 /* distglow *
530 *
531 * A distant glow is an object that sometimes acts as a light source,
532 * but is too far away from the test point to be one in this case.
533 */
534
535 #define distglow(m, r) (m->otype==MAT_GLOW && \
536 r->rot > m->oargs.farg[3])
537
538 /* badambient *
539 *
540 * We must avoid including counting light sources in the ambient calculation,
541 * since the direct component is handled separately. Therefore, any
542 * ambient ray which hits an active light source must be discarded.
543 */
544
545 #define badambient(m, r) ((r->crtype&(AMBIENT|SHADOW))==AMBIENT && \
546 !distglow(m, r))
547
548 /* passillum *
549 *
550 * An illum passes to another material type when we didn't hit it
551 * on purpose (as part of a direct calculation), or it is relaying
552 * a virtual light source.
553 */
554
555 #define passillum(m, r) (m->otype==MAT_ILLUM && \
556 (r->rsrc<0 || source[r->rsrc].so!=r->ro || \
557 source[r->rsrc].sflags&SVIRTUAL))
558
559 /* srcignore *
560 *
561 * The -di flag renders light sources invisible, and here is the test.
562 */
563
564 #define srcignore(m, r) (directinvis && !(r->crtype&SHADOW) && \
565 !distglow(m, r))
566
567
568 m_light(m, r) /* ray hit a light source */
569 register OBJREC *m;
570 register RAY *r;
571 {
572 /* check for over-counting */
573 if (wrongsource(m, r) || badambient(m, r))
574 return;
575 /* check for passed illum */
576 if (passillum(m, r)) {
577 if (m->oargs.nsargs < 1 || !strcmp(m->oargs.sarg[0], VOIDID))
578 raytrans(r);
579 else
580 rayshade(r, modifier(m->oargs.sarg[0]));
581 return;
582 }
583 /* otherwise treat as source */
584 /* check for behind */
585 if (r->rod < 0.0)
586 return;
587 /* check for invisibility */
588 if (srcignore(m, r))
589 return;
590 /* get distribution pattern */
591 raytexture(r, m->omod);
592 /* get source color */
593 setcolor(r->rcol, m->oargs.farg[0],
594 m->oargs.farg[1],
595 m->oargs.farg[2]);
596 /* modify value */
597 multcolor(r->rcol, r->pcol);
598 }