#pragma segment REAL3D /* This file assumes 4 char tabs. Memory related functions are not included in this file. See MacMemory.c */ /* Implementation of Marching Cubes algorithm described in SIGGRAPH 87 proceedings, "Marching Cubes: A High Resolution 3D Surface Contruction Algorithm", by William Lorensen and Harvey Cline. Written by Daniel LaLiberte, Software Development Group, NCSA, September 1989. mcube is a 3D contour program. Call it with a pointer to data, the size of each dimension, the ranges of data to plot for each dimension (0 origin), a pointer to an array of values to find, and the number of values. The data and values are float32 type. The structure of the data is not simply the 3D data. It must include pointers to each X plane and Y line followed by the data. Use td_Malloc3Dfloat32 in td_mem.c to allocate and initialize the pointers for the data memory. Then store your data at &data[0][0][0]. This implementation subdivides cubes into 5 tetrahedrons, which have far fewer cases for how to triangulate then whole cubes. For more implementation notes, see comments before data structures and functions. mcube uses NGL (like Iris GL) to draw line segments of the contour, but you can easily substitute the Nmove and Ndraw in triangTetra and plotLine, to port it to another environment. Also, if you are not on a Mac, change the following defines. */ /* Allocate and dispose of memory. Use malloc and free on Unix. */ #define ALLOCPTR AllocPtr #define DISPOSPTR DisposPtr /* Exit to the calling program. Use exit or abort on Unix. */ #define MacExit ExitToShell /* Test if mouse button is pressed, for exit from mcube. */ #define MacButton Button #include "MIRG.h" #include "Ngl.h" /* the following are all that is needed from the hdf includes */ typedef long int32; typedef float float32; /* 3D array of data, declared globally for access by several functions in this file. */ static float32 ***data; /* indexing and data access - depend on x, y, and z being set */ #define D(x,y,z) (data[x][y][z]) /* test if an integer is even */ #define even(a) (1 & (~a)) /* The coordinates, vertices, and edges of a cube are numbered as follows: xyz 001 4--7--7 101 z /| /| / 8 4 11 | ->x 000 0---3-3..|..100 | | | | 6 y 0 | | | 011..|..5-5|--6 111 | / 2 / |9 |10 010 1--1--2 110 Six more cube edges are the diagonals of each face: edge vertices 12 0123 13 4567 14 0154 15 3267 16 1265 17 0374 The diagonal lines may be oriented in one of two ways. This is because, while marching through the cubes, we want to share diagonal lines between neighboring cubes. The center tetrahedron for one cube might have opposite edges 0-5 and 2-7, and the next cube would have edges 1-4 and 3-6. */ /* Table of coordinates of each vertex, in terms of 0 or 1 increment, x, y, z. This table could be eliminated if coordinates were numbered in binary order. */ typedef int Coordinate[3]; static Coordinate cubeVertices[8] = { {0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {1, 0, 0}, {0, 0, 1}, {0, 1, 1}, {1, 1, 1}, {1, 0, 1}}; /* these are not used */ #define D0 D(x,y,z) #define D1 D(x,y+1,z) #define D2 D(x+1,y+1,z) #define D3 D(x+1,y,z) #define D4 D(x,y,z+1) #define D5 D(x,y+1,z+1) #define D6 D(x+1,y+1,z+1) #define D7 D(x+1,y,z+1) /* Table of vertices for each edge. Vertices of each edge are listed from lowest coordinate to highest, for each dimension, x, y, z, so that interpolation is done the same way each time. Since diagonals can be oriented one of two ways, two tables are needed. */ typedef int Edge[2]; static Edge *cubeEdges; /* pointer to either cubeEdgesEven or cubeEdgesOdd */ static int cubeEdgesEven[18][2] = { {0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6}, {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}, {0, 2}, /* 12 */ {5, 7}, {0, 5}, {7, 2}, {5, 2}, {0, 7} }; static int cubeEdgesOdd[18][2] = { {0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6}, {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}, {1, 3}, /* 12 */ {4, 6}, {4, 1}, {3, 6}, {1, 6}, {4, 3} }; /* Table of vertices of tetrahedron edges, indexed by tetrahedron edge. Vertices of a tetrahedron are numbered as follows: 0-4--2 |. /| 0 ./ 2 | 1. | |/ 3| 1-5--3 where edge 0-3 is in the background of the tetrahedron. Notice: use left-hand rule to orient tetrahedron. */ static int tetraVertices[6][2] = { {0, 1}, /* edge 0 */ {1, 2}, {2, 3}, {0, 3}, {0, 2}, {1, 3}}; /* Table of 5 tetrahedrons in a cube; each row is the list of cube vertices. The order is left-handed, i.e. the first three vertices point to the last. */ typedef TetraVertices[4]; static TetraVertices *cubeTetraVertices; static int cubeTetraVerticesEven[5][4] = { {0, 1, 2, 5}, {2, 5, 6, 7}, {0, 2, 3, 7}, {0, 7, 4, 5}, {0, 5, 2, 7}}; static int cubeTetraVerticesOdd[5][4] = { {0, 1, 3, 4}, {1, 2, 3, 6}, {3, 4, 6, 7}, {1, 4, 5, 6}, {1, 3, 4, 6}}; /* Table of 5 tetrahedrons in a cube; each row is the list of cube edges. The order corresponds to tetrahedron edge order relative to cubeTetraVertices. */ typedef TetraEdges[6]; static TetraEdges *cubeTetraEdges; static int cubeTetraEdgesEven[5][6] = { {0, 1, 16, 14, 12, 9}, {16, 5, 6, 15, 10, 13}, {12, 2, 11, 17, 3, 15}, {17, 7, 4, 14, 8, 13}, {14, 16, 15, 17, 12, 13}}; static int cubeTetraEdgesOdd[5][6] = { {0, 12, 17, 8, 3, 14}, {1, 2, 15, 16, 12, 10}, {17, 13, 6, 11, 15, 7}, {14, 4, 5, 16, 9, 13}, {12, 17, 13, 16, 14, 15}}; /* Table of triangles that intersect a tetrahedron for each of 16 cases. Since a tetrahedron has four vertices, and each vertex is either higher or lower (or equal) to the value being tested, there are 16 cases of how triangles intersect. These reduce (by rotation and inversion) to three basic cases: 0) no intersection, 1) one triangle that intersects three adjacent faces, and 2) two triangles that intersect all faces but miss two opposite edges. An example of opposite edges are 0-3 and 1-2 in the above tetrahedron diagram. Each row in the table has seven elements; the first is the number of triangles that follow, the next two triples are the edges of the tetrahedron that the triangles intersect. The vertices of each triangle are listed in right-hand order (counterclockwise), pointing towards the high value, for use by polygon shading programs. Note that this is not necessarily the same as pointing "out" since we dont know which way is out unless we look at the whole object - the object may be inside out if the high values are on the inside. */ static int triTetra[16][7] = { /* 0123 - vertices of tetrahedron */ {0, 0, 0, 0, 0, 0, 0}, /* 0000 - all vertices are lower than value */ {1, 2, 5, 3, 0, 0, 0}, /* 0001 - only vertex 3 is higher */ {1, 1, 2, 4, 0, 0, 0}, /* 0010 */ {2, 3, 4, 5, 5, 4, 1}, /* 0011 */ {1, 0, 5, 1, 0, 0, 0}, /* 0100 */ {2, 1, 0, 2, 2, 0, 3}, /* 0101 */ {2, 0, 5, 4, 4, 5, 2}, /* 0110 */ {1, 0, 3, 4, 0, 0, 0}, /* 0111 */ {1, 0, 4, 3, 0, 0, 0}, /* 1000 */ {2, 0, 4, 5, 5, 4, 2}, /* 1001 */ {2, 1, 2, 0, 0, 2, 3}, /* 1010 */ {1, 0, 1, 5, 0, 0, 0}, /* 1011 */ {2, 3, 5, 4, 4, 5, 1}, /* 1100 */ {1, 1, 4, 2, 0, 0, 0}, /* 1101 */ {1, 2, 3, 5, 0, 0, 0}, /* 1110 */ {0, 0, 0, 0, 0, 0, 0}}; /* 1111 */ /* Interpolation data for the face of a cube is a triple of three values, one each for the horizontal, vertical, and diagonal lines labeled below as a-b, a-c, a-d. The diagonal line might be b-c instead. (0,0) (0,1) a----b |\ | | \ | | \ | | \| c----d Each interpolated point is stored as the (x,y,z) coordinate */ typedef struct { int good; /* non-zero if value is valid */ float32 x, y, z; } InterpCoor; typedef struct { InterpCoor h; /* horizontal */ InterpCoor v; /* vertical */ InterpCoor d; /* diagonal, either / or \, but not both */ } InterpData; /* Interpolation data is stored to avoid recomputation. */ static InterpData ** zFace, /* Previous z plane of faces */ * yFace, /* Previous y line of faces */ xFace; /* Previous x face */ /* Booleans that indicates if currently indexing on the edge of the data */ static int onXedge, onYedge, onZedge; /* The interpolation function is a macro. Return the fraction of the distance between a and b that c is found. If the result is negative, then c is less than a; if greater than 1, then c is greater than b. */ static float32 aSave, baSave; /* for use by interp macro */ static float32 dataValue; /* the value to interpolate at */ #define interp(a, b) (aSave=a, baSave=b-aSave, \ (baSave ? (((float32)dataValue-aSave) / baSave) : -1 )) /* for debugging, same as interp but no intermediate storage */ #define dinterp(a, b) (((float32)dataValue-a) / (b-a)) static mcubeOneValue( float32 value, int xFirst, int xLast, /* range of data to plot, for each dimension */ int yFirst, int yLast, int zFirst, int zLast); static triangCube(void); static triangTetra( int E[6], /* edges of tetrahedron */ int V[4]); /* vertices .. */ static void (*plotFunc)(); /* function to use to plot a triangle */ static void fillTriangle( InterpCoor p[3] ); static void plotTriangle( InterpCoor p[3] ); extern MacSetGELineHook(void (* hook)()); extern MacSetGEPolyHook(void (* hook)()); extern void fillDepthQue(char *); extern void plotDepthQue(char *); /*==================================================================*/ void mcube(float32 ***d, /* 3D array of data */ int32 dims[3], /* size of each dimension of d */ int xi, int xj, /* range of data to plot, for each dimension */ int yi, int yj, int zi, int zj, float32 values[], /* list of values plot surface of */ int numvalues, /* length of values */ float32 minData, float32 maxData, Boolean fill /* true if polygons should be filled */ ) { int i; float32 denom; data = d; denom = maxData - minData; plotFunc = fill ? fillTriangle : plotTriangle; if (fill) MacSetGEPolyHook(fillDepthQue); else MacSetGELineHook(NULL); /* allocate interpolation data arrays */ yFace = NULL; zFace = NULL; yFace = (InterpData *) ALLOCPTR(dims[0] * sizeof(InterpData)); if (!yFace) goto error; if (MacButton()) goto done; zFace = (InterpData **) Alloc2D(dims[0], dims[1], sizeof(InterpData)); if (!zFace) goto error; PenMode(srcOr); for (i = 0; i < numvalues; i++) { long colorIndex = (1 + (254 * (values[i] - minData)) / denom); /* RGBColor color; PixPatHandle ppat; Index2Color(colorIndex, &color); ppat = NewPixPat(); MakeRGBPat(ppat, &color); PenPixPat(ppat); */ Ncolor((short) colorIndex); if (MacButton()) break; mcubeOneValue(values[i], xi, xj, yi, yj, zi, zj); /* DisposPixPat(ppat); RGBForeColor(&color);*/ } goto done; /* free space of yFace and zFace */ error: ErrorAlert(GetResource('STR ',774)); /* should be a define and a more appropriate message!! */ done: if (yFace) DISPOSPTR((Ptr)yFace); if (zFace) DISPOSPTR((Ptr)zFace); PenPat(qd.black); PenMode(srcCopy); MacSetGELineHook(NULL); MacSetGEPolyHook(NULL); } /* Global vars for index into data */ static int x, y, z; /* Global var to store whether each corner of the cube is higher than the value */ static int cubecase[8]; static mcubeOneValue(value, xFirst, xLast, yFirst, yLast, zFirst, zLast) float32 value; int xFirst, xLast, /* range of data to plot, for each dimension */ yFirst, yLast, zFirst, zLast; { /* set global dataValue for use by interp */ dataValue = value; /* initialize the z=0 plane of interpolated values. This plane has size (xLast-xFirst)+1 by (yLast-yFirst)+1. The edges with # below will be initialized. 4----7 /| /| 0####3 | # | # | # | # | # 5--#-6 #/ #/ 1####2 */ z = 0; for (x = xFirst; x <= xLast; x++) { for (y = yFirst; y <= yLast; y++) { if (MacButton()) break; zFace[x][y].h.good = 0; zFace[x][y].v.good = 0; zFace[x][y].d.good = 0; }}; /* loop through each z plane.*/ for (z = zFirst; z < zLast; z++) { onZedge = (z == zFirst || z == (zLast-1)); if (MacButton()) break; mSpinCursor(); /* initialize the y=0, z=z line. The edges with # below will be initialized. 4####7 #| #| 0----3 | | | | | | | | | | 5--|-6 |/ |/ 1----2 */ y = 0; for (x = xFirst; x <= xLast; x++) { yFace[x].h.good = 0; yFace[x].v.good = 0; yFace[x].d.good = 0; }; /* loop through each y line */ for (y = yFirst; y < yLast; y++) { onYedge = (y == yFirst || y == (yLast-1)); if (MacButton()) break; /* initialize the x=0, y=y, z=z point The edges with # below will be initialized. 4----7 /# /| 0----3 | | # | | | # | | | 5--|-6 |# |/ 1----2 */ x = 0; xFace.h.good = 0; xFace.v.good = 0; xFace.d.good = 0; /* loop through each x value */ for (x = xFirst; x < xLast; x++) { int i; int found; if (MacButton()) break; /* test whether there is any intersection in this cube */ found = 0; cubecase[0] = D(x+cubeVertices[0][0], y+cubeVertices[0][1], z+cubeVertices[0][2]) > dataValue; for (i=1; i<8; i++) { cubecase[i] = D(x+cubeVertices[i][0], y+cubeVertices[i][1], z+cubeVertices[i][2]) > dataValue; found |= (cubecase[i] != cubecase[i-1]); }; if (found) { onXedge = (x == xFirst || x == (xLast-1)); triangCube(); } else { /* skipping a cube, so most values we would save are lost */ xFace.h.good = 0; xFace.v.good = 0; xFace.d.good = 0; yFace[x].h.good = 0; yFace[x].v = xFace.h; yFace[x].d.good = 0; zFace[x][y].h = yFace[x].h; zFace[x][y].v = xFace.v; zFace[x][y].d.good = 0; } } /* reinit some more values after x iteration is done */ yFace[x].v.good = 0; zFace[x][y].v.good = 0; } /* last row of zFace needs to be reinitialized */ y = yLast; for (x = xFirst; x < xLast; x++) zFace[x][y].h.good = 0; } } /* global var to hold interpolation coordinates for current cube. Used by triangCube and triangTetra. */ InterpCoor cubeD[18]; /* 12 edges of cube + 6 diagonal edges */ /* triangCube - given the index of a cube, determine triangulation of the five tetrahedrons that compose it. One tetra is in the center of the cube, and the other four are on the faces of the first. The index points to the vertex closest to the (0,0,0) vertex. The interpolation points for previously visited edges may be found in zFace, yFace, and xFace. */ static triangCube() { /* find old faces */ /* could avoid copying if value isnt valid, ie good=0 */ cubeD[9] = xFace.h; cubeD[4] = xFace.v; cubeD[14] = xFace.d; cubeD[7] = yFace[x].h; cubeD[8] = yFace[x].v; cubeD[17] = yFace[x].d; cubeD[3] = zFace[x][y].h; cubeD[0] = zFace[x][y].v; cubeD[12] = zFace[x][y].d; cubeD[1] = zFace[x][y+1].h; cubeD[2] = zFace[x+1][y].v; cubeD[11] = yFace[x+1].v; /* interpolate the new edges. The edges with # below will be initialized. 4----7 /| /# 0----3 # | | | # -- B | | | # | 5##|#6 -- C |/ |# -- A 1----2 */ cubeD[10].good = 0; /* A */ cubeD[6].good = 0; /* B */ cubeD[5].good = 0; /* C */ cubeD[13].good = 0; cubeD[15].good = 0; cubeD[16].good = 0; if (even(x+y+z)) { cubeEdges = (Edge *)&cubeEdgesEven; cubeTetraEdges = (TetraEdges *)&cubeTetraEdgesEven; cubeTetraVertices = (TetraVertices *)&cubeTetraVerticesEven; } else { cubeEdges = (Edge *)&cubeEdgesOdd; cubeTetraEdges = (TetraEdges *)&cubeTetraEdgesOdd; cubeTetraVertices = (TetraVertices *)&cubeTetraVerticesOdd; }; /* triangulate each tetrahedron in the cube. */ { int i; for (i = 0; i < 5; i++) triangTetra(cubeTetraEdges[i], cubeTetraVertices[i]); } /* save values for future cubes */ xFace.h = cubeD[10]; xFace.v = cubeD[6]; xFace.d = cubeD[15]; yFace[x].h = cubeD[5]; yFace[x].v = cubeD[9]; yFace[x].d = cubeD[16]; zFace[x][y].h = cubeD[7]; zFace[x][y].v = cubeD[4]; zFace[x][y].d = cubeD[13]; }; /* last point moved to */ static InterpCoor lastPt = {0, -1, -1, -1 /* non-existent point */}; static plotLine( InterpCoor p1, InterpCoor p2); static triangTetra(E, V) int E[6]; /* edges of tetrahedron */ int V[4]; /* vertices .. */ { int i, j; Coordinate *c1, *c2; /* pointers to coordinates of cube edge endpoints */ int *tri; /* row of TriTetra */ int numTriangles; InterpCoor p[3]; /* three points of interpolated triangle */ int edge; /* first determine which of 16 cases we have */ { /* register int sum; register int factor; sum = 0; factor = 8; for (i = 0; i < 4; i++) { sum += (cubecase[V[i]]) * factor; factor = factor >> 1; }; */ /* unwrap previous loop for speed */ tri = triTetra[(cubecase[V[0]] << 3) + (cubecase[V[1]] << 2) + (cubecase[V[2]] << 1) + (cubecase[V[3]])]; }; /* for each triangle, find intersection points and draw segments between them */ numTriangles = *tri; tri++; for (i = 0; i < numTriangles; i++) { /* for each intersecting edge */ for (j = 0; j < 3; j++) { register InterpCoor *q; /* where to store interpolation data */ float32 dist; edge = *(tri++); q = &cubeD[E[edge]]; if (!q->good) { q->good = 1; /* get the coordinates of endpoints */ c1 = (Coordinate *) &cubeVertices[cubeEdges[E[edge]][0]]; c2 = (Coordinate *) &cubeVertices[cubeEdges[E[edge]][1]]; dist = interp(D(x+(*c1)[0], y+(*c1)[1], z+(*c1)[2]), D(x+(*c2)[0], y+(*c2)[1], z+(*c2)[2])); /* find intersection point along the edge */ q->x = x + (*c1)[0] + ((*c2)[0] - (*c1)[0]) * dist; q->y = y + (*c1)[1] + ((*c2)[1] - (*c1)[1]) * dist; q->z = z + (*c1)[2] + ((*c2)[2] - (*c1)[2]) * dist; }; p[j] = cubeD[E[edge]]; }; plotFunc(p); } } #define x1 tri[0][0] #define y1 tri[0][1] #define z1 tri[0][2] #define x2 tri[1][0] #define y2 tri[1][1] #define z2 tri[1][2] #define x3 tri[2][0] #define y3 tri[2][1] #define z3 tri[2][2] static void fillTriangle( InterpCoor p[3] /* three points of interpolated triangle */ ) { float tri[3][3]; if (((p[0].x == p[1].x) && (p[0].y == p[1].y) && (p[0].z == p[1].z)) || ((p[0].x == p[2].x) && (p[0].y == p[2].y) && (p[0].z == p[2].z)) || ((p[1].x == p[2].x) && (p[1].y == p[2].y) && (p[1].z == p[2].z))) /* only a single line visible - another triangle will plot it */ return; /* we could avoid this if InterpCoor were the same structure as polf requires. */ x1 = p[0].x; y1 = p[0].y; z1 = p[0].z; x2 = p[1].x; y2 = p[1].y; z2 = p[1].z; x3 = p[2].x; y3 = p[2].y; z3 = p[2].z; Npolf(3, tri); } static void plotTriangle( InterpCoor p[3] /* three points of interpolated triangle */ ) { if (((p[0].x == p[1].x) && (p[0].y == p[1].y) && (p[0].z == p[1].z)) && ((p[0].x == p[2].x) && (p[0].y == p[2].y) && (p[0].z == p[2].z))) /* only a single point visible - another triangle will plot it */ return; if (onXedge || onYedge || onZedge) { /* plot the triangle */ Nmove(p[0].x, p[0].y, p[0].z); Ndraw(p[1].x, p[1].y, p[1].z); Ndraw(p[2].x, p[2].y, p[2].z); Ndraw(p[0].x, p[0].y, p[0].z); lastPt = p[0]; /* this line only needed if plotLine uses lastPt */ } else { /* plot three lines, only if needed */ plotLine(p[0], p[1]); plotLine(p[1], p[2]); plotLine(p[2], p[0]); }; } static plotLine(p1, p2) InterpCoor p1, p2; { /* Dont draw this line if it has been drawn before. */ if ((p1.x < p2.x) || ((p1.x == p2.x) && ((p1.y < p2.y) || ((p1.y == p2.y) && (p1.z < p2.z)) ))) { /* Move to the first point if it is not the same as the last point we visited. */ /* This may not be worth it. */ if ((lastPt.x != p1.x) || (lastPt.y != p1.y) || (lastPt.z != p1.z)) Nmove(p1.x, p1.y, p1.z); Ndraw(p2.x, p2.y, p2.z); lastPt = p2; } };