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.17 by greg, Sat Oct 23 18:13:54 2010 UTC vs.
Revision 2.23 by greg, Fri Dec 13 00:50:55 2024 UTC

# Line 15 | Line 15 | static const char      RCSid[] = "$Id$";
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 < static int  cyl_partit(), flt_partit();
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
101 +        if (r->rsrc < -1 &&     /* ray has custom skip flags? */
102 +                        sskip_chk(sskip_flags(r->rsrc), sn))
103 +                return(1);
104 + #endif
105 +        if ((sp->sflags & (SPROX|SDISTANT)) != SPROX)
106 +                return(0);
107 +
108 +        return(dist2(r->rorg, sp->sloc) >
109 +                        (sp->sl.prox + sp->srad)*(sp->sl.prox + sp->srad));
110 + }
111 +
112   double
113 < nextssamp(r, si)                /* compute sample for source, rtn. distance */
114 < register RAY  *r;               /* origin is read, direction is set */
115 < register SRCINDEX  *si;         /* source index (modified to current) */
113 > nextssamp(                      /* compute sample for source, rtn. distance */
114 >        RAY  *r,                /* origin is read, direction is set */
115 >        SRCINDEX  *si           /* source index (modified to current) */
116 > )
117   {
118          int  cent[3], size[3], parr[2];
119          SRCREC  *srcp;
120 <        FVECT  vpos;
120 >        double  vpos[3];
121          double  d;
122 <        register int  i;
122 >        int  i;
123   nextsample:
124          while (++si->sp >= si->np) {    /* get next sample */
125                  if (++si->sn >= nsources)
126                          return(0.0);    /* no more */
127 <                if (source[si->sn].sflags & SSKIP)
127 >                if (srcskip(si->sn, r))
128                          si->np = 0;
129                  else if (srcsizerat <= FTINY)
130                          nopart(si, r);
# Line 70 | Line 161 | nextsample:
161  
162          VSUM(vpos, vpos, cent, 1.0/MAXSPART);
163                                          /* avoid circular aiming failures */
164 <        if ((srcp->sflags & SCIR) && (si->np > 1 || dstrsrc > 0.7)) {
164 >        if ((srcp->sflags & SCIR) && (si->np > 1) | (dstrsrc > 0.7)) {
165                  FVECT   trim;
166                  if (srcp->sflags & (SFLAT|SDISTANT)) {
167                          d = 1.12837917;         /* correct setflatss() */
# Line 131 | Line 222 | nextsample:
222  
223  
224   int
225 < skipparts(ct, sz, pp, pt)               /* skip to requested partition */
226 < int  ct[3], sz[3];              /* center and size of partition (returned) */
227 < register int  pp[2];            /* current index, number to skip (modified) */
228 < unsigned char  *pt;             /* partition array */
225 > skipparts(                      /* skip to requested partition */
226 >        int  ct[3],
227 >        int  sz[3],             /* center and size of partition (returned) */
228 >        int  pp[2],             /* current index, number to skip (modified) */
229 >        unsigned char  *pt      /* partition array */
230 > )
231   {
232 <        register int  p;
232 >        int  p;
233                                          /* check this partition */
234          p = spart(pt, pp[0]);
235          pp[0]++;
# Line 164 | Line 257 | unsigned char  *pt;            /* partition array */
257  
258  
259   void
260 < nopart(si, r)                   /* single source partition */
261 < register SRCINDEX  *si;
262 < RAY  *r;
260 > nopart(                         /* single source partition */
261 >        SRCINDEX  *si,
262 >        RAY  *r
263 > )
264   {
265          clrpart(si->spt);
266          setpart(si->spt, 0, S0);
# Line 174 | Line 268 | RAY  *r;
268   }
269  
270  
177 – void
178 – cylpart(si, r)                  /* partition a cylinder */
179 – SRCINDEX  *si;
180 – register RAY  *r;
181 – {
182 –        double  dist2, safedist2, dist2cent, rad2;
183 –        FVECT  v;
184 –        register SRCREC  *sp;
185 –        int  pi;
186 –                                        /* first check point location */
187 –        clrpart(si->spt);
188 –        sp = source + si->sn;
189 –        rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
190 –        v[0] = r->rorg[0] - sp->sloc[0];
191 –        v[1] = r->rorg[1] - sp->sloc[1];
192 –        v[2] = r->rorg[2] - sp->sloc[2];
193 –        dist2 = DOT(v,sp->ss[SU]);
194 –        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
195 –        dist2 *= dist2 / safedist2;
196 –        dist2cent = DOT(v,v);
197 –        dist2 = dist2cent - dist2;
198 –        if (dist2 <= rad2) {            /* point inside extended cylinder */
199 –                si->np = 0;
200 –                return;
201 –        }
202 –        safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
203 –        if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
204 –                        dist2cent >= safedist2) {       /* or too far */
205 –                setpart(si->spt, 0, S0);
206 –                si->np = 1;
207 –                return;
208 –        }
209 –        pi = 0;
210 –        si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
211 –                        sp->sloc, sp->ss[SU], safedist2);
212 – }
213 –
214 –
271   static int
272 < cyl_partit(ro, pt, pi, mp, cent, axis, d2)      /* slice a cylinder */
273 < FVECT  ro;
274 < unsigned char  *pt;
275 < register int  *pi;
276 < int  mp;
277 < FVECT  cent, axis;
278 < double  d2;
272 > cyl_partit(                             /* slice a cylinder */
273 >        FVECT  ro,
274 >        unsigned char  *pt,
275 >        int  *pi,
276 >        int  mp,
277 >        FVECT  cent,
278 >        FVECT  axis,
279 >        double  d2
280 > )
281   {
282          FVECT  newct, newax;
283          int  npl, npu;
# Line 252 | Line 310 | double  d2;
310  
311  
312   void
313 < flatpart(si, r)                         /* partition a flat source */
314 < register SRCINDEX  *si;
315 < register RAY  *r;
313 > cylpart(                        /* partition a cylinder */
314 >        SRCINDEX  *si,
315 >        RAY  *r
316 > )
317   {
318 <        register RREAL  *vp;
318 >        double  dist2, safedist2, dist2cent, rad2;
319          FVECT  v;
320 <        double  du2, dv2;
320 >        SRCREC  *sp;
321          int  pi;
322 <
322 >                                        /* first check point location */
323          clrpart(si->spt);
324 <        vp = source[si->sn].sloc;
325 <        v[0] = r->rorg[0] - vp[0];
326 <        v[1] = r->rorg[1] - vp[1];
327 <        v[2] = r->rorg[2] - vp[2];
328 <        vp = source[si->sn].snorm;
329 <        if (DOT(v,vp) <= 0.) {          /* behind source */
324 >        sp = source + si->sn;
325 >        rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
326 >        v[0] = r->rorg[0] - sp->sloc[0];
327 >        v[1] = r->rorg[1] - sp->sloc[1];
328 >        v[2] = r->rorg[2] - sp->sloc[2];
329 >        dist2 = DOT(v,sp->ss[SU]);
330 >        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
331 >        dist2 *= dist2 / safedist2;
332 >        dist2cent = DOT(v,v);
333 >        dist2 = dist2cent - dist2;
334 >        if (dist2 <= rad2) {            /* point inside extended cylinder */
335                  si->np = 0;
336                  return;
337          }
338 <        dv2 = 2.*r->rweight/srcsizerat;
339 <        dv2 *= dv2;
340 <        vp = source[si->sn].ss[SU];
341 <        du2 = dv2 * DOT(vp,vp);
342 <        vp = source[si->sn].ss[SV];
343 <        dv2 *= DOT(vp,vp);
338 >        safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
339 >        if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
340 >                        dist2cent >= safedist2) {       /* or too far */
341 >                setpart(si->spt, 0, S0);
342 >                si->np = 1;
343 >                return;
344 >        }
345          pi = 0;
346 <        si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
347 <                source[si->sn].sloc,
283 <                source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
346 >        si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
347 >                        sp->sloc, sp->ss[SU], safedist2);
348   }
349  
350  
351   static int
352 < flt_partit(ro, pt, pi, mp, cent, u, v, du2, dv2)        /* partition flatty */
353 < FVECT  ro;
354 < unsigned char  *pt;
355 < register int  *pi;
356 < int  mp;
357 < FVECT  cent, u, v;
358 < double  du2, dv2;
352 > flt_partit(                             /* partition flatty */
353 >        FVECT  ro,
354 >        unsigned char  *pt,
355 >        int  *pi,
356 >        int  mp,
357 >        FVECT  cent,
358 >        FVECT  u,
359 >        FVECT  v,
360 >        double  du2,
361 >        double  dv2
362 > )
363   {
364          double  d2;
365          FVECT  newct, newax;
# Line 335 | Line 403 | double  du2, dv2;
403   }
404  
405  
406 + void
407 + flatpart(                               /* partition a flat source */
408 +        SRCINDEX  *si,
409 +        RAY  *r
410 + )
411 + {
412 +        RREAL  *vp;
413 +        FVECT  v;
414 +        double  du2, dv2;
415 +        int  pi;
416 +
417 +        clrpart(si->spt);
418 +        vp = source[si->sn].sloc;
419 +        v[0] = r->rorg[0] - vp[0];
420 +        v[1] = r->rorg[1] - vp[1];
421 +        v[2] = r->rorg[2] - vp[2];
422 +        vp = source[si->sn].snorm;
423 +        if (DOT(v,vp) <= 0.) {          /* behind source */
424 +                si->np = 0;
425 +                return;
426 +        }
427 +        dv2 = 2.*r->rweight/srcsizerat;
428 +        dv2 *= dv2;
429 +        vp = source[si->sn].ss[SU];
430 +        du2 = dv2 * DOT(vp,vp);
431 +        vp = source[si->sn].ss[SV];
432 +        dv2 *= DOT(vp,vp);
433 +        pi = 0;
434 +        si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
435 +                source[si->sn].sloc,
436 +                source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
437 + }
438 +
439 +
440   double
441 < scylform(sn, dir)               /* compute cosine for cylinder's projection */
442 < int  sn;
443 < register FVECT  dir;            /* assume normalized */
441 > scylform(                       /* compute cosine for cylinder's projection */
442 >        int  sn,
443 >        FVECT  dir              /* assume normalized */
444 > )
445   {
446 <        register RREAL  *dv;
446 >        RREAL  *dv;
447          double  d;
448  
449          dv = source[sn].ss[SU];

Diff Legend

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