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;
{
        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->wn[i], hp->xv[(i+1)%3], hp->xv[(i+2)%3]);
                if (normalize(hp->wn[i]) == 0.)
                        error(USER, "degenerate holodeck section");
                hp->wo[i<<1] = DOT(hp->wn[i],hp->orig);
                d = DOT(hp->wn[i],hp->xv[i]);
                hp->wo[i<<1|1] = hp->wo[i<<1] + d;
                hp->wg[i] = (double)hp->grid[i] / 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, 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][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);
}


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

        vt[0] = wp[0] - hp->orig[0];
        vt[1] = wp[1] - hp->orig[1];
        vt[2] = wp[2] - hp->orig[2];
        gp[0] = DOT(vt, hp->wn[0]) * hp->wg[0];
        gp[1] = DOT(vt, hp->wn[1]) * hp->wg[1];
        gp[2] = DOT(vt, hp->wn[2]) * 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 */
        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, ed, hp, ro, rd)  /* compute ray intersection with section */
register GCOORD gc[2];  /* returned */
BYTE    r[2][2];        /* returned */
double  *ed;            /* returned (optional) */
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];
                v = hp->wn[wg0[gc[i].w]];
                d = DOT(vt, v) * hp->wg[wg0[gc[i].w]];
                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->wn[wg1[gc[i].w]];
                d = DOT(vt, v) * hp->wg[wg1[gc[i].w]];
                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]);
        }

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