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];

static double   logstep;

static int      wg0[6] = {1,1,2,2,0,0};
static int      wg1[6] = {2,2,0,0,1,1};


hdcompgrid(hp)                  /* compute derived grid vector and index */
register HOLO   *hp;
{
        FVECT   AxB;
        double  d;
        register FLOAT  *v;
        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(AxB, hp->xv[(i+1)%3], v=hp->xv[(i+2)%3]);
                VCOPY(hp->wn[i], AxB);
                if (normalize(hp->wn[i]) == 0.)
                        error(USER, "degenerate holodeck section");
                hp->wo[i<<1] = DOT(hp->wn[i],hp->orig);
                hp->wo[i<<1|1] = hp->wo[i<<1] + DOT(hp->wn[i],hp->xv[i]);
                fcross(hp->gv[i][0], v, AxB);
                d = DOT(v,v) / DOT(hp->gv[i][0],hp->gv[i][0]) *
                                hp->grid[(i+1)%3];
                for (j = 0; j < 3; j++)
                        hp->gv[i][0][j] *= d;
                fcross(hp->gv[i][1], AxB, v=hp->xv[(i+1)%3]);
                d = DOT(v,v) / DOT(hp->gv[i][1],hp->gv[i][1]) *
                                hp->grid[(i+2)%3];
                for (j = 0; j < 3; j++)
                        hp->gv[i][1][j] *= 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[wg0[j]] * hp->grid[wg1[j]];
                hp->wi[i] *= hp->grid[wg0[i-1]] * hp->grid[wg1[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[wg0[j]] * hp->grid[wg1[j]];
        while (++j < 5) {
                n = n2 * hp->grid[wg0[j]] * hp->grid[wg1[j]];
                if (n > i)
                        break;
                i -= n;
        }
        gc[1].w = j;
                                        /* find position on w0 */
        n2 = hp->grid[wg0[j]] * hp->grid[wg1[j]];
        n = i / n2;
        gc[0].i[1] = n / hp->grid[wg0[gc[0].w]];
        gc[0].i[0] = n - gc[0].i[1]*hp->grid[wg0[gc[0].w]];
        i -= n*n2;
                                        /* find position on w1 */
        gc[1].i[1] = i / hp->grid[wg0[gc[1].w]];
        gc[1].i[0] = i - gc[1].i[1]*hp->grid[wg0[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[wg0[j]] * hp->grid[wg1[j]];
        i *= hp->grid[wg0[gc[0].w]] * hp->grid[wg1[gc[0].w]];
        i += hp->wi[gc[0].w];
        i += (hp->grid[wg0[gc[0].w]]*gc[0].i[1] + gc[0].i[0]) *
                        hp->grid[wg0[gc[1].w]] * hp->grid[wg1[gc[1].w]] ;
        i += hp->grid[wg0[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]) */
HOLO    *hp;
register GCOORD *gc;
{
        register int    i;
        register FLOAT  *v;
        double  d;
                                        /* compute each corner */
        for (i = 0; i < 4; i++) {
                VCOPY(cp[i], hp->orig);
                if (gc->w & 1) {
                        v = hp->xv[gc->w>>1];
                        cp[i][0] += *v++; cp[i][1] += *v++; cp[i][2] += *v;
                }
                d = (double)( gc->i[0] + (i&1) ) / hp->grid[wg0[gc->w]];
                v = hp->xv[wg0[gc->w]];
                cp[i][0] += d * *v++; cp[i][1] += d * *v++; cp[i][2] += d * *v;

                d = (double)( gc->i[1] + (i>>1) ) / hp->grid[wg1[gc->w]];
                v = hp->xv[wg1[gc->w]];
                cp[i][0] += d * *v++; cp[i][1] += d * *v++; cp[i][2] += d * *v;
        }
}


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

        if (!hdbcoord(gc, hp, i))               /* compute grid coordinates */
                return(0);
        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][wg0[gc[k].w]] = gc[k].i[0];
                lseg[k][wg1[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 unsigned       c;

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


double
hdray(ro, rd, hp, gc, r)        /* compute ray within a beam */
FVECT   ro, rd;         /* returned */
register HOLO   *hp;
register GCOORD gc[2];
BYTE    r[2][2];
{
        FVECT   p[2];
        register int    i;
        register FLOAT  *v;
        double  d;
                                        /* compute entry and exit points */
        for (i = 0; i < 2; i++) {
                VCOPY(p[i], hp->orig);
                if (gc[i].w & 1) {
                        v = hp->xv[gc[i].w>>1];
                        p[i][0] += *v++; p[i][1] += *v++; p[i][2] += *v;
                }
                d = ( gc[i].i[0] + (1./256.)*(r[i][0]+.5) ) /
                                hp->grid[wg0[gc[i].w]];
                v = hp->xv[wg0[gc[i].w]];
                p[i][0] += d * *v++; p[i][1] += d * *v++; p[i][2] += d * *v;
                d = ( gc[i].i[1] + (1./256.)*(r[i][1]+.5) ) /
                                hp->grid[wg1[gc[i].w]];
                v = hp->xv[wg1[gc[i].w]];
                p[i][0] += d * *v++; p[i][1] += d * *v++; p[i][2] += d * *v;
        }
        VCOPY(ro, p[0]);                /* assign ray origin and direction */
        rd[0] = p[1][0] - p[0][0];
        rd[1] = p[1][1] - p[0][1];
        rd[2] = p[1][2] - p[0][2];
        return(normalize(rd));          /* return maximum inside distance */
}


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