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.2 by greg, Mon Oct 21 14:27:36 1991 UTC vs.
Revision 2.21 by greg, Sat Nov 9 15:21:32 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  "standard.h"
10 > #include "copyright.h"
11  
12 < #include  "object.h"
12 > #include  "ray.h"
13  
14   #include  "source.h"
15  
16   #include  "random.h"
17  
18  
19 < extern int  dimlist[];          /* dimension list for distribution */
20 < extern int  ndims;              /* number of dimensions so far */
21 < extern int  samplendx;          /* index for this sample */
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  
30 +        if ((sp->sflags & (SPROX|SDISTANT)) != SPROX)
31 +                return(0);
32 +
33 +        return(dist2(r->rorg, sp->sloc) >
34 +                        (sp->sl.prox + sp->srad)*(sp->sl.prox + sp->srad));
35 + }
36 +
37   double
38 < nextssamp(org, dir, si)         /* compute sample for source, rtn. distance */
39 < FVECT  org, dir;                /* origin is read only, direction is set */
40 < register SRCINDEX  *si;         /* source index (modified to current) */
38 > nextssamp(                      /* compute sample for source, rtn. distance */
39 >        RAY  *r,                /* origin is read, direction is set */
40 >        SRCINDEX  *si           /* source index (modified to current) */
41 > )
42   {
43          int  cent[3], size[3], parr[2];
44 <        FVECT  vpos;
44 >        SRCREC  *srcp;
45 >        double  vpos[3];
46          double  d;
47 <        register int  i;
48 < tryagain:
47 >        int  i;
48 > nextsample:
49          while (++si->sp >= si->np) {    /* get next sample */
50                  if (++si->sn >= nsources)
51                          return(0.0);    /* no more */
52 <                if (srcsizerat <= FTINY)
53 <                        nopart(si, org);
52 >                if (srcskip(si->sn, r))
53 >                        si->np = 0;
54 >                else if (srcsizerat <= FTINY)
55 >                        nopart(si, r);
56                  else {
57                          for (i = si->sn; source[i].sflags & SVIRTUAL;
58                                          i = source[i].sa.sv.sn)
59                                  ;               /* partition source */
60 <                        (*sfun[source[i].so->otype].of->partit)(si, org);
60 >                        (*sfun[source[i].so->otype].of->partit)(si, r);
61                  }
62                  si->sp = -1;
63          }
# Line 52 | Line 68 | tryagain:
68          if (!skipparts(cent, size, parr, si->spt))
69                  error(CONSISTENCY, "bad source partition in nextssamp");
70                                          /* compute sample */
71 +        srcp = source + si->sn;
72          if (dstrsrc > FTINY) {                  /* jitter sample */
73                  dimlist[ndims] = si->sn + 8831;
74                  dimlist[ndims+1] = si->sp + 3109;
75                  d = urand(ilhash(dimlist,ndims+2)+samplendx);
76 <                if (source[si->sn].sflags & SFLAT) {
76 >                if (srcp->sflags & SFLAT) {
77                          multisamp(vpos, 2, d);
78 <                        vpos[2] = 0.5;
78 >                        vpos[SW] = 0.5;
79                  } else
80                          multisamp(vpos, 3, d);
81                  for (i = 0; i < 3; i++)
82                          vpos[i] = dstrsrc * (1. - 2.*vpos[i]) *
83 <                                        (double)size[i]/MAXSPART;
83 >                                        (double)size[i]*(1.0/MAXSPART);
84          } else
85                  vpos[0] = vpos[1] = vpos[2] = 0.0;
86  
87 <        for (i = 0; i < 3; i++)
88 <                vpos[i] += (double)cent[i]/MAXSPART;
87 >        VSUM(vpos, vpos, cent, 1.0/MAXSPART);
88 >                                        /* avoid circular aiming failures */
89 >        if ((srcp->sflags & SCIR) && (si->np > 1) | (dstrsrc > 0.7)) {
90 >                FVECT   trim;
91 >                if (srcp->sflags & (SFLAT|SDISTANT)) {
92 >                        d = 1.12837917;         /* correct setflatss() */
93 >                        trim[SU] = d*sqrt(1.0 - 0.5*vpos[SV]*vpos[SV]);
94 >                        trim[SV] = d*sqrt(1.0 - 0.5*vpos[SU]*vpos[SU]);
95 >                        trim[SW] = 0.0;
96 >                } else {
97 >                        trim[SW] = trim[SU] = vpos[SU]*vpos[SU];
98 >                        d = vpos[SV]*vpos[SV];
99 >                        if (d > trim[SW]) trim[SW] = d;
100 >                        trim[SU] += d;
101 >                        d = vpos[SW]*vpos[SW];
102 >                        if (d > trim[SW]) trim[SW] = d;
103 >                        trim[SU] += d;
104 >                        if (trim[SU] > FTINY*FTINY) {
105 >                                d = 1.0/0.7236; /* correct sphsetsrc() */
106 >                                trim[SW] = trim[SV] = trim[SU] =
107 >                                                d*sqrt(trim[SW]/trim[SU]);
108 >                        } else
109 >                                trim[SW] = trim[SV] = trim[SU] = 0.0;
110 >                }
111 >                for (i = 0; i < 3; i++)
112 >                        vpos[i] *= trim[i];
113 >        }
114                                          /* compute direction */
115          for (i = 0; i < 3; i++)
116 <                dir[i] = source[si->sn].sloc[i] +
117 <                                vpos[SU]*source[si->sn].ss[SU][i] +
118 <                                vpos[SV]*source[si->sn].ss[SV][i] +
119 <                                vpos[SW]*source[si->sn].ss[SW][i];
116 >                r->rdir[i] = srcp->sloc[i] +
117 >                                vpos[SU]*srcp->ss[SU][i] +
118 >                                vpos[SV]*srcp->ss[SV][i] +
119 >                                vpos[SW]*srcp->ss[SW][i];
120  
121 <        if (!(source[si->sn].sflags & SDISTANT))
122 <                for (i = 0; i < 3; i++)
81 <                        dir[i] -= org[i];
121 >        if (!(srcp->sflags & SDISTANT))
122 >                VSUB(r->rdir, r->rdir, r->rorg);
123                                          /* compute distance */
124 <        if ((d = normalize(dir)) == 0.0)
125 <                goto tryagain;  /* at source! */
124 >        if ((d = normalize(r->rdir)) == 0.0)
125 >                goto nextsample;                /* at source! */
126  
127                                          /* compute sample size */
128 <        si->dom  = source[si->sn].ss2;
129 <        if (source[si->sn].sflags & SFLAT) {
130 <                si->dom *= sflatform(si->sn, dir);
131 <                if (si->dom <= FTINY) {         /* behind source */
132 <                        si->sp = si->np;
133 <                        goto tryagain;
93 <                }
94 <                si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART);
95 <        } else if (source[si->sn].sflags & SCYL) {
96 <                si->dom *= scylform(si->sn, dir);
97 <                si->dom *= (double)size[SU]/MAXSPART;
128 >        if (srcp->sflags & SFLAT) {
129 >                si->dom = sflatform(si->sn, r->rdir);
130 >                si->dom *= size[SU]*size[SV]*(1.0/MAXSPART/MAXSPART);
131 >        } else if (srcp->sflags & SCYL) {
132 >                si->dom = scylform(si->sn, r->rdir);
133 >                si->dom *= size[SU]*(1.0/MAXSPART);
134          } else {
135 <                si->dom *= (double)(size[SU]*size[SV]*size[SW]) /
136 <                                (MAXSPART*MAXSPART*MAXSPART) ;
135 >                si->dom = size[SU]*size[SV]*(double)size[SW] *
136 >                                (1.0/MAXSPART/MAXSPART/MAXSPART) ;
137          }
138 <        if (source[si->sn].sflags & SDISTANT)
138 >        if (srcp->sflags & SDISTANT) {
139 >                si->dom *= srcp->ss2;
140                  return(FHUGE);
141 <        si->dom /= d*d;
141 >        }
142 >        if (si->dom <= 1e-4)
143 >                goto nextsample;                /* behind source? */
144 >        si->dom *= srcp->ss2/(d*d);
145          return(d);              /* sample OK, return distance */
146   }
147  
148  
149 < skipparts(ct, sz, pp, pt)               /* skip to requested partition */
150 < int  ct[3], sz[3];              /* center and size of partition (returned) */
151 < register int  pp[2];            /* current index, number to skip (modified) */
152 < unsigned char  *pt;             /* partition array */
149 > int
150 > skipparts(                      /* skip to requested partition */
151 >        int  ct[3],
152 >        int  sz[3],             /* center and size of partition (returned) */
153 >        int  pp[2],             /* current index, number to skip (modified) */
154 >        unsigned char  *pt      /* partition array */
155 > )
156   {
157 <        register int  p;
157 >        int  p;
158                                          /* check this partition */
159          p = spart(pt, pp[0]);
160          pp[0]++;
161 <        if (p == S0)                    /* leaf partition */
161 >        if (p == S0) {                  /* leaf partition */
162                  if (pp[1]) {
163                          pp[1]--;
164                          return(0);      /* not there yet */
165                  } else
166                          return(1);      /* we've arrived */
167 +        }
168                                  /* else check lower */
169          sz[p] >>= 1;
170          ct[p] -= sz[p];
# Line 137 | Line 181 | unsigned char  *pt;            /* partition array */
181   }
182  
183  
184 < nopart(si, ro)                  /* single source partition */
185 < register SRCINDEX  *si;
186 < FVECT  ro;
184 > void
185 > nopart(                         /* single source partition */
186 >        SRCINDEX  *si,
187 >        RAY  *r
188 > )
189   {
190          clrpart(si->spt);
191          setpart(si->spt, 0, S0);
# Line 147 | Line 193 | FVECT  ro;
193   }
194  
195  
150 cylpart(si, ro)                 /* partition a cylinder */
151 SRCINDEX  *si;
152 FVECT  ro;
153 {
154        double  dist2, safedist2, dist2cent, rad2;
155        FVECT  v;
156        register SRCREC  *sp;
157        int  pi;
158                                        /* first check point location */
159        clrpart(si->spt);
160        sp = &source[si->sn];
161        rad2 = 1.273 * DOT(sp->ss[SV],sp->ss[SV]);
162        v[0] = ro[0] - sp->sloc[0];
163        v[1] = ro[1] - sp->sloc[1];
164        v[2] = ro[2] - sp->sloc[2];
165        dist2 = DOT(v,sp->ss[SU]);
166        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
167        dist2 *= dist2 / safedist2;
168        dist2cent = DOT(v,v);
169        dist2 = dist2cent - dist2;
170        if (dist2 <= rad2) {            /* point inside extended cylinder */
171                si->np = 0;
172                return;
173        }
174        safedist2 *= 4./(srcsizerat*srcsizerat);
175        if (dist2 <= 4.*rad2 ||         /* point too close to subdivide? */
176                        dist2cent >= safedist2) {
177                setpart(si->spt, 0, S0);
178                si->np = 1;
179                return;
180        }
181        pi = 0;
182        si->np = cyl_partit(ro, si->spt, &pi, MAXSPART,
183                        sp->sloc, sp->ss[SU], safedist2);
184 }
185
186
196   static int
197 < cyl_partit(ro, pt, pi, mp, cent, axis, d2)      /* slice a cylinder */
198 < FVECT  ro;
199 < unsigned char  *pt;
200 < register int  *pi;
201 < int  mp;
202 < FVECT  cent, axis;
203 < double  d2;
197 > cyl_partit(                             /* slice a cylinder */
198 >        FVECT  ro,
199 >        unsigned char  *pt,
200 >        int  *pi,
201 >        int  mp,
202 >        FVECT  cent,
203 >        FVECT  axis,
204 >        double  d2
205 > )
206   {
207          FVECT  newct, newax;
208          int  npl, npu;
# Line 223 | Line 234 | double  d2;
234   }
235  
236  
237 < flatpart(si, ro)                        /* partition a flat source */
238 < register SRCINDEX  *si;
239 < FVECT  ro;
237 > void
238 > cylpart(                        /* partition a cylinder */
239 >        SRCINDEX  *si,
240 >        RAY  *r
241 > )
242   {
243 <        register double  *vp;
244 <        double  du2, dv2;
243 >        double  dist2, safedist2, dist2cent, rad2;
244 >        FVECT  v;
245 >        SRCREC  *sp;
246          int  pi;
247 <
234 <        vp = source[si->sn].ss[SU];
235 <        du2 = 4./(srcsizerat*srcsizerat) * DOT(vp,vp);
236 <        vp = source[si->sn].ss[SV];
237 <        dv2 = 4./(srcsizerat*srcsizerat) * DOT(vp,vp);
247 >                                        /* first check point location */
248          clrpart(si->spt);
249 +        sp = source + si->sn;
250 +        rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
251 +        v[0] = r->rorg[0] - sp->sloc[0];
252 +        v[1] = r->rorg[1] - sp->sloc[1];
253 +        v[2] = r->rorg[2] - sp->sloc[2];
254 +        dist2 = DOT(v,sp->ss[SU]);
255 +        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
256 +        dist2 *= dist2 / safedist2;
257 +        dist2cent = DOT(v,v);
258 +        dist2 = dist2cent - dist2;
259 +        if (dist2 <= rad2) {            /* point inside extended cylinder */
260 +                si->np = 0;
261 +                return;
262 +        }
263 +        safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
264 +        if (dist2 <= 4.*rad2 ||         /* point too close to subdivide */
265 +                        dist2cent >= safedist2) {       /* or too far */
266 +                setpart(si->spt, 0, S0);
267 +                si->np = 1;
268 +                return;
269 +        }
270          pi = 0;
271 <        si->np = flt_partit(ro, si->spt, &pi, MAXSPART, source[si->sn].sloc,
272 <                source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
271 >        si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
272 >                        sp->sloc, sp->ss[SU], safedist2);
273   }
274  
275  
276   static int
277 < flt_partit(ro, pt, pi, mp, cent, u, v, du2, dv2)        /* partition flatty */
278 < FVECT  ro;
279 < unsigned char  *pt;
280 < register int  *pi;
281 < int  mp;
282 < FVECT  cent, u, v;
283 < double  du2, dv2;
277 > flt_partit(                             /* partition flatty */
278 >        FVECT  ro,
279 >        unsigned char  *pt,
280 >        int  *pi,
281 >        int  mp,
282 >        FVECT  cent,
283 >        FVECT  u,
284 >        FVECT  v,
285 >        double  du2,
286 >        double  dv2
287 > )
288   {
289          double  d2;
290          FVECT  newct, newax;
# Line 293 | Line 328 | double  du2, dv2;
328   }
329  
330  
331 + void
332 + flatpart(                               /* partition a flat source */
333 +        SRCINDEX  *si,
334 +        RAY  *r
335 + )
336 + {
337 +        RREAL  *vp;
338 +        FVECT  v;
339 +        double  du2, dv2;
340 +        int  pi;
341 +
342 +        clrpart(si->spt);
343 +        vp = source[si->sn].sloc;
344 +        v[0] = r->rorg[0] - vp[0];
345 +        v[1] = r->rorg[1] - vp[1];
346 +        v[2] = r->rorg[2] - vp[2];
347 +        vp = source[si->sn].snorm;
348 +        if (DOT(v,vp) <= 0.) {          /* behind source */
349 +                si->np = 0;
350 +                return;
351 +        }
352 +        dv2 = 2.*r->rweight/srcsizerat;
353 +        dv2 *= dv2;
354 +        vp = source[si->sn].ss[SU];
355 +        du2 = dv2 * DOT(vp,vp);
356 +        vp = source[si->sn].ss[SV];
357 +        dv2 *= DOT(vp,vp);
358 +        pi = 0;
359 +        si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
360 +                source[si->sn].sloc,
361 +                source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
362 + }
363 +
364 +
365   double
366 < scylform(sn, dir)               /* compute cosine for cylinder's projection */
367 < int  sn;
368 < register FVECT  dir;            /* assume normalized */
366 > scylform(                       /* compute cosine for cylinder's projection */
367 >        int  sn,
368 >        FVECT  dir              /* assume normalized */
369 > )
370   {
371 <        register double  *dv;
371 >        RREAL  *dv;
372          double  d;
373  
374          dv = source[sn].ss[SU];

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines