// // Marching Cubes Example Program // by Cory Bloyd (corysama@yahoo.com) // // A simple, portable and complete implementation of the Marching Cubes // and Marching Tetrahedrons algorithms in a single source file. // There are many ways that this code could be made faster, but the // intent is for the code to be easy to understand. // // For a description of the algorithm go to // http://astronomy.swin.edu.au/pbourke/modelling/polygonise/ // // This code is public domain. // #include "stdio.h" #include "math.h" //This program requires the OpenGL and GLUT libraries // You can obtain them for free from http://www.opengl.org #include "GL/glut.h" struct GLvector { GLfloat fX; GLfloat fY; GLfloat fZ; }; //These tables are used so that everything can be done in little loops that you can look at all at once // rather than in pages and pages of unrolled code. //a2fVertexOffset lists the positions, relative to vertex0, of each of the 8 vertices of a cube static const GLfloat a2fVertexOffset[8][3] = { {0.0, 0.0, 0.0},{1.0, 0.0, 0.0},{1.0, 1.0, 0.0},{0.0, 1.0, 0.0}, {0.0, 0.0, 1.0},{1.0, 0.0, 1.0},{1.0, 1.0, 1.0},{0.0, 1.0, 1.0} }; //a2iEdgeConnection lists the index of the endpoint vertices for each of the 12 edges of the cube static const GLint a2iEdgeConnection[12][2] = { {0,1}, {1,2}, {2,3}, {3,0}, {4,5}, {5,6}, {6,7}, {7,4}, {0,4}, {1,5}, {2,6}, {3,7} }; //a2fEdgeDirection lists the direction vector (vertex1-vertex0) for each edge in the cube static const GLfloat a2fEdgeDirection[12][3] = { {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0}, {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0}, {0.0, 0.0, 1.0},{0.0, 0.0, 1.0},{ 0.0, 0.0, 1.0},{0.0, 0.0, 1.0} }; //a2iTetrahedronEdgeConnection lists the index of the endpoint vertices for each of the 6 edges of the tetrahedron static const GLint a2iTetrahedronEdgeConnection[6][2] = { {0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3} }; //a2iTetrahedronEdgeConnection lists the index of verticies from a cube // that made up each of the six tetrahedrons within the cube static const GLint a2iTetrahedronsInACube[6][4] = { {0,5,1,6}, {0,1,2,6}, {0,2,3,6}, {0,3,7,6}, {0,7,4,6}, {0,4,5,6}, }; static const GLfloat afAmbientWhite [] = {0.25, 0.25, 0.25, 1.00}; static const GLfloat afAmbientRed [] = {0.25, 0.00, 0.00, 1.00}; static const GLfloat afAmbientGreen [] = {0.00, 0.25, 0.00, 1.00}; static const GLfloat afAmbientBlue [] = {0.00, 0.00, 0.25, 1.00}; static const GLfloat afDiffuseWhite [] = {0.75, 0.75, 0.75, 1.00}; static const GLfloat afDiffuseRed [] = {0.75, 0.00, 0.00, 1.00}; static const GLfloat afDiffuseGreen [] = {0.00, 0.75, 0.00, 1.00}; static const GLfloat afDiffuseBlue [] = {0.00, 0.00, 0.75, 1.00}; static const GLfloat afSpecularWhite[] = {1.00, 1.00, 1.00, 1.00}; static const GLfloat afSpecularRed [] = {1.00, 0.25, 0.25, 1.00}; static const GLfloat afSpecularGreen[] = {0.25, 1.00, 0.25, 1.00}; static const GLfloat afSpecularBlue [] = {0.25, 0.25, 1.00, 1.00}; GLenum ePolygonMode = GL_FILL; GLint iDataSetSize = 16; GLfloat fStepSize = 1.0/iDataSetSize; GLfloat fTargetValue = 48.0; GLfloat fTime = 0.0; GLvector sSourcePoint[3]; GLboolean bSpin = true; GLboolean bMove = true; GLboolean bLight = true; void vIdle(); void vDrawScene(); void vResize(GLsizei, GLsizei); void vKeyboard(unsigned char cKey, int iX, int iY); void vSpecial(int iKey, int iX, int iY); GLvoid vPrintHelp(); GLvoid vSetTime(GLfloat fTime); GLfloat fSample1(GLfloat fX, GLfloat fY, GLfloat fZ); GLfloat fSample2(GLfloat fX, GLfloat fY, GLfloat fZ); GLfloat fSample3(GLfloat fX, GLfloat fY, GLfloat fZ); GLfloat (*fSample)(GLfloat fX, GLfloat fY, GLfloat fZ) = fSample1; GLvoid vMarchingCubes(); GLvoid vMarchCube1(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale); GLvoid vMarchCube2(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale); GLvoid (*vMarchCube)(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale) = vMarchCube1; void main(int argc, char **argv) { GLfloat afPropertiesAmbient [] = {0.50, 0.50, 0.50, 1.00}; GLfloat afPropertiesDiffuse [] = {0.75, 0.75, 0.75, 1.00}; GLfloat afPropertiesSpecular[] = {1.00, 1.00, 1.00, 1.00}; GLsizei iWidth = 640.0; GLsizei iHeight = 480.0; glutInit(&argc, argv); glutInitWindowPosition( 0, 0); glutInitWindowSize(iWidth, iHeight); glutInitDisplayMode( GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE ); glutCreateWindow( "Marching Cubes" ); glutDisplayFunc( vDrawScene ); glutIdleFunc( vIdle ); glutReshapeFunc( vResize ); glutKeyboardFunc( vKeyboard ); glutSpecialFunc( vSpecial ); glClearColor( 0.0, 0.0, 0.0, 1.0 ); glClearDepth( 1.0 ); glEnable(GL_DEPTH_TEST); glEnable(GL_LIGHTING); glPolygonMode(GL_FRONT_AND_BACK, ePolygonMode); glLightfv( GL_LIGHT0, GL_AMBIENT, afPropertiesAmbient); glLightfv( GL_LIGHT0, GL_DIFFUSE, afPropertiesDiffuse); glLightfv( GL_LIGHT0, GL_SPECULAR, afPropertiesSpecular); glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, 1.0); glEnable( GL_LIGHT0 ); glMaterialfv(GL_BACK, GL_AMBIENT, afAmbientGreen); glMaterialfv(GL_BACK, GL_DIFFUSE, afDiffuseGreen); glMaterialfv(GL_FRONT, GL_AMBIENT, afAmbientBlue); glMaterialfv(GL_FRONT, GL_DIFFUSE, afDiffuseBlue); glMaterialfv(GL_FRONT, GL_SPECULAR, afSpecularWhite); glMaterialf( GL_FRONT, GL_SHININESS, 25.0); vResize(iWidth, iHeight); vPrintHelp(); glutMainLoop(); } GLvoid vPrintHelp() { printf("Marching Cubes Example by Cory Bloyd (dejaspaminacan@my-deja.com)\n\n"); printf("+/- increase/decrease sample density\n"); printf("PageUp/PageDown increase/decrease surface value\n"); printf("s change sample function\n"); printf("c toggle marching cubes / marching tetrahedrons\n"); printf("w wireframe on/off\n"); printf("l toggle lighting / color-by-normal\n"); printf("Home spin scene on/off\n"); printf("End source point animation on/off\n"); } void vResize( GLsizei iWidth, GLsizei iHeight ) { GLfloat fAspect, fHalfWorldSize = (1.4142135623730950488016887242097/2); glViewport( 0, 0, iWidth, iHeight ); glMatrixMode (GL_PROJECTION); glLoadIdentity (); if(iWidth <= iHeight) { fAspect = (GLfloat)iHeight / (GLfloat)iWidth; glOrtho(-fHalfWorldSize, fHalfWorldSize, -fHalfWorldSize*fAspect, fHalfWorldSize*fAspect, -10*fHalfWorldSize, 10*fHalfWorldSize); } else { fAspect = (GLfloat)iWidth / (GLfloat)iHeight; glOrtho(-fHalfWorldSize*fAspect, fHalfWorldSize*fAspect, -fHalfWorldSize, fHalfWorldSize, -10*fHalfWorldSize, 10*fHalfWorldSize); } glMatrixMode( GL_MODELVIEW ); } void vKeyboard(unsigned char cKey, int iX, int iY) { switch(cKey) { case 'w' : { if(ePolygonMode == GL_LINE) { ePolygonMode = GL_FILL; } else { ePolygonMode = GL_LINE; } glPolygonMode(GL_FRONT_AND_BACK, ePolygonMode); } break; case '+' : case '=' : { ++iDataSetSize; fStepSize = 1.0/iDataSetSize; } break; case '-' : { if(iDataSetSize > 1) { --iDataSetSize; fStepSize = 1.0/iDataSetSize; } } break; case 'c' : { if(vMarchCube == vMarchCube1) { vMarchCube = vMarchCube2;//Use Marching Tetrahedrons } else { vMarchCube = vMarchCube1;//Use Marching Cubes } } break; case 's' : { if(fSample == fSample1) { fSample = fSample2; } else if(fSample == fSample2) { fSample = fSample3; } else { fSample = fSample1; } } break; case 'l' : { if(bLight) { glDisable(GL_LIGHTING);//use vertex colors } else { glEnable(GL_LIGHTING);//use lit material color } bLight = !bLight; }; } } void vSpecial(int iKey, int iX, int iY) { switch(iKey) { case GLUT_KEY_PAGE_UP : { if(fTargetValue < 1000.0) { fTargetValue *= 1.1; } } break; case GLUT_KEY_PAGE_DOWN : { if(fTargetValue > 1.0) { fTargetValue /= 1.1; } } break; case GLUT_KEY_HOME : { bSpin = !bSpin; } break; case GLUT_KEY_END : { bMove = !bMove; } break; } } void vIdle() { glutPostRedisplay(); } void vDrawScene() { static GLfloat fPitch = 0.0; static GLfloat fYaw = 0.0; static GLfloat fTime = 0.0; glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT ); glPushMatrix(); if(bSpin) { fPitch += 4.0; fYaw += 2.5; } if(bMove) { fTime += 0.025; } vSetTime(fTime); glTranslatef(0.0, 0.0, -1.0); glRotatef( -fPitch, 1.0, 0.0, 0.0); glRotatef( 0.0, 0.0, 1.0, 0.0); glRotatef( fYaw, 0.0, 0.0, 1.0); glPushAttrib(GL_LIGHTING_BIT); glDisable(GL_LIGHTING); glColor3f(1.0, 1.0, 1.0); glutWireCube(1.0); glPopAttrib(); glPushMatrix(); glTranslatef(-0.5, -0.5, -0.5); glBegin(GL_TRIANGLES); vMarchingCubes(); glEnd(); glPopMatrix(); glPopMatrix(); glutSwapBuffers(); } //fGetOffset finds the approximate point of intersection of the surface // between two points with the values fValue1 and fValue2 GLfloat fGetOffset(GLfloat fValue1, GLfloat fValue2, GLfloat fValueDesired) { GLdouble fDelta = fValue2 - fValue1; if(fDelta == 0.0) { return 0.5; } return (fValueDesired - fValue1)/fDelta; } //vGetColor generates a color from a given position and normal of a point GLvoid vGetColor(GLvector &rfColor, GLvector &rfPosition, GLvector &rfNormal) { GLfloat fX = rfNormal.fX; GLfloat fY = rfNormal.fY; GLfloat fZ = rfNormal.fZ; rfColor.fX = (fX > 0.0 ? fX : 0.0) + (fY < 0.0 ? -0.5*fY : 0.0) + (fZ < 0.0 ? -0.5*fZ : 0.0); rfColor.fY = (fY > 0.0 ? fY : 0.0) + (fZ < 0.0 ? -0.5*fZ : 0.0) + (fX < 0.0 ? -0.5*fX : 0.0); rfColor.fZ = (fZ > 0.0 ? fZ : 0.0) + (fX < 0.0 ? -0.5*fX : 0.0) + (fY < 0.0 ? -0.5*fY : 0.0); } GLvoid vNormalizeVector(GLvector &rfVectorResult, GLvector &rfVectorSource) { GLfloat fOldLength; GLfloat fScale; fOldLength = sqrtf( (rfVectorSource.fX * rfVectorSource.fX) + (rfVectorSource.fY * rfVectorSource.fY) + (rfVectorSource.fZ * rfVectorSource.fZ) ); if(fOldLength == 0.0) { rfVectorResult.fX = rfVectorSource.fX; rfVectorResult.fY = rfVectorSource.fY; rfVectorResult.fZ = rfVectorSource.fZ; } else { fScale = 1.0/fOldLength; rfVectorResult.fX = rfVectorSource.fX*fScale; rfVectorResult.fY = rfVectorSource.fY*fScale; rfVectorResult.fZ = rfVectorSource.fZ*fScale; } } //Generate a sample data set. fSample1(), fSample2() and fSample3() define three scalar fields whose // values vary by the X,Y and Z coordinates and by the fTime value set by vSetTime() GLvoid vSetTime(GLfloat fNewTime) { GLfloat fOffset; GLint iSourceNum; for(iSourceNum = 0; iSourceNum < 3; iSourceNum++) { sSourcePoint[iSourceNum].fX = 0.5; sSourcePoint[iSourceNum].fY = 0.5; sSourcePoint[iSourceNum].fZ = 0.5; } fTime = fNewTime; fOffset = 1.0 + sinf(fTime); sSourcePoint[0].fX *= fOffset; sSourcePoint[1].fY *= fOffset; sSourcePoint[2].fZ *= fOffset; } //fSample1 finds the distance of (fX, fY, fZ) from three moving points GLfloat fSample1(GLfloat fX, GLfloat fY, GLfloat fZ) { GLdouble fResult = 0.0; GLdouble fDx, fDy, fDz; fDx = fX - sSourcePoint[0].fX; fDy = fY - sSourcePoint[0].fY; fDz = fZ - sSourcePoint[0].fZ; fResult += 0.5/(fDx*fDx + fDy*fDy + fDz*fDz); fDx = fX - sSourcePoint[1].fX; fDy = fY - sSourcePoint[1].fY; fDz = fZ - sSourcePoint[1].fZ; fResult += 1.0/(fDx*fDx + fDy*fDy + fDz*fDz); fDx = fX - sSourcePoint[2].fX; fDy = fY - sSourcePoint[2].fY; fDz = fZ - sSourcePoint[2].fZ; fResult += 1.5/(fDx*fDx + fDy*fDy + fDz*fDz); return fResult; } //fSample2 finds the distance of (fX, fY, fZ) from three moving lines GLfloat fSample2(GLfloat fX, GLfloat fY, GLfloat fZ) { GLdouble fResult = 0.0; GLdouble fDx, fDy, fDz; fDx = fX - sSourcePoint[0].fX; fDy = fY - sSourcePoint[0].fY; fResult += 0.5/(fDx*fDx + fDy*fDy); fDx = fX - sSourcePoint[1].fX; fDz = fZ - sSourcePoint[1].fZ; fResult += 0.75/(fDx*fDx + fDz*fDz); fDy = fY - sSourcePoint[2].fY; fDz = fZ - sSourcePoint[2].fZ; fResult += 1.0/(fDy*fDy + fDz*fDz); return fResult; } //fSample2 defines a height field by plugging the distance from the center into the sin and cos functions GLfloat fSample3(GLfloat fX, GLfloat fY, GLfloat fZ) { GLfloat fHeight = 20.0*(fTime + sqrt((0.5-fX)*(0.5-fX) + (0.5-fY)*(0.5-fY))); fHeight = 1.5 + 0.1*(sinf(fHeight) + cosf(fHeight)); GLdouble fResult = (fHeight - fZ)*50.0; return fResult; } //vGetNormal() finds the gradient of the scalar field at a point //This gradient can be used as a very accurate vertx normal for lighting calculations GLvoid vGetNormal(GLvector &rfNormal, GLfloat fX, GLfloat fY, GLfloat fZ) { rfNormal.fX = fSample(fX-0.01, fY, fZ) - fSample(fX+0.01, fY, fZ); rfNormal.fY = fSample(fX, fY-0.01, fZ) - fSample(fX, fY+0.01, fZ); rfNormal.fZ = fSample(fX, fY, fZ-0.01) - fSample(fX, fY, fZ+0.01); vNormalizeVector(rfNormal, rfNormal); } //vMarchCube1 performs the Marching Cubes algorithm on a single cube GLvoid vMarchCube1(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale) { extern GLint aiCubeEdgeFlags[256]; extern GLint a2iTriangleConnectionTable[256][16]; GLint iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags; GLfloat fOffset; GLvector sColor; GLfloat afCubeValue[8]; GLvector asEdgeVertex[12]; GLvector asEdgeNorm[12]; //Make a local copy of the values at the cube's corners for(iVertex = 0; iVertex < 8; iVertex++) { afCubeValue[iVertex] = fSample(fX + a2fVertexOffset[iVertex][0]*fScale, fY + a2fVertexOffset[iVertex][1]*fScale, fZ + a2fVertexOffset[iVertex][2]*fScale); } //Find which vertices are inside of the surface and which are outside iFlagIndex = 0; for(iVertexTest = 0; iVertexTest < 8; iVertexTest++) { if(afCubeValue[iVertexTest] <= fTargetValue) iFlagIndex |= 1<