--- ray/src/px/warp3d.c 1997/02/04 16:03:48 3.1 +++ ray/src/px/warp3d.c 2004/05/25 22:04:14 3.8 @@ -1,42 +1,49 @@ -/* Copyright (c) 1997 Regents of the University of California */ - #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: warp3d.c,v 3.8 2004/05/25 22:04:14 greg Exp $"; #endif - /* * 3D warping routines. */ #include - +#include #include +#include "rterror.h" +#include "rtio.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 p; + W3VEC v; } KEYDP; /* key/data allocation pair */ -#define fgetvec(f,v) (fgetval(f,'f',v) > 0 && fgetval(f,'f',v+1) > 0 \ - && fgetval(f,'f',v+2) > 0) +#define fgetvec(p,v) (fgetval(p,'f',v) > 0 && fgetval(p,'f',v+1) > 0 \ + && fgetval(p,'f',v+2) > 0) -#define AHUNK 32 /* number of points to allocate at a time */ +#define AHUNK 24 /* number of points to allocate at a time */ -#ifndef malloc -extern char *malloc(), *realloc(); -#endif -extern void free(); +double wpdist2(W3VEC p1, W3VEC p2); +static int gridpoint(GNDX ndx, W3VEC rem, W3VEC ipt, struct grid3d *gp); +static int get3dgpt(W3VEC ov, GNDX ndx, WARP3D *wp); +static int get3dgin(W3VEC ov, GNDX ndx, W3VEC rem, WARP3D *wp); +static void l3interp(W3VEC vo, W3VEC *cl, W3VEC dv, int n); +static int warp3dex(W3VEC ov, W3VEC pi, WARP3D *wp); +//static unsigned long gridhash(const void *gp); +static lut_hashf_t gridhash; +static int new3dgrid(WARP3D *wp); +static void done3dgrid(struct grid3d *gp); + double -wpdist2(p1, p2) /* compute square of distance between points */ -register W3VEC p1, p2; +wpdist2( /* compute square of distance between points */ + register W3VEC p1, + register W3VEC p2 +) { double d, d2; @@ -48,35 +55,39 @@ register W3VEC p1, p2; int -warp3d(po, pi, wp) /* warp 3D point according to the given map */ -W3VEC po, pi; -register WARP3D *wp; +warp3d( /* warp 3D point according to the given map */ + W3VEC po, + W3VEC pi, + register WARP3D *wp +) { int rval = W3OK; GNDX gi; - W3VEC gd; + W3VEC gd, ov; - if (wp->grid.gn[0] == 0 && (rval = new3dwgrid(wp)) != W3OK) + 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(po, pi, wp); - return(rval); - } - if (wp->grid.flags & W3FAST) { - rval |= get3dgpt(po, gi, wp); - return(rval); - } - rval |= get3dgin(po, gi, gd, wp); + 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; +static int +gridpoint( /* compute grid position for ipt */ + GNDX ndx, + W3VEC rem, + W3VEC ipt, + register struct grid3d *gp +) { int rval = W3OK; register int i; @@ -97,11 +108,12 @@ register struct grid3d *gp; } -int -get3dgpt(po, ndx, wp) /* get value for voxel */ -W3VEC po; -GNDX ndx; -register WARP3D *wp; +static int +get3dgpt( /* get value for voxel */ + W3VEC ov, + GNDX ndx, + register WARP3D *wp +) { W3VEC gpt; register LUENT *le; @@ -109,7 +121,7 @@ register WARP3D *wp; int rval = W3OK; register int i; - le = lu_find(&wp->grid.gtab, (char *)ndx); + le = lu_find(&wp->grid.gtab, ndx); if (le == NULL) return(W3ERROR); if (le->data == NULL) { @@ -123,24 +135,26 @@ register WARP3D *wp; if (wp->grid.flags & W3FAST) /* on centers */ gpt[i] += .5*wp->grid.gstep[i]; } - rval |= warp3dex(kd->p, gpt, wp); + rval = warp3dex(kd->v, gpt, wp); le->key = (char *)kd->n; - le->data = (char *)kd->p; + le->data = (char *)kd->v; } - W3VCPY(po, (float *)le->data); + W3VCPY(ov, (float *)le->data); return(rval); } -int -get3dgin(po, ndx, rem, wp) /* interpolate from warp grid */ -W3VEC po, rem; -GNDX ndx; -WARP3D *wp; +static int +get3dgin( /* interpolate from warp grid */ + W3VEC ov, + GNDX ndx, + W3VEC rem, + WARP3D *wp +) { W3VEC cv[8]; - int rval; GNDX gi; + int rval = W3OK; register int i; /* get corner values */ for (i = 0; i < 8; i++) { @@ -151,14 +165,18 @@ WARP3D *wp; } if (rval & W3ERROR) return(rval); - l3interp(po, cv, rem, 3); /* perform trilinear interpolation */ + l3interp(ov, cv, rem, 3); /* perform trilinear interpolation */ return(rval); } -l3interp(vo, cl, dv, n) /* trilinear interpolation (recursive) */ -W3VEC vo, *cl, dv; -int n; +static void +l3interp( /* trilinear interpolation (recursive) */ + W3VEC vo, + W3VEC *cl, + W3VEC dv, + int n +) { W3VEC v0, v1; register int i; @@ -175,38 +193,41 @@ int n; } -int -warp3dex(po, pi, wp) /* compute warp using 1/r^2 weighted avg. */ -W3VEC po, pi; -register WARP3D *wp; +static int +warp3dex( /* compute warp using 1/r^2 weighted avg. */ + W3VEC ov, + W3VEC pi, + register WARP3D *wp +) { double d2, w, wsum; - W3VEC pt; + W3VEC vt; register int i; - pt[0] = pt[1] = pt[2] = 0.; + 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; - pt[0] += w*wp->op[i][0]; - pt[1] += w*wp->op[i][1]; - pt[2] += w*wp->op[i][2]; + 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.) { - po[0] = pt[0]/wsum; - po[1] = pt[1]/wsum; - po[2] = pt[2]/wsum; + 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; +new3dw( /* allocate and initialize WARP3D struct */ + int flgs +) { register WARP3D *wp; @@ -218,7 +239,7 @@ int flgs; eputs("new3dw: no memory\n"); return(NULL); } - wp->ip = wp->op = 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; @@ -227,9 +248,10 @@ int flgs; WARP3D * -load3dw(fn, wp) /* load 3D warp from file */ -char *fn; -WARP3D *wp; +load3dw( /* load 3D warp from file */ + char *fn, + WARP3D *wp +) { FILE *fp; W3VEC inp, outp; @@ -239,7 +261,7 @@ WARP3D *wp; eputs(": cannot open\n"); return(NULL); } - if (wp == NULL && (wp = new3dw()) == NULL) + if (wp == NULL && (wp = new3dw(0)) == NULL) goto memerr; while (fgetvec(fp, inp) && fgetvec(fp, outp)) if (!add3dpt(wp, inp, outp)) @@ -257,10 +279,30 @@ cleanup: } +extern int +set3dwfl( /* 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; +add3dpt( /* add 3D point pair to warp structure */ + register WARP3D *wp, + W3VEC pti, + W3VEC pto +) { double d2; register W3VEC *na; @@ -269,22 +311,22 @@ W3VEC pti, pto; if (wp->npts == 0) { /* initialize */ wp->ip = (W3VEC *)malloc(AHUNK*sizeof(W3VEC)); if (wp->ip == NULL) return(0); - wp->op = (W3VEC *)malloc(AHUNK*sizeof(W3VEC)); - if (wp->op == 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, + na = (W3VEC *)realloc((void *)wp->ip, (wp->npts+AHUNK)*sizeof(W3VEC)); if (na == NULL) return(0); wp->ip = na; - na = (W3VEC *)realloc((char *)wp->op, + na = (W3VEC *)realloc((void *)wp->ov, (wp->npts+AHUNK)*sizeof(W3VEC)); if (na == NULL) return(0); - wp->op = na; + wp->ov = na; } for (i = 0; i < 3; i++) /* check boundaries */ if (pti[i] < wp->llim[i]) @@ -302,35 +344,45 @@ W3VEC pti, pto; } } W3VCPY(wp->ip[wp->npts], pti); /* add new point */ - W3VCPY(wp->op[wp->npts], pto); + 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; +extern void +free3dw( /* free WARP3D data */ + register WARP3D *wp +) { done3dgrid(&wp->grid); - free((char *)wp->ip); - free((char *)wp->op); - free((char *)wp); + free((void *)wp->ip); + free((void *)wp->ov); + free((void *)wp); } -long -gridhash(gp) /* hash a grid point index */ -GNDX gp; +static unsigned long +gridhash( /* hash a grid point index */ + //GNDX gp + const void *gp +) { - return(((unsigned long)gp[0]<grid); @@ -341,23 +393,27 @@ register WARP3D *wp; return(W3BADMAP); /* not enough disparity */ /* compute gamut */ W3VCPY(wp->grid.gmin, wp->llim); - W3VCPY(wp->grid.gmax, wp->ulim); + W3VCPY(gmax, wp->ulim); gridstep = sqrt(wp->d2min); for (i = 0; i < 3; i++) { wp->grid.gmin[i] -= .501*gridstep; - wp->grid.gmax[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 = wp->grid.gmax[i] - wp->grid.gmin[i]; - wp->grid.gn[i] = d/gridstep + .5; - if (wp->grid.gn[i] >= MAXGN) - wp->grid.gn[i] = MAXGN-1; - wp->grid.gstep[i] = d / wp->grid.gn[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; @@ -368,8 +424,10 @@ register WARP3D *wp; } -done3dgrid(gp) /* free interpolating grid for warp */ -register struct grid3d *gp; +static void +done3dgrid( /* free interpolating grid for warp */ + register struct grid3d *gp +) { if (gp->gn[0] == 0) return;