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 1.7 by greg, Mon Oct 28 08:07:44 1991 UTC vs.
Revision 2.24 by greg, Wed Dec 25 17:40:27 2024 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"
15  
16   #include  "random.h"
17  
18 + #ifdef SSKIPOPT
19 + /* The following table is used for skipping sources */
20 + static uby8     *srcskipflags = NULL;           /* source inclusion lookup */
21 + static int      ssf_count = 0;                  /* number of flag entries */
22 + static int      ssf_max = 0;                    /* current array size */
23 + static uby8     *ssf_noskip = NULL;             /* set of zero flags */
24  
25 + uby8            *ssf_select = NULL;             /* sources we may skip */
26 +
27 + /* Find/allocate source skip flag entry (free all if NULL) */
28 + int
29 + sskip_rsi(uby8 *flags)
30 + {
31 +        uby8    *flp;
32 +        int     i;
33 +
34 +        if (flags == NULL) {            /* means clear all */
35 +                efree(srcskipflags); srcskipflags = NULL;
36 +                ssf_count = ssf_max = 0;
37 +                sskip_free(ssf_noskip);
38 +                sskip_free(ssf_select);
39 +                return(0);
40 +        }
41 +        if (ssf_noskip == NULL)         /* first call? */
42 +                ssf_noskip = sskip_new();
43 +
44 +        if (sskip_eq(flags, ssf_noskip))
45 +                return(-1);             /* nothing to skip */
46 +                                        /* search recent entries */
47 +        flp = srcskipflags + ssf_count*SSKIPFLSIZ;
48 +        for (i = ssf_count; i-- > 0; )
49 +                if (sskip_eq(flp -= SSKIPFLSIZ, flags))
50 +                        return(-2-i);   /* found it! */
51 +                                        /* else tack on new entry */
52 +        if (ssf_count >= ssf_max) {     /* need more space? */
53 +                ssf_max = ssf_count + (ssf_count>>2) + 64;
54 +                if (ssf_max <= ssf_count &&
55 +                                (ssf_max = ssf_count+1024) <= ssf_count)
56 +                        error(SYSTEM, "out of space in sskip_rsi()");
57 +
58 +                srcskipflags = (uby8 *)erealloc(srcskipflags,
59 +                                                ssf_max*SSKIPFLSIZ);
60 +        }
61 +        sskip_cpy(srcskipflags + ssf_count*SSKIPFLSIZ, flags);
62 +
63 +        return(-2 - ssf_count++);       /* return index (< -1) */
64 + }
65 +
66 + /* Get skip flags associated with RAY rsrc index (or NULL) */
67 + uby8 *
68 + sskip_flags(int rsi)
69 + {
70 +        if (rsi >= -1)
71 +                return(ssf_noskip);
72 +
73 +        if ((rsi = -2 - rsi) >= ssf_count)
74 +                error(CONSISTENCY, "bad index to sskip_flags()");
75 +
76 +        return(srcskipflags + rsi*SSKIPFLSIZ);
77 + }
78 +
79 + /* OR in a second set of flags into a first */
80 + void
81 + sskip_addflags(uby8 *dfl, const uby8 *sfl)
82 + {
83 +        int     nb = SSKIPFLSIZ;
84 +
85 +        while (nb--)
86 +                *dfl++ |= *sfl++;
87 + }
88 + #endif
89 +
90 + int
91 + srcskip(                        /* pre-emptive test for source to skip */
92 +        int  sn,
93 +        RAY  *r
94 + )
95 + {
96 +        SRCREC  *sp = source + sn;
97 +
98 +        if (sp->sflags & SSKIP)
99 +                return(1);
100 + #ifdef SSKIPOPT                 /* parent ray has custom skip flags? */
101 +        if (r->parent != NULL && r->parent->rsrc < -1 &&
102 +                        sskip_chk(sskip_flags(r->parent->rsrc), sn))
103 +                return(1);
104 + #endif
105 +        if ((sp->sflags & (SPROX|SDISTANT)) == SPROX)
106 +                return(dist2(r->rorg, sp->sloc) >
107 +                        (sp->sl.prox + sp->srad)*(sp->sl.prox + sp->srad));
108 +        return(0);
109 + }
110 +
111   double
112 < nextssamp(r, si)                /* compute sample for source, rtn. distance */
113 < register RAY  *r;               /* origin is read, direction is set */
114 < register SRCINDEX  *si;         /* source index (modified to current) */
112 > nextssamp(                      /* compute sample for source, rtn. distance */
113 >        RAY  *r,                /* origin is read, direction is set */
114 >        SRCINDEX  *si           /* source index (modified to current) */
115 > )
116   {
117          int  cent[3], size[3], parr[2];
118 <        FVECT  vpos;
118 >        SRCREC  *srcp;
119 >        double  vpos[3];
120          double  d;
121 <        register int  i;
122 <
121 >        int  i;
122 > nextsample:
123          while (++si->sp >= si->np) {    /* get next sample */
124                  if (++si->sn >= nsources)
125                          return(0.0);    /* no more */
126 <                if (source[si->sn].sflags & SSKIP)
126 >                if (srcskip(si->sn, r))
127                          si->np = 0;
128                  else if (srcsizerat <= FTINY)
129                          nopart(si, r);
# Line 47 | Line 142 | register SRCINDEX  *si;                /* source index (modified to
142          if (!skipparts(cent, size, parr, si->spt))
143                  error(CONSISTENCY, "bad source partition in nextssamp");
144                                          /* compute sample */
145 +        srcp = source + si->sn;
146          if (dstrsrc > FTINY) {                  /* jitter sample */
147                  dimlist[ndims] = si->sn + 8831;
148                  dimlist[ndims+1] = si->sp + 3109;
149                  d = urand(ilhash(dimlist,ndims+2)+samplendx);
150 <                if (source[si->sn].sflags & SFLAT) {
150 >                if (srcp->sflags & SFLAT) {
151                          multisamp(vpos, 2, d);
152 <                        vpos[2] = 0.5;
152 >                        vpos[SW] = 0.5;
153                  } else
154                          multisamp(vpos, 3, d);
155                  for (i = 0; i < 3; i++)
156                          vpos[i] = dstrsrc * (1. - 2.*vpos[i]) *
157 <                                        (double)size[i]/MAXSPART;
157 >                                        (double)size[i]*(1.0/MAXSPART);
158          } else
159                  vpos[0] = vpos[1] = vpos[2] = 0.0;
160  
161 <        for (i = 0; i < 3; i++)
162 <                vpos[i] += (double)cent[i]/MAXSPART;
161 >        VSUM(vpos, vpos, cent, 1.0/MAXSPART);
162 >                                        /* avoid circular aiming failures */
163 >        if ((srcp->sflags & SCIR) && (si->np > 1) | (dstrsrc > 0.7)) {
164 >                FVECT   trim;
165 >                if (srcp->sflags & (SFLAT|SDISTANT)) {
166 >                        d = 1.12837917;         /* correct setflatss() */
167 >                        trim[SU] = d*sqrt(1.0 - 0.5*vpos[SV]*vpos[SV]);
168 >                        trim[SV] = d*sqrt(1.0 - 0.5*vpos[SU]*vpos[SU]);
169 >                        trim[SW] = 0.0;
170 >                } else {
171 >                        trim[SW] = trim[SU] = vpos[SU]*vpos[SU];
172 >                        d = vpos[SV]*vpos[SV];
173 >                        if (d > trim[SW]) trim[SW] = d;
174 >                        trim[SU] += d;
175 >                        d = vpos[SW]*vpos[SW];
176 >                        if (d > trim[SW]) trim[SW] = d;
177 >                        trim[SU] += d;
178 >                        if (trim[SU] > FTINY*FTINY) {
179 >                                d = 1.0/0.7236; /* correct sphsetsrc() */
180 >                                trim[SW] = trim[SV] = trim[SU] =
181 >                                                d*sqrt(trim[SW]/trim[SU]);
182 >                        } else
183 >                                trim[SW] = trim[SV] = trim[SU] = 0.0;
184 >                }
185 >                for (i = 0; i < 3; i++)
186 >                        vpos[i] *= trim[i];
187 >        }
188                                          /* compute direction */
189          for (i = 0; i < 3; i++)
190 <                r->rdir[i] = source[si->sn].sloc[i] +
191 <                                vpos[SU]*source[si->sn].ss[SU][i] +
192 <                                vpos[SV]*source[si->sn].ss[SV][i] +
193 <                                vpos[SW]*source[si->sn].ss[SW][i];
190 >                r->rdir[i] = srcp->sloc[i] +
191 >                                vpos[SU]*srcp->ss[SU][i] +
192 >                                vpos[SV]*srcp->ss[SV][i] +
193 >                                vpos[SW]*srcp->ss[SW][i];
194  
195 <        if (!(source[si->sn].sflags & SDISTANT))
196 <                for (i = 0; i < 3; i++)
76 <                        r->rdir[i] -= r->rorg[i];
195 >        if (!(srcp->sflags & SDISTANT))
196 >                VSUB(r->rdir, r->rdir, r->rorg);
197                                          /* compute distance */
198          if ((d = normalize(r->rdir)) == 0.0)
199 <                return(nextssamp(r, si));       /* at source! */
199 >                goto nextsample;                /* at source! */
200  
201                                          /* compute sample size */
202 <        si->dom  = source[si->sn].ss2;
203 <        if (source[si->sn].sflags & SFLAT) {
204 <                si->dom *= sflatform(si->sn, r->rdir);
205 <                si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART);
206 <        } else if (source[si->sn].sflags & SCYL) {
207 <                si->dom *= scylform(si->sn, r->rdir);
88 <                si->dom *= (double)size[SU]/MAXSPART;
202 >        if (srcp->sflags & SFLAT) {
203 >                si->dom = sflatform(si->sn, r->rdir);
204 >                si->dom *= size[SU]*size[SV]*(1.0/MAXSPART/MAXSPART);
205 >        } else if (srcp->sflags & SCYL) {
206 >                si->dom = scylform(si->sn, r->rdir);
207 >                si->dom *= size[SU]*(1.0/MAXSPART);
208          } else {
209 <                si->dom *= (double)(size[SU]*size[SV]*size[SW]) /
210 <                                (MAXSPART*MAXSPART*MAXSPART) ;
209 >                si->dom = size[SU]*size[SV]*(double)size[SW] *
210 >                                (1.0/MAXSPART/MAXSPART/MAXSPART) ;
211          }
212 <        if (source[si->sn].sflags & SDISTANT)
212 >        if (srcp->sflags & SDISTANT) {
213 >                si->dom *= srcp->ss2;
214                  return(FHUGE);
215 <        si->dom /= d*d;
215 >        }
216 >        if (si->dom <= 1e-4)
217 >                goto nextsample;                /* behind source? */
218 >        si->dom *= srcp->ss2/(d*d);
219          return(d);              /* sample OK, return distance */
220   }
221  
222  
223 < skipparts(ct, sz, pp, pt)               /* skip to requested partition */
224 < int  ct[3], sz[3];              /* center and size of partition (returned) */
225 < register int  pp[2];            /* current index, number to skip (modified) */
226 < unsigned char  *pt;             /* partition array */
223 > int
224 > skipparts(                      /* skip to requested partition */
225 >        int  ct[3],
226 >        int  sz[3],             /* center and size of partition (returned) */
227 >        int  pp[2],             /* current index, number to skip (modified) */
228 >        unsigned char  *pt      /* partition array */
229 > )
230   {
231 <        register int  p;
231 >        int  p;
232                                          /* check this partition */
233          p = spart(pt, pp[0]);
234          pp[0]++;
235 <        if (p == S0)                    /* leaf partition */
235 >        if (p == S0) {                  /* leaf partition */
236                  if (pp[1]) {
237                          pp[1]--;
238                          return(0);      /* not there yet */
239                  } else
240                          return(1);      /* we've arrived */
241 +        }
242                                  /* else check lower */
243          sz[p] >>= 1;
244          ct[p] -= sz[p];
# Line 128 | Line 255 | unsigned char  *pt;            /* partition array */
255   }
256  
257  
258 < nopart(si, r)                   /* single source partition */
259 < register SRCINDEX  *si;
260 < RAY  *r;
258 > void
259 > nopart(                         /* single source partition */
260 >        SRCINDEX  *si,
261 >        RAY  *r
262 > )
263   {
264          clrpart(si->spt);
265          setpart(si->spt, 0, S0);
# Line 138 | Line 267 | RAY  *r;
267   }
268  
269  
141 – cylpart(si, r)                  /* partition a cylinder */
142 – SRCINDEX  *si;
143 – register RAY  *r;
144 – {
145 –        double  dist2, safedist2, dist2cent, rad2;
146 –        FVECT  v;
147 –        register SRCREC  *sp;
148 –        int  pi;
149 –                                        /* first check point location */
150 –        clrpart(si->spt);
151 –        sp = source + si->sn;
152 –        rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
153 –        v[0] = r->rorg[0] - sp->sloc[0];
154 –        v[1] = r->rorg[1] - sp->sloc[1];
155 –        v[2] = r->rorg[2] - sp->sloc[2];
156 –        dist2 = DOT(v,sp->ss[SU]);
157 –        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
158 –        dist2 *= dist2 / safedist2;
159 –        dist2cent = DOT(v,v);
160 –        dist2 = dist2cent - dist2;
161 –        if (dist2 <= rad2) {            /* point inside extended cylinder */
162 –                si->np = 0;
163 –                return;
164 –        }
165 –        safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
166 –        if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
167 –                        dist2cent >= safedist2) {       /* or too far */
168 –                setpart(si->spt, 0, S0);
169 –                si->np = 1;
170 –                return;
171 –        }
172 –        pi = 0;
173 –        si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
174 –                        sp->sloc, sp->ss[SU], safedist2);
175 – }
176 –
177 –
270   static int
271 < cyl_partit(ro, pt, pi, mp, cent, axis, d2)      /* slice a cylinder */
272 < FVECT  ro;
273 < unsigned char  *pt;
274 < register int  *pi;
275 < int  mp;
276 < FVECT  cent, axis;
277 < double  d2;
271 > cyl_partit(                             /* slice a cylinder */
272 >        FVECT  ro,
273 >        unsigned char  *pt,
274 >        int  *pi,
275 >        int  mp,
276 >        FVECT  cent,
277 >        FVECT  axis,
278 >        double  d2
279 > )
280   {
281          FVECT  newct, newax;
282          int  npl, npu;
# Line 214 | Line 308 | double  d2;
308   }
309  
310  
311 < flatpart(si, r)                         /* partition a flat source */
312 < register SRCINDEX  *si;
313 < register RAY  *r;
311 > void
312 > cylpart(                        /* partition a cylinder */
313 >        SRCINDEX  *si,
314 >        RAY  *r
315 > )
316   {
317 <        register FLOAT  *vp;
317 >        double  dist2, safedist2, dist2cent, rad2;
318          FVECT  v;
319 <        double  du2, dv2;
319 >        SRCREC  *sp;
320          int  pi;
321 <
321 >                                        /* first check point location */
322          clrpart(si->spt);
323 <        vp = source[si->sn].sloc;
324 <        v[0] = r->rorg[0] - vp[0];
325 <        v[1] = r->rorg[1] - vp[1];
326 <        v[2] = r->rorg[2] - vp[2];
327 <        vp = source[si->sn].snorm;
328 <        if (DOT(v,vp) <= FTINY) {       /* behind source */
323 >        sp = source + si->sn;
324 >        rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
325 >        v[0] = r->rorg[0] - sp->sloc[0];
326 >        v[1] = r->rorg[1] - sp->sloc[1];
327 >        v[2] = r->rorg[2] - sp->sloc[2];
328 >        dist2 = DOT(v,sp->ss[SU]);
329 >        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
330 >        dist2 *= dist2 / safedist2;
331 >        dist2cent = DOT(v,v);
332 >        dist2 = dist2cent - dist2;
333 >        if (dist2 <= rad2) {            /* point inside extended cylinder */
334                  si->np = 0;
335                  return;
336          }
337 <        dv2 = 2.*r->rweight/srcsizerat;
338 <        dv2 *= dv2;
339 <        vp = source[si->sn].ss[SU];
340 <        du2 = dv2 * DOT(vp,vp);
341 <        vp = source[si->sn].ss[SV];
342 <        dv2 *= DOT(vp,vp);
337 >        safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
338 >        if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
339 >                        dist2cent >= safedist2) {       /* or too far */
340 >                setpart(si->spt, 0, S0);
341 >                si->np = 1;
342 >                return;
343 >        }
344          pi = 0;
345 <        si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
346 <                source[si->sn].sloc,
245 <                source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
345 >        si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
346 >                        sp->sloc, sp->ss[SU], safedist2);
347   }
348  
349  
350   static int
351 < flt_partit(ro, pt, pi, mp, cent, u, v, du2, dv2)        /* partition flatty */
352 < FVECT  ro;
353 < unsigned char  *pt;
354 < register int  *pi;
355 < int  mp;
356 < FVECT  cent, u, v;
357 < double  du2, dv2;
351 > flt_partit(                             /* partition flatty */
352 >        FVECT  ro,
353 >        unsigned char  *pt,
354 >        int  *pi,
355 >        int  mp,
356 >        FVECT  cent,
357 >        FVECT  u,
358 >        FVECT  v,
359 >        double  du2,
360 >        double  dv2
361 > )
362   {
363          double  d2;
364          FVECT  newct, newax;
# Line 297 | Line 402 | double  du2, dv2;
402   }
403  
404  
405 + void
406 + flatpart(                               /* partition a flat source */
407 +        SRCINDEX  *si,
408 +        RAY  *r
409 + )
410 + {
411 +        RREAL  *vp;
412 +        FVECT  v;
413 +        double  du2, dv2;
414 +        int  pi;
415 +
416 +        clrpart(si->spt);
417 +        vp = source[si->sn].sloc;
418 +        v[0] = r->rorg[0] - vp[0];
419 +        v[1] = r->rorg[1] - vp[1];
420 +        v[2] = r->rorg[2] - vp[2];
421 +        vp = source[si->sn].snorm;
422 +        if (DOT(v,vp) <= 0.) {          /* behind source */
423 +                si->np = 0;
424 +                return;
425 +        }
426 +        dv2 = 2.*r->rweight/srcsizerat;
427 +        dv2 *= dv2;
428 +        vp = source[si->sn].ss[SU];
429 +        du2 = dv2 * DOT(vp,vp);
430 +        vp = source[si->sn].ss[SV];
431 +        dv2 *= DOT(vp,vp);
432 +        pi = 0;
433 +        si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
434 +                source[si->sn].sloc,
435 +                source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
436 + }
437 +
438 +
439   double
440 < scylform(sn, dir)               /* compute cosine for cylinder's projection */
441 < int  sn;
442 < register FVECT  dir;            /* assume normalized */
440 > scylform(                       /* compute cosine for cylinder's projection */
441 >        int  sn,
442 >        FVECT  dir              /* assume normalized */
443 > )
444   {
445 <        register FLOAT  *dv;
445 >        RREAL  *dv;
446          double  d;
447  
448          dv = source[sn].ss[SU];

Diff Legend

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