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

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: srcsamp.c,v 2.18 2011/12/28 18:39:36 greg Exp $";
3 #endif
4 /*
5 * Source sampling routines
6 *
7 * External symbols declared in source.h
8 */
9
10 #include "copyright.h"
11
12 #include "ray.h"
13
14 #include "source.h"
15
16 #include "random.h"
17
18
19 static int
20 srcskip( /* pre-emptive test for out-of-range glow */
21 SRCREC *sp,
22 FVECT orig
23 )
24 {
25 if (sp->sflags & SSKIP)
26 return(1);
27
28 if ((sp->sflags & (SPROX|SDISTANT)) != SPROX)
29 return(0);
30
31 return(dist2(orig, sp->sloc) >
32 (sp->sl.prox + sp->srad)*(sp->sl.prox + sp->srad));
33 }
34
35 double
36 nextssamp( /* compute sample for source, rtn. distance */
37 RAY *r, /* origin is read, direction is set */
38 SRCINDEX *si /* source index (modified to current) */\
39 )
40 {
41 int cent[3], size[3], parr[2];
42 SRCREC *srcp;
43 double vpos[3];
44 double d;
45 int i;
46 nextsample:
47 while (++si->sp >= si->np) { /* get next sample */
48 if (++si->sn >= nsources)
49 return(0.0); /* no more */
50 if (srcskip(source+si->sn, r->rorg))
51 si->np = 0;
52 else if (srcsizerat <= FTINY)
53 nopart(si, r);
54 else {
55 for (i = si->sn; source[i].sflags & SVIRTUAL;
56 i = source[i].sa.sv.sn)
57 ; /* partition source */
58 (*sfun[source[i].so->otype].of->partit)(si, r);
59 }
60 si->sp = -1;
61 }
62 /* get partition */
63 cent[0] = cent[1] = cent[2] = 0;
64 size[0] = size[1] = size[2] = MAXSPART;
65 parr[0] = 0; parr[1] = si->sp;
66 if (!skipparts(cent, size, parr, si->spt))
67 error(CONSISTENCY, "bad source partition in nextssamp");
68 /* compute sample */
69 srcp = source + si->sn;
70 if (dstrsrc > FTINY) { /* jitter sample */
71 dimlist[ndims] = si->sn + 8831;
72 dimlist[ndims+1] = si->sp + 3109;
73 d = urand(ilhash(dimlist,ndims+2)+samplendx);
74 if (srcp->sflags & SFLAT) {
75 multisamp(vpos, 2, d);
76 vpos[SW] = 0.5;
77 } else
78 multisamp(vpos, 3, d);
79 for (i = 0; i < 3; i++)
80 vpos[i] = dstrsrc * (1. - 2.*vpos[i]) *
81 (double)size[i]*(1.0/MAXSPART);
82 } else
83 vpos[0] = vpos[1] = vpos[2] = 0.0;
84
85 VSUM(vpos, vpos, cent, 1.0/MAXSPART);
86 /* avoid circular aiming failures */
87 if ((srcp->sflags & SCIR) && (si->np > 1 || dstrsrc > 0.7)) {
88 FVECT trim;
89 if (srcp->sflags & (SFLAT|SDISTANT)) {
90 d = 1.12837917; /* correct setflatss() */
91 trim[SU] = d*sqrt(1.0 - 0.5*vpos[SV]*vpos[SV]);
92 trim[SV] = d*sqrt(1.0 - 0.5*vpos[SU]*vpos[SU]);
93 trim[SW] = 0.0;
94 } else {
95 trim[SW] = trim[SU] = vpos[SU]*vpos[SU];
96 d = vpos[SV]*vpos[SV];
97 if (d > trim[SW]) trim[SW] = d;
98 trim[SU] += d;
99 d = vpos[SW]*vpos[SW];
100 if (d > trim[SW]) trim[SW] = d;
101 trim[SU] += d;
102 if (trim[SU] > FTINY*FTINY) {
103 d = 1.0/0.7236; /* correct sphsetsrc() */
104 trim[SW] = trim[SV] = trim[SU] =
105 d*sqrt(trim[SW]/trim[SU]);
106 } else
107 trim[SW] = trim[SV] = trim[SU] = 0.0;
108 }
109 for (i = 0; i < 3; i++)
110 vpos[i] *= trim[i];
111 }
112 /* compute direction */
113 for (i = 0; i < 3; i++)
114 r->rdir[i] = srcp->sloc[i] +
115 vpos[SU]*srcp->ss[SU][i] +
116 vpos[SV]*srcp->ss[SV][i] +
117 vpos[SW]*srcp->ss[SW][i];
118
119 if (!(srcp->sflags & SDISTANT))
120 VSUB(r->rdir, r->rdir, r->rorg);
121 /* compute distance */
122 if ((d = normalize(r->rdir)) == 0.0)
123 goto nextsample; /* at source! */
124
125 /* compute sample size */
126 if (srcp->sflags & SFLAT) {
127 si->dom = sflatform(si->sn, r->rdir);
128 si->dom *= size[SU]*size[SV]*(1.0/MAXSPART/MAXSPART);
129 } else if (srcp->sflags & SCYL) {
130 si->dom = scylform(si->sn, r->rdir);
131 si->dom *= size[SU]*(1.0/MAXSPART);
132 } else {
133 si->dom = size[SU]*size[SV]*(double)size[SW] *
134 (1.0/MAXSPART/MAXSPART/MAXSPART) ;
135 }
136 if (srcp->sflags & SDISTANT) {
137 si->dom *= srcp->ss2;
138 return(FHUGE);
139 }
140 if (si->dom <= 1e-4)
141 goto nextsample; /* behind source? */
142 si->dom *= srcp->ss2/(d*d);
143 return(d); /* sample OK, return distance */
144 }
145
146
147 int
148 skipparts( /* skip to requested partition */
149 int ct[3],
150 int sz[3], /* center and size of partition (returned) */
151 int pp[2], /* current index, number to skip (modified) */
152 unsigned char *pt /* partition array */
153 )
154 {
155 int p;
156 /* check this partition */
157 p = spart(pt, pp[0]);
158 pp[0]++;
159 if (p == S0) { /* leaf partition */
160 if (pp[1]) {
161 pp[1]--;
162 return(0); /* not there yet */
163 } else
164 return(1); /* we've arrived */
165 }
166 /* else check lower */
167 sz[p] >>= 1;
168 ct[p] -= sz[p];
169 if (skipparts(ct, sz, pp, pt))
170 return(1); /* return hit */
171 /* else check upper */
172 ct[p] += sz[p] << 1;
173 if (skipparts(ct, sz, pp, pt))
174 return(1); /* return hit */
175 /* else return to starting position */
176 ct[p] -= sz[p];
177 sz[p] <<= 1;
178 return(0); /* return miss */
179 }
180
181
182 void
183 nopart( /* single source partition */
184 SRCINDEX *si,
185 RAY *r
186 )
187 {
188 clrpart(si->spt);
189 setpart(si->spt, 0, S0);
190 si->np = 1;
191 }
192
193
194 static int
195 cyl_partit( /* slice a cylinder */
196 FVECT ro,
197 unsigned char *pt,
198 int *pi,
199 int mp,
200 FVECT cent,
201 FVECT axis,
202 double d2
203 )
204 {
205 FVECT newct, newax;
206 int npl, npu;
207
208 if (mp < 2 || dist2(ro, cent) >= d2) { /* hit limit? */
209 setpart(pt, *pi, S0);
210 (*pi)++;
211 return(1);
212 }
213 /* subdivide */
214 setpart(pt, *pi, SU);
215 (*pi)++;
216 newax[0] = .5*axis[0];
217 newax[1] = .5*axis[1];
218 newax[2] = .5*axis[2];
219 d2 *= 0.25;
220 /* lower half */
221 newct[0] = cent[0] - newax[0];
222 newct[1] = cent[1] - newax[1];
223 newct[2] = cent[2] - newax[2];
224 npl = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
225 /* upper half */
226 newct[0] = cent[0] + newax[0];
227 newct[1] = cent[1] + newax[1];
228 newct[2] = cent[2] + newax[2];
229 npu = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
230 /* return total */
231 return(npl + npu);
232 }
233
234
235 void
236 cylpart( /* partition a cylinder */
237 SRCINDEX *si,
238 RAY *r
239 )
240 {
241 double dist2, safedist2, dist2cent, rad2;
242 FVECT v;
243 SRCREC *sp;
244 int pi;
245 /* first check point location */
246 clrpart(si->spt);
247 sp = source + si->sn;
248 rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
249 v[0] = r->rorg[0] - sp->sloc[0];
250 v[1] = r->rorg[1] - sp->sloc[1];
251 v[2] = r->rorg[2] - sp->sloc[2];
252 dist2 = DOT(v,sp->ss[SU]);
253 safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
254 dist2 *= dist2 / safedist2;
255 dist2cent = DOT(v,v);
256 dist2 = dist2cent - dist2;
257 if (dist2 <= rad2) { /* point inside extended cylinder */
258 si->np = 0;
259 return;
260 }
261 safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
262 if (dist2 <= 4.*rad2 || /* point too close to subdivide */
263 dist2cent >= safedist2) { /* or too far */
264 setpart(si->spt, 0, S0);
265 si->np = 1;
266 return;
267 }
268 pi = 0;
269 si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
270 sp->sloc, sp->ss[SU], safedist2);
271 }
272
273
274 static int
275 flt_partit( /* partition flatty */
276 FVECT ro,
277 unsigned char *pt,
278 int *pi,
279 int mp,
280 FVECT cent,
281 FVECT u,
282 FVECT v,
283 double du2,
284 double dv2
285 )
286 {
287 double d2;
288 FVECT newct, newax;
289 int npl, npu;
290
291 if (mp < 2 || ((d2 = dist2(ro, cent)) >= du2
292 && d2 >= dv2)) { /* hit limit? */
293 setpart(pt, *pi, S0);
294 (*pi)++;
295 return(1);
296 }
297 if (du2 > dv2) { /* subdivide in U */
298 setpart(pt, *pi, SU);
299 (*pi)++;
300 newax[0] = .5*u[0];
301 newax[1] = .5*u[1];
302 newax[2] = .5*u[2];
303 u = newax;
304 du2 *= 0.25;
305 } else { /* subdivide in V */
306 setpart(pt, *pi, SV);
307 (*pi)++;
308 newax[0] = .5*v[0];
309 newax[1] = .5*v[1];
310 newax[2] = .5*v[2];
311 v = newax;
312 dv2 *= 0.25;
313 }
314 /* lower half */
315 newct[0] = cent[0] - newax[0];
316 newct[1] = cent[1] - newax[1];
317 newct[2] = cent[2] - newax[2];
318 npl = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
319 /* upper half */
320 newct[0] = cent[0] + newax[0];
321 newct[1] = cent[1] + newax[1];
322 newct[2] = cent[2] + newax[2];
323 npu = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
324 /* return total */
325 return(npl + npu);
326 }
327
328
329 void
330 flatpart( /* partition a flat source */
331 SRCINDEX *si,
332 RAY *r
333 )
334 {
335 RREAL *vp;
336 FVECT v;
337 double du2, dv2;
338 int pi;
339
340 clrpart(si->spt);
341 vp = source[si->sn].sloc;
342 v[0] = r->rorg[0] - vp[0];
343 v[1] = r->rorg[1] - vp[1];
344 v[2] = r->rorg[2] - vp[2];
345 vp = source[si->sn].snorm;
346 if (DOT(v,vp) <= 0.) { /* behind source */
347 si->np = 0;
348 return;
349 }
350 dv2 = 2.*r->rweight/srcsizerat;
351 dv2 *= dv2;
352 vp = source[si->sn].ss[SU];
353 du2 = dv2 * DOT(vp,vp);
354 vp = source[si->sn].ss[SV];
355 dv2 *= DOT(vp,vp);
356 pi = 0;
357 si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
358 source[si->sn].sloc,
359 source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
360 }
361
362
363 double
364 scylform( /* compute cosine for cylinder's projection */
365 int sn,
366 FVECT dir /* assume normalized */
367 )
368 {
369 RREAL *dv;
370 double d;
371
372 dv = source[sn].ss[SU];
373 d = DOT(dir, dv);
374 d *= d / DOT(dv,dv);
375 return(sqrt(1. - d));
376 }