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.20 by gwlarson, Thu Dec 3 15:21:05 1998 UTC vs.
Revision 3.21 by gwlarson, Mon Dec 7 16:56:08 1998 UTC

# Line 48 | Line 48 | GCOORD *gc;
48  
49  
50   static
51 < groweyelim(gcl, gc, r0, r1)     /* grow grid limits about eye point */
51 > groweyelim(gcl, gc, r0, r1, tight)      /* grow grid limits about eye point */
52   register struct gclim   *gcl;
53   GCOORD  *gc;
54   double  r0, r1;
55 + int     tight;
56   {
57          FVECT   gp, ab;
58 <        double  vlen, plen, dv0, dv1;
59 <        double  rd2, dwall, gpos;
60 <        int     eyeout;
61 <        register int    i, g0, g1;
62 <
58 >        double  ab2, od, cfact;
59 >        double  sqcoef[3], ctcoef[3], licoef[3], cnst;
60 >        int     gw, gi[2];
61 >        double  wallpos, a, b, c, d, e, f;
62 >        double  root[2], yex;
63 >        int     n, i, j, nex;
64 >                                                /* point/view cone */
65          i = gc->w>>1;
66 <        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];
66 >        gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
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 <        rd2 = DOT(ab,ab);
71 <        if (rd2 <= gcl->erg2) {
70 >        ab2 = DOT(ab, ab);
71 >        gw = gcl->gc.w>>1;
72 >        if ((i==gw ? ab[gw]*ab[gw] : ab2)  <= gcl->erg2 + FTINY) {
73                  gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
74                  gcl->gmax[0] = gcl->gmax[1] = FHUGE;
75 <                return;
75 >                return;                 /* too close (to wall) */
76          }
77 <        rd2 = gcl->erg2 / rd2;
78 <        vlen = 1. - rd2;
79 <        plen = sqrt(rd2 * vlen);
80 <        g0 = gcl->gc.w>>1;
81 <        dwall = (gcl->gc.w&1 ? gcl->hp->grid[g0] : 0) - gp[g0];
82 <        for (i = 0; i < 4; i++) {
83 <                if (i == 2)
84 <                        plen = -plen;
85 <                g1 = (g0+(i&1)+1)%3;
86 <                dv0 = vlen*ab[g0] + plen*ab[g1];
87 <                dv1 = vlen*ab[g1] - plen*ab[g0];
88 <                if ((dv0 < 0 ^ dwall < 0 ^ eyeout) ||
89 <                                (dv0 <= FTINY && dv0 >= -FTINY)) {
90 <                        if (eyeout)
91 <                                dv1 = -dv1;
92 <                        if (dv1 > FTINY)
93 <                                gcl->gmax[i&1] = FHUGE;
94 <                        else if (dv1 < -FTINY)
95 <                                gcl->gmin[i&1] = -FHUGE;
96 <                } else {
97 <                        gpos = gp[g1] + dv1*dwall/dv0;
98 <                        if (gpos < gcl->gmin[i&1])
99 <                                gcl->gmin[i&1] = gpos;
100 <                        if (gpos > gcl->gmax[i&1])
101 <                                gcl->gmax[i&1] = gpos;
77 >        ab2 = 1./ab2;                           /* 1/norm2(ab) */
78 >        od = DOT(gp, ab);                       /* origin dot direction */
79 >        cfact = 1./(1. - ab2*gcl->erg2);        /* tan^2 + 1 of cone angle */
80 >        for (i = 0; i < 3; i++) {               /* compute cone equation */
81 >                sqcoef[i] = ab[i]*ab[i]*cfact*ab2 - 1.;
82 >                ctcoef[i] = 2.*ab[i]*ab[(i+1)%3]*cfact*ab2;
83 >                licoef[i] = 2.*(gp[i] - ab[i]*cfact*od*ab2);
84 >        }
85 >        cnst = cfact*od*od*ab2 - DOT(gp,gp);
86 >        /*
87 >         * CONE:        sqcoef[0]*x*x + sqcoef[1]*y*y + sqcoef[2]*z*z
88 >         *              + ctcoef[0]*x*y + ctcoef[1]*y*z + ctcoef[2]*z*x
89 >         *              + licoef[0]*x + licoef[1]*y + licoef[2]*z + cnst == 0
90 >         */
91 >                                /* equation for conic section in plane */
92 >        gi[0] = hdwg0[gcl->gc.w];
93 >        gi[1] = hdwg1[gcl->gc.w];
94 >        wallpos = gcl->gc.w&1 ? gcl->hp->grid[gw] : 0;
95 >        a = sqcoef[gi[0]];                                      /* x2 */
96 >        b = ctcoef[gi[0]];                                      /* xy */
97 >        c = sqcoef[gi[1]];                                      /* y2 */
98 >        d = ctcoef[gw]*wallpos + licoef[gi[0]];                 /* x */
99 >        e = ctcoef[gi[1]]*wallpos + licoef[gi[1]];              /* y */
100 >        f = wallpos*(wallpos*sqcoef[gw] + licoef[gw]) + cnst;
101 >        for (i = 0; i < 2; i++) {
102 >                if (i) {                /* swap x and y coefficients */
103 >                        register double t;
104 >                        t = a; a = c; c = t;
105 >                        t = d; d = e; e = t;
106                  }
107 +                nex = 0;                /* check global extrema */
108 +                n = quadratic(root, a*(4.*a*c-b*b), 2.*a*(2.*c*d-b*e),
109 +                                d*(c*d-b*e) + f*b*b);
110 +                while (n-- > 0) {
111 +                        if (gc->w>>1 == gi[i] &&
112 +                                        (gc->w&1) ^ root[n] < gp[gc->w>>1]) {
113 +                                if (gc->w&1)
114 +                                        gcl->gmin[i] = -FHUGE;
115 +                                else
116 +                                        gcl->gmax[i] = FHUGE;
117 +                                nex++;
118 +                                continue;               /* hyperbolic */
119 +                        }
120 +                        if (tight) {
121 +                                yex = (-2.*a*root[n] - d)/b;
122 +                                if (yex < gcl->gc.i[1-i] ||
123 +                                                yex > gcl->gc.i[1-i]+1)
124 +                                        continue;       /* outside cell */
125 +                        }
126 +                        if (root[n] < gcl->gmin[i])
127 +                                gcl->gmin[i] = root[n];
128 +                        if (root[n] > gcl->gmax[i])
129 +                                gcl->gmax[i] = root[n];
130 +                        nex++;
131 +                }
132 +                                        /* check local extrema */
133 +                for (j = nex < 2 ? 2 : 0; j--; ) {
134 +                        yex = gcl->gc.i[1-i] + j;
135 +                        n = quadratic(root, a, b*yex+d, yex*(yex*c+e)+f);
136 +                        while (n-- > 0) {
137 +                                if (gc->w>>1 == gi[i] &&
138 +                                        (gc->w&1) ^ root[n] < gp[gc->w>>1])
139 +                                        continue;
140 +                                if (root[n] < gcl->gmin[i])
141 +                                        gcl->gmin[i] = root[n];
142 +                                if (root[n] > gcl->gmax[i])
143 +                                        gcl->gmax[i] = root[n];
144 +                        }
145 +                }
146          }
147   }
148  
# Line 133 | Line 177 | packrays(rod, p)               /* pack ray origins and directions *
177   register float  *rod;
178   register PACKET *p;
179   {
180 + #if 0
181 +        double  dist2sum = 0.;
182 +        FVECT   vt;
183 + #endif
184          int     nretries = p->nr + 2;
185          struct gclim    eyelim;
186          short   rrng0[2][2], rrng1[2][2];
# Line 146 | Line 194 | register PACKET        *p;
194                  error(CONSISTENCY, "bad beam index in packrays");
195          if ((useyelim = myeye.rng > FTINY)) {
196                  initeyelim(&eyelim, hdlist[p->hd], gc);
197 <                groweyelim(&eyelim, gc+1, 0., 0.);
198 <                groweyelim(&eyelim, gc+1, 1., 1.);
199 <                useyelim &= clipeyelim(rrng0, &eyelim);
197 >                groweyelim(&eyelim, gc+1, 0., 0., 0);
198 >                groweyelim(&eyelim, gc+1, 1., 1., 0);
199 >                useyelim = clipeyelim(rrng0, &eyelim);
200 > #ifdef DEBUG
201 >                if (!useyelim)
202 >                        error(WARNING, "no eye overlap in packrays");
203 > #endif
204          }
205          for (i = 0; i < p->nr; i++) {
206          retry:
# Line 160 | Line 212 | register PACKET        *p;
212                                                  + rrng0[1][0];
213                          groweyelim(&eyelim, gc,
214                                          (1./256.)*(p->ra[i].r[0][0]+.5),
215 <                                        (1./256.)*(p->ra[i].r[0][1]+.5));
215 >                                        (1./256.)*(p->ra[i].r[0][1]+.5), 1);
216                          if (!clipeyelim(rrng1, &eyelim)) {
217 <                                useyelim &= nretries-- > 0;
217 >                                useyelim = nretries-- > 0;
218 > #ifdef DEBUG
219 >                                if (!useyelim)
220 >                                        error(WARNING,
221 >                                        "exceeded retry limit in packrays");
222 > #endif
223                                  goto retry;
224                          }
225                          p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
# Line 176 | Line 233 | register PACKET        *p;
233                          p->ra[i].r[1][1] = frandom() * 256.;
234                  }
235                  d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
236 + #if 0
237 +                VSUM(vt, ro, rd, d);
238 +                dist2sum += dist2line(myeye.vpt, ro, vt);
239 + #endif
240                  if (p->offset != NULL) {
241                          if (!vdef(OBSTRUCTIONS))
242                                  d *= frandom();         /* random offset */
# Line 187 | Line 248 | register PACKET        *p;
248                  VCOPY(rod, rd);
249                  rod += 3;
250          }
251 + #if 0
252 +        fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
253 +                        p->nr + 2 - nretries);
254 + #endif
255   }
256  
257  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines