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

Diff Legend

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