Marching Cube 医学图像三维可视化VC++实例

#include <stdlib.h> #include <string.h> #include "stdio.h" #include "math.h" #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 立方体8个顶点 // 立方体8个顶点 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 立方体12个边 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 立方体 12个边的方向向量 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} }; 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; GLfloat fTargetValue = 100.0; GLfloat fTime = 0.0; GLvector sSourcePoint[3]; 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); void spinDisplay(); void fYawDisplay(void); void myMouseCall(int button, int state, int x ,int y); GLvoid vPrintHelp(); GLvoid vMarchingCubes(); GLvoid vMarchCube1(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale); GLvoid (*vMarchCube)(GLfloat fX, GLfloat fY, GLfloat fZ, GLfloat fScale) = vMarchCube1; // #define NX 200 #define NY 160 #define NZ 160 // //#define NX 48 //#define NY 62 //#define NZ 42 short int ***data; GLfloat HoriRotate = 90; GLfloat VertiRotate = 0; GLfloat Xtranform=-90.0; GLfloat Ytranform= 80.0; GLfloat Ztranform=-10.0; GLfloat ZoonAspect=140.0; GLfloat ZoomScale=1.0; short int themax = 0,themin = 255; void main(int argc, char **argv) { /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// int i,j,k,c; short int isolevel = 128; FILE *fptr; argc=2; argv[1]="C:\\Users\\Conan\\Desktop\\MarchingCube_Good\\mri.raw"; // argv[1]="C:\\Users\\Conan\\Desktop\\MarchingCube_Good\\HeadMRVolume.raw"; for (i=1;i<argc;i++) { if (strcmp(argv[i],"-i") == 0) isolevel = atof(argv[i+1]); } // Malloc the volumetric data, hardwired size! data = (short int ***)malloc(NX*sizeof(short int **)); for (i=0;i<NX;i++) data[i] = (short int **)malloc(NY*sizeof(short int *)); for (i=0;i<NX;i++) for (j=0;j<NY;j++) data[i][j] = (short int *)malloc(NZ*sizeof(short int)); // Open and read the raw data fprintf(stderr,"Reading data ...\n"); if ((fptr = fopen(argv[argc-1],"rb")) == NULL) { fprintf(stderr,"File open failed\n"); exit(-1); } for (k=0;k<NZ;k++) { for (j=0;j<NY;j++) { for (i=0;i<NX;i++) { if ((c = fgetc(fptr)) == EOF) { fprintf(stderr,"Unexpected end of file\n"); exit(-1); } data[i][j][k] = c; if (c > themax) themax = c; if (c < themin) themin = c; } } } fclose(fptr); fprintf(stderr,"Volumetric data range: %d -> %d\n",themin,themax); /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// GLfloat afPropertiesAmbient [] = {0.25, 0.25, 0.25, 1.00}; GLfloat afPropertiesDiffuse [] = {0.75, 0.75, 0.75, 1.00}; GLfloat afPropertiesSpecular[] = {1.00, 0.25, 0.25, 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 ); glutMouseFunc(myMouseCall); 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, afAmbientRed); glMaterialfv(GL_BACK, GL_DIFFUSE, afDiffuseRed); glMaterialfv(GL_FRONT, GL_AMBIENT, afAmbientWhite); glMaterialfv(GL_FRONT, GL_DIFFUSE, afDiffuseWhite); glMaterialfv(GL_FRONT, GL_SPECULAR, afSpecularWhite); glMaterialf( GL_FRONT, GL_SHININESS, 10.0); vResize(iWidth, iHeight); vPrintHelp(); glutMainLoop(); } // 帮助内容 GLvoid vPrintHelp() { printf("Marching Cubes Rountines by Niu Kenan (niukenan@sina.com)\n\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*ZoonAspect, fHalfWorldSize*ZoonAspect, -fHalfWorldSize*fAspect*ZoonAspect, fHalfWorldSize*fAspect*ZoonAspect, -10*fHalfWorldSize*ZoonAspect, 10*fHalfWorldSize*ZoonAspect); } else { fAspect = (GLfloat)iWidth / (GLfloat)iHeight; glOrtho(-fHalfWorldSize*fAspect*ZoonAspect, fHalfWorldSize*fAspect*ZoonAspect, -fHalfWorldSize*ZoonAspect, fHalfWorldSize*ZoonAspect, -10*fHalfWorldSize*ZoonAspect, 10*fHalfWorldSize*ZoonAspect); } glMatrixMode( GL_MODELVIEW ); } void spinDisplay(void) { HoriRotate = HoriRotate+10.0; if (HoriRotate >360.0) HoriRotate = HoriRotate -360.0; glutPostRedisplay(); } void fYawDisplay(void) { VertiRotate = VertiRotate+10.0; if (VertiRotate >360.0) VertiRotate = VertiRotate -360.0; glutPostRedisplay(); } // 鼠标响应 void myMouseCall(int button, int state, int x ,int y) { switch (button) { case GLUT_LEFT_BUTTON: if (state == GLUT_DOWN) spinDisplay(); break; case GLUT_RIGHT_BUTTON: