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.25 by greg, Mon Jul 7 17:21:51 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   * Rtrace support routines for holodeck rendering
6   */
# Line 19 | Line 16 | struct gclim {
16          HOLO    *hp;                    /* holodeck pointer */
17          GCOORD  gc;                     /* grid cell */
18          FVECT   egp;                    /* eye grid point */
19 <        double  erg[2];                 /* eye range in wall grid coords */
19 >        double  erg2;                   /* mean square eye grid range */
20          double  gmin[2], gmax[2];       /* grid coordinate limits */
21   };                              /* a grid coordinate range */
22  
23  
24   static
25 < initeyelim(gcl, hd, gc)         /* initialize grid coordinate limits */
25 > initeyelim(gcl, hp, gc)         /* initialize grid coordinate limits */
26   register struct gclim   *gcl;
27 < int     hd;
27 > register HOLO   *hp;
28   GCOORD  *gc;
29   {
30 <        register FLOAT  *v;
30 >        register RREAL  *v;
31          register int    i;
32  
33 <        gcl->hp = hdlist[hd];
34 <        copystruct(&gcl->gc, gc);
35 <        hdgrid(gcl->egp, gcl->hp, myeye.vpt);
36 <        for (i = 0; i < 2; i++) {
37 <                v = gcl->hp->wg[((gcl->gc.w>>1)+i+1)%3];
38 <                gcl->erg[i] = myeye.rng * VLEN(v);
42 <                gcl->gmin[i] = FHUGE; gcl->gmax[i] = -FHUGE;
33 >        if (hp != NULL) {
34 >                hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
35 >                gcl->erg2 = 0;
36 >                for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
37 >                        gcl->erg2 += DOT(v,v);
38 >                gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
39          }
40 +        if (gc != NULL)
41 +                copystruct(&gcl->gc, gc);
42 +        gcl->gmin[0] = gcl->gmin[1] = FHUGE;
43 +        gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
44   }
45  
46  
47   static
48 < groweyelim(gcl, gp)             /* grow grid limits about eye point */
48 > groweyelim(gcl, gc, r0, r1, tight)      /* grow grid limits about eye point */
49   register struct gclim   *gcl;
50 < FVECT   gp;
50 > GCOORD  *gc;
51 > double  r0, r1;
52 > int     tight;
53   {
54 <        FVECT   ab;
55 <        double  l2, d, mult, wg;
56 <        register int    i, g;
57 <
54 >        FVECT   gp, ab;
55 >        double  ab2, od, cfact;
56 >        double  sqcoef[3], ctcoef[3], licoef[3], cnst;
57 >        int     gw, gi[2];
58 >        double  wallpos, a, b, c, d, e, f;
59 >        double  root[2], yex;
60 >        int     n, i, j, nex;
61 >                                                /* point/view cone */
62 >        i = gc->w>>1;
63 >        gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
64 >        gp[hdwg0[gc->w]] = gc->i[0] + r0;
65 >        gp[hdwg1[gc->w]] = gc->i[1] + r1;
66          VSUB(ab, gcl->egp, gp);
67 <        l2 = DOT(ab,ab);
68 <        if (l2 <= gcl->erg[0]*gcl->erg[1]) {
67 >        ab2 = DOT(ab, ab);
68 >        gw = gcl->gc.w>>1;
69 >        if ((i==gw ? ab[gw]*ab[gw] : ab2)  <= gcl->erg2 + FTINY) {
70                  gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
71                  gcl->gmax[0] = gcl->gmax[1] = FHUGE;
72 <                return;
72 >                return;                 /* too close (to wall) */
73          }
74 <        mult = gp[g = gcl->gc.w>>1];
75 <        if (gcl->gc.w&1)
76 <                mult -= gcl->hp->grid[g];
77 <        if (ab[g]*ab[g] > gcl->erg[0]*gcl->erg[1])
78 <                mult /= -ab[g];
79 <        else if (fabs(ab[hdwg0[gcl->gc.w]]) > fabs(ab[hdwg1[gcl->gc.w]]))
80 <                mult = (gcl->gc.i[0] + .5 - gp[hdwg0[gcl->gc.w]]) /
81 <                                ab[hdwg0[gcl->gc.w]];
82 <        else
83 <                mult = (gcl->gc.i[1] + .5 - gp[hdwg1[gcl->gc.w]]) /
84 <                                ab[hdwg1[gcl->gc.w]];
74 >        ab2 = 1./ab2;                           /* 1/norm2(ab) */
75 >        od = DOT(gp, ab);                       /* origin dot direction */
76 >        cfact = 1./(1. - ab2*gcl->erg2);        /* tan^2 + 1 of cone angle */
77 >        for (i = 0; i < 3; i++) {               /* compute cone equation */
78 >                sqcoef[i] = ab[i]*ab[i]*cfact*ab2 - 1.;
79 >                ctcoef[i] = 2.*ab[i]*ab[(i+1)%3]*cfact*ab2;
80 >                licoef[i] = 2.*(gp[i] - ab[i]*cfact*od*ab2);
81 >        }
82 >        cnst = cfact*od*od*ab2 - DOT(gp,gp);
83 >        /*
84 >         * CONE:        sqcoef[0]*x*x + sqcoef[1]*y*y + sqcoef[2]*z*z
85 >         *              + ctcoef[0]*x*y + ctcoef[1]*y*z + ctcoef[2]*z*x
86 >         *              + licoef[0]*x + licoef[1]*y + licoef[2]*z + cnst == 0
87 >         */
88 >                                /* equation for conic section in plane */
89 >        gi[0] = hdwg0[gcl->gc.w];
90 >        gi[1] = hdwg1[gcl->gc.w];
91 >        wallpos = gcl->gc.w&1 ? gcl->hp->grid[gw] : 0;
92 >        a = sqcoef[gi[0]];                                      /* x2 */
93 >        b = ctcoef[gi[0]];                                      /* xy */
94 >        c = sqcoef[gi[1]];                                      /* y2 */
95 >        d = ctcoef[gw]*wallpos + licoef[gi[0]];                 /* x */
96 >        e = ctcoef[gi[1]]*wallpos + licoef[gi[1]];              /* y */
97 >        f = wallpos*(wallpos*sqcoef[gw] + licoef[gw]) + cnst;
98          for (i = 0; i < 2; i++) {
99 <                g = ((gcl->gc.w>>1)+i+1)%3;
100 <                wg = gp[g] + mult*ab[g];
101 <                d = mult*gcl->erg[i];
102 <                if (d < 0.) d = -d;
103 <                if (wg - d < gcl->gmin[i])
104 <                        gcl->gmin[i] = wg - d;
105 <                if (wg + d > gcl->gmax[i])
106 <                        gcl->gmax[i] = wg + d;
99 >                if (i) {                /* swap x and y coefficients */
100 >                        register double t;
101 >                        t = a; a = c; c = t;
102 >                        t = d; d = e; e = t;
103 >                }
104 >                nex = 0;                /* check global extrema */
105 >                n = quadratic(root, a*(4.*a*c-b*b), 2.*a*(2.*c*d-b*e),
106 >                                d*(c*d-b*e) + f*b*b);
107 >                while (n-- > 0) {
108 >                        if (gc->w>>1 == gi[i] &&
109 >                                        (gc->w&1) ^ root[n] < gp[gc->w>>1]) {
110 >                                if (gc->w&1)
111 >                                        gcl->gmin[i] = -FHUGE;
112 >                                else
113 >                                        gcl->gmax[i] = FHUGE;
114 >                                nex++;
115 >                                continue;               /* hyperbolic */
116 >                        }
117 >                        if (tight) {
118 >                                yex = (-2.*a*root[n] - d)/b;
119 >                                if (yex < gcl->gc.i[1-i] ||
120 >                                                yex > gcl->gc.i[1-i]+1)
121 >                                        continue;       /* outside cell */
122 >                        }
123 >                        if (root[n] < gcl->gmin[i])
124 >                                gcl->gmin[i] = root[n];
125 >                        if (root[n] > gcl->gmax[i])
126 >                                gcl->gmax[i] = root[n];
127 >                        nex++;
128 >                }
129 >                                        /* check local extrema */
130 >                for (j = nex < 2 ? 2 : 0; j--; ) {
131 >                        yex = gcl->gc.i[1-i] + j;
132 >                        n = quadratic(root, a, b*yex+d, yex*(yex*c+e)+f);
133 >                        while (n-- > 0) {
134 >                                if (gc->w>>1 == gi[i] &&
135 >                                        (gc->w&1) ^ root[n] < gp[gc->w>>1])
136 >                                        continue;
137 >                                if (root[n] < gcl->gmin[i])
138 >                                        gcl->gmin[i] = root[n];
139 >                                if (root[n] > gcl->gmax[i])
140 >                                        gcl->gmax[i] = root[n];
141 >                        }
142 >                }
143          }
144   }
145  
# Line 97 | Line 157 | register struct gclim  *gcl;
157                          gcl->gmin[i] = gcl->gc.i[i];
158                  if (gcl->gmax[i] > gcl->gc.i[i]+1)
159                          gcl->gmax[i] = gcl->gc.i[i]+1;
160 <                if ((incell &= gcl->gmax[i] > gcl->gmin[i])) {
160 >                if (gcl->gmax[i] > gcl->gmin[i]) {
161                          rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
162                                          (1.-FTINY);
163                          rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
164                                          (1.-FTINY) - rrng[i][0];
165 <                        incell &= rrng[i][1] > 0;
166 <                }
165 >                } else
166 >                        rrng[i][0] = rrng[i][1] = 0;
167 >                incell &= rrng[i][1] > 0;
168          }
169          return(incell);
170   }
# Line 113 | Line 174 | packrays(rod, p)               /* pack ray origins and directions *
174   register float  *rod;
175   register PACKET *p;
176   {
177 < #define gp      ro
178 < #ifdef DEBUG
118 <        double dist2sum = 0.;
177 > #if 0
178 >        double  dist2sum = 0.;
179          FVECT   vt;
180   #endif
181          int     nretries = p->nr + 2;
# Line 130 | Line 190 | register PACKET        *p;
190          if (!hdbcoord(gc, hdlist[p->hd], p->bi))
191                  error(CONSISTENCY, "bad beam index in packrays");
192          if ((useyelim = myeye.rng > FTINY)) {
193 <                initeyelim(&eyelim, p->hd, gc);
194 <                gp[gc[1].w>>1] = gc[1].w&1 ?
195 <                                hdlist[p->hd]->grid[gc[1].w>>1] : 0;
196 <                gp[hdwg0[gc[1].w]] = gc[1].i[0];
197 <                gp[hdwg1[gc[1].w]] = gc[1].i[1];
198 <                groweyelim(&eyelim, gp);
199 <                gp[hdwg0[gc[1].w]]++;
200 <                gp[hdwg1[gc[1].w]]++;
141 <                groweyelim(&eyelim, gp);
142 <                useyelim &= clipeyelim(rrng0, &eyelim);
193 >                initeyelim(&eyelim, hdlist[p->hd], gc);
194 >                groweyelim(&eyelim, gc+1, 0., 0., 0);
195 >                groweyelim(&eyelim, gc+1, 1., 1., 0);
196 >                useyelim = clipeyelim(rrng0, &eyelim);
197 > #ifdef DEBUG
198 >                if (!useyelim)
199 >                        error(WARNING, "no eye overlap in packrays");
200 > #endif
201          }
202          for (i = 0; i < p->nr; i++) {
203          retry:
204                  if (useyelim) {
205 +                        initeyelim(&eyelim, NULL, gc+1);
206                          p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
207                                                  + rrng0[0][0];
208                          p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
209                                                  + rrng0[1][0];
210 <                        initeyelim(&eyelim, p->hd, gc+1);
211 <                        gp[gc[0].w>>1] = gc[0].w&1 ?
212 <                                        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);
210 >                        groweyelim(&eyelim, gc,
211 >                                        (1./256.)*(p->ra[i].r[0][0]+.5),
212 >                                        (1./256.)*(p->ra[i].r[0][1]+.5), 1);
213                          if (!clipeyelim(rrng1, &eyelim)) {
214 <                                useyelim &= nretries-- > 0;
214 >                                useyelim = nretries-- > 0;
215   #ifdef DEBUG
216                                  if (!useyelim)
217 <                                        error(WARNING, "exceeded retry limit in packrays");
217 >                                        error(WARNING,
218 >                                        "exceeded retry limit in packrays");
219   #endif
220                                  goto retry;
221                          }
# Line 175 | Line 230 | register PACKET        *p;
230                          p->ra[i].r[1][1] = frandom() * 256.;
231                  }
232                  d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
233 < #ifdef DEBUG
233 > #if 0
234                  VSUM(vt, ro, rd, d);
235                  dist2sum += dist2line(myeye.vpt, ro, vt);
236   #endif
# Line 190 | Line 245 | register PACKET        *p;
245                  VCOPY(rod, rd);
246                  rod += 3;
247          }
248 < #ifdef DEBUG
249 <        fprintf(stderr, "RMS distance = %f\n", sqrt(dist2sum/p->nr));
248 > #if 0
249 >        fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
250 >                        p->nr + 2 - nretries);
251   #endif
196 – #undef gp
252   }
253  
254  
# Line 222 | Line 277 | done_rtrace()                  /* clean up and close rtrace calculati
277          int     status;
278                                          /* already closed? */
279          if (!nprocs)
280 <                return;
280 >                return(0);
281                                          /* flush beam queue */
282          done_packets(flush_queue());
283                                          /* sync holodeck */
# Line 293 | Line 348 | getradfile()                   /* run rad and get needed variables */
348                          pippt = NULL;
349                  }
350          if (pippt != NULL)
351 <                strcpy(pippt, "> /dev/null");   /* nothing to match */
351 >                strcpy(pippt, "> " NULL_DEVICE);        /* nothing to match */
352          else
353                  sprintf(cp, ")[ \t]*=' > %s", tf2);
354   #ifdef DEBUG

Diff Legend

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