FOX/ObjCryst++  1.10.X (development)
MC.cpp
1 /* This program is free software; you can redistribute it and/or modify
2  it under the terms of the GNU General Public License as published by
3  the Free Software Foundation; either version 2 of the License, or
4  (at your option) any later version.
5 
6  This program is distributed in the hope that it will be useful,
7  but WITHOUT ANY WARRANTY; without even the implied warranty of
8  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
9  GNU General Public License for more details.
10 
11  You should have received a copy of the GNU General Public License
12  along with this program; if not, write to the Free Software
13  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
14 */
16 // FileName: MarchingCubes.cpp
17 // Author : Michael Y. Polyakov
18 // email : myp@andrew.cmu.edu or mikepolyakov@hotmail.com
19 // website : www.angelfire.com/linux/myp
20 // date : July 2002
21 //
22 // Description: 'Straight' and Recursive Marching Cubes Algorithms
23 // Recursive method is faster than the 'straight' one, especially when intersection does not
24 // have to be searched for every time.
25 // Normal vectors are defined for each vertex as a gradients
26 // For function definitions see MarchingCubes.h
27 // For a tutorial on Marching Cubes please visit www.angelfire.com/myp/linux
28 //
29 // Please email me with any suggestion/bugs.
31 
32 #include "MC.h"
33 #include "MCTable.h"
34 #include <math.h>
35 
36 //Linear Interpolation between two points
37 mpVector LinearInterp(mp4Vector p1, mp4Vector p2, float value)
38 {
39  mpVector p;
40  if(fabs(p1.val - p2.val) > 0.00001)
41  p = (mpVector)p1 + ((mpVector)p2 - (mpVector)p1)/(p2.val - p1.val)*(value - p1.val);
42  else
43  p = (mpVector)p1;
44  return p;
45 }
46 
47 
48 //Macros used to compute gradient vector on each vertex of a cube
49 //argument should be the name of array of vertices
50 //can be verts or *verts if done by reference
51 #define CALC_GRAD_VERT_0(verts) mp4Vector(points[ind-YtimeZ].val-(verts[1]).val,points[ind-pointsZ].val-(verts[4]).val,points[ind-1].val-(verts[3]).val, (verts[0]).val);
52 #define CALC_GRAD_VERT_1(verts) mp4Vector((verts[0]).val-points[ind+2*YtimeZ].val,points[ind+YtimeZ-pointsZ].val-(verts[5]).val,points[ind+YtimeZ-1].val-(verts[2]).val, (verts[1]).val);
53 #define CALC_GRAD_VERT_2(verts) mp4Vector((verts[3]).val-points[ind+2*YtimeZ+1].val,points[ind+YtimeZ-ncellsZ].val-(verts[6]).val,(verts[1]).val-points[ind+YtimeZ+2].val, (verts[2]).val);
54 #define CALC_GRAD_VERT_3(verts) mp4Vector(points[ind-YtimeZ+1].val-(verts[2]).val,points[ind-ncellsZ].val-(verts[7]).val,(verts[0]).val-points[ind+2].val, (verts[3]).val);
55 #define CALC_GRAD_VERT_4(verts) mp4Vector(points[ind-YtimeZ+ncellsZ+1].val-(verts[5]).val,(verts[0]).val-points[ind+2*pointsZ].val,points[ind+ncellsZ].val-(verts[7]).val, (verts[4]).val);
56 #define CALC_GRAD_VERT_5(verts) mp4Vector((verts[4]).val-points[ind+2*YtimeZ+ncellsZ+1].val,(verts[1]).val-points[ind+YtimeZ+2*pointsZ].val,points[ind+YtimeZ+ncellsZ].val-(verts[6]).val, (verts[5]).val);
57 #define CALC_GRAD_VERT_6(verts) mp4Vector((verts[7]).val-points[ind+2*YtimeZ+ncellsZ+2].val,(verts[2]).val-points[ind+YtimeZ+2*ncellsZ+3].val,(verts[5]).val-points[ind+YtimeZ+ncellsZ+3].val, (verts[6]).val);
58 #define CALC_GRAD_VERT_7(verts) mp4Vector(points[ind-YtimeZ+ncellsZ+2].val-(verts[6]).val,(verts[3]).val-points[ind+2*ncellsZ+3].val,(verts[4]).val-points[ind+ncellsZ+3].val, (verts[7]).val);
59 
61 // GLOBAL //
62 //Global variables - so they dont have to be passed into functions
63 int pointsZ; //number of points on Z zxis (equal to ncellsZ+1)
64 int YtimeZ; //'plane' of cubes on YZ (equal to (ncellsY+1)*pointsZ )
66 
67 
69 // 'STRAIGHT' MARCHING CUBES ALGORITHM ///////////////////////////////////////////////////////////
71 //for gradients at the edges values 1.0, 1.0, 1.0, 1.0 are given
72 TRIANGLE* MC(int ncellsX, int ncellsY, int ncellsZ,
73  float gradFactorX, float gradFactorY, float gradFactorZ,
74  float minValue, mp4Vector * points, int &numTriangles)
75 {
76  //this should be enough space, if not change 3 to 4
77  TRIANGLE * triangles = new TRIANGLE[3*ncellsX*ncellsY*ncellsZ];
78  numTriangles = int(0);
79 
80  pointsZ = ncellsZ+1; //initialize global variable (for extra speed)
81  YtimeZ = (ncellsY+1)*pointsZ;
82  int lastX = ncellsX; //left from older version
83  int lastY = ncellsY;
84  int lastZ = ncellsZ;
85 
86  mp4Vector *verts[8]; //vertices of a cube (array of pointers for extra speed)
87  mpVector intVerts[12]; //linearly interpolated vertices on each edge
88  int cubeIndex; //shows which vertices are outside/inside
89  int edgeIndex; //index returned by edgeTable[cubeIndex]
90  mp4Vector gradVerts[8]; //gradients at each vertex of a cube
91  mpVector grads[12]; //linearly interpolated gradients on each edge
92  int indGrad; //shows which gradients already have been computed
93  int ind, ni, nj; //ind: index of vertex 0
94  //factor by which corresponding coordinates of gradient vectors are scaled
95  mpVector factor(1.0/(2.0*gradFactorX), 1.0/(2.0*gradFactorY), 1.0/(2.0*gradFactorZ));
96 
97  //MAIN LOOP: goes through all the points
98  for(int i=0; i < lastX; i++) { //x axis
99  ni = i*YtimeZ;
100  for(int j=0; j < lastY; j++) { //y axis
101  nj = j*pointsZ;
102  for(int k=0; k < lastZ; k++, ind++) //z axis
103  {
104  //initialize vertices
105  ind = ni + nj + k;
106  verts[0] = &points[ind];
107  verts[1] = &points[ind + YtimeZ];
108  verts[4] = &points[ind + pointsZ];
109  verts[5] = &points[ind + YtimeZ + pointsZ];
110  verts[2] = &points[ind + YtimeZ + 1];
111  verts[3] = &points[ind + 1];
112  verts[6] = &points[ind + YtimeZ + pointsZ + 1];
113  verts[7] = &points[ind + pointsZ + 1];
114 
115  //get the index
116  cubeIndex = int(0);
117  for(int n=0; n < 8; n++)
118  if(verts[n]->val <= minValue) cubeIndex |= (1 << n);
119 
120  //check if its completely inside or outside
121  if(!edgeTable[cubeIndex]) continue;
122 
123  indGrad = int(0);
124  edgeIndex = edgeTable[cubeIndex];
125 
126  if(edgeIndex & 1) {
127  intVerts[0] = LinearInterp(*verts[0], *verts[1], minValue);
128  if(i != 0 && j != 0 && k != 0) gradVerts[0] = CALC_GRAD_VERT_0(*verts)
129  else gradVerts[0] = mp4Vector(1.0, 1.0, 1.0, 1.0); //for now do not wrap around
130  if(i != lastX-1 && j != 0 && k != 0) gradVerts[1] = CALC_GRAD_VERT_1(*verts)
131  else gradVerts[1] = mp4Vector(1.0, 1.0, 1.0, 1.0);
132  indGrad |= 3;
133  grads[0] = LinearInterp(gradVerts[0], gradVerts[1], minValue);
134  grads[0].x *= factor.x; grads[0].y *= factor.y; grads[0].z *= factor.z;
135  }
136  if(edgeIndex & 2) {
137  intVerts[1] = LinearInterp(*verts[1], *verts[2], minValue);
138  if(! (indGrad & 2)) {
139  if(i != lastX-1 && j != 0 && k != 0) gradVerts[1] = CALC_GRAD_VERT_1(*verts)
140  else gradVerts[1] = mp4Vector(1.0, 1.0, 1.0, 1.0);
141  indGrad |= 2;
142  }
143  if(i != lastX-1 && j != 0 && k != 0) gradVerts[2] = CALC_GRAD_VERT_2(*verts)
144  else gradVerts[2] = mp4Vector(1.0, 1.0, 1.0, 1.0);
145  indGrad |= 4;
146  grads[1] = LinearInterp(gradVerts[1], gradVerts[2], minValue);
147  grads[1].x *= factor.x; grads[1].y *= factor.y; grads[1].z *= factor.z;
148  }
149  if(edgeIndex & 4) {
150  intVerts[2] = LinearInterp(*verts[2], *verts[3], minValue);
151  if(! (indGrad & 4)) {
152  if(i != lastX-1 && j != 0 && k != 0) gradVerts[2] = CALC_GRAD_VERT_2(*verts)
153  else gradVerts[2] = mp4Vector(1.0, 1.0, 1.0, 1.0);
154  indGrad |= 4;
155  }
156  if(i != 0 && j != 0 && k != lastZ-1) gradVerts[3] = CALC_GRAD_VERT_3(*verts)
157  else gradVerts[3] = mp4Vector(1.0, 1.0, 1.0, 1.0);
158  indGrad |= 8;
159  grads[2] = LinearInterp(gradVerts[2], gradVerts[3], minValue);
160  grads[2].x *= factor.x; grads[2].y *= factor.y; grads[2].z *= factor.z;
161  }
162  if(edgeIndex & 8) {
163  intVerts[3] = LinearInterp(*verts[3], *verts[0], minValue);
164  if(! (indGrad & 8)) {
165  if(i != 0 && j != 0 && k != lastZ-1) gradVerts[3] = CALC_GRAD_VERT_3(*verts)
166  else gradVerts[3] = mp4Vector(1.0, 1.0, 1.0, 1.0);
167  indGrad |= 8;
168  }
169  if(! (indGrad & 1)) {
170  if(i != 0 && j != 0 && k != 0) gradVerts[0] = CALC_GRAD_VERT_0(*verts)
171  else gradVerts[0] = mp4Vector(1.0, 1.0, 1.0, 1.0);
172  indGrad |= 1;
173  }
174  grads[3] = LinearInterp(gradVerts[3], gradVerts[0], minValue);
175  grads[3].x *= factor.x; grads[3].y *= factor.y; grads[3].z *= factor.z;
176  }
177  if(edgeIndex & 16) {
178  intVerts[4] = LinearInterp(*verts[4], *verts[5], minValue);
179 
180  if(i != 0 && j != lastY-1 && k != 0) gradVerts[4] = CALC_GRAD_VERT_4(*verts)
181  else gradVerts[4] = mp4Vector(1.0, 1.0, 1.0, 1.0);
182 
183  if(i != lastX-1 && j != lastY-1 && k != 0) gradVerts[5] = CALC_GRAD_VERT_5(*verts)
184  else gradVerts[5] = mp4Vector(1.0, 1.0, 1.0, 1.0);
185 
186  indGrad |= 48;
187  grads[4] = LinearInterp(gradVerts[4], gradVerts[5], minValue);
188  grads[4].x *= factor.x; grads[4].y *= factor.y; grads[4].z *= factor.z;
189  }
190  if(edgeIndex & 32) {
191  intVerts[5] = LinearInterp(*verts[5], *verts[6], minValue);
192  if(! (indGrad & 32)) {
193  if(i != lastX-1 && j != lastY-1 && k != 0) gradVerts[5] = CALC_GRAD_VERT_5(*verts)
194  else gradVerts[5] = mp4Vector(1.0, 1.0, 1.0, 1.0);
195  indGrad |= 32;
196  }
197 
198  if(i != lastX-1 && j != lastY-1 && k != lastZ-1) gradVerts[6] = CALC_GRAD_VERT_6(*verts)
199  else gradVerts[6] = mp4Vector(1.0, 1.0, 1.0, 1.0);
200  indGrad |= 64;
201  grads[5] = LinearInterp(gradVerts[5], gradVerts[6], minValue);
202  grads[5].x *= factor.x; grads[5].y *= factor.y; grads[5].z *= factor.z;
203  }
204  if(edgeIndex & 64) {
205  intVerts[6] = LinearInterp(*verts[6], *verts[7], minValue);
206  if(! (indGrad & 64)) {
207  if(i != lastX-1 && j != lastY-1 && k != lastZ-1) gradVerts[6] = CALC_GRAD_VERT_6(*verts)
208  else gradVerts[6] = mp4Vector(1.0, 1.0, 1.0, 1.0);
209  indGrad |= 64;
210  }
211 
212  if(i != 0 && j != lastY-1 && k != lastZ-1) gradVerts[7] = CALC_GRAD_VERT_7(*verts)
213  else gradVerts[7] = mp4Vector(1.0, 1.0, 1.0, 1.0);
214  indGrad |= 128;
215  grads[6] = LinearInterp(gradVerts[6], gradVerts[7], minValue);
216  grads[6].x *= factor.x; grads[6].y *= factor.y; grads[6].z *= factor.z;
217  }
218  if(edgeIndex & 128) {
219  intVerts[7] = LinearInterp(*verts[7], *verts[4], minValue);
220  if(! (indGrad & 128)) {
221  if(i != 0 && j != lastY-1 && k != lastZ-1) gradVerts[7] = CALC_GRAD_VERT_7(*verts)
222  else gradVerts[7] = mp4Vector(1.0, 1.0, 1.0, 1.0);
223  indGrad |= 128;
224  }
225  if(! (indGrad & 16)) {
226  if(i != 0 && j != lastY-1 && k != 0) gradVerts[4] = CALC_GRAD_VERT_4(*verts)
227  else gradVerts[4] = mp4Vector(1.0, 1.0, 1.0, 1.0);
228  indGrad |= 16;
229  }
230  grads[7] = LinearInterp(gradVerts[7], gradVerts[4], minValue);
231  grads[7].x *= factor.x; grads[7].y *= factor.y; grads[7].z *= factor.z;
232  }
233  if(edgeIndex & 256) {
234  intVerts[8] = LinearInterp(*verts[0], *verts[4], minValue);
235  if(! (indGrad & 1)) {
236  if(i != 0 && j != 0 && k != 0) gradVerts[0] = CALC_GRAD_VERT_0(*verts)
237  else gradVerts[0] = mp4Vector(1.0, 1.0, 1.0, 1.0);
238  indGrad |= 1;
239  }
240  if(! (indGrad & 16)) {
241  if(i != 0 && j != lastY-1 && k != 0) gradVerts[4] = CALC_GRAD_VERT_4(*verts)
242  else gradVerts[4] = mp4Vector(1.0, 1.0, 1.0, 1.0);
243  indGrad |= 16;
244  }
245  grads[8] = LinearInterp(gradVerts[0], gradVerts[4], minValue);
246  grads[8].x *= factor.x; grads[8].y *= factor.y; grads[8].z *= factor.z;
247  }
248  if(edgeIndex & 512) {
249  intVerts[9] = LinearInterp(*verts[1], *verts[5], minValue);
250  if(! (indGrad & 2)) {
251  if(i != lastX-1 && j != 0 && k != 0) gradVerts[1] = CALC_GRAD_VERT_1(*verts)
252  else gradVerts[1] = mp4Vector(1.0, 1.0, 1.0, 1.0);
253  indGrad |= 2;
254  }
255  if(! (indGrad & 32)) {
256  if(i != lastX-1 && j != lastY-1 && k != 0) gradVerts[5] = CALC_GRAD_VERT_5(*verts)
257  else gradVerts[5] = mp4Vector(1.0, 1.0, 1.0, 1.0);
258  indGrad |= 32;
259  }
260  grads[9] = LinearInterp(gradVerts[1], gradVerts[5], minValue);
261  grads[9].x *= factor.x; grads[9].y *= factor.y; grads[9].z *= factor.z;
262  }
263  if(edgeIndex & 1024) {
264  intVerts[10] = LinearInterp(*verts[2], *verts[6], minValue);
265  if(! (indGrad & 4)) {
266  if(i != lastX-1 && j != 0 && k != 0) gradVerts[2] = CALC_GRAD_VERT_2(*verts)
267  else gradVerts[5] = mp4Vector(1.0, 1.0, 1.0, 1.0);
268  indGrad |= 4;
269  }
270  if(! (indGrad & 64)) {
271  if(i != lastX-1 && j != lastY-1 && k != lastZ-1) gradVerts[6] = CALC_GRAD_VERT_6(*verts)
272  else gradVerts[6] = mp4Vector(1.0, 1.0, 1.0, 1.0);
273  indGrad |= 64;
274  }
275  grads[10] = LinearInterp(gradVerts[2], gradVerts[6], minValue);
276  grads[10].x *= factor.x; grads[10].y *= factor.y; grads[10].z *= factor.z;
277  }
278  if(edgeIndex & 2048) {
279  intVerts[11] = LinearInterp(*verts[3], *verts[7], minValue);
280  if(! (indGrad & 8)) {
281  if(i != 0 && j != 0 && k != lastZ-1) gradVerts[3] = CALC_GRAD_VERT_3(*verts)
282  else gradVerts[3] = mp4Vector(1.0, 1.0, 1.0, 1.0);
283  indGrad |= 8;
284  }
285  if(! (indGrad & 128)) {
286  if(i != 0 && j != lastY-1 && k != lastZ-1) gradVerts[7] = CALC_GRAD_VERT_7(*verts)
287  else gradVerts[7] = mp4Vector(1.0, 1.0, 1.0, 1.0);
288  indGrad |= 128;
289  }
290  grads[11] = LinearInterp(gradVerts[3], gradVerts[7], minValue);
291  grads[11].x *= factor.x; grads[11].y *= factor.y; grads[11].z *= factor.z;
292  }
293 
294  //now build the triangles using triTable
295  for (int n=0; triTable[cubeIndex][n] != -1; n+=3) {
296  int index[3] = {triTable[cubeIndex][n+2], triTable[cubeIndex][n+1], triTable[cubeIndex][n]};
297  for(int h=0; h < 3; h++) { //copy vertices and normals into triangles array
298  triangles[numTriangles].p[h] = intVerts[index[h]];
299  triangles[numTriangles].norm[h] = grads[index[h]];
300  }
301  numTriangles++; //one more triangle has been added
302  }
303  } //END OF FOR LOOP ON Z AXIS
304  } //END OF FOR LOOP ON Y AXIS
305  } //END OF FOR LOOP ON X AXIS
306 
307  //free all wasted space
308  TRIANGLE * retTriangles = new TRIANGLE[numTriangles];
309  for(int i=0; i < numTriangles; i++)
310  retTriangles[i] = triangles[i];
311  delete [] triangles;
312 
313  return retTriangles;
314 }
317 
318 
319 
321 // RECURSIVE MARCHING CUBES ALGORITHM ////////////////////////////////////////////////////////////////////////
323 
324 // MACROS ///////////////////////////////////////////////////////////////////////////////////////////////////
325 //these macros initialize data and then run marching cubes on the cube with the surface having the specified
326 // number as 'recieving' data (but number of that surface for the current cube is going to be 'opposite').
327 // Each runs the corresponding recursive function
328 // For numbering, to see which indices of prevVerts,... correspong to indices of the current cube, see
329 // my webpage at www.angelfire.com/linux/myp
330 
331 #define MC_FACE0 \
332 { \
333  if(! marchedCubes[ind - 1]) { \
334  passGradIndex = 0; \
335  if(gradIndex & 1) passGradIndex |= 8; \
336  if(gradIndex & 2) passGradIndex |= 4; \
337  if(gradIndex & 32) passGradIndex |= 64; \
338  if(gradIndex & 16) passGradIndex |= 128; \
339  passEdgeIndex = 0; \
340  if(edgeIndex & 1) passEdgeIndex |= 4; \
341  if(edgeIndex & 512) passGradIndex |= 1024; \
342  if(edgeIndex & 16) passEdgeIndex |= 64; \
343  if(edgeIndex & 256) passGradIndex |= 2048; \
344  prevVerts[0] = verts[0]; prevVerts[1] = verts[1]; prevVerts[2] = verts[5]; prevVerts[3] = verts[4]; \
345  prevIntVerts[0] = intVerts[0]; prevIntVerts[1] = intVerts[9]; \
346  prevIntVerts[2] = intVerts[4]; prevIntVerts[3] = intVerts[8]; \
347  prevGradVerts[0] = gradVerts[0]; prevGradVerts[1] = gradVerts[1]; \
348  prevGradVerts[2] = gradVerts[5]; prevGradVerts[3] = gradVerts[4]; \
349  prevGrads[0] = grads[0]; prevGrads[1] = grads[9]; prevGrads[2] = grads[4]; prevGrads[3] = grads[8]; \
350  triangles = MCFace0(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, \
351  ind-1, i, j, k-1, minValue, points, triangles, numTriangles, \
352  prevVerts, prevIntVerts, passEdgeIndex, \
353  prevGradVerts, prevGrads, passGradIndex, marchedCubes); \
354  } \
355 }
356 
357 #define MC_FACE1 \
358 { \
359  if(! marchedCubes[ind + YtimeZ]) { \
360  passGradIndex = 0; \
361  if(gradIndex & 4) passGradIndex |= 8; \
362  if(gradIndex & 2) passGradIndex |= 1; \
363  if(gradIndex & 32) passGradIndex |= 16; \
364  if(gradIndex & 64) passGradIndex |= 128; \
365  passEdgeIndex = 0; \
366  if(edgeIndex & 2) passEdgeIndex |= 8; \
367  if(edgeIndex & 512) passEdgeIndex |= 256; \
368  if(edgeIndex & 32) passEdgeIndex |= 128; \
369  if(edgeIndex & 1024) passEdgeIndex |= 2048; \
370  prevVerts[0] = verts[2]; prevVerts[1] = verts[1]; prevVerts[2] = verts[5]; prevVerts[3] = verts[6]; \
371  prevIntVerts[0] = intVerts[1]; prevIntVerts[1] = intVerts[9]; \
372  prevIntVerts[2] = intVerts[5]; prevIntVerts[3] = intVerts[10]; \
373  prevGradVerts[0] = gradVerts[2]; prevGradVerts[1] = gradVerts[1]; \
374  prevGradVerts[2] = gradVerts[5]; prevGradVerts[3] = gradVerts[6]; \
375  prevGrads[0] = grads[1]; prevGrads[1] = grads[9]; prevGrads[2] = grads[5]; prevGrads[3] = grads[10]; \
376  triangles = MCFace1(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, \
377  ind+YtimeZ, i+1, j, k, minValue, points, triangles, numTriangles, \
378  prevVerts, prevIntVerts, passEdgeIndex, \
379  prevGradVerts, prevGrads, passGradIndex, marchedCubes); \
380  } \
381 }
382 
383 #define MC_FACE2 \
384 { \
385  if(! marchedCubes[ind + 1]) { \
386  passGradIndex = 0; \
387  if(gradIndex & 8) passGradIndex |= 1; \
388  if(gradIndex & 4) passGradIndex |= 2; \
389  if(gradIndex & 64) passGradIndex |= 32; \
390  if(gradIndex & 128) passGradIndex |= 16; \
391  passEdgeIndex = 0; \
392  if(edgeIndex & 4) passEdgeIndex |= 1; \
393  if(edgeIndex & 1024) passEdgeIndex |= 512; \
394  if(edgeIndex & 64) passEdgeIndex |= 16; \
395  if(edgeIndex & 2048) passEdgeIndex |= 256; \
396  prevVerts[0] = verts[3]; prevVerts[1] = verts[2]; prevVerts[2] = verts[6]; prevVerts[3] = verts[7]; \
397  prevIntVerts[0] = intVerts[2]; prevIntVerts[1] = intVerts[10]; \
398  prevIntVerts[2] = intVerts[6]; prevIntVerts[3] = intVerts[11]; \
399  prevGradVerts[0] = gradVerts[3]; prevGradVerts[1] = gradVerts[2]; \
400  prevGradVerts[2] = gradVerts[6]; prevGradVerts[3] = gradVerts[7]; \
401  prevGrads[0] = grads[2]; prevGrads[1] = grads[10]; prevGrads[2] = grads[6]; prevGrads[3] = grads[11]; \
402  triangles = MCFace2(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, \
403  ind+1, i, j, k+1, minValue, points, triangles, numTriangles, \
404  prevVerts, prevIntVerts, passEdgeIndex, \
405  prevGradVerts, prevGrads, passGradIndex, marchedCubes); \
406  } \
407 }
408 
409 #define MC_FACE3 \
410 { \
411  if(! marchedCubes[ind - YtimeZ]) { \
412  passGradIndex = 0; \
413  if(gradIndex & 8) passGradIndex |= 4; \
414  if(gradIndex & 1) passGradIndex |= 2; \
415  if(gradIndex & 128) passGradIndex |= 64; \
416  if(gradIndex & 16) passGradIndex |= 32; \
417  passEdgeIndex = 0; \
418  if(edgeIndex & 8) passEdgeIndex |= 2; \
419  if(edgeIndex & 256) passEdgeIndex |= 512; \
420  if(edgeIndex & 128) passEdgeIndex |= 32; \
421  if(edgeIndex & 2048) passEdgeIndex |= 1024; \
422  prevVerts[0] = verts[3]; prevVerts[1] = verts[0]; prevVerts[2] = verts[4]; prevVerts[3] = verts[7]; \
423  prevIntVerts[0] = intVerts[3]; prevIntVerts[1] = intVerts[8]; \
424  prevIntVerts[2] = intVerts[7]; prevIntVerts[3] = intVerts[11]; \
425  prevGradVerts[0] = gradVerts[3]; prevGradVerts[1] = gradVerts[0]; \
426  prevGradVerts[2] = gradVerts[4]; prevGradVerts[3] = gradVerts[7]; \
427  prevGrads[0] = grads[3]; prevGrads[1] = grads[8]; prevGrads[2] = grads[7]; prevGrads[3] = grads[11]; \
428  triangles = MCFace3(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, \
429  ind-YtimeZ, i-1, j, k, minValue, points, triangles, numTriangles, \
430  prevVerts, prevIntVerts, passEdgeIndex, \
431  prevGradVerts, prevGrads, passGradIndex, marchedCubes); \
432  } \
433 }
434 
435 //done
436 #define MC_FACE4 \
437 { \
438  if(! marchedCubes[ind + pointsZ]) { \
439  passGradIndex = 0; \
440  if(gradIndex & 128) passGradIndex |= 8; \
441  if(gradIndex & 64) passGradIndex |= 4; \
442  if(gradIndex & 32) passGradIndex |= 2; \
443  if(gradIndex & 16) passGradIndex |= 1; \
444  passEdgeIndex = 0; \
445  if(edgeIndex & 128) passEdgeIndex |= 8; \
446  if(edgeIndex & 64) passEdgeIndex |= 4; \
447  if(edgeIndex & 32) passEdgeIndex |= 2; \
448  if(edgeIndex & 16) passEdgeIndex |= 1; \
449  prevVerts[0] = verts[7]; prevVerts[1] = verts[6]; prevVerts[2] = verts[5]; prevVerts[3] = verts[4]; \
450  prevIntVerts[0] = intVerts[6]; prevIntVerts[1] = intVerts[5]; \
451  prevIntVerts[2] = intVerts[4]; prevIntVerts[3] = intVerts[7]; \
452  prevGradVerts[0] = gradVerts[7]; prevGradVerts[1] = gradVerts[6]; \
453  prevGradVerts[2] = gradVerts[5]; prevGradVerts[3] = gradVerts[4]; \
454  prevGrads[0] = grads[6]; prevGrads[1] = grads[5]; prevGrads[2] = grads[4]; prevGrads[3] = grads[7]; \
455  triangles = MCFace4(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, \
456  ind+pointsZ, i, j+1, k, minValue, points, triangles, numTriangles, \
457  prevVerts, prevIntVerts, passEdgeIndex, \
458  prevGradVerts, prevGrads, passGradIndex, marchedCubes); \
459  } \
460 }
461 
462 #define MC_FACE5 \
463 { \
464  if(! marchedCubes[ind - ncellsZ - 1]) { \
465  passGradIndex = (gradIndex << 4) & 240; \
466  passEdgeIndex = (edgeIndex << 4) & 240; \
467  prevVerts[0] = verts[3]; prevVerts[1] = verts[2]; prevVerts[2] = verts[1]; prevVerts[3] = verts[0]; \
468  prevIntVerts[0] = intVerts[2]; prevIntVerts[1] = intVerts[1]; \
469  prevIntVerts[2] = intVerts[0]; prevIntVerts[3] = intVerts[3]; \
470  prevGradVerts[0] = gradVerts[3]; prevGradVerts[1] = gradVerts[2]; \
471  prevGradVerts[2] = gradVerts[1]; prevGradVerts[3] = gradVerts[0]; \
472  prevGrads[0] = grads[2]; prevGrads[1] = grads[1]; prevGrads[2] = grads[0]; prevGrads[3] = grads[3]; \
473  triangles = MCFace5(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, \
474  ind-ncellsZ-1, i, j-1, k, minValue, points, triangles, numTriangles, \
475  prevVerts, prevIntVerts, passEdgeIndex, \
476  prevGradVerts, prevGrads, passGradIndex, marchedCubes); \
477  } \
478 }
479 
481 
482 
483 
484 
485 // RECURSIVE Marching Cubes Function - cubes at indexes ii, jj, kk intersect the surface
486 // Number of intersecting cubes = numCubes
487 TRIANGLE* MarchingCubesRec(int ncellsX, int ncellsY, int ncellsZ,
488  float gradFactorX, float gradFactorY, float gradFactorZ,
489  int numCubes, int *ii, int *jj, int *kk,
490  float minValue, mp4Vector * points, int &numTriangles)
491 {
492  TRIANGLE * triangles = new TRIANGLE[3*ncellsX*ncellsY*ncellsZ];
493  numTriangles = int(0);
494  //check if none of the starting points are on the outside perimeter
495  for(int n=0; n < numCubes; n++) {
496  if(ii[n] <= 0 || ii[n] >= ncellsX-1) return NULL;
497  if(jj[n] <= 0 || jj[n] >= ncellsY-1) return NULL;
498  if(kk[n] <= 0 || kk[n] >= ncellsZ-1) return NULL;
499  }
500  //array stores which cubes have been marched through - initialized to FALSE
501  int all = ncellsX*ncellsY*ncellsZ;
502  bool* marchedCubes = new bool[all];
503  for(int i=0; i < all; i++) marchedCubes[i] = FALSE; //initialize
504 
505  mp4Vector verts[8]; //vertices of a starting cube
506  mpVector intVerts[12]; //linearly interpolated vertices on each edge
507  int edgeIndex; //shows which edges are intersected
508  mp4Vector gradVerts[8]; //gradients at each vertex of a cube
509  mpVector grads[12]; //linearly interpolated gradients on each edge
510  int gradIndex; //show on which vertices gradients have been computed
511  //initialize global variables - for speed - these would be used by all the recursive functions
512  pointsZ = ncellsZ+1;
513  YtimeZ = (ncellsY+1)*pointsZ;
514  //arrays used to pass already computed values from this initial cube to the next ones
515  mp4Vector prevVerts[4]; //passes vertices
516  mpVector prevIntVerts[4]; //passes interpolated vertices on edges
517  mp4Vector prevGradVerts[4]; //passes gradients at vertices
518  mpVector prevGrads[4]; //passes interpolated gradients on edges
519 
520  //two new indexes formed for each face
521  int passEdgeIndex, passGradIndex; //used to tell which vertices and which edges have been initialized
522  //indices
523  int i, j, k;
524  //initialize first cubes and 'launch' the recursion for each
525  for(int n=0; n < numCubes; n++) {
526  //init vertices
527  i = ii[n]; j = jj[n]; k = kk[n];
528  int ind = i*YtimeZ + j*pointsZ + k;
529  verts[0] = points[ind];
530  verts[1] = points[ind + YtimeZ];
531  verts[2] = points[ind + YtimeZ + 1];
532  verts[3] = points[ind + 1];
533  verts[4] = points[ind + pointsZ];
534  verts[5] = points[ind + YtimeZ + pointsZ];
535  verts[6] = points[ind + YtimeZ + pointsZ + 1];
536  verts[7] = points[ind + pointsZ + 1];
537 
538  //first check if this cube wasnt marched in recursive calls of the previous cube
539  if(! marchedCubes[ind]) {
540  //run marching cubes on the initial cube
541  gradIndex = edgeIndex = 0;
542  triangles = MarchOneCube(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ,
543  ind, i, j, k, minValue, points, triangles, numTriangles,
544  verts, intVerts, edgeIndex, gradVerts, grads, gradIndex);
545  //this cube has been done:
546  marchedCubes[ind] = TRUE;
547  //run M.C. on all 6 faces
548  MC_FACE0
549  MC_FACE1
550  MC_FACE2
551  MC_FACE3
552  MC_FACE4
553  MC_FACE5
554  }
555  }
556  //free wasted space
557  TRIANGLE * retTriangles = new TRIANGLE[numTriangles];
558  for(int i=0; i < numTriangles; i++) retTriangles[i] = triangles[i];
559  delete [] triangles;
560  delete [] marchedCubes;
561  return retTriangles; //done
562 }
563 
564 //SURFACE 0 - Cube ran on surface 0 of previous cube. Recieving side: 2.
565 TRIANGLE* MCFace0(int ncellsX, int ncellsY, int ncellsZ,
566  float gradFactorX, float gradFactorY, float gradFactorZ,
567  int ind, int i, int j, int k,
568  float minValue, mp4Vector * points, TRIANGLE *triangles, int &numTriangles,
569  mp4Vector prevVerts[4], mpVector prevIntVerts[4], int edgeIndex,
570  mp4Vector prevGradVerts[4], mpVector prevGrads[4], int gradIndex, bool* marchedCubes)
571 {
572  //first check if not outside the region
573  if(i <= 0 || i >= ncellsX-1 || j <= 0 || j >= ncellsY-1 || k <= 0 || k >= ncellsZ-1)
574  return triangles;
575  //make sure we save that this cube was marched through
576  marchedCubes[ind] = TRUE;
577  //initialize vertices
578  mp4Vector verts[8];
579  verts[0] = points[ind];
580  verts[1] = points[ind + YtimeZ];
581  verts[2] = prevVerts[1];
582  verts[3] = prevVerts[0];
583  verts[4] = points[ind + ncellsZ + 1];
584  verts[5] = points[ind + YtimeZ + ncellsZ + 1];
585  verts[6] = prevVerts[2];
586  verts[7] = prevVerts[3];
587  //initialize edges from the last cube
588  mpVector intVerts[12];
589  intVerts[2] = prevIntVerts[0]; intVerts[10] = prevIntVerts[1];
590  intVerts[6] = prevIntVerts[2]; intVerts[11] = prevIntVerts[3];
591  //initialize gradients on vertices
592  mp4Vector gradVerts[8];
593  gradVerts[3] = prevGradVerts[0]; gradVerts[2] = prevGradVerts[1];
594  gradVerts[6] = prevGradVerts[2]; gradVerts[7] = prevGradVerts[3];
595  //initialize gradients on edges from the last cube
596  mpVector grads[12];
597  grads[2] = prevGrads[0]; grads[10] = prevGrads[1]; grads[6] = prevGrads[2]; grads[11] = prevGrads[3];
598  //for test if this cube is intersected:
599  int oldNumTriangles = numTriangles;
600  //run marching cubes on this cube
601  triangles = MarchOneCube(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, ind, i, j, k,
602  minValue, points, triangles, numTriangles, verts, intVerts, edgeIndex,
603  gradVerts, grads, gradIndex);
604  //check if this cube is intersected
605  if(numTriangles == oldNumTriangles)
606  return triangles;
607  //two new indexes formed for each face to be passed into other recursive functions
608  int passEdgeIndex, passGradIndex;
609  //run recursive functions on each surface
610  MC_FACE0
611  MC_FACE1
612  MC_FACE3
613  MC_FACE4
614  MC_FACE5
615  return triangles;
616 }
617 
618 //SURFACE 1 - Cube ran on surface 1 of previous cube. Recieving side: 3.
619 TRIANGLE* MCFace1(int ncellsX, int ncellsY, int ncellsZ,
620  float gradFactorX, float gradFactorY, float gradFactorZ,
621  int ind, int i, int j, int k,
622  float minValue, mp4Vector * points, TRIANGLE *triangles, int &numTriangles,
623  mp4Vector prevVerts[4], mpVector prevIntVerts[4], int edgeIndex,
624  mp4Vector prevGradVerts[4], mpVector prevGrads[4], int gradIndex, bool* marchedCubes)
625 {
626  //first check if not outside the region
627  if(i <= 0 || i >= ncellsX-1 || j <= 0 || j >= ncellsY-1 || k <= 0 || k >= ncellsZ-1)
628  return triangles;
629  //make sure we save that this cube was marched through
630  marchedCubes[ind] = TRUE;
631  //initialize vertices
632  mp4Vector verts[8];
633  verts[0] = prevVerts[1];
634  verts[1] = points[ind + YtimeZ];
635  verts[2] = points[ind + YtimeZ + 1];
636  verts[3] = prevVerts[0];
637  verts[4] = prevVerts[2];
638  verts[5] = points[ind + YtimeZ + ncellsZ + 1];
639  verts[6] = points[ind + YtimeZ + ncellsZ + 2];
640  verts[7] = prevVerts[3];
641  //initialize edges from the last cube
642  mpVector intVerts[12];
643  intVerts[3] = prevIntVerts[0]; intVerts[8] = prevIntVerts[1];
644  intVerts[7] = prevIntVerts[2]; intVerts[11] = prevIntVerts[3];
645  //initialize gradients on vertices
646  mp4Vector gradVerts[8];
647  gradVerts[3] = prevGradVerts[0]; gradVerts[0] = prevGradVerts[1];
648  gradVerts[4] = prevGradVerts[2]; gradVerts[7] = prevGradVerts[3];
649  //initialize gradients on edges from the last cube
650  mpVector grads[12];
651  grads[3] = prevGrads[0]; grads[8] = prevGrads[1]; grads[7] = prevGrads[2]; grads[11] = prevGrads[3];
652  //for test if this cube is intersected:
653  int oldNumTriangles = numTriangles;
654  //run marching cubes on this cube
655  triangles = MarchOneCube(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, ind, i, j, k,
656  minValue, points, triangles, numTriangles, verts, intVerts, edgeIndex,
657  gradVerts, grads, gradIndex);
658  //check if this cube is intersected
659  if(numTriangles == oldNumTriangles)
660  return triangles;
661  //two new indexes formed for each face to be passed into other recursive functions
662  int passEdgeIndex, passGradIndex;
663  //run recursive functions on each surface
664  MC_FACE0
665  MC_FACE1
666  MC_FACE2
667  MC_FACE4
668  MC_FACE5
669  return triangles;
670 }
671 
672 //SURFACE 2 - Cube ran on surface 2 of previous cube. Recieving side: 0.
673 TRIANGLE* MCFace2(int ncellsX, int ncellsY, int ncellsZ,
674  float gradFactorX, float gradFactorY, float gradFactorZ,
675  int ind, int i, int j, int k,
676  float minValue, mp4Vector * points, TRIANGLE *triangles, int &numTriangles,
677  mp4Vector prevVerts[4], mpVector prevIntVerts[4], int edgeIndex,
678  mp4Vector prevGradVerts[4], mpVector prevGrads[4], int gradIndex, bool* marchedCubes)
679 {
680  //first check if not outside the region
681  if(i <= 0 || i >= ncellsX-1 || j <= 0 || j >= ncellsY-1 || k <= 0 || k >= ncellsZ-1)
682  return triangles;
683  //make sure we save that this cube was marched through
684  marchedCubes[ind] = TRUE;
685  //initialize vertices
686  mp4Vector verts[8];
687  verts[0] = prevVerts[0];
688  verts[1] = prevVerts[1];
689  verts[2] = points[ind + YtimeZ + 1];
690  verts[3] = points[ind + 1];
691  verts[4] = prevVerts[3];
692  verts[5] = prevVerts[2];
693  verts[6] = points[ind + YtimeZ + ncellsZ + 2];
694  verts[7] = points[ind + ncellsZ + 2];
695  //initialize edges from the last cube
696  mpVector intVerts[12];
697  intVerts[0] = prevIntVerts[0]; intVerts[9] = prevIntVerts[1];
698  intVerts[4] = prevIntVerts[2]; intVerts[8] = prevIntVerts[3];
699  //initialize gradients on vertices
700  mp4Vector gradVerts[8];
701  gradVerts[0] = prevGradVerts[0]; gradVerts[1] = prevGradVerts[1];
702  gradVerts[5] = prevGradVerts[2]; gradVerts[4] = prevGradVerts[3];
703  //initialize gradients on edges from the last cube
704  mpVector grads[12];
705  grads[0] = prevGrads[0]; grads[9] = prevGrads[1]; grads[4] = prevGrads[2]; grads[8] = prevGrads[3];
706  //for test if this cube is intersected
707  int oldNumTriangles = numTriangles;
708  //run marching cubes on this cube
709  triangles = MarchOneCube(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, ind, i, j, k,
710  minValue, points, triangles, numTriangles, verts, intVerts, edgeIndex,
711  gradVerts, grads, gradIndex);
712  //check if this cube is intersected
713  if(numTriangles == oldNumTriangles)
714  return triangles;
715  //two new indexes formed for each face
716  int passEdgeIndex, passGradIndex;
717  //run recursive functions on each surface
718  MC_FACE1
719  MC_FACE2
720  MC_FACE3
721  MC_FACE4
722  MC_FACE5
723  return triangles;
724 }
725 
726 //SURFACE 3 - Cube ran on surface 3 of previous cube. Recieving side: 1.
727 TRIANGLE* MCFace3(int ncellsX, int ncellsY, int ncellsZ,
728  float gradFactorX, float gradFactorY, float gradFactorZ,
729  int ind, int i, int j, int k,
730  float minValue, mp4Vector * points, TRIANGLE *triangles, int &numTriangles,
731  mp4Vector prevVerts[4], mpVector prevIntVerts[4], int edgeIndex,
732  mp4Vector prevGradVerts[4], mpVector prevGrads[4], int gradIndex, bool* marchedCubes)
733 {
734  //first check if not outside the region
735  if(i <= 0 || i >= ncellsX-1 || j <= 0 || j >= ncellsY-1 || k <= 0 || k >= ncellsZ-1)
736  return triangles;
737  //make sure we save that this cube was marched through
738  marchedCubes[ind] = TRUE;
739  //initialize vertices
740  mp4Vector verts[8];
741  verts[0] = points[ind];
742  verts[1] = prevVerts[1];
743  verts[2] = prevVerts[0];
744  verts[3] = points[ind + 1];
745  verts[4] = points[ind + ncellsZ + 1];
746  verts[5] = prevVerts[2];
747  verts[6] = prevVerts[3];
748  verts[7] = points[ind + ncellsZ + 2];
749  //initialize edges from the last cube
750  mpVector intVerts[12];
751  intVerts[1] = prevIntVerts[0]; intVerts[9] = prevIntVerts[1];
752  intVerts[5] = prevIntVerts[2]; intVerts[10] = prevIntVerts[3];
753  //initialize gradients on vertices
754  mp4Vector gradVerts[8];
755  gradVerts[2] = prevGradVerts[0]; gradVerts[1] = prevGradVerts[1];
756  gradVerts[5] = prevGradVerts[2]; gradVerts[6] = prevGradVerts[3];
757  //initialize gradients on edges from the last cube
758  mpVector grads[12];
759  grads[1] = prevGrads[0]; grads[9] = prevGrads[1]; grads[5] = prevGrads[2]; grads[10] = prevGrads[3];
760  //for test if this cube is intersected
761  int oldNumTriangles = numTriangles;
762  //run marching cubes on this cube
763  triangles = MarchOneCube(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, ind, i, j, k,
764  minValue, points, triangles, numTriangles, verts, intVerts, edgeIndex,
765  gradVerts, grads, gradIndex);
766  //check if this cube is intersected
767  if(numTriangles == oldNumTriangles)
768  return triangles;
769  //two new indexes formed for each face
770  int passEdgeIndex, passGradIndex;
771  //run recursive functions on each surface
772  MC_FACE0
773  MC_FACE2
774  MC_FACE3
775  MC_FACE4
776  MC_FACE5
777  return triangles;
778 }
779 
780 //SURFACE 4 - Cube ran on surface 4 of previous cube. Recieving side: 5.
781 TRIANGLE* MCFace4(int ncellsX, int ncellsY, int ncellsZ,
782  float gradFactorX, float gradFactorY, float gradFactorZ,
783  int ind, int i, int j, int k,
784  float minValue, mp4Vector * points, TRIANGLE *triangles, int &numTriangles,
785  mp4Vector prevVerts[4], mpVector prevIntVerts[4], int edgeIndex,
786  mp4Vector prevGradVerts[4], mpVector prevGrads[4], int gradIndex, bool* marchedCubes)
787 {
788  //first check if not outside the region
789  if(i <= 0 || i >= ncellsX-1 || j <= 0 || j >= ncellsY-1 || k <= 0 || k >= ncellsZ-1)
790  return triangles;
791  //make sure we save that this cube was marched through
792  marchedCubes[ind] = TRUE;
793  //initialize vertices
794  mp4Vector verts[8];
795  verts[0] = prevVerts[3];
796  verts[1] = prevVerts[2];
797  verts[2] = prevVerts[1];
798  verts[3] = prevVerts[0];
799  verts[4] = points[ind + ncellsZ + 1];
800  verts[5] = points[ind + YtimeZ + ncellsZ + 1];
801  verts[6] = points[ind + YtimeZ + ncellsZ + 2];
802  verts[7] = points[ind + ncellsZ + 2];
803  //initialize edges from the last cube
804  mpVector intVerts[12];
805  intVerts[2] = prevIntVerts[0]; intVerts[1] = prevIntVerts[1];
806  intVerts[0] = prevIntVerts[2]; intVerts[3] = prevIntVerts[3];
807  //initialize gradients on vertices
808  mp4Vector gradVerts[8];
809  gradVerts[3] = prevGradVerts[0]; gradVerts[2] = prevGradVerts[1];
810  gradVerts[1] = prevGradVerts[2]; gradVerts[0] = prevGradVerts[3];
811  //initialize gradients on edges from the last cube
812  mpVector grads[12];
813  grads[2] = prevGrads[0]; grads[1] = prevGrads[1]; grads[0] = prevGrads[2]; grads[3] = prevGrads[3];
814  //for test if this cube is intersected
815  int oldNumTriangles = numTriangles;
816  //run marching cubes on this cube
817  triangles = MarchOneCube(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, ind, i, j, k,
818  minValue, points, triangles, numTriangles, verts, intVerts, edgeIndex,
819  gradVerts, grads, gradIndex);
820  //check if this cube is intersected
821  if(numTriangles == oldNumTriangles)
822  return triangles;
823  //two new indexes formed for each face
824  int passEdgeIndex, passGradIndex;
825  //run recursive functions on each surface
826  MC_FACE0
827  MC_FACE1
828  MC_FACE2
829  MC_FACE3
830  MC_FACE4
831  return triangles;
832 }
833 
834 //SURFACE 5 - Cube ran on surface 5 of previous cube. Recieving side: 4.
835 TRIANGLE* MCFace5(int ncellsX, int ncellsY, int ncellsZ,
836  float gradFactorX, float gradFactorY, float gradFactorZ,
837  int ind, int i, int j, int k,
838  float minValue, mp4Vector * points, TRIANGLE *triangles, int &numTriangles,
839  mp4Vector prevVerts[4], mpVector prevIntVerts[4], int edgeIndex,
840  mp4Vector prevGradVerts[4], mpVector prevGrads[4], int gradIndex, bool* marchedCubes)
841 {
842  //first check if not outside the region
843  if(i <= 0 || i >= ncellsX-1 || j <= 0 || j >= ncellsY-1 || k <= 0 || k >= ncellsZ-1)
844  return triangles;
845  //make sure we save that this cube was marched through
846  marchedCubes[ind] = TRUE;
847  //initialize vertices
848  mp4Vector verts[8];
849  verts[0] = points[ind];
850  verts[1] = points[ind + YtimeZ];
851  verts[2] = points[ind + YtimeZ + 1];
852  verts[3] = points[ind + 1];
853  verts[4] = prevVerts[3];
854  verts[5] = prevVerts[2];
855  verts[6] = prevVerts[1];
856  verts[7] = prevVerts[0];
857  //initialize edges from the last cube
858  mpVector intVerts[12];
859  intVerts[6] = prevIntVerts[0]; intVerts[5] = prevIntVerts[1];
860  intVerts[4] = prevIntVerts[2]; intVerts[7] = prevIntVerts[3];
861  //initialize gradients on vertices
862  mp4Vector gradVerts[8];
863  gradVerts[7] = prevGradVerts[0]; gradVerts[6] = prevGradVerts[1];
864  gradVerts[5] = prevGradVerts[2]; gradVerts[4] = prevGradVerts[3];
865  //initialize gradients on edges from the last cube
866  mpVector grads[12];
867  grads[6] = prevGrads[0]; grads[5] = prevGrads[1]; grads[4] = prevGrads[2]; grads[7] = prevGrads[3];
868  //for test if this cube is intersected
869  int oldNumTriangles = numTriangles;
870  //run marching cubes on this cube
871  triangles = MarchOneCube(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ, ind, i, j, k,
872  minValue, points, triangles, numTriangles, verts, intVerts, edgeIndex,
873  gradVerts, grads, gradIndex);
874  //check if this cube is intersected
875  if(numTriangles == oldNumTriangles)
876  return triangles;
877  //two new indexes formed for each face
878  int passEdgeIndex, passGradIndex;
879  //run recursive functions on each surface
880  MC_FACE0
881  MC_FACE1
882  MC_FACE2
883  MC_FACE3
884  MC_FACE5
885  return triangles;
886 }
887 
888 //Marching Cubes on a single cube (i, j, k)
889 // Verts should be initialized before hand
890 // Global variables YtimeZ and pointsZ should be defined and initialized
891 TRIANGLE* MarchOneCube(int ncellsX, int ncellsY, int ncellsZ,
892  float gradFactorX, float gradFactorY, float gradFactorZ,
893  int ind, int i, int j, int k,
894  float minValue, mp4Vector * points, TRIANGLE *triangles, int &numTriangles,
895  mp4Vector verts[8], mpVector intVerts[12], int &edgeIndex,
896  mp4Vector gradVerts[8], mpVector grads[12], int &indGrad)
897 {
898  //factor by which gradients are scaled
899  mpVector factor(1.0/(2.0*gradFactorX), 1.0/(2.0*gradFactorY), 1.0/(2.0*gradFactorZ));
900  //get the index
901  int cubeIndex = int(0);
902  for(int n=0; n < 8; n++) if(verts[n].val <= minValue) cubeIndex |= (1 << n);
903  //check if its completely inside or outside
904  if(!edgeTable[cubeIndex]) return triangles;
905 
906  int prevEdgeIndex = edgeIndex;
907  edgeIndex = edgeTable[cubeIndex];
908 
909  if((edgeIndex & 1) && !(prevEdgeIndex & 1)) {
910  intVerts[0] = LinearInterp(verts[0], verts[1], minValue);
911  if(! (indGrad & 1)) {
912  gradVerts[0] = CALC_GRAD_VERT_0(verts)
913  indGrad |= 1;
914  }
915  if(! (indGrad & 2)) {
916  gradVerts[1] = CALC_GRAD_VERT_1(verts)
917  indGrad |= 2;
918  }
919  grads[0] = LinearInterp(gradVerts[0], gradVerts[1], minValue);
920  grads[0].x *= factor.x; grads[0].y *= factor.y; grads[0].z *= factor.z;
921  }
922  if((edgeIndex & 2) && !(prevEdgeIndex & 2)) {
923  intVerts[1] = LinearInterp(verts[1], verts[2], minValue);
924  if(! (indGrad & 2)) {
925  gradVerts[1] = CALC_GRAD_VERT_1(verts)
926  indGrad |= 2;
927  }
928  if(! (indGrad & 4)) {
929  gradVerts[2] = CALC_GRAD_VERT_2(verts)
930  indGrad |= 4;
931  }
932  grads[1] = LinearInterp(gradVerts[1], gradVerts[2], minValue);
933  grads[1].x *= factor.x; grads[1].y *= factor.y; grads[1].z *= factor.z;
934  }
935  if(edgeIndex & 4) {
936  intVerts[2] = LinearInterp(verts[2], verts[3], minValue);
937  if(! (indGrad & 4)) {
938  gradVerts[2] = CALC_GRAD_VERT_2(verts)
939  indGrad |= 4;
940  }
941  if(! (indGrad & 8)) {
942  gradVerts[3] = CALC_GRAD_VERT_3(verts)
943  indGrad |= 8;
944  }
945  grads[2] = LinearInterp(gradVerts[2], gradVerts[3], minValue);
946  grads[2].x *= factor.x; grads[2].y *= factor.y; grads[2].z *= factor.z;
947  }
948  if((edgeIndex & 8) && !(prevEdgeIndex & 8)) {
949  intVerts[3] = LinearInterp(verts[3], verts[0], minValue);
950  if(! (indGrad & 8)) {
951  gradVerts[3] = CALC_GRAD_VERT_3(verts)
952  indGrad |= 8;
953  }
954  if(! (indGrad & 1)) {
955  gradVerts[0] = CALC_GRAD_VERT_0(verts)
956  indGrad |= 1;
957  }
958  grads[3] = LinearInterp(gradVerts[3], gradVerts[0], minValue);
959  grads[3].x *= factor.x; grads[3].y *= factor.y; grads[3].z *= factor.z;
960  }
961  if((edgeIndex & 16) && !(prevEdgeIndex & 16)) {
962  intVerts[4] = LinearInterp(verts[4], verts[5], minValue);
963  if(! (indGrad & 16)) {
964  gradVerts[4] = CALC_GRAD_VERT_4(verts)
965  indGrad |= 16;
966  }
967  if(! (indGrad & 32)) {
968  gradVerts[5] = CALC_GRAD_VERT_5(verts)
969  indGrad |= 32;
970  }
971  grads[4] = LinearInterp(gradVerts[4], gradVerts[5], minValue);
972  grads[4].x *= factor.x; grads[4].y *= factor.y; grads[4].z *= factor.z;
973  }
974  if((edgeIndex & 32) && !(prevEdgeIndex & 32)) {
975  intVerts[5] = LinearInterp(verts[5], verts[6], minValue);
976  if(! (indGrad & 32)) {
977  gradVerts[5] = CALC_GRAD_VERT_5(verts)
978  indGrad |= 32;
979  }
980  if(! (indGrad & 64)) {
981  gradVerts[6] = CALC_GRAD_VERT_6(verts)
982  indGrad |= 64;
983  }
984  grads[5] = LinearInterp(gradVerts[5], gradVerts[6], minValue);
985  grads[5].x *= factor.x; grads[5].y *= factor.y; grads[5].z *= factor.z;
986  }
987  if((edgeIndex & 64) && !(prevEdgeIndex & 64)) {
988  intVerts[6] = LinearInterp(verts[6], verts[7], minValue);
989  if(! (indGrad & 64)) {
990  gradVerts[6] = CALC_GRAD_VERT_6(verts)
991  indGrad |= 64;
992  }
993  if(! (indGrad & 128)) {
994  gradVerts[7] = CALC_GRAD_VERT_7(verts)
995  indGrad |= 128;
996  }
997  grads[6] = LinearInterp(gradVerts[6], gradVerts[7], minValue);
998  grads[6].x *= factor.x; grads[6].y *= factor.y; grads[6].z *= factor.z;
999  }
1000  if((edgeIndex & 128) && !(prevEdgeIndex & 128)) {
1001  intVerts[7] = LinearInterp(verts[7], verts[4], minValue);
1002  if(! (indGrad & 128)) {
1003  gradVerts[7] = CALC_GRAD_VERT_7(verts)
1004  indGrad |= 128;
1005  }
1006  if(! (indGrad & 16)) {
1007  gradVerts[4] = CALC_GRAD_VERT_4(verts)
1008  indGrad |= 16;
1009  }
1010  grads[7] = LinearInterp(gradVerts[7], gradVerts[4], minValue);
1011  grads[7].x *= factor.x; grads[7].y *= factor.y; grads[7].z *= factor.z;
1012  }
1013  if((edgeIndex & 256) && !(prevEdgeIndex & 256)) {
1014  intVerts[8] = LinearInterp(verts[0], verts[4], minValue);
1015  if(! (indGrad & 1)) {
1016  gradVerts[0] = CALC_GRAD_VERT_0(verts)
1017  indGrad |= 1;
1018  }
1019  if(! (indGrad & 16)) {
1020  gradVerts[4] = CALC_GRAD_VERT_4(verts)
1021  indGrad |= 16;
1022  }
1023  grads[8] = LinearInterp(gradVerts[0], gradVerts[4], minValue);
1024  grads[8].x *= factor.x; grads[8].y *= factor.y; grads[8].z *= factor.z;
1025  }
1026  if((edgeIndex & 512) && !(prevEdgeIndex & 512)) {
1027  intVerts[9] = LinearInterp(verts[1], verts[5], minValue);
1028  if(! (indGrad & 2)) {
1029  gradVerts[1] = CALC_GRAD_VERT_1(verts)
1030  indGrad |= 2;
1031  }
1032  if(! (indGrad & 32)) {
1033  gradVerts[5] = CALC_GRAD_VERT_5(verts)
1034  indGrad |= 32;
1035  }
1036  grads[9] = LinearInterp(gradVerts[1], gradVerts[5], minValue);
1037  grads[9].x *= factor.x; grads[9].y *= factor.y; grads[9].z *= factor.z;
1038  }
1039  if((edgeIndex & 1024) && !(prevEdgeIndex & 1024)) {
1040  intVerts[10] = LinearInterp(verts[2], verts[6], minValue);
1041  if(! (indGrad & 4)) {
1042  gradVerts[2] = CALC_GRAD_VERT_2(verts)
1043  indGrad |= 4;
1044  }
1045  if(! (indGrad & 64)) {
1046  gradVerts[6] = CALC_GRAD_VERT_6(verts)
1047  indGrad |= 64;
1048  }
1049  grads[10] = LinearInterp(gradVerts[2], gradVerts[6], minValue);
1050  grads[10].x *= factor.x; grads[10].y *= factor.y; grads[10].z *= factor.z;
1051  }
1052  if((edgeIndex & 2048) && !(prevEdgeIndex & 2048)) {
1053  intVerts[11] = LinearInterp(verts[3], verts[7], minValue);
1054  if(! (indGrad & 8)) {
1055  gradVerts[3] = CALC_GRAD_VERT_3(verts)
1056  indGrad |= 8;
1057  }
1058  if(! (indGrad & 128)) {
1059  gradVerts[7] = CALC_GRAD_VERT_7(verts)
1060  indGrad |= 128;
1061  }
1062  grads[11] = LinearInterp(gradVerts[3], gradVerts[7], minValue);
1063  grads[11].x *= factor.x; grads[11].y *= factor.y; grads[11].z *= factor.z;
1064  }
1065  //now build the triangles using triTable and add them to the triangle array
1066  for (int n=0; triTable[cubeIndex][n] != -1; n+=3) {
1067  int index[3] = {triTable[cubeIndex][n+2], triTable[cubeIndex][n+1], triTable[cubeIndex][n]};
1068  for(int h=0; h < 3; h++) {
1069  triangles[numTriangles].p[h] = intVerts[index[h]];
1070  triangles[numTriangles].norm[h] = grads[index[h]];
1071  }
1072  numTriangles++; //one more triangle...
1073  }
1074  return triangles;
1075 }
1076 
1077 
1078 
1079 //Finds first intersecting cube and returns its indices (or -1 if nothing is found)
1080 float* MCFind(int ncellsX, int ncellsY, int ncellsZ, float minValue, mp4Vector * points)
1081 {
1082  pointsZ = ncellsZ+1; //initializes global variables
1083  YtimeZ = (ncellsY+1)*pointsZ;
1084  int lastX = ncellsX - 1, lastY = ncellsY - 1, lastZ = ncellsZ - 1; //for a little extra speed
1085  mp4Vector *verts[8]; //store address of each point rather than copying x,y,z, and val
1086  int cubeIndex, ind, ni;
1087  float *found = new float[3];//returned indices: initialized to -1
1088  found[0] = -1; found[0] = -1; found[0] = -1;
1089 
1090  for(int i=1; i < lastX; i++) {
1091  ni = i*YtimeZ;
1092  for(int j=1; j < lastY; j++) {
1093  //get the first 4 vertices (0, 1, 4, 5) and place them into verts at indices 3, 2, 7, 6
1094  ind = ni + j*pointsZ + 1;
1095  verts[3] = &points[ind];
1096  verts[2] = &points[ind + YtimeZ];
1097  verts[7] = &points[ind + pointsZ];
1098  verts[6] = &points[ind + YtimeZ + pointsZ];
1099  for(int k=1; k < lastZ; k++, ind++) {
1100  //initialize vertices
1101  verts[0] = verts[3]; //these are saved from the last iteration
1102  verts[1] = verts[2];
1103  verts[4] = verts[7];
1104  verts[5] = verts[6];
1105  verts[4] = &points[ind + pointsZ];
1106  verts[5] = &points[ind + YtimeZ + pointsZ];
1107  verts[6] = &points[ind + YtimeZ + ncellsZ + 2];
1108  verts[7] = &points[ind + ncellsZ + 2];
1109  //build index - shows which vertices are outside/inside
1110  cubeIndex = int(0);
1111  for(int n=0; n < 8; n++) if(verts[n]->val <= minValue) cubeIndex |= (1 << n);
1112  if(edgeTable[cubeIndex]) { //if found then return them
1113  found[0] = i; found[1] = j; found[2] = k;
1114  return found;
1115  }
1116  }
1117  }
1118  }
1119  return found;
1120 }
1121 
1122 //First finds intersecting cube and runs recursive Marching Cubes.
1123 TRIANGLE* MCRecFind(int ncellsX, int ncellsY, int ncellsZ,
1124  float gradFactorX, float gradFactorY, float gradFactorZ,
1125  float minValue, mp4Vector * points, int &numTriangles)
1126 {
1127  float *found = MCFind(ncellsX, ncellsY, ncellsZ, minValue, points);
1128  if(found[0] == -1) { //if nothing is found return NULL
1129  numTriangles = 0;
1130  return NULL;
1131  }
1132  int i[1] = {(int)found[0]}, j[1] = {(int)found[1]}, k[1] = {(int)found[2]};
1133  delete [] found;
1134  return MarchingCubesRec(ncellsX, ncellsY, ncellsZ, gradFactorX, gradFactorY, gradFactorZ,
1135  1, i, j, k, minValue, points, numTriangles);
1136 }
1137 
TRIANGLE
Definition: MC.h:22