Lines 908
##### Keywords
marching cubes (1) marching tetrahedrons (1)
##### Permissions
Owner: Stou S.
Viewable by Everyone
Editable by All Siafoo Users
Siafoo is here to make coding less frustrating and to save you time.

# An implementation of Marching Cubes and Marching Tetrahedrons 1

 In Brief A self contained implement of the Marching Cubes and Marching Tetrahedrons algorithms for creating surfaces from volume data.... more
 Language C
# 's
`   1//   2// Marching Cubes Example Program    3// by Cory Bloyd (corysama@yahoo.com)   4//   5// A simple, portable and complete implementation of the Marching Cubes   6// and Marching Tetrahedrons algorithms in a single source file.   7// There are many ways that this code could be made faster, but the    8// intent is for the code to be easy to understand.   9//  10// For a description of the algorithm go to  11// http://astronomy.swin.edu.au/pbourke/modelling/polygonise/  12//  13// This code is public domain.  14//  15  16#include "stdio.h"  17#include "math.h"  18//This program requires the OpenGL and GLUT libraries  19// You can obtain them for free from http://www.opengl.org  20#include "GL/glut.h"  21  22struct GLvector  23{  24        GLfloat fX;  25        GLfloat fY;  26        GLfloat fZ;       27};  28  29//These tables are used so that everything can be done in little loops that you can look at all at once  30// rather than in pages and pages of unrolled code.  31  32//a2fVertexOffset lists the positions, relative to vertex0, of each of the 8 vertices of a cube  33static const GLfloat a2fVertexOffset[8][3] =  34{  35        {0.0, 0.0, 0.0},{1.0, 0.0, 0.0},{1.0, 1.0, 0.0},{0.0, 1.0, 0.0},  36        {0.0, 0.0, 1.0},{1.0, 0.0, 1.0},{1.0, 1.0, 1.0},{0.0, 1.0, 1.0}  37};  38  39//a2iEdgeConnection lists the index of the endpoint vertices for each of the 12 edges of the cube  40static const GLint a2iEdgeConnection[12][2] =   41{  42        {0,1}, {1,2}, {2,3}, {3,0},  43        {4,5}, {5,6}, {6,7}, {7,4},  44        {0,4}, {1,5}, {2,6}, {3,7}  45};  46  47//a2fEdgeDirection lists the direction vector (vertex1-vertex0) for each edge in the cube  48static const GLfloat a2fEdgeDirection[12][3] =  49{  50        {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},  51        {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},  52        {0.0, 0.0, 1.0},{0.0, 0.0, 1.0},{ 0.0, 0.0, 1.0},{0.0,  0.0, 1.0}  53};  54  55//a2iTetrahedronEdgeConnection lists the index of the endpoint vertices for each of the 6 edges of the tetrahedron  56static const GLint a2iTetrahedronEdgeConnection[6][2] =  57{  58        {0,1},  {1,2},  {2,0},  {0,3},  {1,3},  {2,3}  59};  60  61//a2iTetrahedronEdgeConnection lists the index of verticies from a cube   62// that made up each of the six tetrahedrons within the cube  63static const GLint a2iTetrahedronsInACube[6][4] =  64{  65        {0,5,1,6},  66        {0,1,2,6},  67        {0,2,3,6},  68        {0,3,7,6},  69        {0,7,4,6},  70        {0,4,5,6},  71};  72  73static const GLfloat afAmbientWhite [] = {0.25, 0.25, 0.25, 1.00};   74static const GLfloat afAmbientRed   [] = {0.25, 0.00, 0.00, 1.00};   75static const GLfloat afAmbientGreen [] = {0.00, 0.25, 0.00, 1.00};   76static const GLfloat afAmbientBlue  [] = {0.00, 0.00, 0.25, 1.00};   77static const GLfloat afDiffuseWhite [] = {0.75, 0.75, 0.75, 1.00};   78static const GLfloat afDiffuseRed   [] = {0.75, 0.00, 0.00, 1.00};   79static const GLfloat afDiffuseGreen [] = {0.00, 0.75, 0.00, 1.00};   80static const GLfloat afDiffuseBlue  [] = {0.00, 0.00, 0.75, 1.00};   81static const GLfloat afSpecularWhite[] = {1.00, 1.00, 1.00, 1.00};   82static const GLfloat afSpecularRed  [] = {1.00, 0.25, 0.25, 1.00};   83static const GLfloat afSpecularGreen[] = {0.25, 1.00, 0.25, 1.00};   84static const GLfloat afSpecularBlue [] = {0.25, 0.25, 1.00, 1.00};   85  86  87GLenum    ePolygonMode = GL_FILL;  88GLint     iDataSetSize = 16;  89GLfloat   fStepSize = 1.0/iDataSetSize;  90GLfloat   fTargetValue = 48.0;  91GLfloat   fTime = 0.0;  92GLvector  sSourcePoint[3];  93GLboolean bSpin = true;  94GLboolean bMove = true;  95GLboolean bLight = true;  96  97  98void vIdle();  99void vDrawScene();  100void vResize(GLsizei, GLsizei); 101void vKeyboard(unsigned char cKey, int iX, int iY); 102void vSpecial(int iKey, int iX, int iY); 103 104GLvoid vPrintHelp(); 105GLvoid vSetTime(GLfloat fTime); 106GLfloat fSample1(GLfloat fX, GLfloat fY, GLfloat fZ); 107GLfloat fSample2(GLfloat fX, GLfloat fY, GLfloat fZ); 108GLfloat fSample3(GLfloat fX, GLfloat fY, GLfloat fZ); 109GLfloat (*fSample)(GLfloat fX, GLfloat fY, GLfloat fZ) = fSample1; 110 111GLvoid vMarchingCubes(); 112GLvoid vMarchCube1(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale); 113GLvoid vMarchCube2(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale); 114GLvoid (*vMarchCube)(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale) = vMarchCube1; 115 116void main(int argc, char **argv)  117{  118        GLfloat afPropertiesAmbient [] = {0.50, 0.50, 0.50, 1.00};  119        GLfloat afPropertiesDiffuse [] = {0.75, 0.75, 0.75, 1.00};  120        GLfloat afPropertiesSpecular[] = {1.00, 1.00, 1.00, 1.00};  121 122        GLsizei iWidth = 640.0;  123        GLsizei iHeight = 480.0;  124 125        glutInit(&argc, argv); 126        glutInitWindowPosition( 0, 0); 127        glutInitWindowSize(iWidth, iHeight); 128        glutInitDisplayMode( GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE ); 129        glutCreateWindow( "Marching Cubes" ); 130        glutDisplayFunc( vDrawScene ); 131        glutIdleFunc( vIdle ); 132        glutReshapeFunc( vResize ); 133        glutKeyboardFunc( vKeyboard ); 134        glutSpecialFunc( vSpecial ); 135 136        glClearColor( 0.0, 0.0, 0.0, 1.0 );  137        glClearDepth( 1.0 );  138 139        glEnable(GL_DEPTH_TEST);  140        glEnable(GL_LIGHTING); 141        glPolygonMode(GL_FRONT_AND_BACK, ePolygonMode); 142 143        glLightfv( GL_LIGHT0, GL_AMBIENT,  afPropertiesAmbient);  144        glLightfv( GL_LIGHT0, GL_DIFFUSE,  afPropertiesDiffuse);  145        glLightfv( GL_LIGHT0, GL_SPECULAR, afPropertiesSpecular);  146        glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, 1.0);  147 148        glEnable( GL_LIGHT0 );  149 150        glMaterialfv(GL_BACK,  GL_AMBIENT,   afAmbientGreen);  151        glMaterialfv(GL_BACK,  GL_DIFFUSE,   afDiffuseGreen);  152        glMaterialfv(GL_FRONT, GL_AMBIENT,   afAmbientBlue);  153        glMaterialfv(GL_FRONT, GL_DIFFUSE,   afDiffuseBlue);  154        glMaterialfv(GL_FRONT, GL_SPECULAR,  afSpecularWhite);  155        glMaterialf( GL_FRONT, GL_SHININESS, 25.0);  156 157        vResize(iWidth, iHeight);  158 159        vPrintHelp(); 160        glutMainLoop();  161} 162 163GLvoid vPrintHelp() 164{ 165        printf("Marching Cubes Example by Cory Bloyd (dejaspaminacan@my-deja.com)\n\n"); 166 167        printf("+/-  increase/decrease sample density\n"); 168        printf("PageUp/PageDown  increase/decrease surface value\n"); 169        printf("s  change sample function\n"); 170        printf("c  toggle marching cubes / marching tetrahedrons\n"); 171        printf("w  wireframe on/off\n"); 172        printf("l  toggle lighting / color-by-normal\n"); 173        printf("Home  spin scene on/off\n"); 174        printf("End  source point animation on/off\n"); 175} 176 177 178void vResize( GLsizei iWidth, GLsizei iHeight )  179{  180        GLfloat fAspect, fHalfWorldSize = (1.4142135623730950488016887242097/2);  181 182        glViewport( 0, 0, iWidth, iHeight );  183        glMatrixMode (GL_PROJECTION); 184        glLoadIdentity (); 185 186        if(iWidth <= iHeight) 187        { 188                fAspect = (GLfloat)iHeight / (GLfloat)iWidth; 189                glOrtho(-fHalfWorldSize, fHalfWorldSize, -fHalfWorldSize*fAspect, 190                        fHalfWorldSize*fAspect, -10*fHalfWorldSize, 10*fHalfWorldSize); 191        } 192        else 193        { 194                fAspect = (GLfloat)iWidth / (GLfloat)iHeight;  195                glOrtho(-fHalfWorldSize*fAspect, fHalfWorldSize*fAspect, -fHalfWorldSize, 196                        fHalfWorldSize, -10*fHalfWorldSize, 10*fHalfWorldSize); 197        } 198  199        glMatrixMode( GL_MODELVIEW ); 200} 201 202void vKeyboard(unsigned char cKey, int iX, int iY) 203{ 204        switch(cKey) 205        { 206                case 'w' : 207                { 208                        if(ePolygonMode == GL_LINE) 209                        { 210                                ePolygonMode = GL_FILL; 211                        } 212                        else 213                        { 214                                ePolygonMode = GL_LINE; 215                        } 216                        glPolygonMode(GL_FRONT_AND_BACK, ePolygonMode); 217                } break; 218                case '+' : 219                case '=' : 220                { 221                        ++iDataSetSize; 222                        fStepSize = 1.0/iDataSetSize; 223                } break; 224                case '-' : 225                { 226                        if(iDataSetSize > 1) 227                        { 228                                --iDataSetSize; 229                                fStepSize = 1.0/iDataSetSize; 230                        } 231                } break; 232                case 'c' : 233                { 234                        if(vMarchCube == vMarchCube1) 235                        { 236                                vMarchCube = vMarchCube2;//Use Marching Tetrahedrons 237                        } 238                        else 239                        { 240                                vMarchCube = vMarchCube1;//Use Marching Cubes 241                        } 242                } break; 243                case 's' : 244                { 245                        if(fSample == fSample1) 246                        { 247                                fSample = fSample2; 248                        } 249                        else if(fSample == fSample2) 250                        { 251                                fSample = fSample3; 252                        } 253                        else 254                        { 255                                fSample = fSample1; 256                        } 257                } break; 258                case 'l' : 259                { 260                        if(bLight) 261                        { 262                                glDisable(GL_LIGHTING);//use vertex colors 263                        } 264                        else 265                        { 266                                glEnable(GL_LIGHTING);//use lit material color 267                        } 268 269                        bLight = !bLight; 270                }; 271        } 272} 273 274 275void vSpecial(int iKey, int iX, int iY) 276{ 277        switch(iKey) 278        { 279                case GLUT_KEY_PAGE_UP : 280                { 281                        if(fTargetValue < 1000.0) 282                        { 283                                fTargetValue *= 1.1; 284                        } 285                } break; 286                case GLUT_KEY_PAGE_DOWN : 287                { 288                        if(fTargetValue > 1.0) 289                        { 290                                fTargetValue /= 1.1; 291                        } 292                } break; 293                case GLUT_KEY_HOME : 294                { 295                        bSpin = !bSpin; 296                } break; 297                case GLUT_KEY_END : 298                { 299                        bMove = !bMove; 300                } break; 301        } 302} 303 304void vIdle() 305{ 306        glutPostRedisplay(); 307} 308 309void vDrawScene()  310{  311        static GLfloat fPitch = 0.0; 312        static GLfloat fYaw   = 0.0; 313        static GLfloat fTime = 0.0; 314 315        glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );  316 317        glPushMatrix();  318 319        if(bSpin) 320        { 321                fPitch += 4.0; 322                fYaw   += 2.5; 323        } 324        if(bMove) 325        { 326                fTime  += 0.025; 327        } 328 329        vSetTime(fTime); 330 331        glTranslatef(0.0, 0.0, -1.0);   332        glRotatef( -fPitch, 1.0, 0.0, 0.0); 333        glRotatef(     0.0, 0.0, 1.0, 0.0); 334        glRotatef(    fYaw, 0.0, 0.0, 1.0); 335 336        glPushAttrib(GL_LIGHTING_BIT); 337                glDisable(GL_LIGHTING); 338                glColor3f(1.0, 1.0, 1.0); 339                glutWireCube(1.0);  340        glPopAttrib();  341 342 343        glPushMatrix();  344        glTranslatef(-0.5, -0.5, -0.5); 345        glBegin(GL_TRIANGLES); 346                vMarchingCubes(); 347        glEnd(); 348        glPopMatrix();  349 350 351        glPopMatrix();  352 353        glutSwapBuffers();  354} 355 356//fGetOffset finds the approximate point of intersection of the surface 357// between two points with the values fValue1 and fValue2 358GLfloat fGetOffset(GLfloat fValue1, GLfloat fValue2, GLfloat fValueDesired) 359{ 360        GLdouble fDelta = fValue2 - fValue1; 361 362        if(fDelta == 0.0) 363        { 364                return 0.5; 365        } 366        return (fValueDesired - fValue1)/fDelta; 367} 368 369 370//vGetColor generates a color from a given position and normal of a point 371GLvoid vGetColor(GLvector &rfColor, GLvector &rfPosition, GLvector &rfNormal) 372{ 373        GLfloat fX = rfNormal.fX; 374        GLfloat fY = rfNormal.fY; 375        GLfloat fZ = rfNormal.fZ; 376        rfColor.fX = (fX > 0.0 ? fX : 0.0) + (fY < 0.0 ? -0.5*fY : 0.0) + (fZ < 0.0 ? -0.5*fZ : 0.0); 377        rfColor.fY = (fY > 0.0 ? fY : 0.0) + (fZ < 0.0 ? -0.5*fZ : 0.0) + (fX < 0.0 ? -0.5*fX : 0.0); 378        rfColor.fZ = (fZ > 0.0 ? fZ : 0.0) + (fX < 0.0 ? -0.5*fX : 0.0) + (fY < 0.0 ? -0.5*fY : 0.0); 379} 380 381GLvoid vNormalizeVector(GLvector &rfVectorResult, GLvector &rfVectorSource) 382{ 383        GLfloat fOldLength; 384        GLfloat fScale; 385 386        fOldLength = sqrtf( (rfVectorSource.fX * rfVectorSource.fX) + 387                            (rfVectorSource.fY * rfVectorSource.fY) + 388                            (rfVectorSource.fZ * rfVectorSource.fZ) ); 389 390        if(fOldLength == 0.0) 391        { 392                rfVectorResult.fX = rfVectorSource.fX; 393                rfVectorResult.fY = rfVectorSource.fY; 394                rfVectorResult.fZ = rfVectorSource.fZ; 395        } 396        else 397        { 398                fScale = 1.0/fOldLength; 399                rfVectorResult.fX = rfVectorSource.fX*fScale; 400                rfVectorResult.fY = rfVectorSource.fY*fScale; 401                rfVectorResult.fZ = rfVectorSource.fZ*fScale; 402        } 403} 404 405 406//Generate a sample data set.  fSample1(), fSample2() and fSample3() define three scalar fields whose 407// values vary by the X,Y and Z coordinates and by the fTime value set by vSetTime() 408GLvoid vSetTime(GLfloat fNewTime) 409{ 410        GLfloat fOffset; 411        GLint iSourceNum; 412 413        for(iSourceNum = 0; iSourceNum < 3; iSourceNum++) 414        { 415                sSourcePoint[iSourceNum].fX = 0.5; 416                sSourcePoint[iSourceNum].fY = 0.5; 417                sSourcePoint[iSourceNum].fZ = 0.5; 418        } 419 420        fTime = fNewTime; 421        fOffset = 1.0 + sinf(fTime); 422        sSourcePoint[0].fX *= fOffset; 423        sSourcePoint[1].fY *= fOffset; 424        sSourcePoint[2].fZ *= fOffset; 425} 426 427//fSample1 finds the distance of (fX, fY, fZ) from three moving points 428GLfloat fSample1(GLfloat fX, GLfloat fY, GLfloat fZ) 429{ 430        GLdouble fResult = 0.0; 431        GLdouble fDx, fDy, fDz; 432        fDx = fX - sSourcePoint[0].fX; 433        fDy = fY - sSourcePoint[0].fY; 434        fDz = fZ - sSourcePoint[0].fZ; 435        fResult += 0.5/(fDx*fDx + fDy*fDy + fDz*fDz); 436 437        fDx = fX - sSourcePoint[1].fX; 438        fDy = fY - sSourcePoint[1].fY; 439        fDz = fZ - sSourcePoint[1].fZ; 440        fResult += 1.0/(fDx*fDx + fDy*fDy + fDz*fDz); 441 442        fDx = fX - sSourcePoint[2].fX; 443        fDy = fY - sSourcePoint[2].fY; 444        fDz = fZ - sSourcePoint[2].fZ; 445        fResult += 1.5/(fDx*fDx + fDy*fDy + fDz*fDz); 446 447        return fResult; 448} 449 450//fSample2 finds the distance of (fX, fY, fZ) from three moving lines 451GLfloat fSample2(GLfloat fX, GLfloat fY, GLfloat fZ) 452{ 453        GLdouble fResult = 0.0; 454        GLdouble fDx, fDy, fDz; 455        fDx = fX - sSourcePoint[0].fX; 456        fDy = fY - sSourcePoint[0].fY; 457        fResult += 0.5/(fDx*fDx + fDy*fDy); 458 459        fDx = fX - sSourcePoint[1].fX; 460        fDz = fZ - sSourcePoint[1].fZ; 461        fResult += 0.75/(fDx*fDx + fDz*fDz); 462 463        fDy = fY - sSourcePoint[2].fY; 464        fDz = fZ - sSourcePoint[2].fZ; 465        fResult += 1.0/(fDy*fDy + fDz*fDz); 466 467        return fResult; 468} 469 470 471//fSample2 defines a height field by plugging the distance from the center into the sin and cos functions 472GLfloat fSample3(GLfloat fX, GLfloat fY, GLfloat fZ) 473{ 474        GLfloat fHeight = 20.0*(fTime + sqrt((0.5-fX)*(0.5-fX) + (0.5-fY)*(0.5-fY))); 475        fHeight = 1.5 + 0.1*(sinf(fHeight) + cosf(fHeight)); 476        GLdouble fResult = (fHeight - fZ)*50.0; 477 478        return fResult; 479} 480 481 482//vGetNormal() finds the gradient of the scalar field at a point 483//This gradient can be used as a very accurate vertx normal for lighting calculations 484GLvoid vGetNormal(GLvector &rfNormal, GLfloat fX, GLfloat fY, GLfloat fZ) 485{ 486        rfNormal.fX = fSample(fX-0.01, fY, fZ) - fSample(fX+0.01, fY, fZ); 487        rfNormal.fY = fSample(fX, fY-0.01, fZ) - fSample(fX, fY+0.01, fZ); 488        rfNormal.fZ = fSample(fX, fY, fZ-0.01) - fSample(fX, fY, fZ+0.01); 489        vNormalizeVector(rfNormal, rfNormal); 490} 491 492 493//vMarchCube1 performs the Marching Cubes algorithm on a single cube 494GLvoid vMarchCube1(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale) 495{ 496        extern GLint aiCubeEdgeFlags[256]; 497        extern GLint a2iTriangleConnectionTable[256][16]; 498 499        GLint iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags; 500        GLfloat fOffset; 501        GLvector sColor; 502        GLfloat afCubeValue[8]; 503        GLvector asEdgeVertex[12]; 504        GLvector asEdgeNorm[12]; 505 506        //Make a local copy of the values at the cube's corners 507        for(iVertex = 0; iVertex < 8; iVertex++) 508        { 509                afCubeValue[iVertex] = fSample(fX + a2fVertexOffset[iVertex][0]*fScale, 510                                                   fY + a2fVertexOffset[iVertex][1]*fScale, 511                                                   fZ + a2fVertexOffset[iVertex][2]*fScale); 512        } 513 514        //Find which vertices are inside of the surface and which are outside 515        iFlagIndex = 0; 516        for(iVertexTest = 0; iVertexTest < 8; iVertexTest++) 517        { 518                if(afCubeValue[iVertexTest] <= fTargetValue)  519                        iFlagIndex |= 1<<iVertexTest; 520        } 521 522        //Find which edges are intersected by the surface 523        iEdgeFlags = aiCubeEdgeFlags[iFlagIndex]; 524 525        //If the cube is entirely inside or outside of the surface, then there will be no intersections 526        if(iEdgeFlags == 0)  527        { 528                return; 529        } 530 531        //Find the point of intersection of the surface with each edge 532        //Then find the normal to the surface at those points 533        for(iEdge = 0; iEdge < 12; iEdge++) 534        { 535                //if there is an intersection on this edge 536                if(iEdgeFlags & (1<<iEdge)) 537                { 538                        fOffset = fGetOffset(afCubeValue[ a2iEdgeConnection[iEdge][0] ],  539                                                     afCubeValue[ a2iEdgeConnection[iEdge][1] ], fTargetValue); 540 541                        asEdgeVertex[iEdge].fX = fX + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][0]  +  fOffset * a2fEdgeDirection[iEdge][0]) * fScale; 542                        asEdgeVertex[iEdge].fY = fY + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][1]  +  fOffset * a2fEdgeDirection[iEdge][1]) * fScale; 543                        asEdgeVertex[iEdge].fZ = fZ + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][2]  +  fOffset * a2fEdgeDirection[iEdge][2]) * fScale; 544 545                        vGetNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ); 546                } 547        } 548 549 550        //Draw the triangles that were found.  There can be up to five per cube 551        for(iTriangle = 0; iTriangle < 5; iTriangle++) 552        { 553                if(a2iTriangleConnectionTable[iFlagIndex][3*iTriangle] < 0) 554                        break; 555 556                for(iCorner = 0; iCorner < 3; iCorner++) 557                { 558                        iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner]; 559 560                        vGetColor(sColor, asEdgeVertex[iVertex], asEdgeNorm[iVertex]); 561                        glColor3f(sColor.fX, sColor.fY, sColor.fZ); 562                        glNormal3f(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ); 563                        glVertex3f(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ); 564                } 565        } 566} 567 568//vMarchTetrahedron performs the Marching Tetrahedrons algorithm on a single tetrahedron 569GLvoid vMarchTetrahedron(GLvector *pasTetrahedronPosition, GLfloat *pafTetrahedronValue) 570{ 571        extern GLint aiTetrahedronEdgeFlags[16]; 572        extern GLint a2iTetrahedronTriangles[16][7]; 573 574        GLint iEdge, iVert0, iVert1, iEdgeFlags, iTriangle, iCorner, iVertex, iFlagIndex = 0; 575        GLfloat fOffset, fInvOffset, fValue = 0.0; 576        GLvector asEdgeVertex[6]; 577        GLvector asEdgeNorm[6]; 578        GLvector sColor; 579 580        //Find which vertices are inside of the surface and which are outside 581        for(iVertex = 0; iVertex < 4; iVertex++) 582        { 583                if(pafTetrahedronValue[iVertex] <= fTargetValue)  584                        iFlagIndex |= 1<<iVertex; 585        } 586 587        //Find which edges are intersected by the surface 588        iEdgeFlags = aiTetrahedronEdgeFlags[iFlagIndex]; 589 590        //If the tetrahedron is entirely inside or outside of the surface, then there will be no intersections 591        if(iEdgeFlags == 0) 592        { 593                return; 594        } 595        //Find the point of intersection of the surface with each edge 596        // Then find the normal to the surface at those points 597        for(iEdge = 0; iEdge < 6; iEdge++) 598        { 599                //if there is an intersection on this edge 600                if(iEdgeFlags & (1<<iEdge)) 601                { 602                        iVert0 = a2iTetrahedronEdgeConnection[iEdge][0]; 603                        iVert1 = a2iTetrahedronEdgeConnection[iEdge][1]; 604                        fOffset = fGetOffset(pafTetrahedronValue[iVert0], pafTetrahedronValue[iVert1], fTargetValue); 605                        fInvOffset = 1.0 - fOffset; 606 607                        asEdgeVertex[iEdge].fX = fInvOffset*pasTetrahedronPosition[iVert0].fX  +  fOffset*pasTetrahedronPosition[iVert1].fX; 608                        asEdgeVertex[iEdge].fY = fInvOffset*pasTetrahedronPosition[iVert0].fY  +  fOffset*pasTetrahedronPosition[iVert1].fY; 609                        asEdgeVertex[iEdge].fZ = fInvOffset*pasTetrahedronPosition[iVert0].fZ  +  fOffset*pasTetrahedronPosition[iVert1].fZ; 610                         611                        vGetNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ); 612                } 613        } 614        //Draw the triangles that were found.  There can be up to 2 per tetrahedron 615        for(iTriangle = 0; iTriangle < 2; iTriangle++) 616        { 617                if(a2iTetrahedronTriangles[iFlagIndex][3*iTriangle] < 0) 618                        break; 619 620                for(iCorner = 0; iCorner < 3; iCorner++) 621                { 622                        iVertex = a2iTetrahedronTriangles[iFlagIndex][3*iTriangle+iCorner]; 623 624                        vGetColor(sColor, asEdgeVertex[iVertex], asEdgeNorm[iVertex]); 625                        glColor3f(sColor.fX, sColor.fY, sColor.fZ); 626                        glNormal3f(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ); 627                        glVertex3f(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ); 628                } 629        } 630} 631 632 633 634//vMarchCube2 performs the Marching Tetrahedrons algorithm on a single cube by making six calls to vMarchTetrahedron 635GLvoid vMarchCube2(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale) 636{ 637        GLint iVertex, iTetrahedron, iVertexInACube; 638        GLvector asCubePosition[8]; 639        GLfloat  afCubeValue[8]; 640        GLvector asTetrahedronPosition[4]; 641        GLfloat  afTetrahedronValue[4]; 642 643        //Make a local copy of the cube's corner positions 644        for(iVertex = 0; iVertex < 8; iVertex++) 645        { 646                asCubePosition[iVertex].fX = fX + a2fVertexOffset[iVertex][0]*fScale; 647                asCubePosition[iVertex].fY = fY + a2fVertexOffset[iVertex][1]*fScale; 648                asCubePosition[iVertex].fZ = fZ + a2fVertexOffset[iVertex][2]*fScale; 649        } 650 651        //Make a local copy of the cube's corner values 652        for(iVertex = 0; iVertex < 8; iVertex++) 653        { 654                afCubeValue[iVertex] = fSample(asCubePosition[iVertex].fX, 655                                                   asCubePosition[iVertex].fY, 656                                               asCubePosition[iVertex].fZ); 657        } 658 659        for(iTetrahedron = 0; iTetrahedron < 6; iTetrahedron++) 660        { 661                for(iVertex = 0; iVertex < 4; iVertex++) 662                { 663                        iVertexInACube = a2iTetrahedronsInACube[iTetrahedron][iVertex]; 664                        asTetrahedronPosition[iVertex].fX = asCubePosition[iVertexInACube].fX; 665                        asTetrahedronPosition[iVertex].fY = asCubePosition[iVertexInACube].fY; 666                        asTetrahedronPosition[iVertex].fZ = asCubePosition[iVertexInACube].fZ; 667                        afTetrahedronValue[iVertex] = afCubeValue[iVertexInACube]; 668                } 669                vMarchTetrahedron(asTetrahedronPosition, afTetrahedronValue); 670        } 671} 672         673 674//vMarchingCubes iterates over the entire dataset, calling vMarchCube on each cube 675GLvoid vMarchingCubes() 676{ 677        GLint iX, iY, iZ; 678        for(iX = 0; iX < iDataSetSize; iX++) 679        for(iY = 0; iY < iDataSetSize; iY++) 680        for(iZ = 0; iZ < iDataSetSize; iZ++) 681        { 682                vMarchCube(iX*fStepSize, iY*fStepSize, iZ*fStepSize, fStepSize); 683        } 684} 685 686 687// For any edge, if one vertex is inside of the surface and the other is outside of the surface 688//  then the edge intersects the surface 689// For each of the 4 vertices of the tetrahedron can be two possible states : either inside or outside of the surface 690// For any tetrahedron the are 2^4=16 possible sets of vertex states 691// This table lists the edges intersected by the surface for all 16 possible vertex states 692// There are 6 edges.  For each entry in the table, if edge #n is intersected, then bit #n is set to 1 693 694GLint aiTetrahedronEdgeFlags[16]= 695{ 696        0x00, 0x0d, 0x13, 0x1e, 0x26, 0x2b, 0x35, 0x38, 0x38, 0x35, 0x2b, 0x26, 0x1e, 0x13, 0x0d, 0x00,  697}; 698 699 700// For each of the possible vertex states listed in aiTetrahedronEdgeFlags there is a specific triangulation 701// of the edge intersection points.  a2iTetrahedronTriangles lists all of them in the form of 702// 0-2 edge triples with the list terminated by the invalid value -1. 703// 704// I generated this table by hand 705 706GLint a2iTetrahedronTriangles[16][7] = 707{ 708        {-1, -1, -1, -1, -1, -1, -1}, 709        { 0,  3,  2, -1, -1, -1, -1}, 710        { 0,  1,  4, -1, -1, -1, -1}, 711        { 1,  4,  2,  2,  4,  3, -1}, 712 713        { 1,  2,  5, -1, -1, -1, -1}, 714        { 0,  3,  5,  0,  5,  1, -1}, 715        { 0,  2,  5,  0,  5,  4, -1}, 716        { 5,  4,  3, -1, -1, -1, -1}, 717 718        { 3,  4,  5, -1, -1, -1, -1}, 719        { 4,  5,  0,  5,  2,  0, -1}, 720        { 1,  5,  0,  5,  3,  0, -1}, 721        { 5,  2,  1, -1, -1, -1, -1}, 722 723        { 3,  4,  2,  2,  4,  1, -1}, 724        { 4,  1,  0, -1, -1, -1, -1}, 725        { 2,  3,  0, -1, -1, -1, -1}, 726        {-1, -1, -1, -1, -1, -1, -1}, 727}; 728 729// For any edge, if one vertex is inside of the surface and the other is outside of the surface 730//  then the edge intersects the surface 731// For each of the 8 vertices of the cube can be two possible states : either inside or outside of the surface 732// For any cube the are 2^8=256 possible sets of vertex states 733// This table lists the edges intersected by the surface for all 256 possible vertex states 734// There are 12 edges.  For each entry in the table, if edge #n is intersected, then bit #n is set to 1 735 736GLint aiCubeEdgeFlags[256]= 737{ 738        0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,  739        0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,  740        0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,  741        0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,  742        0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,  743        0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,  744        0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,  745        0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,  746        0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,  747        0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,  748        0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,  749        0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,  750        0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,  751        0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,  752        0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,  753        0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000 754}; 755 756//  For each of the possible vertex states listed in aiCubeEdgeFlags there is a specific triangulation 757//  of the edge intersection points.  a2iTriangleConnectionTable lists all of them in the form of 758//  0-5 edge triples with the list terminated by the invalid value -1. 759//  For example: a2iTriangleConnectionTable[3] list the 2 triangles formed when corner[0]  760//  and corner[1] are inside of the surface, but the rest of the cube is not. 761// 762//  I found this table in an example program someone wrote long ago.  It was probably generated by hand 763 764GLint a2iTriangleConnectionTable[256][16] =   765{ 766        {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 767        {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 768        {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 769        {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 770        {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 771        {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 772        {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 773        {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, 774        {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 775        {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 776        {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 777        {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, 778        {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 779        {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, 780        {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, 781        {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 782        {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 783        {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 784        {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 785        {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, 786        {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 787        {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, 788        {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, 789        {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, 790        {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 791        {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, 792        {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, 793        {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, 794        {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, 795        {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, 796        {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, 797        {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, 798        {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 799        {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 800        {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 801        {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, 802        {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 803        {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, 804        {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, 805        {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, 806        {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 807        {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, 808        {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, 809        {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, 810        {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, 811        {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, 812        {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, 813        {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, 814        {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 815        {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, 816        {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, 817        {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 818        {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, 819        {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, 820        {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, 821        {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, 822        {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, 823        {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, 824        {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, 825        {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, 826        {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, 827        {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, 828        {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, 829        {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 830        {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 831        {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 832        {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 833        {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, 834        {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 835        {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, 836        {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, 837        {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, 838        {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 839        {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, 840        {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, 841        {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, 842        {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, 843        {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, 844        {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, 845        {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, 846        {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 847        {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, 848        {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, 849        {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, 850        {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, 851        {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, 852        {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, 853        {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, 854        {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, 855        {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, 856        {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, 857        {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, 858        {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, 859        {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, 860        {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, 861        {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, 862        {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 863        {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, 864        {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, 865        {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, 866        {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, 867        {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, 868        {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 869        {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, 870        {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, 871        {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, 872        {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, 873        {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, 874        {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, 875        {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, 876        {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, 877        {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 878        {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, 879        {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, 880        {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, 881        {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, 882        {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, 883        {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, 884        {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, 885        {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 886        {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, 887        {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, 888        {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, 889        {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, 890        {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, 891        {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 892        {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, 893        {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 894        {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 895        {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 896        {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 897        {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, 898        {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 899        {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, 900        {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, 901        {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, 902        {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 903        {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, 904        {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, 905        {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, 906        {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, 907        {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, 908        {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, 909        {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, 910        {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 911        {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, 912        {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, 913        {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, 914        {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, 915        {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, 916        {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, 917        {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, 918        {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, 919        {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 920        {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, 921        {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, 922        {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, 923        {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, 924        {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, 925        {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 926        {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 927        {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, 928        {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, 929        {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, 930        {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, 931        {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, 932        {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, 933        {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, 934        {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, 935        {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, 936        {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, 937        {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, 938        {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, 939        {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, 940        {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, 941        {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, 942        {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, 943        {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, 944        {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, 945        {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, 946        {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, 947        {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, 948        {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, 949        {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, 950        {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, 951        {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, 952        {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, 953        {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 954        {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, 955        {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, 956        {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 957        {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 958        {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 959        {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, 960        {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, 961        {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, 962        {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, 963        {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, 964        {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, 965        {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, 966        {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, 967        {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, 968        {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, 969        {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, 970        {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 971        {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, 972        {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, 973        {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 974        {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, 975        {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, 976        {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, 977        {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, 978        {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, 979        {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, 980        {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, 981        {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 982        {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, 983        {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, 984        {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, 985        {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, 986        {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, 987        {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 988        {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, 989        {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 990        {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, 991        {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, 992        {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, 993        {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, 994        {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, 995        {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, 996        {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, 997        {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, 998        {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, 999        {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},1000        {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},1001        {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1002        {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},1003        {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},1004        {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1005        {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1006        {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1007        {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},1008        {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},1009        {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1010        {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},1011        {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},1012        {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1013        {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1014        {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},1015        {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1016        {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},1017        {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1018        {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1019        {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1020        {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},1021        {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}1022};`

A self contained implement of the Marching Cubes and Marching Tetrahedrons algorithms for creating surfaces from volume data.

It is fairly easy to integrate this code into your own application as I did once for a volume data visualization class.