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]) */
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[wg0[gc->w]];
        d = (double)gc->i[0] / hp->grid[wg0[gc->w]];
        VSUM(cp[0], cp[0], v, d);
        v = hp->xv[wg1[gc->w]];
        d = (double)gc->i[1] / hp->grid[wg1[gc->w]];
        VSUM(cp[0], cp[0], v, d);
                                        /* compute x1 sums */
        v = hp->xv[wg0[gc->w]];
        d = 1.0 / hp->grid[wg0[gc->w]];
        VSUM(cp[1], cp[0], v, d);
        VSUM(cp[3], cp[0], v, d);
                                        /* compute y1 sums */
        v = hp->xv[wg1[gc->w]];
        d = 1.0 / hp->grid[wg1[gc->w]];
        VSUM(cp[2], cp[0], v, d);
        VSUM(cp[3], cp[3], v, d);
}


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 */
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 */
        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.5
Committed:	Tue Nov 11 15:14:02 1997 UTC (26 years, 5 months ago) by gregl
Content type:	text/plain
Branch:	MAIN
Changes since 3.4:	+33 -35 lines
Log Message:	consolidated code for hdcell() and hdray()
#	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	gregl	3.5	register HOLO *hp;
184	gregl	3.4	register GCOORD *gc;
185			{
186			register FLOAT *v;
187			double d;
188	gregl	3.5	/* compute common component */
189			VCOPY(cp[0], hp->orig);
190			if (gc->w & 1) {
191			v = hp->xv[gc->w>>1];
192			cp[0][0] += v[0]; cp[0][1] += v[1]; cp[0][2] += v[2];
193	gregl	3.4	}
194	gregl	3.5	v = hp->xv[wg0[gc->w]];
195			d = (double)gc->i[0] / hp->grid[wg0[gc->w]];
196			VSUM(cp[0], cp[0], v, d);
197			v = hp->xv[wg1[gc->w]];
198			d = (double)gc->i[1] / hp->grid[wg1[gc->w]];
199			VSUM(cp[0], cp[0], v, d);
200			/* compute x1 sums */
201			v = hp->xv[wg0[gc->w]];
202			d = 1.0 / hp->grid[wg0[gc->w]];
203			VSUM(cp[1], cp[0], v, d);
204			VSUM(cp[3], cp[0], v, d);
205			/* compute y1 sums */
206			v = hp->xv[wg1[gc->w]];
207			d = 1.0 / hp->grid[wg1[gc->w]];
208			VSUM(cp[2], cp[0], v, d);
209			VSUM(cp[3], cp[3], v, d);
210	gregl	3.4	}
211
212
213	gregl	3.1	hdlseg(lseg, hp, i) /* compute line segment for beam */
214	gregl	3.2	register int lseg[2][3];
215	gregl	3.1	register HOLO *hp;
216			int i;
217			{
218	gregl	3.3	GCOORD gc[2];
219	gregl	3.1	register int k;
220
221	gregl	3.2	if (!hdbcoord(gc, hp, i)) /* compute grid coordinates */
222	gregl	3.1	return(0);
223	gregl	3.2	for (k = 0; k < 2; k++) { /* compute end points */
224			lseg[k][gc[k].w>>1] = gc[k].w&1 ? hp->grid[gc[k].w>>1]-1 : 0 ;
225			lseg[k][wg0[gc[k].w]] = gc[k].i[0];
226			lseg[k][wg1[gc[k].w]] = gc[k].i[1];
227			}
228	gregl	3.1	return(1);
229			}
230
231
232			unsigned
233			hdcode(hp, d) /* compute depth code for d */
234			HOLO *hp;
235			double d;
236			{
237			double tl = hp->tlin;
238			register unsigned c;
239
240			if (d <= 0.)
241			return(0);
242			if (d >= .99*FHUGE)
243			return(DCINF);
244			if (d < tl)
245			return((unsigned)(d*DCLIN/tl));
246			c = (unsigned)(log(d/tl)/logstep) + DCLIN;
247			return(c > DCINF ? DCINF : c);
248			}
249
250
251			double
252			hdray(ro, rd, hp, gc, r) /* compute ray within a beam */
253			FVECT ro, rd; /* returned */
254	gregl	3.5	HOLO *hp;
255			GCOORD gc[2];
256	gregl	3.1	BYTE r[2][2];
257			{
258	gregl	3.5	FVECT cp[4], p[2];
259			register int i, j;
260			double d0, d1;
261	gregl	3.1	/* compute entry and exit points */
262			for (i = 0; i < 2; i++) {
263	gregl	3.5	hdcell(cp, hp, gc+i);
264			d0 = (1./256.)*(r[i][0]+.5);
265			d1 = (1./256.)*(r[i][1]+.5);
266			for (j = 0; j < 3; j++)
267			p[i][j] = (1.-d0-d1)*cp[0][j] +
268			d0cp[1][j] + d1cp[2][j];
269	gregl	3.1	}
270			VCOPY(ro, p[0]); /* assign ray origin and direction */
271			rd[0] = p[1][0] - p[0][0];
272			rd[1] = p[1][1] - p[0][1];
273			rd[2] = p[1][2] - p[0][2];
274			return(normalize(rd)); /* return maximum inside distance */
275			}
276
277
278			double
279			hdinter(gc, r, hp, ro, rd) /* compute ray intersection with section */
280	gregl	3.3	register GCOORD gc[2]; /* returned */
281	gregl	3.1	BYTE r[2][2]; /* returned */
282			register HOLO *hp;
283			FVECT ro, rd; /* rd should be normalized */
284			{
285			FVECT p[2], vt;
286			double d, t0, t1, d0, d1;
287			register FLOAT *v;
288			register int i;
289			/* first, intersect walls */
290			gc[0].w = gc[1].w = -1;
291			t0 = -FHUGE; t1 = FHUGE;
292			for (i = 0; i < 3; i++) { /* for each wall pair */
293			d = -DOT(rd, hp->wn[i]); /* plane distance */
294			if (d <= FTINY && d >= -FTINY) /* check for parallel */
295			continue;
296			d1 = DOT(ro, hp->wn[i]); /* ray distances */
297			d0 = (d1 - hp->wo[i<<1]) / d;
298			d1 = (d1 - hp->wo[i<<1\|1]) / d;
299			if (d0 < d1) { /* check against best */
300			if (d0 > t0) {
301			t0 = d0;
302			gc[0].w = i<<1;
303			}
304			if (d1 < t1) {
305			t1 = d1;
306			gc[1].w = i<<1 \| 1;
307			}
308			} else {
309			if (d1 > t0) {
310			t0 = d1;
311			gc[0].w = i<<1 \| 1;
312			}
313			if (d0 < t1) {
314			t1 = d0;
315			gc[1].w = i<<1;
316			}
317			}
318			}
319			if (gc[0].w < 0 \| gc[1].w < 0) /* paranoid check */
320			return(FHUGE);
321			/* compute intersections */
322			for (i = 0; i < 3; i++) {
323			p[0][i] = ro[i] + rd[i]*t0;
324			p[1][i] = ro[i] + rd[i]*t1;
325			}
326			/* now, compute grid coordinates */
327			for (i = 0; i < 2; i++) {
328			vt[0] = p[i][0] - hp->orig[0];
329			vt[1] = p[i][1] - hp->orig[1];
330			vt[2] = p[i][2] - hp->orig[2];
331			if (gc[i].w & 1) {
332			v = hp->xv[gc[i].w>>1];
333			vt[0] -= v++; vt[1] -= v++; vt[2] -= *v;
334			}
335			v = hp->gv[gc[i].w>>1][0];
336			d = DOT(vt, v);
337	gregl	3.2	if (d < 0. \|\| (gc[i].i[0] = d) >= hp->grid[wg0[gc[i].w]])
338	gregl	3.1	return(FHUGE); /* outside wall */
339			r[i][0] = 256. * (d - gc[i].i[0]);
340			v = hp->gv[gc[i].w>>1][1];
341			d = DOT(vt, v);
342	gregl	3.2	if (d < 0. \|\| (gc[i].i[1] = d) >= hp->grid[wg1[gc[i].w]])
343	gregl	3.1	return(FHUGE); /* outside wall */
344			r[i][1] = 256. * (d - gc[i].i[1]);
345			}
346			/* return distance from entry point */
347			vt[0] = ro[0] - p[0][0];
348			vt[1] = ro[1] - p[0][1];
349			vt[2] = ro[2] - p[0][2];
350			return(DOT(vt,rd));
351			}