ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/srcsamp.c
(Generate patch)

Comparing ray/src/rt/srcsamp.c (file contents):
Revision 2.2 by greg, Mon Nov 25 09:52:49 1991 UTC vs.
Revision 2.25 by greg, Sat Jan 18 03:49:00 2025 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1991 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
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"
# Line 15 | Line 16 | static char SCCSid[] = "$SunId$ LBL";
16   #include  "random.h"
17  
18  
19 + int
20 + srcskip(                        /* pre-emptive test for source to skip */
21 +        int  sn,
22 +        RAY  *r
23 + )
24 + {
25 +        SRCREC  *sp = source + sn;
26 +
27 +        if (sp->sflags & SSKIP)
28 +                return(1);
29 +        if ((sp->sflags & (SPROX|SDISTANT)) == SPROX)
30 +                return(dist2(r->rorg, sp->sloc) >
31 +                        (sp->sl.prox + sp->srad)*(sp->sl.prox + sp->srad));
32 +        return(0);
33 + }
34 +
35   double
36 < nextssamp(r, si)                /* compute sample for source, rtn. distance */
37 < register RAY  *r;               /* origin is read, direction is set */
38 < register SRCINDEX  *si;         /* source index (modified to current) */
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 <        FVECT  vpos;
42 >        SRCREC  *srcp;
43 >        double  vpos[3];
44          double  d;
45 <        register int  i;
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 (source[si->sn].sflags & SSKIP)
50 >                if (srcskip(si->sn, r))
51                          si->np = 0;
52                  else if (srcsizerat <= FTINY)
53                          nopart(si, r);
# Line 47 | Line 66 | nextsample:
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 (source[si->sn].sflags & SFLAT) {
74 >                if (srcp->sflags & SFLAT) {
75                          multisamp(vpos, 2, d);
76 <                        vpos[2] = 0.5;
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]/MAXSPART;
81 >                                        (double)size[i]*(1.0/MAXSPART);
82          } else
83                  vpos[0] = vpos[1] = vpos[2] = 0.0;
84  
85 <        for (i = 0; i < 3; i++)
86 <                vpos[i] += (double)cent[i]/MAXSPART;
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] = source[si->sn].sloc[i] +
115 <                                vpos[SU]*source[si->sn].ss[SU][i] +
116 <                                vpos[SV]*source[si->sn].ss[SV][i] +
117 <                                vpos[SW]*source[si->sn].ss[SW][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 (!(source[si->sn].sflags & SDISTANT))
120 <                for (i = 0; i < 3; i++)
76 <                        r->rdir[i] -= r->rorg[i];
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 <        si->dom  = source[si->sn].ss2;
127 <        if (source[si->sn].sflags & SFLAT) {
128 <                si->dom *= sflatform(si->sn, r->rdir);
129 <                si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART);
130 <        } else if (source[si->sn].sflags & SCYL) {
131 <                si->dom *= scylform(si->sn, r->rdir);
88 <                si->dom *= (double)size[SU]/MAXSPART;
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 *= (double)(size[SU]*size[SV]*size[SW]) /
134 <                                (MAXSPART*MAXSPART*MAXSPART) ;
133 >                si->dom = size[SU]*size[SV]*(double)size[SW] *
134 >                                (1.0/MAXSPART/MAXSPART/MAXSPART) ;
135          }
136 <        if (source[si->sn].sflags & SDISTANT)
136 >        if (srcp->sflags & SDISTANT) {
137 >                si->dom *= srcp->ss2;
138                  return(FHUGE);
139 <        if (si->dom <= FTINY)
139 >        }
140 >        if (si->dom <= 1e-4)
141                  goto nextsample;                /* behind source? */
142 <        si->dom /= d*d;
142 >        si->dom *= srcp->ss2/(d*d);
143          return(d);              /* sample OK, return distance */
144   }
145  
146  
147 < skipparts(ct, sz, pp, pt)               /* skip to requested partition */
148 < int  ct[3], sz[3];              /* center and size of partition (returned) */
149 < register int  pp[2];            /* current index, number to skip (modified) */
150 < unsigned char  *pt;             /* partition array */
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 <        register int  p;
155 >        int  p;
156                                          /* check this partition */
157          p = spart(pt, pp[0]);
158          pp[0]++;
159 <        if (p == S0)                    /* leaf partition */
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];
# Line 130 | Line 179 | unsigned char  *pt;            /* partition array */
179   }
180  
181  
182 < nopart(si, r)                   /* single source partition */
183 < register SRCINDEX  *si;
184 < RAY  *r;
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);
# Line 140 | Line 191 | RAY  *r;
191   }
192  
193  
143 – cylpart(si, r)                  /* partition a cylinder */
144 – SRCINDEX  *si;
145 – register RAY  *r;
146 – {
147 –        double  dist2, safedist2, dist2cent, rad2;
148 –        FVECT  v;
149 –        register SRCREC  *sp;
150 –        int  pi;
151 –                                        /* first check point location */
152 –        clrpart(si->spt);
153 –        sp = source + si->sn;
154 –        rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
155 –        v[0] = r->rorg[0] - sp->sloc[0];
156 –        v[1] = r->rorg[1] - sp->sloc[1];
157 –        v[2] = r->rorg[2] - sp->sloc[2];
158 –        dist2 = DOT(v,sp->ss[SU]);
159 –        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
160 –        dist2 *= dist2 / safedist2;
161 –        dist2cent = DOT(v,v);
162 –        dist2 = dist2cent - dist2;
163 –        if (dist2 <= rad2) {            /* point inside extended cylinder */
164 –                si->np = 0;
165 –                return;
166 –        }
167 –        safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
168 –        if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
169 –                        dist2cent >= safedist2) {       /* or too far */
170 –                setpart(si->spt, 0, S0);
171 –                si->np = 1;
172 –                return;
173 –        }
174 –        pi = 0;
175 –        si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
176 –                        sp->sloc, sp->ss[SU], safedist2);
177 – }
178 –
179 –
194   static int
195 < cyl_partit(ro, pt, pi, mp, cent, axis, d2)      /* slice a cylinder */
196 < FVECT  ro;
197 < unsigned char  *pt;
198 < register int  *pi;
199 < int  mp;
200 < FVECT  cent, axis;
201 < double  d2;
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;
# Line 216 | Line 232 | double  d2;
232   }
233  
234  
235 < flatpart(si, r)                         /* partition a flat source */
236 < register SRCINDEX  *si;
237 < register RAY  *r;
235 > void
236 > cylpart(                        /* partition a cylinder */
237 >        SRCINDEX  *si,
238 >        RAY  *r
239 > )
240   {
241 <        register FLOAT  *vp;
241 >        double  dist2, safedist2, dist2cent, rad2;
242          FVECT  v;
243 <        double  du2, dv2;
243 >        SRCREC  *sp;
244          int  pi;
245 <
245 >                                        /* first check point location */
246          clrpart(si->spt);
247 <        vp = source[si->sn].sloc;
248 <        v[0] = r->rorg[0] - vp[0];
249 <        v[1] = r->rorg[1] - vp[1];
250 <        v[2] = r->rorg[2] - vp[2];
251 <        vp = source[si->sn].snorm;
252 <        if (DOT(v,vp) <= FTINY) {       /* behind source */
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 <        dv2 = 2.*r->rweight/srcsizerat;
262 <        dv2 *= dv2;
263 <        vp = source[si->sn].ss[SU];
264 <        du2 = dv2 * DOT(vp,vp);
265 <        vp = source[si->sn].ss[SV];
266 <        dv2 *= DOT(vp,vp);
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 = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
270 <                source[si->sn].sloc,
247 <                source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
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(ro, pt, pi, mp, cent, u, v, du2, dv2)        /* partition flatty */
276 < FVECT  ro;
277 < unsigned char  *pt;
278 < register int  *pi;
279 < int  mp;
280 < FVECT  cent, u, v;
281 < double  du2, dv2;
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;
# Line 299 | Line 326 | double  du2, dv2;
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(sn, dir)               /* compute cosine for cylinder's projection */
365 < int  sn;
366 < register FVECT  dir;            /* assume normalized */
364 > scylform(                       /* compute cosine for cylinder's projection */
365 >        int  sn,
366 >        FVECT  dir              /* assume normalized */
367 > )
368   {
369 <        register FLOAT  *dv;
369 >        RREAL  *dv;
370          double  d;
371  
372          dv = source[sn].ss[SU];

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines