ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/virtuals.c
Revision: 1.6
Committed: Fri Jun 21 13:25:42 1991 UTC (32 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.5: +79 -73 lines
Log Message:
added test for relay object on same side as virtual source

File Contents

# User Rev Content
1 greg 1.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     * Routines for simulating virtual light sources
9     * Thus far, we only support planar mirrors.
10     */
11    
12     #include "ray.h"
13    
14 greg 1.3 #include "otypes.h"
15    
16 greg 1.1 #include "source.h"
17    
18    
19 greg 1.6 double intercircle(), getdisk();
20 greg 1.1
21     static OBJECT *vobject; /* virtual source objects */
22     static int nvobjects = 0; /* number of virtual source objects */
23    
24    
25     markvirtuals() /* find and mark virtual sources */
26     {
27     register OBJREC *o;
28     register int i;
29     /* check number of direct relays */
30     if (directrelay <= 0)
31     return;
32     /* find virtual source objects */
33     for (i = 0; i < nobjects; i++) {
34     o = objptr(i);
35 greg 1.3 if (!issurface(o->otype) || o->omod == OVOID)
36 greg 1.1 continue;
37     if (!isvlight(objptr(o->omod)->otype))
38     continue;
39 greg 1.3 if (sfun[o->otype].of == NULL ||
40     sfun[o->otype].of->getpleq == NULL)
41     objerror(o, USER, "illegal material");
42 greg 1.1 if (nvobjects == 0)
43     vobject = (OBJECT *)malloc(sizeof(OBJECT));
44     else
45     vobject = (OBJECT *)realloc((char *)vobject,
46     (unsigned)(nvobjects+1)*sizeof(OBJECT));
47     if (vobject == NULL)
48     error(SYSTEM, "out of memory in addvirtuals");
49     vobject[nvobjects++] = i;
50     }
51     if (nvobjects == 0)
52     return;
53 greg 1.4 #ifdef DEBUG
54     fprintf(stderr, "found %d virtual source objects\n", nvobjects);
55     #endif
56 greg 1.1 /* append virtual sources */
57     for (i = nsources; i-- > 0; )
58     if (!(source[i].sflags & SSKIP))
59 greg 1.4 addvirtuals(i, directrelay);
60 greg 1.1 /* done with our object list */
61     free((char *)vobject);
62     nvobjects = 0;
63     }
64    
65    
66 greg 1.4 addvirtuals(sn, nr) /* add virtuals associated with source */
67     int sn;
68 greg 1.1 int nr;
69     {
70     register int i;
71     /* check relay limit first */
72     if (nr <= 0)
73     return;
74     /* check each virtual object for projection */
75     for (i = 0; i < nvobjects; i++)
76 greg 1.3 /* vproject() calls us recursively */
77 greg 1.4 vproject(objptr(vobject[i]), sn, nr-1);
78 greg 1.1 }
79    
80    
81 greg 1.4 vproject(o, sn, n) /* create projected source(s) if they exist */
82 greg 1.3 OBJREC *o;
83 greg 1.4 int sn;
84 greg 1.3 int n;
85     {
86     register int i;
87     register VSMATERIAL *vsmat;
88     MAT4 proj;
89 greg 1.4 int ns;
90    
91     if (o == source[sn].so) /* objects cannot project themselves */
92     return;
93 greg 1.3 /* get virtual source material */
94     vsmat = sfun[objptr(o->omod)->otype].mf;
95     /* project virtual sources */
96     for (i = 0; i < vsmat->nproj; i++)
97 greg 1.4 if ((*vsmat->vproj)(proj, o, &source[sn], i))
98     if ((ns = makevsrc(o, sn, proj)) >= 0) {
99     #ifdef DEBUG
100 greg 1.6 virtverb(ns, stderr);
101 greg 1.4 #endif
102 greg 1.3 addvirtuals(ns, n);
103 greg 1.4 }
104 greg 1.3 }
105    
106    
107 greg 1.4 int
108     makevsrc(op, sn, pm) /* make virtual source if reasonable */
109 greg 1.1 OBJREC *op;
110 greg 1.4 register int sn;
111 greg 1.1 MAT4 pm;
112     {
113 greg 1.6 FVECT nsloc, nsnorm, ocent;
114     double maxrad2;
115 greg 1.3 int nsflags;
116     double d1;
117 greg 1.1 SPOT theirspot, ourspot;
118     register int i;
119 greg 1.3
120 greg 1.6 nsflags = source[sn].sflags | (SVIRTUAL|SSPOT|SFOLLOW);
121 greg 1.1 /* get object center and max. radius */
122 greg 1.6 maxrad2 = getdisk(ocent, op, sn);
123     if (maxrad2 <= FTINY) /* too small? */
124     return(-1);
125 greg 1.1 /* get location and spot */
126 greg 1.4 if (source[sn].sflags & SDISTANT) { /* distant source */
127     if (source[sn].sflags & SPROX)
128 greg 1.5 return(-1); /* should never get here! */
129 greg 1.4 multv3(nsloc, source[sn].sloc, pm);
130 greg 1.6 VCOPY(ourspot.aim, ocent);
131     ourspot.siz = PI*maxrad2;
132     ourspot.flen = 0.;
133 greg 1.4 if (source[sn].sflags & SSPOT) {
134     copystruct(&theirspot, source[sn].sl.s);
135     multp3(theirspot.aim, source[sn].sl.s->aim, pm);
136 greg 1.6 if (!commonbeam(&ourspot, &theirspot, nsloc))
137 greg 1.5 return(-1); /* no overlap */
138 greg 1.1 }
139     } else { /* local source */
140 greg 1.4 multp3(nsloc, source[sn].sloc, pm);
141 greg 1.6 for (i = 0; i < 3; i++)
142     ourspot.aim[i] = ocent[i] - nsloc[i];
143     if ((d1 = normalize(ourspot.aim)) == 0.)
144     return(-1); /* at source!! */
145     if (source[sn].sflags & SPROX && d1 > source[sn].sl.prox)
146     return(-1); /* too far away */
147     ourspot.siz = 2.*PI*(1. - d1/sqrt(d1*d1+maxrad2));
148     ourspot.flen = 0.;
149 greg 1.4 if (source[sn].sflags & SSPOT) {
150     copystruct(&theirspot, source[sn].sl.s);
151     multv3(theirspot.aim, source[sn].sl.s->aim, pm);
152 greg 1.6 if (!commonspot(&ourspot, &theirspot, nsloc))
153     return(-1); /* no overlap */
154     ourspot.flen = theirspot.flen;
155 greg 1.1 }
156 greg 1.4 if (source[sn].sflags & SFLAT) { /* behind source? */
157     multv3(nsnorm, source[sn].snorm, pm);
158 greg 1.6 if (checkspot(&ourspot, nsnorm) < 0)
159 greg 1.5 return(-1);
160 greg 1.1 }
161     }
162     /* everything is OK, make source */
163 greg 1.6 if ((i = newsource()) < 0)
164 greg 1.1 goto memerr;
165 greg 1.6 source[i].sflags = nsflags;
166     VCOPY(source[i].sloc, nsloc);
167 greg 1.3 if (nsflags & SFLAT)
168 greg 1.6 VCOPY(source[i].snorm, nsnorm);
169     source[i].ss = source[sn].ss; source[i].ss2 = source[sn].ss2;
170     if ((source[i].sl.s = (SPOT *)malloc(sizeof(SPOT))) == NULL)
171     goto memerr;
172     copystruct(source[i].sl.s, &ourspot);
173 greg 1.3 if (nsflags & SPROX)
174 greg 1.6 source[i].sl.prox = source[sn].sl.prox;
175     source[i].sa.svnext = sn;
176     source[i].so = op;
177     return(i);
178 greg 1.1 memerr:
179     error(SYSTEM, "out of memory in makevsrc");
180     }
181    
182    
183 greg 1.6 double
184     getdisk(oc, op, sn) /* get visible object disk */
185     FVECT oc;
186     OBJREC *op;
187     register int sn;
188     {
189     double rad2, roffs, offs, d, rd, rdoto;
190     FVECT rnrm, nrm;
191     /* first, use object getdisk function */
192     rad2 = (*sfun[op->otype].of->getdisk)(oc, op);
193     if (!(source[sn].sflags & SVIRTUAL))
194     return(rad2); /* all done for normal source */
195     /* check for correct side of relay surface */
196     roffs = (*sfun[source[sn].so->otype].of->getpleq)(rnrm, source[sn].so);
197     rd = DOT(rnrm, source[sn].sloc); /* source projection */
198     if (!(source[sn].sflags & SDISTANT))
199     rd -= roffs;
200     d = DOT(rnrm, oc) - roffs; /* disk distance to relay plane */
201     if ((d > 0.) ^ (rd > 0.))
202     return(rad2); /* OK if opposite sides */
203     if (d*d >= rad2)
204     return(.0); /* no relay is possible */
205     /* we need a closer look */
206     offs = (*sfun[op->otype].of->getpleq)(nrm, op);
207     rdoto = DOT(rnrm, nrm);
208     if (d*d >= rad2*(1.-rdoto*rdoto))
209     return(0.); /* disk entirely on projection side */
210     /* should shrink disk but I'm lazy */
211     return(rad2);
212     }
213    
214    
215 greg 1.1 commonspot(sp1, sp2, org) /* set sp1 to intersection of sp1 and sp2 */
216     register SPOT *sp1, *sp2;
217     FVECT org;
218     {
219     FVECT cent;
220 greg 1.2 double rad2, cos1, cos2;
221 greg 1.1
222 greg 1.2 cos1 = 1. - sp1->siz/(2.*PI);
223     cos2 = 1. - sp2->siz/(2.*PI);
224 greg 1.1 if (sp2->siz >= 2.*PI-FTINY) /* BIG, just check overlap */
225 greg 1.2 return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
226     sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
227 greg 1.1 /* compute and check disks */
228 greg 1.2 rad2 = intercircle(cent, sp1->aim, sp2->aim,
229     1./(cos1*cos1) - 1., 1./(cos2*cos2) - 1.);
230 greg 1.1 if (rad2 <= FTINY || normalize(cent) == 0.)
231     return(0);
232     VCOPY(sp1->aim, cent);
233     sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
234     return(1);
235     }
236    
237    
238     commonbeam(sp1, sp2, dir) /* set sp1 to intersection of sp1 and sp2 */
239     register SPOT *sp1, *sp2;
240     FVECT dir;
241     {
242     FVECT cent, c1, c2;
243     double rad2, d;
244     register int i;
245     /* move centers to common plane */
246     d = DOT(sp1->aim, dir);
247     for (i = 0; i < 3; i++)
248 greg 1.2 c1[i] = sp1->aim[i] - d*dir[i];
249 greg 1.1 d = DOT(sp2->aim, dir);
250     for (i = 0; i < 3; i++)
251     c2[i] = sp2->aim[i] - d*dir[i];
252     /* compute overlap */
253     rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
254     if (rad2 <= FTINY)
255     return(0);
256     VCOPY(sp1->aim, cent);
257     sp1->siz = PI*rad2;
258     return(1);
259     }
260    
261    
262     checkspot(sp, nrm) /* check spotlight for behind source */
263     register SPOT *sp;
264     FVECT nrm;
265     {
266     double d, d1;
267    
268     d = DOT(sp->aim, nrm);
269     if (d > FTINY) /* center in front? */
270     return(0);
271     /* else check horizon */
272     d1 = 1. - sp->siz/(2.*PI);
273     return(1.-FTINY-d*d > d1*d1);
274     }
275    
276    
277     double
278     intercircle(cc, c1, c2, r1s, r2s) /* intersect two circles */
279     FVECT cc; /* midpoint (return value) */
280     FVECT c1, c2; /* circle centers */
281     double r1s, r2s; /* radii squared */
282     {
283     double a2, d2, l;
284     FVECT disp;
285     register int i;
286    
287     for (i = 0; i < 3; i++)
288     disp[i] = c2[i] - c1[i];
289     d2 = DOT(disp,disp);
290     /* circle within overlap? */
291     if (r1s < r2s) {
292     if (r2s >= r1s + d2) {
293     VCOPY(cc, c1);
294     return(r1s);
295     }
296     } else {
297     if (r1s >= r2s + d2) {
298     VCOPY(cc, c2);
299     return(r2s);
300     }
301     }
302     a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
303     /* no overlap? */
304     if (a2 <= 0.)
305     return(0.);
306 greg 1.2 /* overlap, compute center */
307 greg 1.1 l = sqrt((r1s - a2)/d2);
308     for (i = 0; i < 3; i++)
309     cc[i] = c1[i] + l*disp[i];
310     return(a2);
311     }
312 greg 1.4
313    
314     #ifdef DEBUG
315 greg 1.6 virtverb(sn, fp) /* print verbose description of virtual source */
316     register int sn;
317 greg 1.4 FILE *fp;
318     {
319     register int i;
320    
321     fprintf(fp, "%s virtual source %d in %s %s\n",
322 greg 1.6 source[sn].sflags & SDISTANT ? "distant" : "local",
323     sn, ofun[source[sn].so->otype].funame,
324     source[sn].so->oname);
325 greg 1.4 fprintf(fp, "\tat (%f,%f,%f)\n",
326 greg 1.6 source[sn].sloc[0], source[sn].sloc[1], source[sn].sloc[2]);
327 greg 1.4 fprintf(fp, "\tlinked to source %d (%s)\n",
328 greg 1.6 source[sn].sa.svnext, source[source[sn].sa.svnext].so->oname);
329     if (source[sn].sflags & SFOLLOW)
330 greg 1.4 fprintf(fp, "\talways followed\n");
331     else
332     fprintf(fp, "\tnever followed\n");
333 greg 1.6 if (!(source[sn].sflags & SSPOT))
334 greg 1.4 return;
335     fprintf(fp, "\twith spot aim (%f,%f,%f) and size %f\n",
336 greg 1.6 source[sn].sl.s->aim[0], source[sn].sl.s->aim[1],
337     source[sn].sl.s->aim[2], source[sn].sl.s->siz);
338 greg 1.4 }
339     #endif