面绘制经典算法:MarchingCube实现(控制台篇)
来源:互联网 发布:核弹 知乎 编辑:程序博客网 时间:2024/05/22 06:25
1.MarchingCube
marchingcube是一个比较经典而古老的算法,也是面绘制中应用比较多的算法,Marchingcube发展到今天也遇到了几何拓扑、一致性的问题仍待改善。本文研究的就是经典的marchingcube(Paul,1997).
关于MarchingCube原理(含代码),重点推荐如下:
1.《医学图像三维重建和可视化:VC++实现实例》(张俊华)
2.Bourke P. Polygonising a scalar field[J]. 1997.
3.论文+源代码+交互界面
2.控制台上实现与讨论
算法设计步骤:
1.二维断层图像顺序地读入内存并构成三位数据场;
2.依次扫描相邻两层数据,逐个构造立方体;
3.将立方体每个顶点灰度和给定等值面阈值进行比较,计算索引值;
4.对于含有等值面的立方体,通过灰度差分计算立方体各顶点的梯度;
5.根据索引值查找边索引表,获得和等值面有交点的当前立方体的相交边;
6.根据相交边的两个顶点坐标及其法向量,利用差值(中值)计算等值点的坐标与法向量;
7.根据索引值查找三角片索引表,确定当前立方体内构成三角片的等值点的组合方式;
8.由各个立方体内的三角面片构成等值面。
#include <stdio.h>#include <stdlib.h>#include <math.h>#include <iostream>using namespace std;static int edgeTable[256]={ }; //根据索引值确定立方体和等值面交点情况static int triTable[256][16] ={ };//根据索引值确定三角面片的个数以及位置typedef struct {double x,y,z;} XYZ;typedef struct{XYZ p[8]; //顶点坐标XYZ n[8]; //法向量double val[8]; //顶点对应的灰度值} GRIDCELL; //立方体typedef struct {XYZ p[3]; //顶点坐标XYZ c; //几何中心XYZ n[3]; //法向量} TRIANGLE;#define ABS(x) (x < 0 ? -(x) : (x))//声明int PolygoniseCube(GRIDCELL , double , TRIANGLE * );XYZ VertexInterp(double , XYZ , XYZ , double, double );//3D数据场大小#define NX 200#define NY 160#define NZ 160int main( int argc, char **argv){short int*** data = NULL;short int isolevel = 128, themax = 0, themin = 255; //等值面 数据场的范围GRIDCELL grid;TRIANGLE triangles[10];TRIANGLE* tri = NULL; //保存的是所有三角形面片的顶点坐标int ntri = 0;FILE* fptr = NULL;//打开并阅读原始文件fprintf(stderr, "Reading data...\n");errno_t err;err = fopen_s(&fptr,"mri.raw", "rb");if( err != NULL ) {fprintf(stderr, "File open failed\n");exit(-1);}//为三维数据场分配空间data = (short int***) malloc( NX * sizeof(short int **) );for (int i=0;i<NX;i++)data[i] = (short int**) malloc( NY * sizeof(short int *));for (int i=0;i<NX;i++)for (int j=0;j<NY;j++)data[i][j] = (short int*) malloc(NZ * sizeof(short int));int str; //从文件中读取字符for (int k=0; k<NZ; k++){for(int j=0; j<NY; j++){for(int i=0; i<NX; i++){if ((str = fgetc(fptr)) == EOF) { //fgetc顺序读取数据,直到遇见-文件结束符EOFfprintf(stderr,"Unexpected end of file\n");exit(-1);}data[i][j][k] = str;if( str > themax )themax = str;if( str < themin )themin = str;}}}fclose(fptr);fprintf(stderr,"Volumetric data range: %d -> %d\n",themin,themax);//3D数据场处理fprintf(stderr, "polygonising data...\n");for (int i=0; i<NX-1; i++){if(i % (NX/10) == 0)fprintf(stderr, " Slice %d of %d \n", i, NX);for (int j=0; j<NY-1; j++){for (int k=0; k<NZ-1; k++){grid.p[0].x = i; grid.p[0].y = j; grid.p[0].z = k; grid.val[0] = data[i] [j] [k];grid.p[1].x = i+1;grid.p[1].y = j; grid.p[1].z = k; grid.val[1] = data[i+1][j] [k];grid.p[2].x = i+1;grid.p[2].y = j+1;grid.p[2].z = k; grid.val[2] = data[i+1][j+1][k];grid.p[3].x = i; grid.p[3].y = j+1;grid.p[3].z = k; grid.val[3] = data[i] [j+1][k];grid.p[4].x = i; grid.p[4].y = j; grid.p[4].z = k+1;grid.val[4] = data[i] [j] [k+1];grid.p[5].x = i+1;grid.p[5].y = j; grid.p[5].z = k+1;grid.val[5] = data[i+1][j] [k+1];grid.p[6].x = i+1;grid.p[6].y = j+1;grid.p[6].z = k+1;grid.val[6] = data[i+1][j+1][k+1];grid.p[7].x = i; grid.p[7].y = j+1;grid.p[7].z = k+1;grid.val[7] = data[i] [j+1][k+1];//计算当前立方体中的三角面片数目int numOftriangle = PolygoniseCube(grid, isolevel, triangles);tri = (TRIANGLE*) realloc(tri , (ntri + numOftriangle) * sizeof(TRIANGLE) ); for (int i=0; i<numOftriangle; i++) //将当前立方体中包含的三角面片保存tri[ntri+i] = triangles[i];ntri += numOftriangle; }}}fprintf(stderr,"Total of %d triangles\n",ntri);// Now do something with the triangles .... Here I just write them to a geom filefprintf(stderr,"Writing triangles ...\n");if ((fopen_s (&fptr,"output.geom","w")) != NULL) {fprintf(stderr,"Failed to open output file\n");exit(-1);}for (int i=0;i<ntri;i++) { fprintf(fptr,"f3 "); for (int k=0;k<3;k++) { fprintf(fptr,"%g %g %g ",tri[i].p[k].x,tri[i].p[k].y,tri[i].p[k].z); } fprintf(fptr,"0.5 0.5 0.5\n"); // colour}fclose(fptr); free(data); //释放堆上的内存exit(0);}/*///////////////////////////////////////////////////////////////////////// 4--------5 *---4----* /| /| /| /| / | / | 7 | 5 | / | / | / 8 / 9 7--------6 | *----6---* | | | | | | | | | | 0----|---1 | *---0|---* | / | / 11 / 10 / | / | / | 3 | 1 |/ |/ |/ |/ 3--------2 *---2----*////////////////////////////////////////////////////////////////////////////*/int PolygoniseCube(GRIDCELL grid, double isolevel, TRIANGLE* tri){int numOftriangle = 0;int cubeIndex = 0;XYZ vertlist[12]; //保存等值面与立方体各边相交的坐标//确定那个顶点位于等值面内部if (grid.val[0] < isolevel) cubeIndex |= 1;if (grid.val[1] < isolevel) cubeIndex |= 2;if (grid.val[2] < isolevel) cubeIndex |= 4;if (grid.val[3] < isolevel) cubeIndex |= 8;if (grid.val[4] < isolevel) cubeIndex |= 16;if (grid.val[5] < isolevel) cubeIndex |= 32;if (grid.val[6] < isolevel) cubeIndex |= 64;if (grid.val[7] < isolevel) cubeIndex |= 128;//异常:立方体所有顶点都在或者都不在等值面内部if ( edgeTable[cubeIndex] == 0 )return 0;//确定等值面与立方体交点坐标if (edgeTable[cubeIndex] & 1) {vertlist[0] = VertexInterp(isolevel,grid.p[0],grid.p[1],grid.val[0],grid.val[1]); }if (edgeTable[cubeIndex] & 2) {vertlist[1] = VertexInterp(isolevel,grid.p[1],grid.p[2],grid.val[1],grid.val[2]);}if (edgeTable[cubeIndex] & 4) {vertlist[2] = VertexInterp(isolevel,grid.p[2],grid.p[3],grid.val[2],grid.val[3]);}if (edgeTable[cubeIndex] & 8) {vertlist[3] = VertexInterp(isolevel,grid.p[3],grid.p[0],grid.val[3],grid.val[0]);}if (edgeTable[cubeIndex] & 16) {vertlist[4] = VertexInterp(isolevel,grid.p[4],grid.p[5],grid.val[4],grid.val[5]);}if (edgeTable[cubeIndex] & 32) {vertlist[5] = VertexInterp(isolevel,grid.p[5],grid.p[6],grid.val[5],grid.val[6]);}if (edgeTable[cubeIndex] & 64) {vertlist[6] = VertexInterp(isolevel,grid.p[6],grid.p[7],grid.val[6],grid.val[7]);}if (edgeTable[cubeIndex] & 128) {vertlist[7] = VertexInterp(isolevel,grid.p[7],grid.p[4],grid.val[7],grid.val[4]);}if (edgeTable[cubeIndex] & 256) {vertlist[8] = VertexInterp(isolevel,grid.p[0],grid.p[4],grid.val[0],grid.val[4]);}if (edgeTable[cubeIndex] & 512) {vertlist[9] = VertexInterp(isolevel,grid.p[1],grid.p[5],grid.val[1],grid.val[5]);}if (edgeTable[cubeIndex] & 1024) {vertlist[10] = VertexInterp(isolevel,grid.p[2],grid.p[6],grid.val[2],grid.val[6]);}if (edgeTable[cubeIndex] & 2048) {vertlist[11] = VertexInterp(isolevel,grid.p[3],grid.p[7],grid.val[3],grid.val[7]);}//根据交点坐标确定三角形面片,并进行保存for (int i=0; triTable[cubeIndex][i] != -1; i+=3){ tri[numOftriangle].p[0] = vertlist[ triTable[cubeIndex][i ] ]; tri[numOftriangle].p[1] = vertlist[ triTable[cubeIndex][i+1] ]; tri[numOftriangle].p[2] = vertlist[ triTable[cubeIndex][i+2] ]; numOftriangle++;}return (numOftriangle);}//等值面与立方体交点坐标XYZ VertexInterp(double isolevel, XYZ p1, XYZ p2, double valp1, double valp2){XYZ p;if (ABS(isolevel-valp1) < 0.00001)return(p1); if (ABS(isolevel-valp2) < 0.00001)return(p2); if (ABS(valp1-valp2) < 0.00001)return(p1);double coef = (isolevel - valp1 ) / (valp2 - valp1);p.x = p1.x + coef * (p2.x - p1.x); p.y = p1.y + coef * (p2.y - p1.y);p.z = p1.z + coef * (p2.z - p1.z); return (p);}
3.学习心得
1.学会定义结构体,例如三维空间的立方体,需要先定义一个坐标XYZ结构体,然后再定义一个立方体结构体,而他们之间恰恰是相互依赖的。
typedef struct {double x,y,z;} XYZ;typedef struct{XYZ p[8]; //顶点坐标XYZ n[8]; //法向量double val[8]; //顶点对应的灰度值} GRIDCELL; //立方体2.如何动态为三位数据场分配空间?data = (short int***) malloc( NX * sizeof(short int **) );for (int i=0;i<NX;i++)data[i] = (short int**) malloc( NY * sizeof(short int *));for (int i=0;i<NX;i++)for (int j=0;j<NY;j++)data[i][j] = (short int*) malloc(NZ * sizeof(short int));我们要注意到,利用malloc()函数进行内存申请之后返回来的是void*型的指针,我们需要进行指针的强制类型转换。该段代码也给我一个很深的感悟就是"指针确实是一门艺术"。3.读取文件以及数据元素读取
errno_t err;err = fopen_s(&fptr,"mri.raw", "rb");if( err != NULL ) {fprintf(stderr, "File open failed\n");exit(-1);}读取文件可以采用fopen或fopen_s两个函数,不同之处在于前者需要提供两个参数,并返回FILE*类型;而后者需要提供三个参数,返回的是errno_t(实质就是int的别名)类型。相比较而言,fopen_s更安全。int str; //从文件中读取字符for (int k=0; k<NZ; k++){for(int j=0; j<NY; j++){for(int i=0; i<NX; i++){if ((str = fgetc(fptr)) == EOF) { fprintf(stderr,"Unexpected end of file\n");exit(-1);}data[i][j][k] = str;}}}fgetc()顺序读取数据。
0 0
- 面绘制经典算法:MarchingCube实现(控制台篇)
- 面绘制经典算法:MarchingCube实现(C++ OpenGl代码篇)
- MarchingCube实现(C++ OpenGl代码篇)
- 控制台经典算法
- 算法经典面试题
- 经典算法面试题
- 经典算法面试题(-)
- 经典算法面试题
- 算法经典面试题整理(java实现)
- 算法经典面试题整理(java实现)
- vtk面绘制的实现
- 经典算法面试题一览
- 经典算法面试题(二)
- 经典算法面试题整理
- 连连看算法(控制台实现)
- 控制台绘制登陆框(四) 实现简单的字符输入
- c#经典算法实现
- c#经典算法实现
- struts的工作原理
- FFmpeg滤镜使用指南
- github团队协作教程
- docker 运行实例报错
- 不同厂家网络设备添加/删除用户
- 面绘制经典算法:MarchingCube实现(控制台篇)
- 编译器错误 C3254:centos能编译正常,window编译报错
- Socket实现网络聊天室和联系人私聊功能
- 处理 Maven 项目名称红色感叹号的问题
- css3实现字体渐变效果
- html拼接在方法中添加参数
- 共享单车原理大揭秘:小编亲自示范如何“撬锁”
- C++第4次实验(基础班)—循环结构程序设计(参考答案)-项目7-2:年龄几何
- 开发WebRTC使用什么语言?