ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcsupp.c
Revision: 1.14
Committed: Mon Oct 21 12:58:16 1991 UTC (32 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.13: +95 -13 lines
Log Message:
added source sampling (-ds option)

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] = .886227 * 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 for (i = 0; i < 3; i++)
216 src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
217 src->srad = co->al;
218 src->ss2 = 2.*CO_R0(co)*co->al;
219 /* set sampling vectors */
220 for (i = 0; i < 3; i++)
221 src->ss[SU][i] = .5 * co->al * co->ad[i];
222 src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
223 for (i = 0; i < 3; i++)
224 if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
225 break;
226 src->ss[SV][i] = 1.0;
227 fcross(src->ss[SW], src->ss[SV], co->ad);
228 normalize(src->ss[SW]);
229 for (i = 0; i < 3; i++)
230 src->ss[SW][i] *= .886227 * CO_R0(co);
231 fcross(src->ss[SV], src->ss[SW], co->ad);
232 }
233
234
235 SPOT *
236 makespot(m) /* make a spotlight */
237 register OBJREC *m;
238 {
239 register SPOT *ns;
240
241 if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
242 return(NULL);
243 ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
244 VCOPY(ns->aim, m->oargs.farg+4);
245 if ((ns->flen = normalize(ns->aim)) == 0.0)
246 objerror(m, USER, "zero focus vector");
247 return(ns);
248 }
249
250
251 double
252 fgetmaxdisk(ocent, op) /* get center and squared radius of face */
253 FVECT ocent;
254 OBJREC *op;
255 {
256 double maxrad2;
257 double d;
258 register int i, j;
259 register FACE *f;
260
261 f = getface(op);
262 if (f->area == 0.)
263 return(0.);
264 for (i = 0; i < 3; i++) {
265 ocent[i] = 0.;
266 for (j = 0; j < f->nv; j++)
267 ocent[i] += VERTEX(f,j)[i];
268 ocent[i] /= (double)f->nv;
269 }
270 d = DOT(ocent,f->norm);
271 for (i = 0; i < 3; i++)
272 ocent[i] += (f->offset - d)*f->norm[i];
273 maxrad2 = 0.;
274 for (j = 0; j < f->nv; j++) {
275 d = dist2(VERTEX(f,j), ocent);
276 if (d > maxrad2)
277 maxrad2 = d;
278 }
279 return(maxrad2);
280 }
281
282
283 double
284 rgetmaxdisk(ocent, op) /* get center and squared radius of ring */
285 FVECT ocent;
286 OBJREC *op;
287 {
288 register CONE *co;
289
290 co = getcone(op, 0);
291 VCOPY(ocent, CO_P0(co));
292 return(CO_R1(co)*CO_R1(co));
293 }
294
295
296 double
297 fgetplaneq(nvec, op) /* get plane equation for face */
298 FVECT nvec;
299 OBJREC *op;
300 {
301 register FACE *fo;
302
303 fo = getface(op);
304 VCOPY(nvec, fo->norm);
305 return(fo->offset);
306 }
307
308
309 double
310 rgetplaneq(nvec, op) /* get plane equation for ring */
311 FVECT nvec;
312 OBJREC *op;
313 {
314 register CONE *co;
315
316 co = getcone(op, 0);
317 VCOPY(nvec, co->ad);
318 return(DOT(nvec, CO_P0(co)));
319 }
320
321
322 commonspot(sp1, sp2, org) /* set sp1 to intersection of sp1 and sp2 */
323 register SPOT *sp1, *sp2;
324 FVECT org;
325 {
326 FVECT cent;
327 double rad2, cos1, cos2;
328
329 cos1 = 1. - sp1->siz/(2.*PI);
330 cos2 = 1. - sp2->siz/(2.*PI);
331 if (sp2->siz >= 2.*PI-FTINY) /* BIG, just check overlap */
332 return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
333 sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
334 /* compute and check disks */
335 rad2 = intercircle(cent, sp1->aim, sp2->aim,
336 1./(cos1*cos1) - 1., 1./(cos2*cos2) - 1.);
337 if (rad2 <= FTINY || normalize(cent) == 0.)
338 return(0);
339 VCOPY(sp1->aim, cent);
340 sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
341 return(1);
342 }
343
344
345 commonbeam(sp1, sp2, dir) /* set sp1 to intersection of sp1 and sp2 */
346 register SPOT *sp1, *sp2;
347 FVECT dir;
348 {
349 FVECT cent, c1, c2;
350 double rad2, d;
351 register int i;
352 /* move centers to common plane */
353 d = DOT(sp1->aim, dir);
354 for (i = 0; i < 3; i++)
355 c1[i] = sp1->aim[i] - d*dir[i];
356 d = DOT(sp2->aim, dir);
357 for (i = 0; i < 3; i++)
358 c2[i] = sp2->aim[i] - d*dir[i];
359 /* compute overlap */
360 rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
361 if (rad2 <= FTINY)
362 return(0);
363 VCOPY(sp1->aim, cent);
364 sp1->siz = PI*rad2;
365 return(1);
366 }
367
368
369 checkspot(sp, nrm) /* check spotlight for behind source */
370 register SPOT *sp; /* spotlight */
371 FVECT nrm; /* source surface normal */
372 {
373 double d, d1;
374
375 d = DOT(sp->aim, nrm);
376 if (d > FTINY) /* center in front? */
377 return(1);
378 /* else check horizon */
379 d1 = 1. - sp->siz/(2.*PI);
380 return(1.-FTINY-d*d < d1*d1);
381 }
382
383
384 double
385 spotdisk(oc, op, sp, pos) /* intersect spot with object op */
386 FVECT oc;
387 OBJREC *op;
388 register SPOT *sp;
389 FVECT pos;
390 {
391 FVECT onorm;
392 double offs, d, dist;
393 register int i;
394
395 offs = getplaneq(onorm, op);
396 d = -DOT(onorm, sp->aim);
397 if (d >= -FTINY && d <= FTINY)
398 return(0.);
399 dist = (DOT(pos, onorm) - offs)/d;
400 if (dist < 0.)
401 return(0.);
402 for (i = 0; i < 3; i++)
403 oc[i] = pos[i] + dist*sp->aim[i];
404 return(sp->siz*dist*dist/PI/(d*d));
405 }
406
407
408 double
409 beamdisk(oc, op, sp, dir) /* intersect beam with object op */
410 FVECT oc;
411 OBJREC *op;
412 register SPOT *sp;
413 FVECT dir;
414 {
415 FVECT onorm;
416 double offs, d, dist;
417 register int i;
418
419 offs = getplaneq(onorm, op);
420 d = -DOT(onorm, dir);
421 if (d >= -FTINY && d <= FTINY)
422 return(0.);
423 dist = (DOT(sp->aim, onorm) - offs)/d;
424 for (i = 0; i < 3; i++)
425 oc[i] = sp->aim[i] + dist*dir[i];
426 return(sp->siz/PI/(d*d));
427 }
428
429
430 double
431 intercircle(cc, c1, c2, r1s, r2s) /* intersect two circles */
432 FVECT cc; /* midpoint (return value) */
433 FVECT c1, c2; /* circle centers */
434 double r1s, r2s; /* radii squared */
435 {
436 double a2, d2, l;
437 FVECT disp;
438 register int i;
439
440 for (i = 0; i < 3; i++)
441 disp[i] = c2[i] - c1[i];
442 d2 = DOT(disp,disp);
443 /* circle within overlap? */
444 if (r1s < r2s) {
445 if (r2s >= r1s + d2) {
446 VCOPY(cc, c1);
447 return(r1s);
448 }
449 } else {
450 if (r1s >= r2s + d2) {
451 VCOPY(cc, c2);
452 return(r2s);
453 }
454 }
455 a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
456 /* no overlap? */
457 if (a2 <= 0.)
458 return(0.);
459 /* overlap, compute center */
460 l = sqrt((r1s - a2)/d2);
461 for (i = 0; i < 3; i++)
462 cc[i] = c1[i] + l*disp[i];
463 return(a2);
464 }
465
466
467 sourcehit(r) /* check to see if ray hit distant source */
468 register RAY *r;
469 {
470 int first, last;
471 register int i;
472
473 if (r->rsrc >= 0) { /* check only one if aimed */
474 first = last = r->rsrc;
475 } else { /* otherwise check all */
476 first = 0; last = nsources-1;
477 }
478 for (i = first; i <= last; i++)
479 if ((source[i].sflags & (SDISTANT|SVIRTUAL)) == SDISTANT)
480 /*
481 * Check to see if ray is within
482 * solid angle of source.
483 */
484 if (2.0*PI * (1.0 - DOT(source[i].sloc,r->rdir))
485 <= source[i].ss2) {
486 r->ro = source[i].so;
487 if (!(source[i].sflags & SSKIP))
488 break;
489 }
490
491 if (r->ro != NULL) {
492 for (i = 0; i < 3; i++)
493 r->ron[i] = -r->rdir[i];
494 r->rod = 1.0;
495 r->rox = NULL;
496 return(1);
497 }
498 return(0);
499 }
500
501
502 #define wrongsource(m, r) (r->rsrc>=0 && \
503 source[r->rsrc].so!=r->ro && \
504 (m->otype!=MAT_ILLUM || \
505 objptr(source[r->rsrc].so->omod)->otype==MAT_ILLUM))
506
507 #define distglow(m, r) (m->otype==MAT_GLOW && \
508 r->rot > m->oargs.farg[3])
509
510 #define badambient(m, r) ((r->crtype&(AMBIENT|SHADOW))==AMBIENT && \
511 !distglow(m, r))
512
513 #define passillum(m, r) (m->otype==MAT_ILLUM && \
514 !(r->rsrc>=0&&source[r->rsrc].so==r->ro))
515
516 #define srcignore(m, r) (directinvis && !(r->crtype&SHADOW) && \
517 !distglow(m, r))
518
519
520 m_light(m, r) /* ray hit a light source */
521 register OBJREC *m;
522 register RAY *r;
523 {
524 /* check for over-counting */
525 if (wrongsource(m, r) || badambient(m, r))
526 return;
527 /* check for passed illum */
528 if (passillum(m, r)) {
529 if (m->oargs.nsargs < 1 || !strcmp(m->oargs.sarg[0], VOIDID))
530 raytrans(r);
531 else
532 rayshade(r, modifier(m->oargs.sarg[0]));
533 return;
534 }
535 /* otherwise treat as source */
536 /* check for behind */
537 if (r->rod < 0.0)
538 return;
539 /* check for invisibility */
540 if (srcignore(m, r))
541 return;
542 /* get distribution pattern */
543 raytexture(r, m->omod);
544 /* get source color */
545 setcolor(r->rcol, m->oargs.farg[0],
546 m->oargs.farg[1],
547 m->oargs.farg[2]);
548 /* modify value */
549 multcolor(r->rcol, r->pcol);
550 }