ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/virtuals.c
Revision: 2.17
Committed: Thu Sep 7 05:20:54 2006 UTC (17 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0, rad3R8, rad3R9
Changes since 2.16: +14 -5 lines
Log Message:
Fixed bug in virtual sources with intervening transmitters

File Contents

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