ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/virtuals.c
Revision: 1.3
Committed: Thu Jun 20 13:43:38 1991 UTC (32 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.2: +80 -216 lines
Log Message:
major reorganization using dispatch table for source types

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 "otypes.h"
15
16 #include "source.h"
17
18 #include "cone.h"
19
20 #include "face.h"
21
22
23 double intercircle();
24 SRCREC *makevsrc();
25
26 static OBJECT *vobject; /* virtual source objects */
27 static int nvobjects = 0; /* number of virtual source objects */
28
29
30 markvirtuals() /* 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 < nobjects; i++) {
39 o = objptr(i);
40 if (!issurface(o->otype) || o->omod == OVOID)
41 continue;
42 if (!isvlight(objptr(o->omod)->otype))
43 continue;
44 if (sfun[o->otype].of == NULL ||
45 sfun[o->otype].of->getpleq == NULL)
46 objerror(o, USER, "illegal material");
47 if (nvobjects == 0)
48 vobject = (OBJECT *)malloc(sizeof(OBJECT));
49 else
50 vobject = (OBJECT *)realloc((char *)vobject,
51 (unsigned)(nvobjects+1)*sizeof(OBJECT));
52 if (vobject == NULL)
53 error(SYSTEM, "out of memory in addvirtuals");
54 vobject[nvobjects++] = i;
55 }
56 if (nvobjects == 0)
57 return;
58 /* append virtual sources */
59 for (i = nsources; i-- > 0; )
60 if (!(source[i].sflags & SSKIP))
61 addvirtuals(&source[i], directrelay);
62 /* done with our object list */
63 free((char *)vobject);
64 nvobjects = 0;
65 }
66
67
68 addvirtuals(sr, nr) /* add virtual sources associated with sr */
69 SRCREC *sr;
70 int nr;
71 {
72 register int i;
73 /* check relay limit first */
74 if (nr <= 0)
75 return;
76 /* check each virtual object for projection */
77 for (i = 0; i < nvobjects; i++)
78 /* vproject() calls us recursively */
79 vproject(objptr(i), sr, nr-1);
80 }
81
82
83 vproject(o, s, n) /* create projected source(s) if they exist */
84 OBJREC *o;
85 SRCREC *s;
86 int n;
87 {
88 register int i;
89 register VSMATERIAL *vsmat;
90 MAT4 proj;
91 SRCREC *ns;
92 /* get virtual source material */
93 vsmat = sfun[objptr(o->omod)->otype].mf;
94 /* project virtual sources */
95 for (i = 0; i < vsmat->nproj; i++)
96 if ((*vsmat->vproj)(proj, o, s, i))
97 if ((ns = makevsrc(o, s, proj)) != NULL)
98 addvirtuals(ns, n);
99 }
100
101
102 SRCREC *
103 makevsrc(op, sp, pm) /* make virtual source if reasonable */
104 OBJREC *op;
105 register SRCREC *sp;
106 MAT4 pm;
107 {
108 register SRCREC *newsrc;
109 FVECT nsloc, ocent, nsnorm;
110 int nsflags;
111 double maxrad2;
112 double d1;
113 SPOT theirspot, ourspot;
114 register int i;
115
116 nsflags = (sp->sflags|(SVIRTUAL|SFOLLOW)) & ~SSPOT;
117 /* get object center and max. radius */
118 if (sfun[op->otype].of->getdisk != NULL) {
119 maxrad2 = (*sfun[op->otype].of->getdisk)(ocent, op);
120 if (maxrad2 <= FTINY) /* too small? */
121 return(NULL);
122 nsflags |= SSPOT;
123 }
124 /* get location and spot */
125 if (sp->sflags & SDISTANT) { /* distant source */
126 if (sp->sflags & SPROX)
127 return(NULL); /* should never get here! */
128 multv3(nsloc, sp->sloc, pm);
129 if (nsflags & SSPOT) {
130 VCOPY(ourspot.aim, ocent);
131 ourspot.siz = PI*maxrad2;
132 ourspot.flen = 0.;
133 }
134 if (sp->sflags & SSPOT) {
135 copystruct(&theirspot, sp->sl.s);
136 multp3(theirspot.aim, sp->sl.s->aim, pm);
137 if (nsflags & SSPOT &&
138 !commonbeam(&ourspot, &theirspot, nsloc))
139 return(NULL); /* no overlap */
140 }
141 } else { /* local source */
142 multp3(nsloc, sp->sloc, pm);
143 if (nsflags & SSPOT) {
144 for (i = 0; i < 3; i++)
145 ourspot.aim[i] = ocent[i] - nsloc[i];
146 if ((d1 = normalize(ourspot.aim)) == 0.)
147 return(NULL); /* at source!! */
148 if (sp->sflags & SPROX && d1 > sp->sl.prox)
149 return(NULL); /* too far away */
150 ourspot.siz = 2.*PI*(1. - d1/sqrt(d1*d1+maxrad2));
151 ourspot.flen = 0.;
152 } else if (sp->sflags & SPROX) {
153 FVECT norm;
154 double offs;
155 /* use distance from plane */
156 offs = (*sfun[op->otype].of->getpleq)(norm, op);
157 d1 = DOT(norm, nsloc) - offs;
158 if (d1 > sp->sl.prox || d1 < -sp->sl.prox)
159 return(NULL); /* too far away */
160 }
161 if (sp->sflags & SSPOT) {
162 copystruct(&theirspot, sp->sl.s);
163 multv3(theirspot.aim, sp->sl.s->aim, pm);
164 if (nsflags & SSPOT) {
165 if (!commonspot(&ourspot, &theirspot, nsloc))
166 return(NULL); /* no overlap */
167 ourspot.flen = theirspot.flen;
168 }
169 }
170 if (sp->sflags & SFLAT) { /* check for behind source */
171 multv3(nsnorm, sp->snorm, pm);
172 if (nsflags & SSPOT && checkspot(&ourspot, nsnorm) < 0)
173 return(NULL);
174 }
175 }
176 /* everything is OK, make source */
177 if ((newsrc = newsource()) == NULL)
178 goto memerr;
179 newsrc->sflags = nsflags;
180 VCOPY(newsrc->sloc, nsloc);
181 if (nsflags & SFLAT)
182 VCOPY(newsrc->snorm, nsnorm);
183 newsrc->ss = sp->ss; newsrc->ss2 = sp->ss2;
184 if ((nsflags | sp->sflags) & SSPOT) {
185 if ((newsrc->sl.s = (SPOT *)malloc(sizeof(SPOT))) == NULL)
186 goto memerr;
187 if (nsflags & SSPOT)
188 copystruct(newsrc->sl.s, &ourspot);
189 else
190 copystruct(newsrc->sl.s, &theirspot);
191 newsrc->sflags |= SSPOT;
192 }
193 if (nsflags & SPROX)
194 newsrc->sl.prox = sp->sl.prox;
195 newsrc->sa.svnext = sp - source;
196 return(newsrc);
197 memerr:
198 error(SYSTEM, "out of memory in makevsrc");
199 }
200
201
202 commonspot(sp1, sp2, org) /* set sp1 to intersection of sp1 and sp2 */
203 register SPOT *sp1, *sp2;
204 FVECT org;
205 {
206 FVECT cent;
207 double rad2, cos1, cos2;
208
209 cos1 = 1. - sp1->siz/(2.*PI);
210 cos2 = 1. - sp2->siz/(2.*PI);
211 if (sp2->siz >= 2.*PI-FTINY) /* BIG, just check overlap */
212 return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
213 sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
214 /* compute and check disks */
215 rad2 = intercircle(cent, sp1->aim, sp2->aim,
216 1./(cos1*cos1) - 1., 1./(cos2*cos2) - 1.);
217 if (rad2 <= FTINY || normalize(cent) == 0.)
218 return(0);
219 VCOPY(sp1->aim, cent);
220 sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
221 return(1);
222 }
223
224
225 commonbeam(sp1, sp2, dir) /* set sp1 to intersection of sp1 and sp2 */
226 register SPOT *sp1, *sp2;
227 FVECT dir;
228 {
229 FVECT cent, c1, c2;
230 double rad2, d;
231 register int i;
232 /* move centers to common plane */
233 d = DOT(sp1->aim, dir);
234 for (i = 0; i < 3; i++)
235 c1[i] = sp1->aim[i] - d*dir[i];
236 d = DOT(sp2->aim, dir);
237 for (i = 0; i < 3; i++)
238 c2[i] = sp2->aim[i] - d*dir[i];
239 /* compute overlap */
240 rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
241 if (rad2 <= FTINY)
242 return(0);
243 VCOPY(sp1->aim, cent);
244 sp1->siz = PI*rad2;
245 return(1);
246 }
247
248
249 checkspot(sp, nrm) /* check spotlight for behind source */
250 register SPOT *sp;
251 FVECT nrm;
252 {
253 double d, d1;
254
255 d = DOT(sp->aim, nrm);
256 if (d > FTINY) /* center in front? */
257 return(0);
258 /* else check horizon */
259 d1 = 1. - sp->siz/(2.*PI);
260 return(1.-FTINY-d*d > d1*d1);
261 }
262
263
264 double
265 intercircle(cc, c1, c2, r1s, r2s) /* intersect two circles */
266 FVECT cc; /* midpoint (return value) */
267 FVECT c1, c2; /* circle centers */
268 double r1s, r2s; /* radii squared */
269 {
270 double a2, d2, l;
271 FVECT disp;
272 register int i;
273
274 for (i = 0; i < 3; i++)
275 disp[i] = c2[i] - c1[i];
276 d2 = DOT(disp,disp);
277 /* circle within overlap? */
278 if (r1s < r2s) {
279 if (r2s >= r1s + d2) {
280 VCOPY(cc, c1);
281 return(r1s);
282 }
283 } else {
284 if (r1s >= r2s + d2) {
285 VCOPY(cc, c2);
286 return(r2s);
287 }
288 }
289 a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
290 /* no overlap? */
291 if (a2 <= 0.)
292 return(0.);
293 /* overlap, compute center */
294 l = sqrt((r1s - a2)/d2);
295 for (i = 0; i < 3; i++)
296 cc[i] = c1[i] + l*disp[i];
297 return(a2);
298 }