ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/virtuals.c
Revision: 1.20
Committed: Fri Aug 2 10:30:27 1991 UTC (32 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.19: +6 -2 lines
Log Message:
made computation of number of sample rays more accurate
and fixed calculation of solid angle for virtuals sources spots

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 * Routines for simulating virtual light sources
9 * Thus far, we only support planar mirrors.
10 */
11
12 #include "ray.h"
13
14 #include "octree.h"
15
16 #include "otypes.h"
17
18 #include "source.h"
19
20 #include "random.h"
21
22 #define MINSAMPLES 5 /* minimum number of pretest samples */
23 #define STESTMAX 30 /* maximum seeks per sample */
24
25
26 double getdisk();
27
28 static OBJECT *vobject; /* virtual source objects */
29 static int nvobjects = 0; /* number of virtual source objects */
30
31
32 markvirtuals() /* find and mark virtual sources */
33 {
34 register OBJREC *o;
35 register int i;
36 /* check number of direct relays */
37 if (directrelay <= 0)
38 return;
39 /* find virtual source objects */
40 for (i = 0; i < nobjects; i++) {
41 o = objptr(i);
42 if (!issurface(o->otype) || o->omod == OVOID)
43 continue;
44 if (!isvlight(objptr(o->omod)->otype))
45 continue;
46 if (sfun[o->otype].of == NULL ||
47 sfun[o->otype].of->getpleq == NULL)
48 objerror(o, USER, "illegal material");
49 if (nvobjects == 0)
50 vobject = (OBJECT *)malloc(sizeof(OBJECT));
51 else
52 vobject = (OBJECT *)realloc((char *)vobject,
53 (unsigned)(nvobjects+1)*sizeof(OBJECT));
54 if (vobject == NULL)
55 error(SYSTEM, "out of memory in addvirtuals");
56 vobject[nvobjects++] = i;
57 }
58 if (nvobjects == 0)
59 return;
60 #ifdef DEBUG
61 fprintf(stderr, "found %d virtual source objects\n", nvobjects);
62 #endif
63 /* append virtual sources */
64 for (i = nsources; i-- > 0; )
65 addvirtuals(i, directrelay);
66 /* done with our object list */
67 free((char *)vobject);
68 nvobjects = 0;
69 }
70
71
72 addvirtuals(sn, nr) /* add virtuals associated with source */
73 int sn;
74 int nr;
75 {
76 register int i;
77 /* check relay limit first */
78 if (nr <= 0)
79 return;
80 if (source[sn].sflags & SSKIP)
81 return;
82 /* check each virtual object for projection */
83 for (i = 0; i < nvobjects; i++)
84 /* vproject() calls us recursively */
85 vproject(objptr(vobject[i]), sn, nr-1);
86 }
87
88
89 vproject(o, sn, n) /* create projected source(s) if they exist */
90 OBJREC *o;
91 int sn;
92 int n;
93 {
94 register int i;
95 register VSMATERIAL *vsmat;
96 MAT4 proj;
97 int ns;
98
99 if (o == source[sn].so) /* objects cannot project themselves */
100 return;
101 /* get virtual source material */
102 vsmat = sfun[objptr(o->omod)->otype].mf;
103 /* project virtual sources */
104 for (i = 0; i < vsmat->nproj; i++)
105 if ((*vsmat->vproj)(proj, o, &source[sn], i))
106 if ((ns = makevsrc(o, sn, proj)) >= 0) {
107 source[ns].sa.sv.pn = i;
108 #ifdef DEBUG
109 virtverb(ns, stderr);
110 #endif
111 addvirtuals(ns, n);
112 }
113 }
114
115
116 int
117 makevsrc(op, sn, pm) /* make virtual source if reasonable */
118 OBJREC *op;
119 register int sn;
120 MAT4 pm;
121 {
122 FVECT nsloc, nsnorm, ocent, v;
123 double maxrad2, d;
124 int nsflags;
125 SPOT theirspot, ourspot;
126 register int i;
127
128 nsflags = source[sn].sflags | (SVIRTUAL|SSPOT|SFOLLOW);
129 /* get object center and max. radius */
130 maxrad2 = getdisk(ocent, op, sn);
131 if (maxrad2 <= FTINY) /* too small? */
132 return(-1);
133 /* get location and spot */
134 if (source[sn].sflags & SDISTANT) { /* distant source */
135 if (source[sn].sflags & SPROX)
136 return(-1); /* should never get here! */
137 multv3(nsloc, source[sn].sloc, pm);
138 normalize(nsloc);
139 VCOPY(ourspot.aim, ocent);
140 ourspot.siz = PI*maxrad2;
141 ourspot.flen = 0.;
142 if (source[sn].sflags & SSPOT) {
143 multp3(theirspot.aim, source[sn].sl.s->aim, pm);
144 d = sqrt(dist2(ourspot.aim, theirspot.aim));
145 d = sqrt(source[sn].sl.s->siz/PI) + d*source[sn].ss;
146 theirspot.siz = PI*d*d;
147 ourspot.flen = theirspot.flen = source[sn].sl.s->flen;
148 d = ourspot.siz;
149 if (!commonbeam(&ourspot, &theirspot, nsloc))
150 return(-1); /* no overlap */
151 if (ourspot.siz < d-FTINY) { /* it shrunk */
152 d = beamdisk(v, op, &ourspot, nsloc);
153 if (d <= FTINY)
154 return(-1);
155 if (d < maxrad2) {
156 maxrad2 = d;
157 VCOPY(ocent, v);
158 }
159 }
160 }
161 } else { /* local source */
162 multp3(nsloc, source[sn].sloc, pm);
163 for (i = 0; i < 3; i++)
164 ourspot.aim[i] = ocent[i] - nsloc[i];
165 if ((d = normalize(ourspot.aim)) == 0.)
166 return(-1); /* at source!! */
167 if (source[sn].sflags & SPROX && d > source[sn].sl.prox)
168 return(-1); /* too far away */
169 ourspot.flen = 0.;
170 d = (sqrt(maxrad2) + source[sn].ss) / d;
171 if (d < 1.-FTINY)
172 ourspot.siz = 2.*PI*(1. - sqrt(1.-d*d));
173 else
174 nsflags &= ~SSPOT;
175 if (source[sn].sflags & SSPOT) {
176 copystruct(&theirspot, source[sn].sl.s);
177 multv3(theirspot.aim, source[sn].sl.s->aim, pm);
178 normalize(theirspot.aim);
179 if (nsflags & SSPOT) {
180 ourspot.flen = theirspot.flen;
181 d = ourspot.siz;
182 if (!commonspot(&ourspot, &theirspot, nsloc))
183 return(-1); /* no overlap */
184 } else {
185 nsflags |= SSPOT;
186 copystruct(&ourspot, &theirspot);
187 d = 2.*ourspot.siz;
188 }
189 if (ourspot.siz < d-FTINY) { /* it shrunk */
190 d = spotdisk(v, op, &ourspot, nsloc);
191 if (d <= FTINY)
192 return(-1);
193 if (d < maxrad2) {
194 maxrad2 = d;
195 VCOPY(ocent, v);
196 }
197 }
198 }
199 if (source[sn].sflags & SFLAT) { /* behind source? */
200 multv3(nsnorm, source[sn].snorm, pm);
201 normalize(nsnorm);
202 if (nsflags & SSPOT && !checkspot(&ourspot, nsnorm))
203 return(-1);
204 }
205 }
206 /* pretest visibility */
207 nsflags = vstestvis(nsflags, op, ocent, maxrad2, sn);
208 if (nsflags & SSKIP)
209 return(-1); /* obstructed */
210 /* it all checks out, so make it */
211 if ((i = newsource()) < 0)
212 goto memerr;
213 source[i].sflags = nsflags;
214 VCOPY(source[i].sloc, nsloc);
215 if (nsflags & SFLAT)
216 VCOPY(source[i].snorm, nsnorm);
217 source[i].ss = source[sn].ss; source[i].ss2 = source[sn].ss2;
218 if (nsflags & SSPOT) {
219 if ((source[i].sl.s = (SPOT *)malloc(sizeof(SPOT))) == NULL)
220 goto memerr;
221 copystruct(source[i].sl.s, &ourspot);
222 }
223 if (nsflags & SPROX)
224 source[i].sl.prox = source[sn].sl.prox;
225 source[i].sa.sv.sn = sn;
226 source[i].so = op;
227 return(i);
228 memerr:
229 error(SYSTEM, "out of memory in makevsrc");
230 }
231
232
233 double
234 getdisk(oc, op, sn) /* get visible object disk */
235 FVECT oc;
236 OBJREC *op;
237 register int sn;
238 {
239 double rad2, roffs, offs, d, rd, rdoto;
240 FVECT rnrm, nrm;
241 /* first, use object getdisk function */
242 rad2 = getmaxdisk(oc, op);
243 if (!(source[sn].sflags & SVIRTUAL))
244 return(rad2); /* all done for normal source */
245 /* check for correct side of relay surface */
246 roffs = getplaneq(rnrm, source[sn].so);
247 rd = DOT(rnrm, source[sn].sloc); /* source projection */
248 if (!(source[sn].sflags & SDISTANT))
249 rd -= roffs;
250 d = DOT(rnrm, oc) - roffs; /* disk distance to relay plane */
251 if ((d > 0.) ^ (rd > 0.))
252 return(rad2); /* OK if opposite sides */
253 if (d*d >= rad2)
254 return(0.); /* no relay is possible */
255 /* we need a closer look */
256 offs = getplaneq(nrm, op);
257 rdoto = DOT(rnrm, nrm);
258 if (d*d >= rad2*(1.-rdoto*rdoto))
259 return(0.); /* disk entirely on projection side */
260 /* should shrink disk but I'm lazy */
261 return(rad2);
262 }
263
264
265 int
266 vstestvis(f, o, oc, or2, sn) /* pretest source visibility */
267 int f; /* virtual source flags */
268 OBJREC *o; /* relay object */
269 FVECT oc; /* relay object center */
270 double or2; /* relay object radius squared */
271 register int sn; /* target source number */
272 {
273 RAY sr;
274 FVECT onorm;
275 FVECT offsdir;
276 double or, d;
277 int infront;
278 int stestlim, ssn;
279 int nhit, nok;
280 register int i, n;
281 /* return if pretesting disabled */
282 if (vspretest <= 0)
283 return(f);
284 /* get surface normal */
285 getplaneq(onorm, o);
286 /* set number of rays to sample */
287 if (source[sn].sflags & SDISTANT) {
288 n = (2./3.*PI*PI)*or2/(thescene.cusize*thescene.cusize)*
289 vspretest + .5;
290 infront = DOT(onorm, source[sn].sloc) > 0.;
291 } else {
292 for (i = 0; i < 3; i++)
293 offsdir[i] = source[sn].sloc[i] - oc[i];
294 d = DOT(offsdir,offsdir);
295 if (d <= FTINY)
296 n = 2.*PI * vspretest + .5;
297 else
298 n = 2.*PI * (1.-sqrt(1./(1.+or2/d)))*vspretest + .5;
299 infront = DOT(onorm, offsdir) > 0.;
300 }
301 if (n < MINSAMPLES) n = MINSAMPLES;
302 #ifdef DEBUG
303 fprintf(stderr, "pretesting source %d in object %s with %d rays\n",
304 sn, o->oname, n);
305 #endif
306 /* sample */
307 or = sqrt(or2);
308 stestlim = n*STESTMAX;
309 ssn = 0;
310 nhit = nok = 0;
311 while (n-- > 0) {
312 /* get sample point */
313 do {
314 if (ssn >= stestlim) {
315 #ifdef DEBUG
316 fprintf(stderr, "\ttoo hard to hit\n");
317 #endif
318 return(f); /* too small a target! */
319 }
320 for (i = 0; i < 3; i++)
321 offsdir[i] = or*(1. -
322 2.*urand(urind(931*i+5827,ssn)));
323 ssn++;
324 for (i = 0; i < 3; i++)
325 sr.rorg[i] = oc[i] + offsdir[i];
326 d = DOT(offsdir,onorm);
327 if (infront)
328 for (i = 0; i < 3; i++) {
329 sr.rorg[i] -= (d-.0001)*onorm[i];
330 sr.rdir[i] = -onorm[i];
331 }
332 else
333 for (i = 0; i < 3; i++) {
334 sr.rorg[i] -= (d+.0001)*onorm[i];
335 sr.rdir[i] = onorm[i];
336 }
337 rayorigin(&sr, NULL, PRIMARY, 1.0);
338 } while (!(*ofun[o->otype].funp)(o, &sr));
339 /* check against source */
340 samplendx++;
341 if (srcray(&sr, NULL, sn) == 0.)
342 continue;
343 sr.revf = srcvalue;
344 rayvalue(&sr);
345 if (bright(sr.rcol) <= FTINY)
346 continue;
347 nok++;
348 /* check against obstructions */
349 rayclear(&sr);
350 sr.revf = raytrace;
351 rayvalue(&sr);
352 if (bright(sr.rcol) > FTINY)
353 nhit++;
354 if (nhit > 0 && nhit < nok) {
355 #ifdef DEBUG
356 fprintf(stderr, "\tpartially occluded\n");
357 #endif
358 return(f); /* need to shadow test */
359 }
360 }
361 if (nhit == 0) {
362 #ifdef DEBUG
363 fprintf(stderr, "\t0%% hit rate\n");
364 #endif
365 return(f | SSKIP); /* 0% hit rate: totally occluded */
366 }
367 #ifdef DEBUG
368 fprintf(stderr, "\t100%% hit rate\n");
369 #endif
370 return(f & ~SFOLLOW); /* 100% hit rate: no occlusion */
371 }
372
373
374 #ifdef DEBUG
375 virtverb(sn, fp) /* print verbose description of virtual source */
376 register int sn;
377 FILE *fp;
378 {
379 register int i;
380
381 fprintf(fp, "%s virtual source %d in %s %s\n",
382 source[sn].sflags & SDISTANT ? "distant" : "local",
383 sn, ofun[source[sn].so->otype].funame,
384 source[sn].so->oname);
385 fprintf(fp, "\tat (%f,%f,%f)\n",
386 source[sn].sloc[0], source[sn].sloc[1], source[sn].sloc[2]);
387 fprintf(fp, "\tlinked to source %d (%s)\n",
388 source[sn].sa.sv.sn, source[source[sn].sa.sv.sn].so->oname);
389 if (source[sn].sflags & SFOLLOW)
390 fprintf(fp, "\talways followed\n");
391 else
392 fprintf(fp, "\tnever followed\n");
393 if (!(source[sn].sflags & SSPOT))
394 return;
395 fprintf(fp, "\twith spot aim (%f,%f,%f) and size %f\n",
396 source[sn].sl.s->aim[0], source[sn].sl.s->aim[1],
397 source[sn].sl.s->aim[2], source[sn].sl.s->siz);
398 }
399 #endif