ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcsupp.c
Revision: 2.9
Committed: Sat Feb 22 02:07:29 2003 UTC (21 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.8: +71 -8 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id$";
3 #endif
4 /*
5 * Support routines for source objects and materials
6 *
7 * External symbols declared in source.h
8 */
9
10 /* ====================================================================
11 * The Radiance Software License, Version 1.0
12 *
13 * Copyright (c) 1990 - 2002 The Regents of the University of California,
14 * through Lawrence Berkeley National Laboratory. All rights reserved.
15 *
16 * Redistribution and use in source and binary forms, with or without
17 * modification, are permitted provided that the following conditions
18 * are met:
19 *
20 * 1. Redistributions of source code must retain the above copyright
21 * notice, this list of conditions and the following disclaimer.
22 *
23 * 2. Redistributions in binary form must reproduce the above copyright
24 * notice, this list of conditions and the following disclaimer in
25 * the documentation and/or other materials provided with the
26 * distribution.
27 *
28 * 3. The end-user documentation included with the redistribution,
29 * if any, must include the following acknowledgment:
30 * "This product includes Radiance software
31 * (http://radsite.lbl.gov/)
32 * developed by the Lawrence Berkeley National Laboratory
33 * (http://www.lbl.gov/)."
34 * Alternately, this acknowledgment may appear in the software itself,
35 * if and wherever such third-party acknowledgments normally appear.
36 *
37 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
38 * and "The Regents of the University of California" must
39 * not be used to endorse or promote products derived from this
40 * software without prior written permission. For written
41 * permission, please contact [email protected].
42 *
43 * 5. Products derived from this software may not be called "Radiance",
44 * nor may "Radiance" appear in their name, without prior written
45 * permission of Lawrence Berkeley National Laboratory.
46 *
47 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
48 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
49 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
50 * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
51 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
52 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
53 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
54 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
55 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
56 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
57 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
58 * SUCH DAMAGE.
59 * ====================================================================
60 *
61 * This software consists of voluntary contributions made by many
62 * individuals on behalf of Lawrence Berkeley National Laboratory. For more
63 * information on Lawrence Berkeley National Laboratory, please see
64 * <http://www.lbl.gov/>.
65 */
66
67 #include "ray.h"
68
69 #include "otypes.h"
70
71 #include "source.h"
72
73 #include "cone.h"
74
75 #include "face.h"
76
77 #define SRCINC 4 /* realloc increment for array */
78
79 SRCREC *source = NULL; /* our list of sources */
80 int nsources = 0; /* the number of sources */
81
82 SRCFUNC sfun[NUMOTYPE]; /* source dispatch table */
83
84
85 void
86 initstypes() /* initialize source dispatch table */
87 {
88 extern VSMATERIAL mirror_vs, direct1_vs, direct2_vs;
89 static SOBJECT fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
90 static SOBJECT ssobj = {ssetsrc, nopart};
91 static SOBJECT sphsobj = {sphsetsrc, nopart};
92 static SOBJECT cylsobj = {cylsetsrc, cylpart};
93 static SOBJECT rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};
94
95 sfun[MAT_MIRROR].mf = &mirror_vs;
96 sfun[MAT_DIRECT1].mf = &direct1_vs;
97 sfun[MAT_DIRECT2].mf = &direct2_vs;
98 sfun[OBJ_FACE].of = &fsobj;
99 sfun[OBJ_SOURCE].of = &ssobj;
100 sfun[OBJ_SPHERE].of = &sphsobj;
101 sfun[OBJ_CYLINDER].of = &cylsobj;
102 sfun[OBJ_RING].of = &rsobj;
103 }
104
105
106 int
107 newsource() /* allocate new source in our array */
108 {
109 if (nsources == 0)
110 source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
111 else if (nsources%SRCINC == 0)
112 source = (SRCREC *)realloc((char *)source,
113 (unsigned)(nsources+SRCINC)*sizeof(SRCREC));
114 if (source == NULL)
115 return(-1);
116 source[nsources].sflags = 0;
117 source[nsources].nhits = 1;
118 source[nsources].ntests = 2; /* initial hit probability = 1/2 */
119 return(nsources++);
120 }
121
122
123 void
124 setflatss(src) /* set sampling for a flat source */
125 register SRCREC *src;
126 {
127 double mult;
128 register int i;
129
130 src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
131 for (i = 0; i < 3; i++)
132 if (src->snorm[i] < 0.6 && src->snorm[i] > -0.6)
133 break;
134 src->ss[SV][i] = 1.0;
135 fcross(src->ss[SU], src->ss[SV], src->snorm);
136 mult = .5 * sqrt( src->ss2 / DOT(src->ss[SU],src->ss[SU]) );
137 for (i = 0; i < 3; i++)
138 src->ss[SU][i] *= mult;
139 fcross(src->ss[SV], src->snorm, src->ss[SU]);
140 }
141
142
143 void
144 fsetsrc(src, so) /* set a face as a source */
145 register SRCREC *src;
146 OBJREC *so;
147 {
148 register FACE *f;
149 register int i, j;
150 double d;
151
152 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
153 src->so = so;
154 /* get the face */
155 f = getface(so);
156 /* find the center */
157 for (j = 0; j < 3; j++) {
158 src->sloc[j] = 0.0;
159 for (i = 0; i < f->nv; i++)
160 src->sloc[j] += VERTEX(f,i)[j];
161 src->sloc[j] /= (double)f->nv;
162 }
163 if (!inface(src->sloc, f))
164 objerror(so, USER, "cannot hit center");
165 src->sflags |= SFLAT;
166 VCOPY(src->snorm, f->norm);
167 src->ss2 = f->area;
168 /* find maximum radius */
169 src->srad = 0.;
170 for (i = 0; i < f->nv; i++) {
171 d = dist2(VERTEX(f,i), src->sloc);
172 if (d > src->srad)
173 src->srad = d;
174 }
175 src->srad = sqrt(src->srad);
176 /* compute size vectors */
177 if (f->nv == 4) /* parallelogram case */
178 for (j = 0; j < 3; j++) {
179 src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
180 src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
181 }
182 else
183 setflatss(src);
184 }
185
186
187 void
188 ssetsrc(src, so) /* set a source as a source */
189 register SRCREC *src;
190 register OBJREC *so;
191 {
192 double theta;
193
194 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
195 src->so = so;
196 if (so->oargs.nfargs != 4)
197 objerror(so, USER, "bad arguments");
198 src->sflags |= SDISTANT;
199 VCOPY(src->sloc, so->oargs.farg);
200 if (normalize(src->sloc) == 0.0)
201 objerror(so, USER, "zero direction");
202 theta = PI/180.0/2.0 * so->oargs.farg[3];
203 if (theta <= FTINY)
204 objerror(so, USER, "zero size");
205 src->ss2 = 2.0*PI * (1.0 - cos(theta));
206 /* the following is approximate */
207 src->srad = sqrt(src->ss2/PI);
208 VCOPY(src->snorm, src->sloc);
209 setflatss(src); /* hey, whatever works */
210 src->ss[SW][0] = src->ss[SW][1] = src->ss[SW][2] = 0.0;
211 }
212
213
214 void
215 sphsetsrc(src, so) /* set a sphere as a source */
216 register SRCREC *src;
217 register OBJREC *so;
218 {
219 register int i;
220
221 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
222 src->so = so;
223 if (so->oargs.nfargs != 4)
224 objerror(so, USER, "bad # arguments");
225 if (so->oargs.farg[3] <= FTINY)
226 objerror(so, USER, "illegal radius");
227 VCOPY(src->sloc, so->oargs.farg);
228 src->srad = so->oargs.farg[3];
229 src->ss2 = PI * src->srad * src->srad;
230 for (i = 0; i < 3; i++)
231 src->ss[SU][i] = src->ss[SV][i] = src->ss[SW][i] = 0.0;
232 for (i = 0; i < 3; i++)
233 src->ss[i][i] = .7236 * so->oargs.farg[3];
234 }
235
236
237 void
238 rsetsrc(src, so) /* set a ring (disk) as a source */
239 register SRCREC *src;
240 OBJREC *so;
241 {
242 register CONE *co;
243
244 src->sa.success = 2*AIMREQT-1; /* bitch on second failure */
245 src->so = so;
246 /* get the ring */
247 co = getcone(so, 0);
248 VCOPY(src->sloc, CO_P0(co));
249 if (CO_R0(co) > 0.0)
250 objerror(so, USER, "cannot hit center");
251 src->sflags |= SFLAT;
252 VCOPY(src->snorm, co->ad);
253 src->srad = CO_R1(co);
254 src->ss2 = PI * src->srad * src->srad;
255 setflatss(src);
256 }
257
258
259 void
260 cylsetsrc(src, so) /* set a cylinder as a source */
261 register SRCREC *src;
262 OBJREC *so;
263 {
264 register CONE *co;
265 register int i;
266
267 src->sa.success = 4*AIMREQT-1; /* bitch on fourth failure */
268 src->so = so;
269 /* get the cylinder */
270 co = getcone(so, 0);
271 if (CO_R0(co) > .2*co->al) /* heuristic constraint */
272 objerror(so, WARNING, "source aspect too small");
273 src->sflags |= SCYL;
274 for (i = 0; i < 3; i++)
275 src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
276 src->srad = .5*co->al;
277 src->ss2 = 2.*CO_R0(co)*co->al;
278 /* set sampling vectors */
279 for (i = 0; i < 3; i++)
280 src->ss[SU][i] = .5 * co->al * co->ad[i];
281 src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
282 for (i = 0; i < 3; i++)
283 if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
284 break;
285 src->ss[SV][i] = 1.0;
286 fcross(src->ss[SW], src->ss[SV], co->ad);
287 normalize(src->ss[SW]);
288 for (i = 0; i < 3; i++)
289 src->ss[SW][i] *= .8559 * CO_R0(co);
290 fcross(src->ss[SV], src->ss[SW], co->ad);
291 }
292
293
294 SPOT *
295 makespot(m) /* make a spotlight */
296 register OBJREC *m;
297 {
298 register SPOT *ns;
299
300 if ((ns = (SPOT *)m->os) != NULL)
301 return(ns);
302 if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
303 return(NULL);
304 ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
305 VCOPY(ns->aim, m->oargs.farg+4);
306 if ((ns->flen = normalize(ns->aim)) == 0.0)
307 objerror(m, USER, "zero focus vector");
308 m->os = (char *)ns;
309 return(ns);
310 }
311
312
313 int
314 spotout(r, s) /* check if we're outside spot region */
315 register RAY *r;
316 register SPOT *s;
317 {
318 double d;
319 FVECT vd;
320
321 if (s == NULL)
322 return(0);
323 if (s->flen < -FTINY) { /* distant source */
324 vd[0] = s->aim[0] - r->rorg[0];
325 vd[1] = s->aim[1] - r->rorg[1];
326 vd[2] = s->aim[2] - r->rorg[2];
327 d = DOT(r->rdir,vd);
328 /* wrong side?
329 if (d <= FTINY)
330 return(1); */
331 d = DOT(vd,vd) - d*d;
332 if (PI*d > s->siz)
333 return(1); /* out */
334 return(0); /* OK */
335 }
336 /* local source */
337 if (s->siz < 2.0*PI * (1.0 + DOT(s->aim,r->rdir)))
338 return(1); /* out */
339 return(0); /* OK */
340 }
341
342
343 double
344 fgetmaxdisk(ocent, op) /* get center and squared radius of face */
345 FVECT ocent;
346 OBJREC *op;
347 {
348 double maxrad2;
349 double d;
350 register int i, j;
351 register FACE *f;
352
353 f = getface(op);
354 if (f->area == 0.)
355 return(0.);
356 for (i = 0; i < 3; i++) {
357 ocent[i] = 0.;
358 for (j = 0; j < f->nv; j++)
359 ocent[i] += VERTEX(f,j)[i];
360 ocent[i] /= (double)f->nv;
361 }
362 d = DOT(ocent,f->norm);
363 for (i = 0; i < 3; i++)
364 ocent[i] += (f->offset - d)*f->norm[i];
365 maxrad2 = 0.;
366 for (j = 0; j < f->nv; j++) {
367 d = dist2(VERTEX(f,j), ocent);
368 if (d > maxrad2)
369 maxrad2 = d;
370 }
371 return(maxrad2);
372 }
373
374
375 double
376 rgetmaxdisk(ocent, op) /* get center and squared radius of ring */
377 FVECT ocent;
378 OBJREC *op;
379 {
380 register CONE *co;
381
382 co = getcone(op, 0);
383 VCOPY(ocent, CO_P0(co));
384 return(CO_R1(co)*CO_R1(co));
385 }
386
387
388 double
389 fgetplaneq(nvec, op) /* get plane equation for face */
390 FVECT nvec;
391 OBJREC *op;
392 {
393 register FACE *fo;
394
395 fo = getface(op);
396 VCOPY(nvec, fo->norm);
397 return(fo->offset);
398 }
399
400
401 double
402 rgetplaneq(nvec, op) /* get plane equation for ring */
403 FVECT nvec;
404 OBJREC *op;
405 {
406 register CONE *co;
407
408 co = getcone(op, 0);
409 VCOPY(nvec, co->ad);
410 return(DOT(nvec, CO_P0(co)));
411 }
412
413
414 int
415 commonspot(sp1, sp2, org) /* set sp1 to intersection of sp1 and sp2 */
416 register SPOT *sp1, *sp2;
417 FVECT org;
418 {
419 FVECT cent;
420 double rad2, cos1, cos2;
421
422 cos1 = 1. - sp1->siz/(2.*PI);
423 cos2 = 1. - sp2->siz/(2.*PI);
424 if (sp2->siz >= 2.*PI-FTINY) /* BIG, just check overlap */
425 return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
426 sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
427 /* compute and check disks */
428 rad2 = intercircle(cent, sp1->aim, sp2->aim,
429 1./(cos1*cos1) - 1., 1./(cos2*cos2) - 1.);
430 if (rad2 <= FTINY || normalize(cent) == 0.)
431 return(0);
432 VCOPY(sp1->aim, cent);
433 sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
434 return(1);
435 }
436
437
438 int
439 commonbeam(sp1, sp2, dir) /* set sp1 to intersection of sp1 and sp2 */
440 register SPOT *sp1, *sp2;
441 FVECT dir;
442 {
443 FVECT cent, c1, c2;
444 double rad2, d;
445 register int i;
446 /* move centers to common plane */
447 d = DOT(sp1->aim, dir);
448 for (i = 0; i < 3; i++)
449 c1[i] = sp1->aim[i] - d*dir[i];
450 d = DOT(sp2->aim, dir);
451 for (i = 0; i < 3; i++)
452 c2[i] = sp2->aim[i] - d*dir[i];
453 /* compute overlap */
454 rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
455 if (rad2 <= FTINY)
456 return(0);
457 VCOPY(sp1->aim, cent);
458 sp1->siz = PI*rad2;
459 return(1);
460 }
461
462
463 int
464 checkspot(sp, nrm) /* check spotlight for behind source */
465 register SPOT *sp; /* spotlight */
466 FVECT nrm; /* source surface normal */
467 {
468 double d, d1;
469
470 d = DOT(sp->aim, nrm);
471 if (d > FTINY) /* center in front? */
472 return(1);
473 /* else check horizon */
474 d1 = 1. - sp->siz/(2.*PI);
475 return(1.-FTINY-d*d < d1*d1);
476 }
477
478
479 double
480 spotdisk(oc, op, sp, pos) /* intersect spot with object op */
481 FVECT oc;
482 OBJREC *op;
483 register SPOT *sp;
484 FVECT pos;
485 {
486 FVECT onorm;
487 double offs, d, dist;
488 register int i;
489
490 offs = getplaneq(onorm, op);
491 d = -DOT(onorm, sp->aim);
492 if (d >= -FTINY && d <= FTINY)
493 return(0.);
494 dist = (DOT(pos, onorm) - offs)/d;
495 if (dist < 0.)
496 return(0.);
497 for (i = 0; i < 3; i++)
498 oc[i] = pos[i] + dist*sp->aim[i];
499 return(sp->siz*dist*dist/PI/(d*d));
500 }
501
502
503 double
504 beamdisk(oc, op, sp, dir) /* intersect beam with object op */
505 FVECT oc;
506 OBJREC *op;
507 register SPOT *sp;
508 FVECT dir;
509 {
510 FVECT onorm;
511 double offs, d, dist;
512 register int i;
513
514 offs = getplaneq(onorm, op);
515 d = -DOT(onorm, dir);
516 if (d >= -FTINY && d <= FTINY)
517 return(0.);
518 dist = (DOT(sp->aim, onorm) - offs)/d;
519 for (i = 0; i < 3; i++)
520 oc[i] = sp->aim[i] + dist*dir[i];
521 return(sp->siz/PI/(d*d));
522 }
523
524
525 double
526 intercircle(cc, c1, c2, r1s, r2s) /* intersect two circles */
527 FVECT cc; /* midpoint (return value) */
528 FVECT c1, c2; /* circle centers */
529 double r1s, r2s; /* radii squared */
530 {
531 double a2, d2, l;
532 FVECT disp;
533 register int i;
534
535 for (i = 0; i < 3; i++)
536 disp[i] = c2[i] - c1[i];
537 d2 = DOT(disp,disp);
538 /* circle within overlap? */
539 if (r1s < r2s) {
540 if (r2s >= r1s + d2) {
541 VCOPY(cc, c1);
542 return(r1s);
543 }
544 } else {
545 if (r1s >= r2s + d2) {
546 VCOPY(cc, c2);
547 return(r2s);
548 }
549 }
550 a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
551 /* no overlap? */
552 if (a2 <= 0.)
553 return(0.);
554 /* overlap, compute center */
555 l = sqrt((r1s - a2)/d2);
556 for (i = 0; i < 3; i++)
557 cc[i] = c1[i] + l*disp[i];
558 return(a2);
559 }