1.MarchingCube

marchingcube是一個(gè)比較經(jīng)典而古老的算法，也是面繪制中應(yīng)用比較多的算法，Marchingcube發(fā)展到今天也遇到了幾何拓?fù)?、一致性的?wèn)題仍待改善。本文研究的就是經(jīng)典的marchingcube(Paul,1997).

關(guān)于MarchingCube原理（含代碼），重點(diǎn)推薦如下：

1.《醫(yī)學(xué)圖像三維重建和可視化:VC++實(shí)現(xiàn)實(shí)例》(張俊華)

2.Bourke P. Polygonising a scalar field[J]. 1997.

3.論文+源代碼+交互界面

2.控制臺(tái)上實(shí)現(xiàn)與討論

算法設(shè)計(jì)步驟：

1.二維斷層圖像順序地讀入內(nèi)存并構(gòu)成三位數(shù)據(jù)場(chǎng)；

2.依次掃描相鄰兩層數(shù)據(jù)，逐個(gè)構(gòu)造立方體；

3.將立方體每個(gè)頂點(diǎn)灰度和給定等值面閾值進(jìn)行比較，計(jì)算索引值；

4.對(duì)于含有等值面的立方體，通過(guò)灰度差分計(jì)算立方體各頂點(diǎn)的梯度；

5.根據(jù)索引值查找邊索引表，獲得和等值面有交點(diǎn)的當(dāng)前立方體的相交邊；

6.根據(jù)相交邊的兩個(gè)頂點(diǎn)坐標(biāo)及其法向量，利用差值（中值）計(jì)算等值點(diǎn)的坐標(biāo)與法向量；

7.根據(jù)索引值查找三角片索引表，確定當(dāng)前立方體內(nèi)構(gòu)成三角片的等值點(diǎn)的組合方式；

8.由各個(gè)立方體內(nèi)的三角面片構(gòu)成等值面。

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <iostream> using namespace std;static int edgeTable[256]={ }; //根據(jù)索引值確定立方體和等值面交點(diǎn)情況 static int triTable[256][16] ={ };//根據(jù)索引值確定三角面片的個(gè)數(shù)以及位置typedef struct {double x,y,z; } XYZ;typedef struct {XYZ p[8]; //頂點(diǎn)坐標(biāo)XYZ n[8]; //法向量double val[8]; //頂點(diǎn)對(duì)應(yīng)的灰度值 } GRIDCELL; //立方體typedef struct {XYZ p[3]; //頂點(diǎn)坐標(biāo)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數(shù)據(jù)場(chǎng)大小 #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; //等值面 數(shù)據(jù)場(chǎng)的范圍GRIDCELL grid;TRIANGLE triangles[10];TRIANGLE* tri = NULL; //保存的是所有三角形面片的頂點(diǎn)坐標(biāo)int ntri = 0;FILE* fptr = NULL;//打開(kāi)并閱讀原始文件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);}//為三維數(shù)據(jù)場(chǎng)分配空間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順序讀取數(shù)據(jù)，直到遇見(jiàn)-文件結(jié)束符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數(shù)據(jù)場(chǎng)處理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];//計(jì)算當(dāng)前立方體中的三角面片數(shù)目int numOftriangle = PolygoniseCube(grid, isolevel, triangles);tri = (TRIANGLE*) realloc(tri , (ntri + numOftriangle) * sizeof(TRIANGLE) ); for (int i=0; i<numOftriangle; i++) //將當(dāng)前立方體中包含的三角面片保存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); //釋放堆上的內(nèi)存exit(0); } /*/ 4--------5 *---4----*/| /| /| /|/ | / | 7 | 5 |/ | / | / 8 / 97--------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]; //保存等值面與立方體各邊相交的坐標(biāo)//確定那個(gè)頂點(diǎn)位于等值面內(nèi)部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;//異常：立方體所有頂點(diǎn)都在或者都不在等值面內(nèi)部if ( edgeTable[cubeIndex] == 0 )return 0;//確定等值面與立方體交點(diǎn)坐標(biāo)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]);}//根據(jù)交點(diǎn)坐標(biāo)確定三角形面片，并進(jìn)行保存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); } //等值面與立方體交點(diǎn)坐標(biāo) 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.學(xué)習(xí)心得

1.學(xué)會(huì)定義結(jié)構(gòu)體，例如三維空間的立方體，需要先定義一個(gè)坐標(biāo)XYZ結(jié)構(gòu)體，然后再定義一個(gè)立方體結(jié)構(gòu)體，而他們之間恰恰是相互依賴(lài)的。

typedef struct {double x,y,z; } XYZ;typedef struct {XYZ p[8]; //頂點(diǎn)坐標(biāo)XYZ n[8]; //法向量double val[8]; //頂點(diǎn)對(duì)應(yīng)的灰度值 } GRIDCELL; //立方體2.如何動(dòng)態(tài)為三位數(shù)據(jù)場(chǎng)分配空間？

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()函數(shù)進(jìn)行內(nèi)存申請(qǐng)之后返回來(lái)的是void*型的指針，我們需要進(jìn)行指針的強(qiáng)制類(lèi)型轉(zhuǎn)換。該段代碼也給我一個(gè)很深的感悟就是"指針確實(shí)是一門(mén)藝術(shù)"。

3.讀取文件以及數(shù)據(jù)元素讀取

errno_t err; err = fopen_s(&fptr,"mri.raw", "rb"); if( err != NULL ) {fprintf(stderr, "File open failed\n");exit(-1); }讀取文件可以采用fopen或fopen_s兩個(gè)函數(shù)，不同之處在于前者需要提供兩個(gè)參數(shù)，并返回FILE*類(lèi)型；而后者需要提供三個(gè)參數(shù)，返回的是errno_t(實(shí)質(zhì)就是int的別名)類(lèi)型。相比較而言，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()順序讀取數(shù)據(jù)。

總結(jié)

以上是生活随笔為你收集整理的面绘制经典算法：MarchingCube实现(控制台篇)的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。