/* * C code from the article * "Tri-linear Interpolation" * by Steve Hill, sah@ukc.ac.uk * in "Graphics Gems IV", Academic Press, 1994 * * compile with "cc -DMAIN ..." to create a test program */ #include #include "GraphicsGems.h" double trilinear(p, d, xsize, ysize, zsize, def) Point3 *p; double *d; int xsize, ysize, zsize; double def; { # define DENS(X, Y, Z) d[(X)+xsize*((Y)+ysize*(Z))] int x0, y0, z0, x1, y1, z1; double *dp, fx, fy, fz, d000, d001, d010, d011, d100, d101, d110, d111, dx00, dx01, dx10, dx11, dxy0, dxy1, dxyz; x0 = floor(p->x); fx = p->x - x0; y0 = floor(p->y); fy = p->y - y0; z0 = floor(p->z); fz = p->z - z0; x1 = x0 + 1; y1 = y0 + 1; z1 = z0 + 1; if (x0 >= 0 && x1 < xsize && y0 >= 0 && y1 < ysize && z0 >= 0 && z1 < zsize) { dp = &DENS(x0, y0, z0); d000 = dp[0]; d100 = dp[1]; dp += xsize; d010 = dp[0]; d110 = dp[1]; dp += xsize*ysize; d011 = dp[0]; d111 = dp[1]; dp -= xsize; d001 = dp[0]; d101 = dp[1]; } else { # define INRANGE(X, Y, Z) \ ((X) >= 0 && (X) < xsize && \ (Y) >= 0 && (Y) < ysize && \ (Z) >= 0 && (Z) < zsize) d000 = INRANGE(x0, y0, z0) ? DENS(x0, y0, z0) : def; d001 = INRANGE(x0, y0, z1) ? DENS(x0, y0, z1) : def; d010 = INRANGE(x0, y1, z0) ? DENS(x0, y1, z0) : def; d011 = INRANGE(x0, y1, z1) ? DENS(x0, y1, z1) : def; d100 = INRANGE(x1, y0, z0) ? DENS(x1, y0, z0) : def; d101 = INRANGE(x1, y0, z1) ? DENS(x1, y0, z1) : def; d110 = INRANGE(x1, y1, z0) ? DENS(x1, y1, z0) : def; d111 = INRANGE(x1, y1, z1) ? DENS(x1, y1, z1) : def; } dx00 = LERP(fx, d000, d100); dx01 = LERP(fx, d001, d101); dx10 = LERP(fx, d010, d110); dx11 = LERP(fx, d011, d111); dxy0 = LERP(fy, dx00, dx10); dxy1 = LERP(fy, dx01, dx11); dxyz = LERP(fz, dxy0, dxy1); return dxyz; } #ifdef MAIN /* test program for trilinear interpolation */ #include #define XSIZE 2 #define YSIZE 2 #define ZSIZE 2 main() { #define TDENS(X,Y,Z) d[(X)+XSIZE*((Y)+YSIZE*(Z))] extern double trilinear(); double *d, def; int x, y, z; Point3 p; d = (double *) malloc(sizeof(double) * XSIZE * YSIZE * ZSIZE); printf("Test for trilinear interpolation\n"); printf("Enter the densities for the corners of a cube\n"); for (x = 0; x < XSIZE; x++) for (y = 0; y < YSIZE; y++) for (z = 0; z < ZSIZE; z++) { printf("Point (%d, %d, %d) is: ", x, y, z); scanf("%lf", &TDENS(x, y, z)); } printf("Enter the default density: "); scanf("%lf", &def); printf("Enter point of interest: "); while (scanf("%lf %lf %lf", &p.x, &p.y, &p.z) == 3) { double n; n = trilinear(&p, d, XSIZE, YSIZE, ZSIZE, def); printf("Density at (%lf, %lf, %lf) is %lf\n", p.x, p.y, p.z, n); printf("Enter point of interest: "); } } #endif