ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/virtuals.c
Revision: 2.21
Committed: Tue Jul 15 23:44:53 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad4R2, rad4R2P1
Changes since 2.20: +3 -4 lines
Log Message:
Fixed potential type conflict

File Contents

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