ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/virtuals.c
Revision: 1.29
Committed: Tue Oct 22 11:23:51 1991 UTC (32 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.28: +3 -0 lines
Log Message:
bug fixes in source sampling routines

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