--- ray/src/px/warp3d.c 1997/02/04 16:03:48 3.1 +++ ray/src/px/warp3d.c 2024/02/22 17:45:27 3.12 @@ -1,42 +1,47 @@ -/* 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.12 2024/02/22 17:45:27 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 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 */ + W3VEC p1, + W3VEC p2 +) { double d, d2; @@ -48,38 +53,42 @@ 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, + 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, + struct grid3d *gp +) { int rval = W3OK; - register int i; + int i; for (i = 0; i < 3; i++) { rem[i] = (ipt[i] - gp->gmin[i])/gp->gstep[i]; @@ -97,19 +106,20 @@ 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, + WARP3D *wp +) { W3VEC gpt; - register LUENT *le; + LUENT *le; KEYDP *kd; int rval = W3OK; - register int i; + 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,25 +133,27 @@ 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; - register int i; + int rval = W3OK; + int i; /* get corner values */ for (i = 0; i < 8; i++) { gi[0] = ndx[0] + (i & 1); @@ -151,17 +163,21 @@ 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; + int i; if (--n) { l3interp(v0, cl, dv, n); @@ -175,40 +191,43 @@ 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, + WARP3D *wp +) { double d2, w, wsum; - W3VEC pt; - register int i; + W3VEC vt; + 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; + WARP3D *wp; if ((flgs & (W3EXACT|W3FAST)) == (W3EXACT|W3FAST)) { eputs("new3dw: only one of W3EXACT or W3FAST\n"); @@ -218,7 +237,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 +246,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 +259,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,34 +277,54 @@ cleanup: } +extern int +set3dwfl( /* reset warp flags */ + 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 */ + WARP3D *wp, + W3VEC pti, + W3VEC pto +) { double d2; - register W3VEC *na; - register int i; + W3VEC *na; + int i; 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,36 +342,44 @@ 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 */ + WARP3D *wp +) { done3dgrid(&wp->grid); - free((char *)wp->ip); - free((char *)wp->op); - free((char *)wp); + free(wp->ip); + free(wp->ov); + free(wp); } -long -gridhash(gp) /* hash a grid point index */ -GNDX gp; +static unsigned long +gridhash( /* hash a grid point index */ + const char *gp +) { - return(((unsigned long)gp[0]<grid); /* check parameters */ @@ -341,23 +389,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 +420,10 @@ register WARP3D *wp; } -done3dgrid(gp) /* free interpolating grid for warp */ -register struct grid3d *gp; +static void +done3dgrid( /* free interpolating grid for warp */ + struct grid3d *gp +) { if (gp->gn[0] == 0) return;