src/px/warp3d.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 * 3D warping routines.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fvect.h"
#include "warp3d.h"

#define MIND    1e-5            /* minimum distance between input points */

typedef struct {
        GNDX    n;              /* index must be first */
        W3VEC   v;
} KEYDP;                /* key/data allocation pair */

#define fgetvec(p,v)    (fgetval(p,'f',v) > 0 && fgetval(p,'f',v+1) > 0 \
                                && fgetval(p,'f',v+2) > 0)

#define AHUNK   24              /* number of points to allocate at a time */


double
wpdist2(p1, p2)                 /* compute square of distance between points */
register W3VEC  p1, p2;
{
        double  d, d2;

        d = p1[0] - p2[0]; d2  = d*d;
        d = p1[1] - p2[1]; d2 += d*d;
        d = p1[2] - p2[2]; d2 += d*d;
        return(d2);
}


int
warp3d(po, pi, wp)              /* warp 3D point according to the given map */
W3VEC   po, pi;
register WARP3D *wp;
{
        int     rval = W3OK;
        GNDX    gi;
        W3VEC   gd, ov;

        if (wp->grid.gn[0] == 0 && (rval = new3dgrid(wp)) != W3OK)
                return(rval);
        rval |= gridpoint(gi, gd, pi, &wp->grid);
        if (wp->grid.flags & W3EXACT)
                rval |= warp3dex(ov, pi, wp);
        else if (wp->grid.flags & W3FAST)
                rval |= get3dgpt(ov, gi, wp);
        else
                rval |= get3dgin(ov, gi, gd, wp);
        po[0] = pi[0] + ov[0];
        po[1] = pi[1] + ov[1];
        po[2] = pi[2] + ov[2];
        return(rval);
}


int
gridpoint(ndx, rem, ipt, gp)    /* compute grid position for ipt */
GNDX    ndx;
W3VEC   rem, ipt;
register struct grid3d  *gp;
{
        int     rval = W3OK;
        register int    i;

        for (i = 0; i < 3; i++) {
                rem[i] = (ipt[i] - gp->gmin[i])/gp->gstep[i];
                if (rem[i] < 0.) {
                        ndx[i] = 0;
                        rval = W3GAMUT;
                } else if ((int)rem[i] >= gp->gn[i]) {
                        ndx[i] = gp->gn[i] - 1;
                        rval = W3GAMUT;
                } else
                        ndx[i] = (int)rem[i];
                rem[i] -= (double)ndx[i];
        }
        return(rval);
}


int
get3dgpt(ov, ndx, wp)           /* get value for voxel */
W3VEC   ov;
GNDX    ndx;
register WARP3D *wp;
{
        W3VEC   gpt;
        register LUENT  *le;
        KEYDP   *kd;
        int     rval = W3OK;
        register int    i;

        le = lu_find(&wp->grid.gtab, ndx);
        if (le == NULL)
                return(W3ERROR);
        if (le->data == NULL) {
                if (le->key != NULL)
                        kd = (KEYDP *)le->key;
                else if ((kd = (KEYDP *)malloc(sizeof(KEYDP))) == NULL)
                        return(W3ERROR);
                for (i = 0; i < 3; i++) {
                        kd->n[i] = ndx[i];
                        gpt[i] = wp->grid.gmin[i] + ndx[i]*wp->grid.gstep[i];
                        if (wp->grid.flags & W3FAST)    /* on centers */
                                gpt[i] += .5*wp->grid.gstep[i];
                }
                rval = warp3dex(kd->v, gpt, wp);
                le->key = (char *)kd->n;
                le->data = (char *)kd->v;
        }
        W3VCPY(ov, (float *)le->data);
        return(rval);
}


int
get3dgin(ov, ndx, rem, wp)      /* interpolate from warp grid */
W3VEC   ov, rem;
GNDX    ndx;
WARP3D  *wp;
{
        W3VEC   cv[8];
        GNDX    gi;
        int     rval = W3OK;
        register int    i;
                                        /* get corner values */
        for (i = 0; i < 8; i++) {
                gi[0] = ndx[0] + (i & 1);
                gi[1] = ndx[1] + (i>>1 & 1);
                gi[2] = ndx[2] + (i>>2);
                rval |= get3dgpt(cv[i], gi, wp);
        }
        if (rval & W3ERROR)
                return(rval);
        l3interp(ov, cv, rem, 3);       /* perform trilinear interpolation */
        return(rval);
}


l3interp(vo, cl, dv, n)         /* trilinear interpolation (recursive) */
W3VEC   vo, *cl, dv;
int     n;
{
        W3VEC   v0, v1;
        register int    i;

        if (--n) {
                l3interp(v0, cl, dv, n);
                l3interp(v1, cl+(1<<n), dv, n);
        } else {
                W3VCPY(v0, cl[0]);
                W3VCPY(v1, cl[1]);
        }
        for (i = 0; i < 3; i++)
                vo[i] = (1.-dv[n])*v0[i] + dv[n]*v1[i];
}


int
warp3dex(ov, pi, wp)            /* compute warp using 1/r^2 weighted avg. */
W3VEC   ov, pi;
register WARP3D *wp;
{
        double  d2, w, wsum;
        W3VEC   vt;
        register int    i;

        vt[0] = vt[1] = vt[2] = 0.;
        wsum = 0.;
        for (i = wp->npts; i--; ) {
                d2 = wpdist2(pi, wp->ip[i]);
                if (d2 <= MIND*MIND) w = 1./(MIND*MIND);
                else w = 1./d2;
                vt[0] += w*wp->ov[i][0];
                vt[1] += w*wp->ov[i][1];
                vt[2] += w*wp->ov[i][2];
                wsum += w;
        }
        if (wsum > 0.) {
                ov[0] = vt[0]/wsum;
                ov[1] = vt[1]/wsum;
                ov[2] = vt[2]/wsum;
        }
        return(W3OK);
}


WARP3D *
new3dw(flgs)                    /* allocate and initialize WARP3D struct */
int     flgs;
{
        register WARP3D  *wp;

        if ((flgs & (W3EXACT|W3FAST)) == (W3EXACT|W3FAST)) {
                eputs("new3dw: only one of W3EXACT or W3FAST\n");
                return(NULL);
        }
        if ((wp = (WARP3D *)malloc(sizeof(WARP3D))) == NULL) {
                eputs("new3dw: no memory\n");
                return(NULL);
        }
        wp->ip = wp->ov = NULL;
        wp->npts = 0;
        wp->grid.flags = flgs;
        wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 0;
        return(wp);
}


WARP3D *
load3dw(fn, wp)                 /* load 3D warp from file */
char    *fn;
WARP3D  *wp;
{
        FILE    *fp;
        W3VEC   inp, outp;

        if ((fp = fopen(fn, "r")) == NULL) {
                eputs(fn);
                eputs(": cannot open\n");
                return(NULL);
        }
        if (wp == NULL && (wp = new3dw(0)) == NULL)
                goto memerr;
        while (fgetvec(fp, inp) && fgetvec(fp, outp))
                if (!add3dpt(wp, inp, outp))
                        goto memerr;
        if (!feof(fp)) {
                wputs(fn);
                wputs(": non-number in 3D warp file\n");
        }
        goto cleanup;
memerr:
        eputs("load3dw: out of memory\n");
cleanup:
        fclose(fp);
        return(wp);
}


int
set3dwfl(wp, flgs)              /* reset warp flags */
register WARP3D *wp;
int     flgs;
{
        if (flgs == wp->grid.flags)
                return(0);
        if ((flgs & (W3EXACT|W3FAST)) == (W3EXACT|W3FAST)) {
                eputs("set3dwfl: only one of W3EXACT or W3FAST\n");
                return(-1);
        }
        wp->grid.flags = flgs;
        done3dgrid(&wp->grid);          /* old grid is invalid */
        return(0);
}


int
add3dpt(wp, pti, pto)           /* add 3D point pair to warp structure */
register WARP3D *wp;
W3VEC   pti, pto;
{
        double  d2;
        register W3VEC  *na;
        register int    i;

        if (wp->npts == 0) {                    /* initialize */
                wp->ip = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
                if (wp->ip == NULL) return(0);
                wp->ov = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
                if (wp->ov == NULL) return(0);
                wp->d2min = 1e10;
                wp->d2max = 0.;
                W3VCPY(wp->llim, pti);
                W3VCPY(wp->ulim, pti);
        } else {
                if (wp->npts % AHUNK == 0) {    /* allocate another hunk */
                        na = (W3VEC *)realloc((char *)wp->ip,
                                        (wp->npts+AHUNK)*sizeof(W3VEC));
                        if (na == NULL) return(0);
                        wp->ip = na;
                        na = (W3VEC *)realloc((char *)wp->ov,
                                        (wp->npts+AHUNK)*sizeof(W3VEC));
                        if (na == NULL) return(0);
                        wp->ov = na;
                }
                for (i = 0; i < 3; i++)         /* check boundaries */
                        if (pti[i] < wp->llim[i])
                                wp->llim[i] = pti[i];
                        else if (pti[i] > wp->ulim[i])
                                wp->ulim[i] = pti[i];
                for (i = wp->npts; i--; ) {     /* check distances */
                        d2 = wpdist2(pti, wp->ip[i]);
                        if (d2 < MIND*MIND)     /* value is too close */
                                return(wp->npts);
                        if (d2 < wp->d2min)
                                wp->d2min = d2;
                        if (d2 > wp->d2max)
                                wp->d2max = d2;
                }
        }
        W3VCPY(wp->ip[wp->npts], pti);  /* add new point */
        wp->ov[wp->npts][0] = pto[0] - pti[0];
        wp->ov[wp->npts][1] = pto[1] - pti[1];
        wp->ov[wp->npts][2] = pto[2] - pti[2];
        done3dgrid(&wp->grid);          /* old grid is invalid */
        return(++wp->npts);
}


free3dw(wp)                     /* free WARP3D data */
register WARP3D *wp;
{
        done3dgrid(&wp->grid);
        free((void *)wp->ip);
        free((void *)wp->ov);
        free((void *)wp);
}


unsigned long
gridhash(gp)                    /* hash a grid point index */
GNDX    gp;
{
        return(((unsigned long)gp[0]<<GNBITS | gp[1])<<GNBITS | gp[2]);
}


int
new3dgrid(wp)                   /* initialize interpolating grid for warp */
register WARP3D *wp;
{
        W3VEC   gmax;
        double  gridstep, d;
        int     n;
        register int    i;
                                /* free old grid (if any) */
        done3dgrid(&wp->grid);
                                /* check parameters */
        if (wp->npts < 2)
                return(W3BADMAP);       /* undefined for less than 2 points */
        if (wp->d2max < MIND)
                return(W3BADMAP);       /* not enough disparity */
                                /* compute gamut */
        W3VCPY(wp->grid.gmin, wp->llim);
        W3VCPY(gmax, wp->ulim);
        gridstep = sqrt(wp->d2min);
        for (i = 0; i < 3; i++) {
                wp->grid.gmin[i] -= .501*gridstep;
                gmax[i] += .501*gridstep;
        }
        if (wp->grid.flags & W3EXACT) {
                wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 1;
                wp->grid.gstep[0] = gmax[0] - wp->grid.gmin[0];
                wp->grid.gstep[1] = gmax[1] - wp->grid.gmin[1];
                wp->grid.gstep[2] = gmax[2] - wp->grid.gmin[2];
                return(W3OK);           /* no interpolation, so no grid */
        }
                                /* create grid */
        for (i = 0; i < 3; i++) {
                d = gmax[i] - wp->grid.gmin[i];
                n = d/gridstep + .5;
                if (n >= MAXGN-1)
                        n = MAXGN-2;
                wp->grid.gstep[i] = d / n;
                wp->grid.gn[i] = n;
        }
                                /* initialize lookup table */
        wp->grid.gtab.hashf = gridhash;
        wp->grid.gtab.keycmp = NULL;
        wp->grid.gtab.freek = free;
        wp->grid.gtab.freed = NULL;
        return(lu_init(&wp->grid.gtab, 1024) ? W3OK : W3ERROR);
}


done3dgrid(gp)                  /* free interpolating grid for warp */
register struct grid3d  *gp;
{
        if (gp->gn[0] == 0)
                return;
        lu_done(&gp->gtab);
        gp->gn[0] = gp->gn[1] = gp->gn[2] = 0;
}
Revision:	3.5
Committed:	Sat Feb 22 02:07:28 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R5
Changes since 3.4:	+5 -12 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	greg	3.1	#ifndef lint
2	greg	3.5	static const char RCSid[] = "$Id$";
3	greg	3.1	#endif
4			/*
5			* 3D warping routines.
6			*/
7
8			#include <stdio.h>
9	greg	3.5	#include <stdlib.h>
10	greg	3.1	#include <math.h>
11			#include "fvect.h"
12			#include "warp3d.h"
13
14			#define MIND 1e-5 /* minimum distance between input points */
15
16			typedef struct {
17			GNDX n; /* index must be first */
18	greg	3.2	W3VEC v;
19	greg	3.1	} KEYDP; /* key/data allocation pair */
20
21	greg	3.2	#define fgetvec(p,v) (fgetval(p,'f',v) > 0 && fgetval(p,'f',v+1) > 0 \
22			&& fgetval(p,'f',v+2) > 0)
23	greg	3.1
24	greg	3.2	#define AHUNK 24 /* number of points to allocate at a time */
25	greg	3.1
26
27			double
28			wpdist2(p1, p2) /* compute square of distance between points */
29			register W3VEC p1, p2;
30			{
31			double d, d2;
32
33			d = p1[0] - p2[0]; d2 = d*d;
34			d = p1[1] - p2[1]; d2 += d*d;
35			d = p1[2] - p2[2]; d2 += d*d;
36			return(d2);
37			}
38
39
40			int
41			warp3d(po, pi, wp) /* warp 3D point according to the given map */
42			W3VEC po, pi;
43			register WARP3D *wp;
44			{
45			int rval = W3OK;
46			GNDX gi;
47	greg	3.2	W3VEC gd, ov;
48	greg	3.1
49	greg	3.2	if (wp->grid.gn[0] == 0 && (rval = new3dgrid(wp)) != W3OK)
50	greg	3.1	return(rval);
51			rval \|= gridpoint(gi, gd, pi, &wp->grid);
52	greg	3.2	if (wp->grid.flags & W3EXACT)
53			rval \|= warp3dex(ov, pi, wp);
54			else if (wp->grid.flags & W3FAST)
55			rval \|= get3dgpt(ov, gi, wp);
56			else
57			rval \|= get3dgin(ov, gi, gd, wp);
58			po[0] = pi[0] + ov[0];
59			po[1] = pi[1] + ov[1];
60			po[2] = pi[2] + ov[2];
61	greg	3.1	return(rval);
62			}
63
64
65			int
66			gridpoint(ndx, rem, ipt, gp) /* compute grid position for ipt */
67			GNDX ndx;
68			W3VEC rem, ipt;
69			register struct grid3d *gp;
70			{
71			int rval = W3OK;
72			register int i;
73
74			for (i = 0; i < 3; i++) {
75			rem[i] = (ipt[i] - gp->gmin[i])/gp->gstep[i];
76			if (rem[i] < 0.) {
77			ndx[i] = 0;
78			rval = W3GAMUT;
79			} else if ((int)rem[i] >= gp->gn[i]) {
80			ndx[i] = gp->gn[i] - 1;
81			rval = W3GAMUT;
82			} else
83			ndx[i] = (int)rem[i];
84			rem[i] -= (double)ndx[i];
85			}
86			return(rval);
87			}
88
89
90			int
91	greg	3.2	get3dgpt(ov, ndx, wp) /* get value for voxel */
92			W3VEC ov;
93	greg	3.1	GNDX ndx;
94			register WARP3D *wp;
95			{
96			W3VEC gpt;
97			register LUENT *le;
98			KEYDP *kd;
99			int rval = W3OK;
100			register int i;
101
102	greg	3.2	le = lu_find(&wp->grid.gtab, ndx);
103	greg	3.1	if (le == NULL)
104			return(W3ERROR);
105			if (le->data == NULL) {
106			if (le->key != NULL)
107			kd = (KEYDP *)le->key;
108			else if ((kd = (KEYDP *)malloc(sizeof(KEYDP))) == NULL)
109			return(W3ERROR);
110			for (i = 0; i < 3; i++) {
111			kd->n[i] = ndx[i];
112			gpt[i] = wp->grid.gmin[i] + ndx[i]*wp->grid.gstep[i];
113			if (wp->grid.flags & W3FAST) /* on centers */
114			gpt[i] += .5*wp->grid.gstep[i];
115			}
116	greg	3.2	rval = warp3dex(kd->v, gpt, wp);
117	greg	3.1	le->key = (char *)kd->n;
118	greg	3.2	le->data = (char *)kd->v;
119	greg	3.1	}
120	greg	3.2	W3VCPY(ov, (float *)le->data);
121	greg	3.1	return(rval);
122			}
123
124
125			int
126	greg	3.2	get3dgin(ov, ndx, rem, wp) /* interpolate from warp grid */
127			W3VEC ov, rem;
128	greg	3.1	GNDX ndx;
129			WARP3D *wp;
130			{
131			W3VEC cv[8];
132			GNDX gi;
133	greg	3.3	int rval = W3OK;
134	greg	3.1	register int i;
135			/* get corner values */
136			for (i = 0; i < 8; i++) {
137			gi[0] = ndx[0] + (i & 1);
138			gi[1] = ndx[1] + (i>>1 & 1);
139			gi[2] = ndx[2] + (i>>2);
140			rval \|= get3dgpt(cv[i], gi, wp);
141			}
142			if (rval & W3ERROR)
143			return(rval);
144	greg	3.2	l3interp(ov, cv, rem, 3); /* perform trilinear interpolation */
145	greg	3.1	return(rval);
146			}
147
148
149			l3interp(vo, cl, dv, n) /* trilinear interpolation (recursive) */
150			W3VEC vo, *cl, dv;
151			int n;
152			{
153			W3VEC v0, v1;
154			register int i;
155
156			if (--n) {
157			l3interp(v0, cl, dv, n);
158			l3interp(v1, cl+(1<<n), dv, n);
159			} else {
160			W3VCPY(v0, cl[0]);
161			W3VCPY(v1, cl[1]);
162			}
163			for (i = 0; i < 3; i++)
164			vo[i] = (1.-dv[n])v0[i] + dv[n]v1[i];
165			}
166
167
168			int
169	greg	3.2	warp3dex(ov, pi, wp) /* compute warp using 1/r^2 weighted avg. */
170			W3VEC ov, pi;
171	greg	3.1	register WARP3D *wp;
172			{
173			double d2, w, wsum;
174	greg	3.2	W3VEC vt;
175	greg	3.1	register int i;
176
177	greg	3.2	vt[0] = vt[1] = vt[2] = 0.;
178	greg	3.1	wsum = 0.;
179			for (i = wp->npts; i--; ) {
180			d2 = wpdist2(pi, wp->ip[i]);
181			if (d2 <= MINDMIND) w = 1./(MINDMIND);
182			else w = 1./d2;
183	greg	3.2	vt[0] += w*wp->ov[i][0];
184			vt[1] += w*wp->ov[i][1];
185			vt[2] += w*wp->ov[i][2];
186	greg	3.1	wsum += w;
187			}
188			if (wsum > 0.) {
189	greg	3.2	ov[0] = vt[0]/wsum;
190			ov[1] = vt[1]/wsum;
191			ov[2] = vt[2]/wsum;
192	greg	3.1	}
193			return(W3OK);
194			}
195
196
197			WARP3D *
198			new3dw(flgs) /* allocate and initialize WARP3D struct */
199			int flgs;
200			{
201			register WARP3D *wp;
202
203			if ((flgs & (W3EXACT\|W3FAST)) == (W3EXACT\|W3FAST)) {
204			eputs("new3dw: only one of W3EXACT or W3FAST\n");
205			return(NULL);
206			}
207			if ((wp = (WARP3D *)malloc(sizeof(WARP3D))) == NULL) {
208			eputs("new3dw: no memory\n");
209			return(NULL);
210			}
211	greg	3.2	wp->ip = wp->ov = NULL;
212	greg	3.1	wp->npts = 0;
213			wp->grid.flags = flgs;
214			wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 0;
215			return(wp);
216			}
217
218
219			WARP3D *
220			load3dw(fn, wp) /* load 3D warp from file */
221			char *fn;
222			WARP3D *wp;
223			{
224			FILE *fp;
225			W3VEC inp, outp;
226
227			if ((fp = fopen(fn, "r")) == NULL) {
228			eputs(fn);
229			eputs(": cannot open\n");
230			return(NULL);
231			}
232	greg	3.2	if (wp == NULL && (wp = new3dw(0)) == NULL)
233	greg	3.1	goto memerr;
234			while (fgetvec(fp, inp) && fgetvec(fp, outp))
235			if (!add3dpt(wp, inp, outp))
236			goto memerr;
237			if (!feof(fp)) {
238			wputs(fn);
239			wputs(": non-number in 3D warp file\n");
240			}
241			goto cleanup;
242			memerr:
243			eputs("load3dw: out of memory\n");
244			cleanup:
245			fclose(fp);
246			return(wp);
247			}
248
249
250			int
251	greg	3.2	set3dwfl(wp, flgs) /* reset warp flags */
252			register WARP3D *wp;
253			int flgs;
254			{
255			if (flgs == wp->grid.flags)
256			return(0);
257			if ((flgs & (W3EXACT\|W3FAST)) == (W3EXACT\|W3FAST)) {
258			eputs("set3dwfl: only one of W3EXACT or W3FAST\n");
259			return(-1);
260			}
261			wp->grid.flags = flgs;
262			done3dgrid(&wp->grid); /* old grid is invalid */
263			return(0);
264			}
265
266
267			int
268	greg	3.1	add3dpt(wp, pti, pto) /* add 3D point pair to warp structure */
269			register WARP3D *wp;
270			W3VEC pti, pto;
271			{
272			double d2;
273			register W3VEC *na;
274			register int i;
275
276			if (wp->npts == 0) { /* initialize */
277			wp->ip = (W3VEC )malloc(AHUNKsizeof(W3VEC));
278			if (wp->ip == NULL) return(0);
279	greg	3.2	wp->ov = (W3VEC )malloc(AHUNKsizeof(W3VEC));
280			if (wp->ov == NULL) return(0);
281	greg	3.1	wp->d2min = 1e10;
282			wp->d2max = 0.;
283			W3VCPY(wp->llim, pti);
284			W3VCPY(wp->ulim, pti);
285			} else {
286			if (wp->npts % AHUNK == 0) { /* allocate another hunk */
287			na = (W3VEC )realloc((char )wp->ip,
288			(wp->npts+AHUNK)*sizeof(W3VEC));
289			if (na == NULL) return(0);
290			wp->ip = na;
291	greg	3.2	na = (W3VEC )realloc((char )wp->ov,
292	greg	3.1	(wp->npts+AHUNK)*sizeof(W3VEC));
293			if (na == NULL) return(0);
294	greg	3.2	wp->ov = na;
295	greg	3.1	}
296			for (i = 0; i < 3; i++) /* check boundaries */
297			if (pti[i] < wp->llim[i])
298			wp->llim[i] = pti[i];
299			else if (pti[i] > wp->ulim[i])
300			wp->ulim[i] = pti[i];
301			for (i = wp->npts; i--; ) { /* check distances */
302			d2 = wpdist2(pti, wp->ip[i]);
303			if (d2 < MINDMIND) / value is too close */
304			return(wp->npts);
305			if (d2 < wp->d2min)
306			wp->d2min = d2;
307			if (d2 > wp->d2max)
308			wp->d2max = d2;
309			}
310			}
311			W3VCPY(wp->ip[wp->npts], pti); /* add new point */
312	greg	3.2	wp->ov[wp->npts][0] = pto[0] - pti[0];
313			wp->ov[wp->npts][1] = pto[1] - pti[1];
314			wp->ov[wp->npts][2] = pto[2] - pti[2];
315	greg	3.1	done3dgrid(&wp->grid); /* old grid is invalid */
316			return(++wp->npts);
317			}
318
319
320			free3dw(wp) /* free WARP3D data */
321			register WARP3D *wp;
322			{
323			done3dgrid(&wp->grid);
324	greg	3.5	free((void *)wp->ip);
325			free((void *)wp->ov);
326			free((void *)wp);
327	greg	3.1	}
328
329
330	gwlarson	3.4	unsigned long
331	greg	3.1	gridhash(gp) /* hash a grid point index */
332			GNDX gp;
333			{
334			return(((unsigned long)gp[0]<<GNBITS \| gp[1])<<GNBITS \| gp[2]);
335			}
336
337
338			int
339	greg	3.2	new3dgrid(wp) /* initialize interpolating grid for warp */
340	greg	3.1	register WARP3D *wp;
341			{
342	greg	3.2	W3VEC gmax;
343	greg	3.1	double gridstep, d;
344	greg	3.2	int n;
345	greg	3.1	register int i;
346			/* free old grid (if any) */
347			done3dgrid(&wp->grid);
348			/* check parameters */
349			if (wp->npts < 2)
350			return(W3BADMAP); /* undefined for less than 2 points */
351			if (wp->d2max < MIND)
352			return(W3BADMAP); /* not enough disparity */
353			/* compute gamut */
354			W3VCPY(wp->grid.gmin, wp->llim);
355	greg	3.2	W3VCPY(gmax, wp->ulim);
356	greg	3.1	gridstep = sqrt(wp->d2min);
357			for (i = 0; i < 3; i++) {
358			wp->grid.gmin[i] -= .501*gridstep;
359	greg	3.2	gmax[i] += .501*gridstep;
360	greg	3.1	}
361			if (wp->grid.flags & W3EXACT) {
362			wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 1;
363	greg	3.3	wp->grid.gstep[0] = gmax[0] - wp->grid.gmin[0];
364			wp->grid.gstep[1] = gmax[1] - wp->grid.gmin[1];
365			wp->grid.gstep[2] = gmax[2] - wp->grid.gmin[2];
366	greg	3.1	return(W3OK); /* no interpolation, so no grid */
367			}
368			/* create grid */
369			for (i = 0; i < 3; i++) {
370	greg	3.2	d = gmax[i] - wp->grid.gmin[i];
371			n = d/gridstep + .5;
372			if (n >= MAXGN-1)
373			n = MAXGN-2;
374			wp->grid.gstep[i] = d / n;
375			wp->grid.gn[i] = n;
376	greg	3.1	}
377			/* initialize lookup table */
378			wp->grid.gtab.hashf = gridhash;
379			wp->grid.gtab.keycmp = NULL;
380			wp->grid.gtab.freek = free;
381			wp->grid.gtab.freed = NULL;
382			return(lu_init(&wp->grid.gtab, 1024) ? W3OK : W3ERROR);
383			}
384
385
386			done3dgrid(gp) /* free interpolating grid for warp */
387			register struct grid3d *gp;
388			{
389			if (gp->gn[0] == 0)
390			return;
391			lu_done(&gp->gtab);
392			gp->gn[0] = gp->gn[1] = gp->gn[2] = 0;
393			}