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

Comparing ray/src/hd/holo.c (file contents):
Revision 3.6 by gregl, Tue Nov 11 17:03:28 1997 UTC vs.
Revision 3.21 by schorsch, Thu Jan 1 11:21:55 2004 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1997 Silicon Graphics, Inc. */
2
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6
4   /*
5   * Routines for converting holodeck coordinates, etc.
6   *
# Line 14 | Line 11 | static char SCCSid[] = "$SunId$ SGI";
11  
12   float   hd_depthmap[DCINF-DCLIN];
13  
14 + int     hdwg0[6] = {1,1,2,2,0,0};
15 + int     hdwg1[6] = {2,2,0,0,1,1};
16 +
17   static double   logstep;
18  
19 static int      wg0[6] = {1,1,2,2,0,0};
20 static int      wg1[6] = {2,2,0,0,1,1};
19  
20 <
21 < hdcompgrid(hp)                  /* compute derived grid vector and index */
22 < register HOLO   *hp;
20 > extern void
21 > hdcompgrid(                     /* compute derived grid vector and index */
22 >        register HOLO   *hp
23 > )
24   {
25          double  d;
26          register int    i, j;
# Line 36 | Line 35 | register HOLO  *hp;
35          }
36                                  /* compute grid coordinate vectors */
37          for (i = 0; i < 3; i++) {
38 <                fcross(hp->wn[i], hp->xv[(i+1)%3], hp->xv[(i+2)%3]);
39 <                if (normalize(hp->wn[i]) == 0.)
38 >                fcross(hp->wg[i], hp->xv[(i+1)%3], hp->xv[(i+2)%3]);
39 >                d = DOT(hp->wg[i],hp->xv[i]);
40 >                if ((d <= FTINY) & (d >= -FTINY))
41                          error(USER, "degenerate holodeck section");
42 <                hp->wo[i<<1] = DOT(hp->wn[i],hp->orig);
43 <                d = DOT(hp->wn[i],hp->xv[i]);
44 <                hp->wo[i<<1|1] = hp->wo[i<<1] + d;
45 <                hp->wg[i] = (double)hp->grid[i] / d;
42 >                d = hp->grid[i] / d;
43 >                hp->wg[i][0] *= d; hp->wg[i][1] *= d; hp->wg[i][2] *= d;
44          }
45                                  /* compute linear depth range */
46          hp->tlin = VLEN(hp->xv[0]) + VLEN(hp->xv[1]) + VLEN(hp->xv[2]);
# Line 51 | Line 49 | register HOLO  *hp;
49          for (i = 1; i < 6; i++) {
50                  hp->wi[i] = 0;
51                  for (j = i; j < 6; j++)
52 <                        hp->wi[i] += hp->grid[wg0[j]] * hp->grid[wg1[j]];
53 <                hp->wi[i] *= hp->grid[wg0[i-1]] * hp->grid[wg1[i-1]];
52 >                        hp->wi[i] += hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
53 >                hp->wi[i] *= hp->grid[hdwg0[i-1]] * hp->grid[hdwg1[i-1]];
54                  hp->wi[i] += hp->wi[i-1];
55          }
56   }
57  
58  
59 < HOLO *
60 < hdalloc(hproto)         /* allocate and set holodeck section based on grid */
61 < HDGRID  *hproto;
59 > extern int
60 > hdbcoord(               /* compute beam coordinates from index */
61 >        GCOORD  gc[2],          /* returned */
62 >        register HOLO   *hp,
63 >        register int    i
64 > )
65   {
65        HOLO    hdhead;
66        register HOLO   *hp;
67        int     n;
68                                /* copy grid to temporary header */
69        bcopy((char *)hproto, (char *)&hdhead, sizeof(HDGRID));
70                                /* compute grid vectors and sizes */
71        hdcompgrid(&hdhead);
72                                /* allocate header with directory */
73        n = sizeof(HOLO)+nbeams(&hdhead)*sizeof(BEAMI);
74        if ((hp = (HOLO *)malloc(n)) == NULL)
75                return(NULL);
76                                /* copy header information */
77        copystruct(hp, &hdhead);
78                                /* allocate and clear beam list */
79        hp->bl = (BEAM **)malloc((nbeams(hp)+1)*sizeof(BEAM *)+sizeof(BEAM));
80        if (hp->bl == NULL) {
81                free((char *)hp);
82                return(NULL);
83        }
84        bzero((char *)hp->bl, (nbeams(hp)+1)*sizeof(BEAM *)+sizeof(BEAM));
85        hp->bl[0] = (BEAM *)(hp->bl+nbeams(hp)+1);      /* set blglob(hp) */
86        hp->fd = -1;
87        hp->dirty = 0;
88        hp->priv = NULL;
89                                /* clear beam directory */
90        bzero((char *)hp->bi, (nbeams(hp)+1)*sizeof(BEAMI));
91        return(hp);             /* all is well */
92 }
93
94
95 hdbcoord(gc, hp, i)             /* compute beam coordinates from index */
96 GCOORD  gc[2];          /* returned */
97 register HOLO   *hp;
98 register int    i;
99 {
66          register int    j, n;
67          int     n2, reverse;
68          GCOORD  g2[2];
69                                          /* check range */
70 <        if (i < 1 | i > nbeams(hp))
70 >        if ((i < 1) | (i > nbeams(hp)))
71                  return(0);
72 <        if (reverse = i >= hp->wi[5])
72 >        if ( (reverse = i >= hp->wi[5]) )
73                  i -= hp->wi[5] - 1;
74          for (j = 0; j < 5; j++)         /* find w0 */
75                  if (hp->wi[j+1] > i)
76                          break;
77          i -= hp->wi[gc[0].w=j];
78                                          /* find w1 */
79 <        n2 = hp->grid[wg0[j]] * hp->grid[wg1[j]];
79 >        n2 = hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
80          while (++j < 5) {
81 <                n = n2 * hp->grid[wg0[j]] * hp->grid[wg1[j]];
81 >                n = n2 * hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
82                  if (n > i)
83                          break;
84                  i -= n;
85          }
86          gc[1].w = j;
87                                          /* find position on w0 */
88 <        n2 = hp->grid[wg0[j]] * hp->grid[wg1[j]];
88 >        n2 = hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
89          n = i / n2;
90 <        gc[0].i[1] = n / hp->grid[wg0[gc[0].w]];
91 <        gc[0].i[0] = n - gc[0].i[1]*hp->grid[wg0[gc[0].w]];
90 >        gc[0].i[1] = n / hp->grid[hdwg0[gc[0].w]];
91 >        gc[0].i[0] = n - gc[0].i[1]*hp->grid[hdwg0[gc[0].w]];
92          i -= n*n2;
93                                          /* find position on w1 */
94 <        gc[1].i[1] = i / hp->grid[wg0[gc[1].w]];
95 <        gc[1].i[0] = i - gc[1].i[1]*hp->grid[wg0[gc[1].w]];
94 >        gc[1].i[1] = i / hp->grid[hdwg0[gc[1].w]];
95 >        gc[1].i[0] = i - gc[1].i[1]*hp->grid[hdwg0[gc[1].w]];
96          if (reverse) {
97 <                copystruct(g2, gc+1);
98 <                copystruct(gc+1, gc);
99 <                copystruct(gc, g2);
97 >                *g2 = *(gc+1);
98 >                *(gc+1) = *gc;
99 >                *gc = *g2;
100          }
101          return(1);                      /* we're done */
102   }
103  
104  
105 < int
106 < hdbindex(hp, gc)                /* compute index from beam coordinates */
107 < register HOLO   *hp;
108 < register GCOORD gc[2];
105 > extern int
106 > hdbindex(               /* compute index from beam coordinates */
107 >        register HOLO   *hp,
108 >        register GCOORD gc[2]
109 > )
110   {
111          GCOORD  g2[2];
112          int     reverse;
113          register int    i, j;
114                                          /* check ordering and limits */
115 <        if (reverse = gc[0].w > gc[1].w) {
116 <                copystruct(g2, gc+1);
117 <                copystruct(g2+1, gc);
115 >        if ( (reverse = gc[0].w > gc[1].w) ) {
116 >                *g2 = *(gc+1);
117 >                *(g2+1) = *gc;
118                  gc = g2;
119          } else if (gc[0].w == gc[1].w)
120                  return(0);
121 <        if (gc[0].w < 0 | gc[1].w > 5)
121 >        if ((gc[0].w < 0) | (gc[1].w > 5))
122                  return(0);
123          i = 0;                          /* compute index */
124          for (j = gc[0].w+1; j < gc[1].w; j++)
125 <                i += hp->grid[wg0[j]] * hp->grid[wg1[j]];
126 <        i *= hp->grid[wg0[gc[0].w]] * hp->grid[wg1[gc[0].w]];
125 >                i += hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
126 >        i *= hp->grid[hdwg0[gc[0].w]] * hp->grid[hdwg1[gc[0].w]];
127          i += hp->wi[gc[0].w];
128 <        i += (hp->grid[wg0[gc[0].w]]*gc[0].i[1] + gc[0].i[0]) *
129 <                        hp->grid[wg0[gc[1].w]] * hp->grid[wg1[gc[1].w]] ;
130 <        i += hp->grid[wg0[gc[1].w]]*gc[1].i[1] + gc[1].i[0];
128 >        i += (hp->grid[hdwg0[gc[0].w]]*gc[0].i[1] + gc[0].i[0]) *
129 >                        hp->grid[hdwg0[gc[1].w]] * hp->grid[hdwg1[gc[1].w]] ;
130 >        i += hp->grid[hdwg0[gc[1].w]]*gc[1].i[1] + gc[1].i[0];
131          if (reverse)
132                  i += hp->wi[5] - 1;
133          return(i);
134   }
135  
136  
137 < hdcell(cp, hp, gc)              /* compute cell coordinates */
138 < register FVECT  cp[4];  /* returned (may be passed as FVECT cp[2][2]) */
139 < register HOLO   *hp;
140 < register GCOORD *gc;
137 > extern void
138 > hdcell(         /* compute cell coordinates */
139 >        register FVECT  cp[4],  /* returned (may be passed as FVECT cp[2][2]) */
140 >        register HOLO   *hp,
141 >        register GCOORD *gc
142 > )
143   {
144 <        register FLOAT  *v;
144 >        register RREAL  *v;
145          double  d;
146                                          /* compute common component */
147          VCOPY(cp[0], hp->orig);
# Line 180 | Line 149 | register GCOORD        *gc;
149                  v = hp->xv[gc->w>>1];
150                  cp[0][0] += v[0]; cp[0][1] += v[1]; cp[0][2] += v[2];
151          }
152 <        v = hp->xv[wg0[gc->w]];
153 <        d = (double)gc->i[0] / hp->grid[wg0[gc->w]];
152 >        v = hp->xv[hdwg0[gc->w]];
153 >        d = (double)gc->i[0] / hp->grid[hdwg0[gc->w]];
154          VSUM(cp[0], cp[0], v, d);
155 <        v = hp->xv[wg1[gc->w]];
156 <        d = (double)gc->i[1] / hp->grid[wg1[gc->w]];
155 >        v = hp->xv[hdwg1[gc->w]];
156 >        d = (double)gc->i[1] / hp->grid[hdwg1[gc->w]];
157          VSUM(cp[0], cp[0], v, d);
158                                          /* compute x1 sums */
159 <        v = hp->xv[wg0[gc->w]];
160 <        d = 1.0 / hp->grid[wg0[gc->w]];
159 >        v = hp->xv[hdwg0[gc->w]];
160 >        d = 1.0 / hp->grid[hdwg0[gc->w]];
161          VSUM(cp[1], cp[0], v, d);
162          VSUM(cp[3], cp[0], v, d);
163                                          /* compute y1 sums */
164 <        v = hp->xv[wg1[gc->w]];
165 <        d = 1.0 / hp->grid[wg1[gc->w]];
164 >        v = hp->xv[hdwg1[gc->w]];
165 >        d = 1.0 / hp->grid[hdwg1[gc->w]];
166          VSUM(cp[2], cp[0], v, d);
167          VSUM(cp[3], cp[3], v, d);
168   }
169  
170  
171 < hdlseg(lseg, hp, i)                     /* compute line segment for beam */
172 < register int    lseg[2][3];
173 < register HOLO   *hp;
174 < int     i;
171 > extern int
172 > hdlseg(                 /* compute line segment for beam */
173 >        register int    lseg[2][3],
174 >        register HOLO   *hp,
175 >        GCOORD  gc[2]
176 > )
177   {
207        GCOORD  gc[2];
178          register int    k;
179  
210        if (!hdbcoord(gc, hp, i))               /* compute grid coordinates */
211                return(0);
180          for (k = 0; k < 2; k++) {               /* compute end points */
181                  lseg[k][gc[k].w>>1] = gc[k].w&1 ? hp->grid[gc[k].w>>1]-1 : 0 ;
182 <                lseg[k][wg0[gc[k].w]] = gc[k].i[0];
183 <                lseg[k][wg1[gc[k].w]] = gc[k].i[1];
182 >                lseg[k][hdwg0[gc[k].w]] = gc[k].i[0];
183 >                lseg[k][hdwg1[gc[k].w]] = gc[k].i[1];
184          }
185          return(1);
186   }
187  
188  
189 < unsigned
190 < hdcode(hp, d)                   /* compute depth code for d */
191 < HOLO    *hp;
192 < double  d;
189 > extern unsigned int
190 > hdcode(                 /* compute depth code for d */
191 >        HOLO    *hp,
192 >        double  d
193 > )
194   {
195          double  tl = hp->tlin;
196 <        register unsigned       c;
196 >        register long   c;
197  
198          if (d <= 0.)
199                  return(0);
# Line 232 | Line 201 | double d;
201                  return(DCINF);
202          if (d < tl)
203                  return((unsigned)(d*DCLIN/tl));
204 <        c = (unsigned)(log(d/tl)/logstep) + DCLIN;
205 <        return(c > DCINF ? DCINF : c);
204 >        c = (long)(log(d/tl)/logstep) + DCLIN;
205 >        return(c > DCINF ? (unsigned)DCINF : (unsigned)c);
206   }
207  
208  
209 < hdgrid(gp, hp, wp)              /* compute grid coordinates */
210 < FVECT   gp;             /* returned */
211 < register HOLO   *hp;
212 < FVECT   wp;
209 > extern void
210 > hdgrid(         /* compute grid coordinates */
211 >        FVECT   gp,             /* returned */
212 >        register HOLO   *hp,
213 >        FVECT   wp
214 > )
215   {
216          FVECT   vt;
217  
218 <        vt[0] = wp[0] - hp->orig[0];
219 <        vt[1] = wp[1] - hp->orig[1];
220 <        vt[2] = wp[2] - hp->orig[2];
221 <        gp[0] = DOT(vt, hp->wn[0]) * hp->wg[0];
251 <        gp[1] = DOT(vt, hp->wn[1]) * hp->wg[1];
252 <        gp[2] = DOT(vt, hp->wn[2]) * hp->wg[2];
218 >        VSUB(vt, wp, hp->orig);
219 >        gp[0] = DOT(vt, hp->wg[0]);
220 >        gp[1] = DOT(vt, hp->wg[1]);
221 >        gp[2] = DOT(vt, hp->wg[2]);
222   }
223  
224  
225 < double
226 < hdray(ro, rd, hp, gc, r)        /* compute ray within a beam */
227 < FVECT   ro, rd;         /* returned */
228 < HOLO    *hp;
229 < GCOORD  gc[2];
230 < BYTE    r[2][2];
225 > extern void
226 > hdworld(                /* compute world coordinates */
227 >        register FVECT  wp,
228 >        register HOLO   *hp,
229 >        FVECT   gp
230 > )
231   {
232 +        register double d;
233 +
234 +        d = gp[0]/hp->grid[0];
235 +        VSUM(wp, hp->orig, hp->xv[0], d);
236 +
237 +        d = gp[1]/hp->grid[1];
238 +        VSUM(wp, wp, hp->xv[1], d);
239 +
240 +        d = gp[2]/hp->grid[2];
241 +        VSUM(wp, wp, hp->xv[2], d);
242 + }
243 +
244 +
245 + extern double
246 + hdray(  /* compute ray within a beam */
247 +        FVECT   ro,
248 +        FVECT   rd,             /* returned */
249 +        HOLO    *hp,
250 +        GCOORD  gc[2],
251 +        BYTE    r[2][2]
252 + )
253 + {
254          FVECT   cp[4], p[2];
255          register int    i, j;
256          double  d0, d1;
# Line 273 | Line 264 | BYTE   r[2][2];
264                                          d0*cp[1][j] + d1*cp[2][j];
265          }
266          VCOPY(ro, p[0]);                /* assign ray origin and direction */
267 <        rd[0] = p[1][0] - p[0][0];
277 <        rd[1] = p[1][1] - p[0][1];
278 <        rd[2] = p[1][2] - p[0][2];
267 >        VSUB(rd, p[1], p[0]);
268          return(normalize(rd));          /* return maximum inside distance */
269   }
270  
271  
272 < double
273 < hdinter(gc, r, hp, ro, rd)      /* compute ray intersection with section */
274 < register GCOORD gc[2];  /* returned */
275 < BYTE    r[2][2];        /* returned */
276 < register HOLO   *hp;
277 < FVECT   ro, rd;         /* rd should be normalized */
272 > extern double
273 > hdinter(        /* compute ray intersection with section */
274 >        register GCOORD gc[2],  /* returned */
275 >        BYTE    r[2][2],        /* returned (optional) */
276 >        double  *ed,            /* returned (optional) */
277 >        register HOLO   *hp,
278 >        FVECT   ro,
279 >        FVECT   rd              /* normalization of rd affects distances */
280 > )
281   {
282          FVECT   p[2], vt;
283          double  d, t0, t1, d0, d1;
284 <        register FLOAT  *v;
284 >        register RREAL  *v;
285          register int    i;
286                                          /* first, intersect walls */
287          gc[0].w = gc[1].w = -1;
288          t0 = -FHUGE; t1 = FHUGE;
289 +        VSUB(vt, ro, hp->orig);
290          for (i = 0; i < 3; i++) {               /* for each wall pair */
291 <                d = -DOT(rd, hp->wn[i]);        /* plane distance */
291 >                d = -DOT(rd, hp->wg[i]);        /* plane distance */
292                  if (d <= FTINY && d >= -FTINY)  /* check for parallel */
293                          continue;
294 <                d1 = DOT(ro, hp->wn[i]);        /* ray distances */
295 <                d0 = (d1 - hp->wo[i<<1]) / d;
296 <                d1 = (d1 - hp->wo[i<<1|1]) / d;
297 <                if (d0 < d1) {          /* check against best */
294 >                d1 = DOT(vt, hp->wg[i]);        /* ray distances */
295 >                d0 = d1 / d;
296 >                d1 = (d1 - hp->grid[i]) / d;
297 >                if (d < 0) {                    /* check against best */
298                          if (d0 > t0) {
299                                  t0 = d0;
300                                  gc[0].w = i<<1;
# Line 321 | Line 314 | FVECT  ro, rd;         /* rd should be normalized */
314                          }
315                  }
316          }
317 <        if (gc[0].w < 0 | gc[1].w < 0)          /* paranoid check */
317 >        if ((gc[0].w < 0) | (gc[1].w < 0))              /* paranoid check */
318                  return(FHUGE);
319                                                  /* compute intersections */
320 <        for (i = 0; i < 3; i++) {
321 <                p[0][i] = ro[i] + rd[i]*t0;
329 <                p[1][i] = ro[i] + rd[i]*t1;
330 <        }
320 >        VSUM(p[0], ro, rd, t0);
321 >        VSUM(p[1], ro, rd, t1);
322                                          /* now, compute grid coordinates */
323          for (i = 0; i < 2; i++) {
324 <                vt[0] = p[i][0] - hp->orig[0];
325 <                vt[1] = p[i][1] - hp->orig[1];
326 <                vt[2] = p[i][2] - hp->orig[2];
327 <                v = hp->wn[wg0[gc[i].w]];
337 <                d = DOT(vt, v) * hp->wg[wg0[gc[i].w]];
338 <                if (d < 0. || (gc[i].i[0] = d) >= hp->grid[wg0[gc[i].w]])
324 >                VSUB(vt, p[i], hp->orig);
325 >                v = hp->wg[hdwg0[gc[i].w]];
326 >                d = DOT(vt, v);
327 >                if (d < 0 || d >= hp->grid[hdwg0[gc[i].w]])
328                          return(FHUGE);          /* outside wall */
329 <                r[i][0] = 256. * (d - gc[i].i[0]);
330 <                v = hp->wn[wg1[gc[i].w]];
331 <                d = DOT(vt, v) * hp->wg[wg1[gc[i].w]];
332 <                if (d < 0. || (gc[i].i[1] = d) >= hp->grid[wg1[gc[i].w]])
329 >                gc[i].i[0] = d;
330 >                if (r != NULL)
331 >                        r[i][0] = 256. * (d - gc[i].i[0]);
332 >                v = hp->wg[hdwg1[gc[i].w]];
333 >                d = DOT(vt, v);
334 >                if (d < 0 || d >= hp->grid[hdwg1[gc[i].w]])
335                          return(FHUGE);          /* outside wall */
336 <                r[i][1] = 256. * (d - gc[i].i[1]);
336 >                gc[i].i[1] = d;
337 >                if (r != NULL)
338 >                        r[i][1] = 256. * (d - gc[i].i[1]);
339          }
340 <                                        /* return distance from entry point */
341 <        vt[0] = ro[0] - p[0][0];
342 <        vt[1] = ro[1] - p[0][1];
350 <        vt[2] = ro[2] - p[0][2];
351 <        return(DOT(vt,rd));
340 >        if (ed != NULL)                 /* assign distance to exit point */
341 >                *ed = t1;
342 >        return(t0);                     /* return distance to entry point */
343   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines