src/hd/holo.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#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 = hp->grid[i] / d;
                hp->wg[i][0] *= d; hp->wg[i][1] *= d; hp->wg[i][2] *= 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];
        }
}


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