ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcsamp.c
Revision: 2.11
Committed: Sat Sep 13 17:31:35 2003 UTC (20 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 2.10: +2 -2 lines
Log Message:
Minor fix in source sampling calculation

File Contents

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