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

Comparing ray/src/hd/rholo2.c (file contents):
Revision 3.19 by gwlarson, Tue Dec 1 15:47:05 1998 UTC vs.
Revision 3.20 by gwlarson, Thu Dec 3 15:21:05 1998 UTC

# Line 19 | Line 19 | struct gclim {
19          HOLO    *hp;                    /* holodeck pointer */
20          GCOORD  gc;                     /* grid cell */
21          FVECT   egp;                    /* eye grid point */
22 <        double  erg[2];                 /* eye range in wall grid coords */
22 >        double  erg2;                   /* mean square eye grid range */
23          double  gmin[2], gmax[2];       /* grid coordinate limits */
24   };                              /* a grid coordinate range */
25  
26  
27   static
28 < initeyelim(gcl, hd, gc)         /* initialize grid coordinate limits */
28 > initeyelim(gcl, hp, gc)         /* initialize grid coordinate limits */
29   register struct gclim   *gcl;
30 < int     hd;
30 > register HOLO   *hp;
31   GCOORD  *gc;
32   {
33          register FLOAT  *v;
34          register int    i;
35  
36 <        gcl->hp = hdlist[hd];
37 <        copystruct(&gcl->gc, gc);
38 <        hdgrid(gcl->egp, gcl->hp, myeye.vpt);
39 <        for (i = 0; i < 2; i++) {
40 <                v = gcl->hp->wg[((gcl->gc.w>>1)+i+1)%3];
41 <                gcl->erg[i] = myeye.rng * VLEN(v);
42 <                gcl->gmin[i] = FHUGE; gcl->gmax[i] = -FHUGE;
36 >        if (hp != NULL) {
37 >                hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
38 >                gcl->erg2 = 0;
39 >                for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
40 >                        gcl->erg2 += DOT(v,v);
41 >                gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
42          }
43 +        if (gc != NULL)
44 +                copystruct(&gcl->gc, gc);
45 +        gcl->gmin[0] = gcl->gmin[1] = FHUGE;
46 +        gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
47   }
48  
49  
50   static
51 < groweyelim(gcl, gp)             /* grow grid limits about eye point */
51 > groweyelim(gcl, gc, r0, r1)     /* grow grid limits about eye point */
52   register struct gclim   *gcl;
53 < FVECT   gp;
53 > GCOORD  *gc;
54 > double  r0, r1;
55   {
56 <        FVECT   ab;
57 <        double  l2, d, mult, wg;
58 <        register int    i, g;
56 >        FVECT   gp, ab;
57 >        double  vlen, plen, dv0, dv1;
58 >        double  rd2, dwall, gpos;
59 >        int     eyeout;
60 >        register int    i, g0, g1;
61  
62 +        i = gc->w>>1;
63 +        if (gc->w&1)
64 +                eyeout = (gp[i] = gcl->hp->grid[i]) < gcl->egp[i];
65 +        else
66 +                eyeout = (gp[i] = 0) > gcl->egp[i];
67 +        gp[hdwg0[gc->w]] = gc->i[0] + r0;
68 +        gp[hdwg1[gc->w]] = gc->i[1] + r1;
69          VSUB(ab, gcl->egp, gp);
70 <        l2 = DOT(ab,ab);
71 <        if (l2 <= gcl->erg[0]*gcl->erg[1]) {
70 >        rd2 = DOT(ab,ab);
71 >        if (rd2 <= gcl->erg2) {
72                  gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
73                  gcl->gmax[0] = gcl->gmax[1] = FHUGE;
74                  return;
75          }
76 <        mult = gp[g = gcl->gc.w>>1];
77 <        if (gcl->gc.w&1)
78 <                mult -= gcl->hp->grid[g];
79 <        if (ab[g]*ab[g] > gcl->erg[0]*gcl->erg[1])
80 <                mult /= -ab[g];
81 <        else if (fabs(ab[hdwg0[gcl->gc.w]]) > fabs(ab[hdwg1[gcl->gc.w]]))
82 <                mult = (gcl->gc.i[0] + .5 - gp[hdwg0[gcl->gc.w]]) /
83 <                                ab[hdwg0[gcl->gc.w]];
84 <        else
85 <                mult = (gcl->gc.i[1] + .5 - gp[hdwg1[gcl->gc.w]]) /
86 <                                ab[hdwg1[gcl->gc.w]];
87 <        for (i = 0; i < 2; i++) {
88 <                g = ((gcl->gc.w>>1)+i+1)%3;
89 <                wg = gp[g] + mult*ab[g];
90 <                d = mult*gcl->erg[i];
91 <                if (d < 0.) d = -d;
92 <                if (wg - d < gcl->gmin[i])
93 <                        gcl->gmin[i] = wg - d;
94 <                if (wg + d > gcl->gmax[i])
95 <                        gcl->gmax[i] = wg + d;
76 >        rd2 = gcl->erg2 / rd2;
77 >        vlen = 1. - rd2;
78 >        plen = sqrt(rd2 * vlen);
79 >        g0 = gcl->gc.w>>1;
80 >        dwall = (gcl->gc.w&1 ? gcl->hp->grid[g0] : 0) - gp[g0];
81 >        for (i = 0; i < 4; i++) {
82 >                if (i == 2)
83 >                        plen = -plen;
84 >                g1 = (g0+(i&1)+1)%3;
85 >                dv0 = vlen*ab[g0] + plen*ab[g1];
86 >                dv1 = vlen*ab[g1] - plen*ab[g0];
87 >                if ((dv0 < 0 ^ dwall < 0 ^ eyeout) ||
88 >                                (dv0 <= FTINY && dv0 >= -FTINY)) {
89 >                        if (eyeout)
90 >                                dv1 = -dv1;
91 >                        if (dv1 > FTINY)
92 >                                gcl->gmax[i&1] = FHUGE;
93 >                        else if (dv1 < -FTINY)
94 >                                gcl->gmin[i&1] = -FHUGE;
95 >                } else {
96 >                        gpos = gp[g1] + dv1*dwall/dv0;
97 >                        if (gpos < gcl->gmin[i&1])
98 >                                gcl->gmin[i&1] = gpos;
99 >                        if (gpos > gcl->gmax[i&1])
100 >                                gcl->gmax[i&1] = gpos;
101 >                }
102          }
103   }
104  
# Line 97 | Line 116 | register struct gclim  *gcl;
116                          gcl->gmin[i] = gcl->gc.i[i];
117                  if (gcl->gmax[i] > gcl->gc.i[i]+1)
118                          gcl->gmax[i] = gcl->gc.i[i]+1;
119 <                if ((incell &= gcl->gmax[i] > gcl->gmin[i])) {
119 >                if (gcl->gmax[i] > gcl->gmin[i]) {
120                          rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
121                                          (1.-FTINY);
122                          rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
123                                          (1.-FTINY) - rrng[i][0];
124 <                        incell &= rrng[i][1] > 0;
125 <                }
124 >                } else
125 >                        rrng[i][0] = rrng[i][1] = 0;
126 >                incell &= rrng[i][1] > 0;
127          }
128          return(incell);
129   }
# Line 113 | Line 133 | packrays(rod, p)               /* pack ray origins and directions *
133   register float  *rod;
134   register PACKET *p;
135   {
116 #define gp      ro
117 #ifdef DEBUG
118        double dist2sum = 0.;
119        FVECT   vt;
120 #endif
136          int     nretries = p->nr + 2;
137          struct gclim    eyelim;
138          short   rrng0[2][2], rrng1[2][2];
# Line 130 | Line 145 | register PACKET        *p;
145          if (!hdbcoord(gc, hdlist[p->hd], p->bi))
146                  error(CONSISTENCY, "bad beam index in packrays");
147          if ((useyelim = myeye.rng > FTINY)) {
148 <                initeyelim(&eyelim, p->hd, gc);
149 <                gp[gc[1].w>>1] = gc[1].w&1 ?
150 <                                hdlist[p->hd]->grid[gc[1].w>>1] : 0;
136 <                gp[hdwg0[gc[1].w]] = gc[1].i[0];
137 <                gp[hdwg1[gc[1].w]] = gc[1].i[1];
138 <                groweyelim(&eyelim, gp);
139 <                gp[hdwg0[gc[1].w]]++;
140 <                gp[hdwg1[gc[1].w]]++;
141 <                groweyelim(&eyelim, gp);
148 >                initeyelim(&eyelim, hdlist[p->hd], gc);
149 >                groweyelim(&eyelim, gc+1, 0., 0.);
150 >                groweyelim(&eyelim, gc+1, 1., 1.);
151                  useyelim &= clipeyelim(rrng0, &eyelim);
152          }
153          for (i = 0; i < p->nr; i++) {
154          retry:
155                  if (useyelim) {
156 +                        initeyelim(&eyelim, NULL, gc+1);
157                          p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
158                                                  + rrng0[0][0];
159                          p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
160                                                  + rrng0[1][0];
161 <                        initeyelim(&eyelim, p->hd, gc+1);
162 <                        gp[gc[0].w>>1] = gc[0].w&1 ?
163 <                                        hdlist[p->hd]->grid[gc[0].w>>1] : 0;
154 <                        gp[hdwg0[gc[0].w]] = gc[0].i[0] +
155 <                                        (1./256.)*(p->ra[i].r[0][0]+.5);
156 <                        gp[hdwg1[gc[0].w]] = gc[0].i[1] +
157 <                                        (1./256.)*(p->ra[i].r[0][1]+.5);
158 <                        groweyelim(&eyelim, gp);
161 >                        groweyelim(&eyelim, gc,
162 >                                        (1./256.)*(p->ra[i].r[0][0]+.5),
163 >                                        (1./256.)*(p->ra[i].r[0][1]+.5));
164                          if (!clipeyelim(rrng1, &eyelim)) {
165                                  useyelim &= nretries-- > 0;
161 #ifdef DEBUG
162                                if (!useyelim)
163                                        error(WARNING, "exceeded retry limit in packrays");
164 #endif
166                                  goto retry;
167                          }
168                          p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
# Line 175 | Line 176 | register PACKET        *p;
176                          p->ra[i].r[1][1] = frandom() * 256.;
177                  }
178                  d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
178 #ifdef DEBUG
179                VSUM(vt, ro, rd, d);
180                dist2sum += dist2line(myeye.vpt, ro, vt);
181 #endif
179                  if (p->offset != NULL) {
180                          if (!vdef(OBSTRUCTIONS))
181                                  d *= frandom();         /* random offset */
# Line 190 | Line 187 | register PACKET        *p;
187                  VCOPY(rod, rd);
188                  rod += 3;
189          }
193 #ifdef DEBUG
194        fprintf(stderr, "RMS distance = %f\n", sqrt(dist2sum/p->nr));
195 #endif
196 #undef gp
190   }
191  
192  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines