src/hd/holo.c

/* Copyright (c) 1997 Silicon Graphics, Inc. */

#ifndef lint
static char SCCSid[] = "$SunId$ SGI";
#endif

/*
 * Routines for converting holodeck coordinates, etc.
 *
 *      10/22/97        GWLarson
 */

#include "holo.h"

float   hd_depthmap[DCINF-DCLIN];

int     hdwg0[6] = {1,1,2,2,0,0};
int     hdwg1[6] = {2,2,0,0,1,1};

static double   logstep;


hdcompgrid(hp)                  /* compute derived grid vector and index */
register HOLO   *hp;
{
        double  d;
        register int    i, j;
                                /* initialize depth map */
        if (hd_depthmap[0] < 1.) {
                d = 1. + .5/DCLIN;
                for (i = 0; i < DCINF-DCLIN; i++) {
                        hd_depthmap[i] = d;
                        d *= 1. + 1./DCLIN;
                }
                logstep = log(1. + 1./DCLIN);
        }
                                /* compute grid coordinate vectors */
        for (i = 0; i < 3; i++) {
                fcross(hp->wg[i], hp->xv[(i+1)%3], hp->xv[(i+2)%3]);
                d = DOT(hp->wg[i],hp->xv[i]);
                if (d <= FTINY & d >= -FTINY)
                        error(USER, "degenerate holodeck section");
                d = (double)hp->grid[i] / d;
                hp->wg[i][0] *= d; hp->wg[i][1] *= d; hp->wg[i][2] *= d;
                hp->wo[i<<1] = DOT(hp->wg[i],hp->orig);
                d = DOT(hp->wg[i],hp->xv[i]);
                hp->wo[i<<1|1] = hp->wo[i<<1] + d;
        }
                                /* compute linear depth range */
        hp->tlin = VLEN(hp->xv[0]) + VLEN(hp->xv[1]) + VLEN(hp->xv[2]);
                                /* compute wall super-indices from grid */
        hp->wi[0] = 1;          /**** index values begin at 1 ****/
        for (i = 1; i < 6; i++) {
                hp->wi[i] = 0;
                for (j = i; j < 6; j++)
                        hp->wi[i] += hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
                hp->wi[i] *= hp->grid[hdwg0[i-1]] * hp->grid[hdwg1[i-1]];
                hp->wi[i] += hp->wi[i-1];
        }
}


HOLO *
hdalloc(hproto)         /* allocate and set holodeck section based on grid */
HDGRID  *hproto;
{
        HOLO    hdhead;
        register HOLO   *hp;
        int     n;
                                /* copy grid to temporary header */
        bcopy((char *)hproto, (char *)&hdhead, sizeof(HDGRID));
                                /* compute grid vectors and sizes */
        hdcompgrid(&hdhead);
                                /* allocate header with directory */
        n = sizeof(HOLO)+nbeams(&hdhead)*sizeof(BEAMI);
        if ((hp = (HOLO *)malloc(n)) == NULL)
                return(NULL);
                                /* copy header information */
        copystruct(hp, &hdhead);
                                /* allocate and clear beam list */
        hp->bl = (BEAM **)malloc((nbeams(hp)+1)*sizeof(BEAM *)+sizeof(BEAM));
        if (hp->bl == NULL) {
                free((char *)hp);
                return(NULL);
        }
        bzero((char *)hp->bl, (nbeams(hp)+1)*sizeof(BEAM *)+sizeof(BEAM));
        hp->bl[0] = (BEAM *)(hp->bl+nbeams(hp)+1);      /* set blglob(hp) */
        hp->fd = -1;
        hp->dirty = 0;
        hp->priv = NULL;
                                /* clear beam directory */
        bzero((char *)hp->bi, (nbeams(hp)+1)*sizeof(BEAMI));
        return(hp);             /* all is well */
}


hdbcoord(gc, hp, i)             /* compute beam coordinates from index */
GCOORD  gc[2];          /* returned */
register HOLO   *hp;
register int    i;
{
        register int    j, n;
        int     n2, reverse;
        GCOORD  g2[2];
                                        /* check range */
        if (i < 1 | i > nbeams(hp))
                return(0);
        if (reverse = i >= hp->wi[5])
                i -= hp->wi[5] - 1;
        for (j = 0; j < 5; j++)         /* find w0 */
                if (hp->wi[j+1] > i)
                        break;
        i -= hp->wi[gc[0].w=j];
                                        /* find w1 */
        n2 = hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
        while (++j < 5) {
                n = n2 * hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
                if (n > i)
                        break;
                i -= n;
        }
        gc[1].w = j;
                                        /* find position on w0 */
        n2 = hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
        n = i / n2;
        gc[0].i[1] = n / hp->grid[hdwg0[gc[0].w]];
        gc[0].i[0] = n - gc[0].i[1]*hp->grid[hdwg0[gc[0].w]];
        i -= n*n2;
                                        /* find position on w1 */
        gc[1].i[1] = i / hp->grid[hdwg0[gc[1].w]];
        gc[1].i[0] = i - gc[1].i[1]*hp->grid[hdwg0[gc[1].w]];
        if (reverse) {
                copystruct(g2, gc+1);
                copystruct(gc+1, gc);
                copystruct(gc, g2);
        }
        return(1);                      /* we're done */
}


int
hdbindex(hp, gc)                /* compute index from beam coordinates */
register HOLO   *hp;
register GCOORD gc[2];
{
        GCOORD  g2[2];
        int     reverse;
        register int    i, j;
                                        /* check ordering and limits */
        if (reverse = gc[0].w > gc[1].w) {
                copystruct(g2, gc+1);
                copystruct(g2+1, gc);
                gc = g2;
        } else if (gc[0].w == gc[1].w)
                return(0);
        if (gc[0].w < 0 | gc[1].w > 5)
                return(0);
        i = 0;                          /* compute index */
        for (j = gc[0].w+1; j < gc[1].w; j++)
                i += hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
        i *= hp->grid[hdwg0[gc[0].w]] * hp->grid[hdwg1[gc[0].w]];
        i += hp->wi[gc[0].w];
        i += (hp->grid[hdwg0[gc[0].w]]*gc[0].i[1] + gc[0].i[0]) *
                        hp->grid[hdwg0[gc[1].w]] * hp->grid[hdwg1[gc[1].w]] ;
        i += hp->grid[hdwg0[gc[1].w]]*gc[1].i[1] + gc[1].i[0];
        if (reverse)
                i += hp->wi[5] - 1;
        return(i);
}


hdcell(cp, hp, gc)              /* compute cell coordinates */
register FVECT  cp[4];  /* returned (may be passed as FVECT cp[2][2]) */
register HOLO   *hp;
register GCOORD *gc;
{
        register FLOAT  *v;
        double  d;
                                        /* compute common component */
        VCOPY(cp[0], hp->orig);
        if (gc->w & 1) {
                v = hp->xv[gc->w>>1];
                cp[0][0] += v[0]; cp[0][1] += v[1]; cp[0][2] += v[2];
        }
        v = hp->xv[hdwg0[gc->w]];
        d = (double)gc->i[0] / hp->grid[hdwg0[gc->w]];
        VSUM(cp[0], cp[0], v, d);
        v = hp->xv[hdwg1[gc->w]];
        d = (double)gc->i[1] / hp->grid[hdwg1[gc->w]];
        VSUM(cp[0], cp[0], v, d);
                                        /* compute x1 sums */
        v = hp->xv[hdwg0[gc->w]];
        d = 1.0 / hp->grid[hdwg0[gc->w]];
        VSUM(cp[1], cp[0], v, d);
        VSUM(cp[3], cp[0], v, d);
                                        /* compute y1 sums */
        v = hp->xv[hdwg1[gc->w]];
        d = 1.0 / hp->grid[hdwg1[gc->w]];
        VSUM(cp[2], cp[0], v, d);
        VSUM(cp[3], cp[3], v, d);
}


hdlseg(lseg, hp, gc)                    /* compute line segment for beam */
register int    lseg[2][3];
register HOLO   *hp;
GCOORD  gc[2];
{
        register int    k;

        for (k = 0; k < 2; k++) {               /* compute end points */
                lseg[k][gc[k].w>>1] = gc[k].w&1 ? hp->grid[gc[k].w>>1]-1 : 0 ;
                lseg[k][hdwg0[gc[k].w]] = gc[k].i[0];
                lseg[k][hdwg1[gc[k].w]] = gc[k].i[1];
        }
        return(1);
}


unsigned
hdcode(hp, d)                   /* compute depth code for d */
HOLO    *hp;
double  d;
{
        double  tl = hp->tlin;
        register long   c;

        if (d <= 0.)
                return(0);
        if (d >= .99*FHUGE)
                return(DCINF);
        if (d < tl)
                return((unsigned)(d*DCLIN/tl));
        c = (long)(log(d/tl)/logstep) + DCLIN;
        return(c > DCINF ? (unsigned)DCINF : (unsigned)c);
}


hdgrid(gp, hp, wp)              /* compute grid coordinates */
FVECT   gp;             /* returned */
register HOLO   *hp;
FVECT   wp;
{
        FVECT   vt;

        VSUB(vt, wp, hp->orig);
        gp[0] = DOT(vt, hp->wg[0]);
        gp[1] = DOT(vt, hp->wg[1]);
        gp[2] = DOT(vt, hp->wg[2]);
}


hdworld(wp, hp, gp)             /* compute world coordinates */
register FVECT  wp;
register HOLO   *hp;
FVECT   gp;
{
        register double d;

        d = gp[0]/hp->grid[0];
        VSUM(wp, hp->orig, hp->xv[0], d);

        d = gp[1]/hp->grid[1];
        VSUM(wp, wp, hp->xv[1], d);

        d = gp[2]/hp->grid[2];
        VSUM(wp, wp, hp->xv[2], d);
}


double
hdray(ro, rd, hp, gc, r)        /* compute ray within a beam */
FVECT   ro, rd;         /* returned */
HOLO    *hp;
GCOORD  gc[2];
BYTE    r[2][2];
{
        FVECT   cp[4], p[2];
        register int    i, j;
        double  d0, d1;
                                        /* compute entry and exit points */
        for (i = 0; i < 2; i++) {
                hdcell(cp, hp, gc+i);
                d0 = (1./256.)*(r[i][0]+.5);
                d1 = (1./256.)*(r[i][1]+.5);
                for (j = 0; j < 3; j++)
                        p[i][j] = (1.-d0-d1)*cp[0][j] +
                                        d0*cp[1][j] + d1*cp[2][j];
        }
        VCOPY(ro, p[0]);                /* assign ray origin and direction */
        VSUB(rd, p[1], p[0]);
        return(normalize(rd));          /* return maximum inside distance */
}


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