ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcsupp.c
Revision: 2.5
Committed: Fri Aug 28 14:07:45 1992 UTC (31 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +33 -108 lines
Log Message:
fixed bug with virtual spotlight sources
moved m_light() back to source.c

File Contents

# Content
1 /* Copyright (c) 1992 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 *)m->os) != NULL)
243 return(ns);
244 if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
245 return(NULL);
246 ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
247 VCOPY(ns->aim, m->oargs.farg+4);
248 if ((ns->flen = normalize(ns->aim)) == 0.0)
249 objerror(m, USER, "zero focus vector");
250 m->os = (char *)ns;
251 return(ns);
252 }
253
254
255 spotout(r, s, dist) /* check if we're outside spot region */
256 register RAY *r;
257 register SPOT *s;
258 int dist;
259 {
260 double d;
261 FVECT vd;
262
263 if (s == NULL)
264 return(0);
265 if (dist) { /* distant source */
266 vd[0] = s->aim[0] - r->rorg[0];
267 vd[1] = s->aim[1] - r->rorg[1];
268 vd[2] = s->aim[2] - r->rorg[2];
269 d = DOT(r->rdir,vd);
270 /* wrong side?
271 if (d <= FTINY)
272 return(1); */
273 d = DOT(vd,vd) - d*d;
274 if (PI*d > s->siz)
275 return(1); /* out */
276 return(0); /* OK */
277 }
278 /* local source */
279 if (s->siz < 2.0*PI * (1.0 + DOT(s->aim,r->rdir)))
280 return(1); /* out */
281 return(0); /* OK */
282 }
283
284
285 double
286 fgetmaxdisk(ocent, op) /* get center and squared radius of face */
287 FVECT ocent;
288 OBJREC *op;
289 {
290 double maxrad2;
291 double d;
292 register int i, j;
293 register FACE *f;
294
295 f = getface(op);
296 if (f->area == 0.)
297 return(0.);
298 for (i = 0; i < 3; i++) {
299 ocent[i] = 0.;
300 for (j = 0; j < f->nv; j++)
301 ocent[i] += VERTEX(f,j)[i];
302 ocent[i] /= (double)f->nv;
303 }
304 d = DOT(ocent,f->norm);
305 for (i = 0; i < 3; i++)
306 ocent[i] += (f->offset - d)*f->norm[i];
307 maxrad2 = 0.;
308 for (j = 0; j < f->nv; j++) {
309 d = dist2(VERTEX(f,j), ocent);
310 if (d > maxrad2)
311 maxrad2 = d;
312 }
313 return(maxrad2);
314 }
315
316
317 double
318 rgetmaxdisk(ocent, op) /* get center and squared radius of ring */
319 FVECT ocent;
320 OBJREC *op;
321 {
322 register CONE *co;
323
324 co = getcone(op, 0);
325 VCOPY(ocent, CO_P0(co));
326 return(CO_R1(co)*CO_R1(co));
327 }
328
329
330 double
331 fgetplaneq(nvec, op) /* get plane equation for face */
332 FVECT nvec;
333 OBJREC *op;
334 {
335 register FACE *fo;
336
337 fo = getface(op);
338 VCOPY(nvec, fo->norm);
339 return(fo->offset);
340 }
341
342
343 double
344 rgetplaneq(nvec, op) /* get plane equation for ring */
345 FVECT nvec;
346 OBJREC *op;
347 {
348 register CONE *co;
349
350 co = getcone(op, 0);
351 VCOPY(nvec, co->ad);
352 return(DOT(nvec, CO_P0(co)));
353 }
354
355
356 commonspot(sp1, sp2, org) /* set sp1 to intersection of sp1 and sp2 */
357 register SPOT *sp1, *sp2;
358 FVECT org;
359 {
360 FVECT cent;
361 double rad2, cos1, cos2;
362
363 cos1 = 1. - sp1->siz/(2.*PI);
364 cos2 = 1. - sp2->siz/(2.*PI);
365 if (sp2->siz >= 2.*PI-FTINY) /* BIG, just check overlap */
366 return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
367 sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
368 /* compute and check disks */
369 rad2 = intercircle(cent, sp1->aim, sp2->aim,
370 1./(cos1*cos1) - 1., 1./(cos2*cos2) - 1.);
371 if (rad2 <= FTINY || normalize(cent) == 0.)
372 return(0);
373 VCOPY(sp1->aim, cent);
374 sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
375 return(1);
376 }
377
378
379 commonbeam(sp1, sp2, dir) /* set sp1 to intersection of sp1 and sp2 */
380 register SPOT *sp1, *sp2;
381 FVECT dir;
382 {
383 FVECT cent, c1, c2;
384 double rad2, d;
385 register int i;
386 /* move centers to common plane */
387 d = DOT(sp1->aim, dir);
388 for (i = 0; i < 3; i++)
389 c1[i] = sp1->aim[i] - d*dir[i];
390 d = DOT(sp2->aim, dir);
391 for (i = 0; i < 3; i++)
392 c2[i] = sp2->aim[i] - d*dir[i];
393 /* compute overlap */
394 rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
395 if (rad2 <= FTINY)
396 return(0);
397 VCOPY(sp1->aim, cent);
398 sp1->siz = PI*rad2;
399 return(1);
400 }
401
402
403 checkspot(sp, nrm) /* check spotlight for behind source */
404 register SPOT *sp; /* spotlight */
405 FVECT nrm; /* source surface normal */
406 {
407 double d, d1;
408
409 d = DOT(sp->aim, nrm);
410 if (d > FTINY) /* center in front? */
411 return(1);
412 /* else check horizon */
413 d1 = 1. - sp->siz/(2.*PI);
414 return(1.-FTINY-d*d < d1*d1);
415 }
416
417
418 double
419 spotdisk(oc, op, sp, pos) /* intersect spot with object op */
420 FVECT oc;
421 OBJREC *op;
422 register SPOT *sp;
423 FVECT pos;
424 {
425 FVECT onorm;
426 double offs, d, dist;
427 register int i;
428
429 offs = getplaneq(onorm, op);
430 d = -DOT(onorm, sp->aim);
431 if (d >= -FTINY && d <= FTINY)
432 return(0.);
433 dist = (DOT(pos, onorm) - offs)/d;
434 if (dist < 0.)
435 return(0.);
436 for (i = 0; i < 3; i++)
437 oc[i] = pos[i] + dist*sp->aim[i];
438 return(sp->siz*dist*dist/PI/(d*d));
439 }
440
441
442 double
443 beamdisk(oc, op, sp, dir) /* intersect beam with object op */
444 FVECT oc;
445 OBJREC *op;
446 register SPOT *sp;
447 FVECT dir;
448 {
449 FVECT onorm;
450 double offs, d, dist;
451 register int i;
452
453 offs = getplaneq(onorm, op);
454 d = -DOT(onorm, dir);
455 if (d >= -FTINY && d <= FTINY)
456 return(0.);
457 dist = (DOT(sp->aim, onorm) - offs)/d;
458 for (i = 0; i < 3; i++)
459 oc[i] = sp->aim[i] + dist*dir[i];
460 return(sp->siz/PI/(d*d));
461 }
462
463
464 double
465 intercircle(cc, c1, c2, r1s, r2s) /* intersect two circles */
466 FVECT cc; /* midpoint (return value) */
467 FVECT c1, c2; /* circle centers */
468 double r1s, r2s; /* radii squared */
469 {
470 double a2, d2, l;
471 FVECT disp;
472 register int i;
473
474 for (i = 0; i < 3; i++)
475 disp[i] = c2[i] - c1[i];
476 d2 = DOT(disp,disp);
477 /* circle within overlap? */
478 if (r1s < r2s) {
479 if (r2s >= r1s + d2) {
480 VCOPY(cc, c1);
481 return(r1s);
482 }
483 } else {
484 if (r1s >= r2s + d2) {
485 VCOPY(cc, c2);
486 return(r2s);
487 }
488 }
489 a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
490 /* no overlap? */
491 if (a2 <= 0.)
492 return(0.);
493 /* overlap, compute center */
494 l = sqrt((r1s - a2)/d2);
495 for (i = 0; i < 3; i++)
496 cc[i] = c1[i] + l*disp[i];
497 return(a2);
498 }
499
500
501 sourcehit(r) /* check to see if ray hit distant source */
502 register RAY *r;
503 {
504 int first, last;
505 register int i;
506
507 if (r->rsrc >= 0) { /* check only one if aimed */
508 first = last = r->rsrc;
509 } else { /* otherwise check all */
510 first = 0; last = nsources-1;
511 }
512 for (i = first; i <= last; i++)
513 if ((source[i].sflags & (SDISTANT|SVIRTUAL)) == SDISTANT)
514 /*
515 * Check to see if ray is within
516 * solid angle of source.
517 */
518 if (2.0*PI * (1.0 - DOT(source[i].sloc,r->rdir))
519 <= source[i].ss2) {
520 r->ro = source[i].so;
521 if (!(source[i].sflags & SSKIP))
522 break;
523 }
524
525 if (r->ro != NULL) {
526 for (i = 0; i < 3; i++)
527 r->ron[i] = -r->rdir[i];
528 r->rod = 1.0;
529 r->rox = NULL;
530 return(1);
531 }
532 return(0);
533 }