src/hd/rhdisp3.c

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

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

/*
 * Holodeck beam support for display process
 */

#include "rholo.h"
#include "rhdisp.h"
#include "view.h"

struct cellist {
        GCOORD  *cl;
        int     n;
};


int
npixels(vp, hr, vr, hp, bi)     /* compute appropriate number to evaluate */
register VIEW   *vp;
int     hr, vr;
HOLO    *hp;
int     bi;
{
        VIEW    vrev;
        GCOORD  gc[2];
        FVECT   cp[4], ip[4];
        double  af, ab;
        register int    i;
                                        /* compute cell corners in image */
        if (!hdbcoord(gc, hp, bi))
                error(CONSISTENCY, "bad beam index in npixels");
        hdcell(cp, hp, gc+1);           /* find cell on front image */
        for (i = 0; i < 4; i++) {
                viewloc(ip[i], vp, cp[i]);
                if (ip[i][2] < 0.) {
                        af = 0;
                        goto getback;
                }
                ip[i][0] *= (double)hr; /* scale by resolution */
                ip[i][1] *= (double)vr;
        }
                                        /* compute front area */
        af = (ip[1][0]-ip[0][0])*(ip[2][1]-ip[0][1]) -
                (ip[2][0]-ip[0][0])*(ip[1][1]-ip[0][1]);
        af += (ip[2][0]-ip[3][0])*(ip[1][1]-ip[3][1]) -
                (ip[1][0]-ip[3][0])*(ip[2][1]-ip[3][1]);
        if (af >= 0) af *= 0.5;
        else af *= -0.5;
getback:
        copystruct(&vrev, vp);          /* compute reverse view */
        for (i = 0; i < 3; i++) {
                vrev.vdir[i] = -vp->vdir[i];
                vrev.vup[i] = -vp->vup[i];
                vrev.hvec[i] = -vp->hvec[i];
                vrev.vvec[i] = -vp->vvec[i];
        }
        hdcell(cp, hp, gc);             /* find cell on back image */
        for (i = 0; i < 4; i++) {
                viewloc(ip[i], &vrev, cp[i]);
                if (ip[i][2] < 0.)
                        return((int)(af + 0.5));
                ip[i][0] *= (double)hr; /* scale by resolution */
                ip[i][1] *= (double)vr;
        }
                                        /* compute back area */
        ab = (ip[1][0]-ip[0][0])*(ip[2][1]-ip[0][1]) -
                (ip[2][0]-ip[0][0])*(ip[1][1]-ip[0][1]);
        ab += (ip[2][0]-ip[3][0])*(ip[1][1]-ip[3][1]) -
                (ip[1][0]-ip[3][0])*(ip[2][1]-ip[3][1]);
        if (ab >= 0) ab *= 0.5;
        else ab *= -0.5;
                                        /* round off smaller area */
        if (af <= ab)
                return((int)(af + 0.5));
        return((int)(ab + 0.5));
}


/*
 * The ray directions that define the pyramid in visit_cells() needn't
 * be normalized, but they must be given in clockwise order as seen
 * from the pyramid's apex (origin).
 */
int
visit_cells(orig, pyrd, hp, vf, dp)     /* visit cells within a pyramid */
FVECT   orig, pyrd[4];          /* pyramid ray directions in clockwise order */
HOLO    *hp;
int     (*vf)();
char    *dp;
{
        int     n = 0;
        int     inflags = 0;
        FVECT   gp, pn[4], lo, ld;
        double  po[4], lbeg, lend, d, t;
        GCOORD  gc;
        register int    i;
                                        /* figure out whose side we're on */
        hdgrid(gp, hp, orig);
        for (i = 0; i < 3; i++) {
                inflags |= (gp[i] > FTINY) << (i<<1);
                inflags |= (gp[i] < hp->grid[i]-FTINY) << (i<<1 | 1);
        }
                                        /* compute pyramid planes */
        for (i = 0; i < 4; i++) {
                fcross(pn[i], pyrd[i], pyrd[(i+1)&03]);
                po[i] = DOT(pn[i], orig);
        }
                                        /* traverse each wall */
        for (gc.w = 0; gc.w < 6; gc.w++) {
                if (!(inflags & 1<<gc.w))       /* origin on wrong side */
                        continue;
                                        /* scanline algorithm */
                for (gc.i[1] = hp->grid[((gc.w>>1)+2)%3]; gc.i[1]--; ) {
                                                /* compute scanline */
                        gp[gc.w>>1] = gc.w&1 ? hp->grid[gc.w>>1] : 0;
                        gp[((gc.w>>1)+1)%3] = 0;
                        gp[((gc.w>>1)+2)%3] = gc.i[1] + 0.5;
                        hdworld(lo, hp, gp);
                        gp[((gc.w>>1)+1)%3] = 1;
                        hdworld(ld, hp, gp);
                        ld[0] -= lo[0]; ld[1] -= lo[1]; ld[2] -= lo[2];
                                                /* find scanline limits */
                        lbeg = 0; lend = hp->grid[((gc.w>>1)+1)%3];
                        for (i = 0; i < 4; i++) {
                                t = DOT(pn[i], lo) - po[i];
                                d = -DOT(pn[i], ld);
                                if (d > FTINY) {                /* <- plane */
                                        if ((t /= d) < lend)
                                                lend = t;
                                } else if (d < -FTINY) {        /* plane -> */
                                        if ((t /= d) > lbeg)
                                                lbeg = t;
                                } else if (t < 0) {             /* outside */
                                        lend = -1;
                                        break;
                                }
                        }
                        if (lbeg >= lend)
                                continue;
                        i = lend + .5;          /* visit cells on this scan */
                        for (gc.i[0] = lbeg + .5; gc.i[0] < i; gc.i[0]++)
                                n += (*vf)(&gc, dp);
                }
        }
        return(n);
}


sect_behind(hp, vp)             /* check if section is "behind" viewpoint */
register HOLO   *hp;
register VIEW   *vp;
{
        FVECT   hcent;
                                        /* compute holodeck section center */
        VSUM(hcent, hp->orig, hp->xv[0], 0.5);
        VSUM(hcent, hcent, hp->xv[1], 0.5);
        VSUM(hcent, hcent, hp->xv[2], 0.5);
                                        /* behind if center is behind */
        return(DOT(vp->vdir,hcent) < DOT(vp->vdir,vp->vp));
}


viewpyramid(org, dir, hp, vp)   /* compute view pyramid */
FVECT   org, dir[4];
HOLO    *hp;
VIEW    *vp;
{
        register int    i;
                                        /* check view type */
        if (vp->type == VT_PAR)
                return(0);
                                        /* in front or behind? */
        if (!sect_behind(hp, vp)) {
                if (viewray(org, dir[0], vp, 0., 0.) < -FTINY)
                        return(0);
                if (viewray(org, dir[1], vp, 0., 1.) < -FTINY)
                        return(0);
                if (viewray(org, dir[2], vp, 1., 1.) < -FTINY)
                        return(0);
                if (viewray(org, dir[3], vp, 1., 0.) < -FTINY)
                        return(0);
                return(1);
        }                               /* reverse pyramid */
        if (viewray(org, dir[3], vp, 0., 0.) < -FTINY)
                return(0);
        if (viewray(org, dir[2], vp, 0., 1.) < -FTINY)
                return(0);
        if (viewray(org, dir[1], vp, 1., 1.) < -FTINY)
                return(0);
        if (viewray(org, dir[0], vp, 1., 0.) < -FTINY)
                return(0);
        for (i = 0; i < 3; i++) {
                dir[0][i] = -dir[0][i];
                dir[1][i] = -dir[1][i];
                dir[2][i] = -dir[2][i];
                dir[3][i] = -dir[3][i];
        }
        return(-1);
}


int
addcell(gcp, cl)                /* add a cell to a list */
GCOORD  *gcp;
register struct cellist *cl;
{
        copystruct(cl->cl+cl->n, gcp);
        cl->n++;
        return(1);
}


int
cellcmp(gcp1, gcp2)             /* visit_cells() cell ordering */
register GCOORD *gcp1, *gcp2;
{
        register int    c;

        if ((c = gcp1->w - gcp2->w))
                return(c);
        if ((c = gcp2->i[1] - gcp1->i[1]))      /* wg1 is reverse-ordered */
                return(c);
        return(gcp1->i[0] - gcp2->i[0]);
}


GCOORD *
getviewcells(np, hp, vp)        /* get ordered cell list for section view */
int     *np;            /* returned number of cells (negative if reversed) */
register HOLO   *hp;
VIEW    *vp;
{
        FVECT   org, dir[4];
        int     orient;
        struct cellist  cl;
                                        /* compute view pyramid */
        *np = 0;
        orient = viewpyramid(org, dir, hp, vp);
        if (!orient)
                return(NULL);
                                        /* allocate enough list space */
        cl.n = 2*(      hp->grid[0]*hp->grid[1] +
                        hp->grid[0]*hp->grid[2] +
                        hp->grid[1]*hp->grid[2] );
        cl.cl = (GCOORD *)malloc(cl.n*sizeof(GCOORD));
        if (cl.cl == NULL)
                goto memerr;
        cl.n = 0;                       /* add cells within pyramid */
        visit_cells(org, dir, hp, addcell, &cl);
        if (!cl.n) {
                free((char *)cl.cl);
                return(NULL);
        }
        *np = cl.n * orient;
#if 0
        /* We're just going to free this memory in a moment, and list is
         * sorted automatically by visit_cells(), so we don't need this.
         */
                                        /* optimize memory use */
        cl.cl = (GCOORD *)realloc((char *)cl.cl, cl.n*sizeof(GCOORD));
        if (cl.cl == NULL)
                goto memerr;
                                        /* sort the list */
        qsort((char *)cl.cl, cl.n, sizeof(GCOORD), cellcmp);
#endif
        return(cl.cl);
memerr:
        error(SYSTEM, "out of memory in getviewcells");
}


gridlines(f)                    /* run through holodeck section grid lines */
int     (*f)();
{
        register int    hd, w, i;
        int     g0, g1;
        FVECT   wp[2], mov;
        double  d;
                                        /* do each wall on each section */
        for (hd = 0; hdlist[hd] != NULL; hd++)
                for (w = 0; w < 6; w++) {
                        g0 = ((w>>1)+1)%3;
                        g1 = ((w>>1)+2)%3;
                        d = 1.0/hdlist[hd]->grid[g0];
                        mov[0] = d * hdlist[hd]->xv[g0][0];
                        mov[1] = d * hdlist[hd]->xv[g0][1];
                        mov[2] = d * hdlist[hd]->xv[g0][2];
                        if (w & 1) {
                                VSUM(wp[0], hdlist[hd]->orig,
                                                hdlist[hd]->xv[w>>1], 1.);
                                VSUM(wp[0], wp[0], mov, 1.);
                        } else
                                VCOPY(wp[0], hdlist[hd]->orig);
                        VSUM(wp[1], wp[0], hdlist[hd]->xv[g1], 1.);
                        for (i = hdlist[hd]->grid[g0]; ; ) {    /* g0 lines */
                                (*f)(wp);
                                if (!--i) break;
                                wp[0][0] += mov[0]; wp[0][1] += mov[1];
                                wp[0][2] += mov[2]; wp[1][0] += mov[0];
                                wp[1][1] += mov[1]; wp[1][2] += mov[2];
                        }
                        d = 1.0/hdlist[hd]->grid[g1];
                        mov[0] = d * hdlist[hd]->xv[g1][0];
                        mov[1] = d * hdlist[hd]->xv[g1][1];
                        mov[2] = d * hdlist[hd]->xv[g1][2];
                        if (w & 1)
                                VSUM(wp[0], hdlist[hd]->orig,
                                                hdlist[hd]->xv[w>>1], 1.);
                        else
                                VSUM(wp[0], hdlist[hd]->orig, mov, 1.);
                        VSUM(wp[1], wp[0], hdlist[hd]->xv[g0], 1.);
                        for (i = hdlist[hd]->grid[g1]; ; ) {    /* g1 lines */
                                (*f)(wp);
                                if (!--i) break;
                                wp[0][0] += mov[0]; wp[0][1] += mov[1];
                                wp[0][2] += mov[2]; wp[1][0] += mov[0];
                                wp[1][1] += mov[1]; wp[1][2] += mov[2];
                        }
                }
}
Revision:	3.7
Committed:	Wed Nov 26 20:12:19 1997 UTC (27 years, 11 months ago) by gregl
Content type:	text/plain
Branch:	MAIN
Changes since 3.6:	+28 -13 lines
Log Message:	fixed bug in gridlines()
#	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	gregl	3.2	* Holodeck beam support for display process
9	gregl	3.1	*/
10
11			#include "rholo.h"
12			#include "rhdisp.h"
13			#include "view.h"
14
15	gregl	3.4	struct cellist {
16			GCOORD *cl;
17			int n;
18			};
19	gregl	3.1
20	gregl	3.4
21	gregl	3.1	int
22			npixels(vp, hr, vr, hp, bi) /* compute appropriate number to evaluate */
23	gregl	3.5	register VIEW *vp;
24	gregl	3.1	int hr, vr;
25			HOLO *hp;
26			int bi;
27			{
28	gregl	3.5	VIEW vrev;
29	gregl	3.1	GCOORD gc[2];
30	gregl	3.5	FVECT cp[4], ip[4];
31			double af, ab;
32	gregl	3.1	register int i;
33			/* compute cell corners in image */
34			if (!hdbcoord(gc, hp, bi))
35			error(CONSISTENCY, "bad beam index in npixels");
36	gregl	3.5	hdcell(cp, hp, gc+1); /* find cell on front image */
37			for (i = 0; i < 4; i++) {
38			viewloc(ip[i], vp, cp[i]);
39			if (ip[i][2] < 0.) {
40			af = 0;
41			goto getback;
42	gregl	3.4	}
43	gregl	3.5	ip[i][0] = (double)hr; / scale by resolution */
44			ip[i][1] *= (double)vr;
45	gregl	3.4	}
46	gregl	3.5	/* compute front area */
47			af = (ip[1][0]-ip[0][0])*(ip[2][1]-ip[0][1]) -
48			(ip[2][0]-ip[0][0])*(ip[1][1]-ip[0][1]);
49			af += (ip[2][0]-ip[3][0])*(ip[1][1]-ip[3][1]) -
50			(ip[1][0]-ip[3][0])*(ip[2][1]-ip[3][1]);
51			if (af >= 0) af *= 0.5;
52			else af *= -0.5;
53			getback:
54			copystruct(&vrev, vp); /* compute reverse view */
55			for (i = 0; i < 3; i++) {
56			vrev.vdir[i] = -vp->vdir[i];
57			vrev.vup[i] = -vp->vup[i];
58			vrev.hvec[i] = -vp->hvec[i];
59			vrev.vvec[i] = -vp->vvec[i];
60			}
61			hdcell(cp, hp, gc); /* find cell on back image */
62	gregl	3.1	for (i = 0; i < 4; i++) {
63	gregl	3.5	viewloc(ip[i], &vrev, cp[i]);
64	gregl	3.1	if (ip[i][2] < 0.)
65	gregl	3.5	return((int)(af + 0.5));
66	gregl	3.1	ip[i][0] = (double)hr; / scale by resolution */
67			ip[i][1] *= (double)vr;
68			}
69	gregl	3.5	/* compute back area */
70			ab = (ip[1][0]-ip[0][0])*(ip[2][1]-ip[0][1]) -
71	gregl	3.1	(ip[2][0]-ip[0][0])*(ip[1][1]-ip[0][1]);
72	gregl	3.5	ab += (ip[2][0]-ip[3][0])*(ip[1][1]-ip[3][1]) -
73	gregl	3.1	(ip[1][0]-ip[3][0])*(ip[2][1]-ip[3][1]);
74	gregl	3.5	if (ab >= 0) ab *= 0.5;
75			else ab *= -0.5;
76			/* round off smaller area */
77			if (af <= ab)
78			return((int)(af + 0.5));
79			return((int)(ab + 0.5));
80	gregl	3.1	}
81
82
83			/*
84			* The ray directions that define the pyramid in visit_cells() needn't
85			* be normalized, but they must be given in clockwise order as seen
86			* from the pyramid's apex (origin).
87			*/
88			int
89			visit_cells(orig, pyrd, hp, vf, dp) /* visit cells within a pyramid */
90			FVECT orig, pyrd[4]; /* pyramid ray directions in clockwise order */
91			HOLO *hp;
92			int (*vf)();
93			char *dp;
94			{
95			int n = 0;
96			int inflags = 0;
97			FVECT gp, pn[4], lo, ld;
98			double po[4], lbeg, lend, d, t;
99			GCOORD gc;
100			register int i;
101			/* figure out whose side we're on */
102			hdgrid(gp, hp, orig);
103			for (i = 0; i < 3; i++) {
104			inflags \|= (gp[i] > FTINY) << (i<<1);
105			inflags \|= (gp[i] < hp->grid[i]-FTINY) << (i<<1 \| 1);
106			}
107			/* compute pyramid planes */
108			for (i = 0; i < 4; i++) {
109			fcross(pn[i], pyrd[i], pyrd[(i+1)&03]);
110			po[i] = DOT(pn[i], orig);
111			}
112			/* traverse each wall */
113			for (gc.w = 0; gc.w < 6; gc.w++) {
114			if (!(inflags & 1<<gc.w)) /* origin on wrong side */
115			continue;
116			/* scanline algorithm */
117			for (gc.i[1] = hp->grid[((gc.w>>1)+2)%3]; gc.i[1]--; ) {
118			/* compute scanline */
119			gp[gc.w>>1] = gc.w&1 ? hp->grid[gc.w>>1] : 0;
120			gp[((gc.w>>1)+1)%3] = 0;
121			gp[((gc.w>>1)+2)%3] = gc.i[1] + 0.5;
122			hdworld(lo, hp, gp);
123			gp[((gc.w>>1)+1)%3] = 1;
124			hdworld(ld, hp, gp);
125	gregl	3.2	ld[0] -= lo[0]; ld[1] -= lo[1]; ld[2] -= lo[2];
126	gregl	3.1	/* find scanline limits */
127			lbeg = 0; lend = hp->grid[((gc.w>>1)+1)%3];
128			for (i = 0; i < 4; i++) {
129			t = DOT(pn[i], lo) - po[i];
130			d = -DOT(pn[i], ld);
131	gregl	3.2	if (d > FTINY) { /* <- plane */
132	gregl	3.1	if ((t /= d) < lend)
133			lend = t;
134	gregl	3.2	} else if (d < -FTINY) { /* plane -> */
135	gregl	3.1	if ((t /= d) > lbeg)
136			lbeg = t;
137	gregl	3.3	} else if (t < 0) { /* outside */
138			lend = -1;
139			break;
140			}
141	gregl	3.1	}
142	gregl	3.3	if (lbeg >= lend)
143			continue;
144	gregl	3.1	i = lend + .5; /* visit cells on this scan */
145			for (gc.i[0] = lbeg + .5; gc.i[0] < i; gc.i[0]++)
146			n += (*vf)(&gc, dp);
147			}
148			}
149			return(n);
150			}
151
152
153	gregl	3.4	sect_behind(hp, vp) /* check if section is "behind" viewpoint */
154			register HOLO *hp;
155			register VIEW *vp;
156			{
157			FVECT hcent;
158			/* compute holodeck section center */
159			VSUM(hcent, hp->orig, hp->xv[0], 0.5);
160			VSUM(hcent, hcent, hp->xv[1], 0.5);
161			VSUM(hcent, hcent, hp->xv[2], 0.5);
162			/* behind if center is behind */
163			return(DOT(vp->vdir,hcent) < DOT(vp->vdir,vp->vp));
164			}
165
166
167			viewpyramid(org, dir, hp, vp) /* compute view pyramid */
168			FVECT org, dir[4];
169			HOLO *hp;
170			VIEW *vp;
171			{
172			register int i;
173			/* check view type */
174			if (vp->type == VT_PAR)
175			return(0);
176			/* in front or behind? */
177			if (!sect_behind(hp, vp)) {
178			if (viewray(org, dir[0], vp, 0., 0.) < -FTINY)
179			return(0);
180			if (viewray(org, dir[1], vp, 0., 1.) < -FTINY)
181			return(0);
182			if (viewray(org, dir[2], vp, 1., 1.) < -FTINY)
183			return(0);
184			if (viewray(org, dir[3], vp, 1., 0.) < -FTINY)
185			return(0);
186			return(1);
187			} /* reverse pyramid */
188			if (viewray(org, dir[3], vp, 0., 0.) < -FTINY)
189			return(0);
190			if (viewray(org, dir[2], vp, 0., 1.) < -FTINY)
191			return(0);
192			if (viewray(org, dir[1], vp, 1., 1.) < -FTINY)
193			return(0);
194			if (viewray(org, dir[0], vp, 1., 0.) < -FTINY)
195			return(0);
196			for (i = 0; i < 3; i++) {
197			dir[0][i] = -dir[0][i];
198			dir[1][i] = -dir[1][i];
199			dir[2][i] = -dir[2][i];
200			dir[3][i] = -dir[3][i];
201			}
202			return(-1);
203			}
204
205
206	gregl	3.1	int
207			addcell(gcp, cl) /* add a cell to a list */
208			GCOORD *gcp;
209	gregl	3.4	register struct cellist *cl;
210	gregl	3.1	{
211	gregl	3.4	copystruct(cl->cl+cl->n, gcp);
212			cl->n++;
213	gregl	3.1	return(1);
214			}
215
216
217			int
218			cellcmp(gcp1, gcp2) /* visit_cells() cell ordering */
219			register GCOORD gcp1, gcp2;
220			{
221			register int c;
222
223			if ((c = gcp1->w - gcp2->w))
224			return(c);
225			if ((c = gcp2->i[1] - gcp1->i[1])) /* wg1 is reverse-ordered */
226			return(c);
227			return(gcp1->i[0] - gcp2->i[0]);
228			}
229
230
231	gregl	3.4	GCOORD *
232			getviewcells(np, hp, vp) /* get ordered cell list for section view */
233			int np; / returned number of cells (negative if reversed) */
234	gregl	3.1	register HOLO *hp;
235			VIEW *vp;
236			{
237			FVECT org, dir[4];
238	gregl	3.4	int orient;
239			struct cellist cl;
240	gregl	3.1	/* compute view pyramid */
241	gregl	3.4	*np = 0;
242			orient = viewpyramid(org, dir, hp, vp);
243			if (!orient)
244			return(NULL);
245	gregl	3.1	/* allocate enough list space */
246	gregl	3.4	cl.n = 2( hp->grid[0]hp->grid[1] +
247			hp->grid[0]*hp->grid[2] +
248			hp->grid[1]*hp->grid[2] );
249			cl.cl = (GCOORD )malloc(cl.nsizeof(GCOORD));
250			if (cl.cl == NULL)
251	gregl	3.1	goto memerr;
252	gregl	3.4	cl.n = 0; /* add cells within pyramid */
253			visit_cells(org, dir, hp, addcell, &cl);
254			if (!cl.n) {
255			free((char *)cl.cl);
256	gregl	3.1	return(NULL);
257			}
258	gregl	3.4	np = cl.n orient;
259	gregl	3.1	#if 0
260	gregl	3.2	/* We're just going to free this memory in a moment, and list is
261			* sorted automatically by visit_cells(), so we don't need this.
262			*/
263	gregl	3.4	/* optimize memory use */
264			cl.cl = (GCOORD )realloc((char )cl.cl, cl.n*sizeof(GCOORD));
265			if (cl.cl == NULL)
266			goto memerr;
267	gregl	3.1	/* sort the list */
268	gregl	3.4	qsort((char *)cl.cl, cl.n, sizeof(GCOORD), cellcmp);
269	gregl	3.1	#endif
270	gregl	3.4	return(cl.cl);
271	gregl	3.1	memerr:
272			error(SYSTEM, "out of memory in getviewcells");
273			}
274	gregl	3.6
275
276			gridlines(f) /* run through holodeck section grid lines */
277			int (*f)();
278			{
279			register int hd, w, i;
280			int g0, g1;
281	gregl	3.7	FVECT wp[2], mov;
282	gregl	3.6	double d;
283			/* do each wall on each section */
284			for (hd = 0; hdlist[hd] != NULL; hd++)
285			for (w = 0; w < 6; w++) {
286			g0 = ((w>>1)+1)%3;
287			g1 = ((w>>1)+2)%3;
288	gregl	3.7	d = 1.0/hdlist[hd]->grid[g0];
289			mov[0] = d * hdlist[hd]->xv[g0][0];
290			mov[1] = d * hdlist[hd]->xv[g0][1];
291			mov[2] = d * hdlist[hd]->xv[g0][2];
292			if (w & 1) {
293	gregl	3.6	VSUM(wp[0], hdlist[hd]->orig,
294			hdlist[hd]->xv[w>>1], 1.);
295	gregl	3.7	VSUM(wp[0], wp[0], mov, 1.);
296			} else
297			VCOPY(wp[0], hdlist[hd]->orig);
298			VSUM(wp[1], wp[0], hdlist[hd]->xv[g1], 1.);
299			for (i = hdlist[hd]->grid[g0]; ; ) { /* g0 lines */
300	gregl	3.6	(*f)(wp);
301	gregl	3.7	if (!--i) break;
302			wp[0][0] += mov[0]; wp[0][1] += mov[1];
303			wp[0][2] += mov[2]; wp[1][0] += mov[0];
304			wp[1][1] += mov[1]; wp[1][2] += mov[2];
305	gregl	3.6	}
306	gregl	3.7	d = 1.0/hdlist[hd]->grid[g1];
307			mov[0] = d * hdlist[hd]->xv[g1][0];
308			mov[1] = d * hdlist[hd]->xv[g1][1];
309			mov[2] = d * hdlist[hd]->xv[g1][2];
310			if (w & 1)
311	gregl	3.6	VSUM(wp[0], hdlist[hd]->orig,
312			hdlist[hd]->xv[w>>1], 1.);
313	gregl	3.7	else
314			VSUM(wp[0], hdlist[hd]->orig, mov, 1.);
315			VSUM(wp[1], wp[0], hdlist[hd]->xv[g0], 1.);
316			for (i = hdlist[hd]->grid[g1]; ; ) { /* g1 lines */
317	gregl	3.6	(*f)(wp);
318	gregl	3.7	if (!--i) break;
319			wp[0][0] += mov[0]; wp[0][1] += mov[1];
320			wp[0][2] += mov[2]; wp[1][0] += mov[0];
321			wp[1][1] += mov[1]; wp[1][2] += mov[2];
322	gregl	3.6	}
323			}
324			}