Marching Cubes算法工作原理

MarchingCubes算法工作原理

MarchingCubes算法是三維數(shù)據(jù)場等值面生成的經(jīng)典算法，是體素單元內(nèi)等值面抽取技術(shù)的代表。

等值面是空間中所有具有某個相同值的點的集合。它可以表示為，，c是常數(shù)。則稱F(f)為體數(shù)據(jù)f中的等值面。

在MC算法中，假定原始數(shù)據(jù)是離散的三維空間規(guī)則數(shù)據(jù)場。用于醫(yī)療診斷的斷層掃描(CT)及核磁共振成像(MRI)等產(chǎn)生的圖像均屬于這一類型。MC算法的基本思想是逐個處理數(shù)據(jù)場中的體素，分類出與等值面相交的體素，采用插值計算出等值面與體素棱邊的交點。根據(jù)體素中每一頂點與等值面的相對位置，將等值面與立方體邊的交點按一定方式連接生成等值面，作為等值面在該立方體內(nèi)的一個逼近表示。在計算出關(guān)于體數(shù)據(jù)場內(nèi)等值面的有關(guān)參數(shù)后山常用的圖形軟件包或硬件提供的面繪制功能繪制出等值面。

圖2.1離散的三維空間規(guī)則數(shù)據(jù)場中的一個體素

2.1.1MC算法的主要步驟

1．確定包含等值面的體素

離散的三維空間規(guī)則數(shù)據(jù)場中的一個體素可以用圖2.1表示。8個數(shù)據(jù)點分別位于該體素的8個角點上。MC算法的基本假設是:沿著體素的邊其數(shù)據(jù)場呈局部連續(xù)線性變化，根據(jù)這個假設，可認為，如果兩個相鄰采樣點一個為正點，一個為負點，則它們連成的邊上一定存在且僅有一個等值點(設等值面值為c)。如果得到了體素各條邊上的等值點，就可以以這些點為頂點，用一系列的三角形擬合出該體素中的等值面。因此確定立方體體素中等值面的分布是該算法的基礎。

首先對體素的8個頂點進行分類，以判斷其頂點是位于等值面之外，還是位于等值面之內(nèi)。再根據(jù)8個頂點的狀態(tài)，確定等值面的剖分模式。頂點分類規(guī)則為:

1.如體素頂點的數(shù)據(jù)值大于或等于等值面的值，則定義該頂點位于等值面之外，記為正點，即“1“

2.如體素頂點的數(shù)據(jù)值小于等值面的值，則定義該頂點位于等值面之內(nèi)，記為負點，即“0"

由于每個體素共有8個頂點，且每個頂點有正負兩種狀態(tài)，所以等值面可能以=256種方式與一個體素相交。通過列舉出這256種情況，就能創(chuàng)建一張表格，利用它可以查出任意體素中的等值面的三角面片表示。如果考慮互補對稱性，將等值面的值和8個角點的函數(shù)值的大小關(guān)系顛倒過來，即將體素的頂點標記置反(0變?yōu)?,1變?yōu)?)，這樣做不會影響該體素的8個角點和等值面之間的拓撲結(jié)構(gòu)，可將256種方式簡化成128種。其次，再利用旋轉(zhuǎn)對稱性，可將這128種構(gòu)型進一步簡化成15種。圖3.2給出了這15種基本構(gòu)型[131其中黑點標記為“1”的角點。

圖2.2分布狀態(tài)表

圖2.3體素角點分布不同情況

基于上面的分析，MC算法中用一個字節(jié)的空間構(gòu)造了一個體素狀態(tài)表，如圖2.2所示，該狀態(tài)表中的每一位可表示出該體元中的一個角點的0或1的狀態(tài)。根據(jù)這一狀態(tài)表，就可知道當前體素屬于圖2.3中哪一種情況，以及等值面將與哪一條邊相交。

2.求等值面與體元邊界的交點

在確定體素內(nèi)三角剖分模式后，就要計算三角片頂點位置。當三維離散數(shù)據(jù)場的密度較高時，即當體素很小時，可以假定函數(shù)值沿體素邊界呈線性變化，這就是MC算法的基本假設。因此，根據(jù)這一基本假設，可以直接用線性插值計算等值面與體素邊的交點。

對于當前被處理體素的某一條邊，如果其兩頂點，的標記值不同，那么等值面一定與此邊相交，且僅有一個交點。交點為其中P代表等值點坐標,,代表兩個端點的坐標，，代表兩個端點的灰度值，v代表域值。求出等值面與體素棱邊的交點以后，就可以將這些交點連接成三角形或多邊形，形成等值面的一部分。

3.等值面的法向量計算

為了利用圖形硬件顯示等值面圖象，必須給出形成等值面的各三角面片的法向分量，選擇適當?shù)木植棵婀庹漳Ｐ瓦M行光照計算，以生成真實感圖形。

對于等值面上的每一點，其沿面的切線方向的梯度分量應該是零，因此，該點的梯度矢量的方向也就代表了等值面在該點的法向量，當梯度值非零。所幸的是等值面往往是由兩種具有不同密度的物質(zhì)的分解面，因此其上的每點的梯度矢量均不為零，即

Mc算法采用中心差分方法求采樣點p〔m,n,k)處的梯度矢量，公式如下:

Gx=〔g(i+1,j,k)-g(i-1,j,k)〕/2dx

Gy=〔g(i,j+1,k)-g(i,j-1,k)〕/2dy

Gz=〔g(i,j,k+1)-g(i,j,k-1)〕/2dz

其中D（i，j,k）是切片k在像素（i，j）的密度，,，是立方體邊的長度。對g進行歸一化，得到(gx/｜g｜，gy/｜g｜，gz/｜g｜）作為（i,j,k)上的單位法向量。然后，對體素八個頂點上法向量進行線性插值就可得到位于體素棱邊上的三角片的各個頂點上的法向量。設計算得到的某個三角片的三個頂點上的單位法向量分別為(,和)，這個三角片的幾何重心為，則該三角片的法向量起始于，終止于。代入Gourand光照模型公式，就可計算出小三角片表面的光強(灰度)。將其投影在某個特定的二維平面上進行顯示，從而顯示出物體富有光感的整個表面形態(tài)。其中我們在內(nèi)存中保留四個切片來計算立方體中所有頂點梯度。

2.1.2MC算法流程

1、將三維離散規(guī)則數(shù)據(jù)場分層讀入內(nèi)存;

2、掃描兩層數(shù)據(jù)，逐個構(gòu)造體素，每個體素中的8個角點取自相鄰的兩層;

3、將體素每個角點的函數(shù)值與給定的等值面值c做比較，根據(jù)比較結(jié)果，構(gòu)造

該體素的狀態(tài)表;

4、根據(jù)狀態(tài)表，得出將與等值面有交點的邊界體素;

5、通過線性插值方法計算出體素棱邊與等值面的交點;

6、利用中心差分方法，求出體素各角點處的法向量，再通過線性插值方法，求出三角面片各頂點處的法向;

7,根據(jù)各三角面片上各頂點的坐標及法向量繪制等值面圖像。

========================================================

MC代碼

MarchingCubes(floatlowThreshold,floathighThreshold,floatXSpace,floatYSpace,floatZSpace)

{

//記錄生成的頂點數(shù)和面數(shù)初始時應該為0

m_vNumber=0;

m_fNumber=0;

//當前Cube中生成的頂點和面數(shù)

intvertPos,facePos;

//包圍盒的尺寸用于繪制程序計算整個場景的包圍盒，用于調(diào)整觀察位置，以使整個場景盡可能占滿整個窗口。

floatmin[3],max[3];

min[0]=min[1]=min[2]=max[0]=max[1]=max[2]=0;//初始化

//當前掃描層的切片數(shù)據(jù)和一個臨時的切片數(shù)據(jù)

short*pSliceA,*pSliceB,*pSliceC,*pSliceD,*tempSlice;

pSliceA=pSliceB=pSliceC=tempSlice=NULL;

intimageWidth,imageHeight,imageSize,sliceNumber;

imageWidth=imageHeight=512;//我們是512×512的數(shù)據(jù)

imageSize=imageWidth*imageHeight;

sliceNumber=m_FileCount-1;

if((highThreshold*lowThreshold)==0)

{

return0;

}

pSliceD=newshort[imageSize];

//因為等值面是每相鄰兩層切片為單位進行提取的，所以在處理后兩層切片時難免生成前兩層切片已經(jīng)生成的頂點，這時候就用下面的數(shù)組記錄哪些邊上的頂點已經(jīng)生成了，如果遇到已經(jīng)生成的頂點就不再重復計算而是直接使用記錄的索引，否則就生成新的頂點。

long*bottomXEdge=newlong[imageSize];

long*bottomYEdge=newlong[imageSize];

long*topXEdge=newlong[imageSize];

long*topYEdge=newlong[imageSize];

long*zEdge=newlong[imageSize];

tempSlice=newshort[imageSize];

if(bottomXEdge==NULL||bottomYEdge==NULL||

topXEdge==NULL||topYEdge==NULL||

zEdge==NULL||tempSlice==NULL)

{

return0;//錯誤

}

//初始化數(shù)據(jù)

memset(bottomXEdge,-1,sizeof(long)*imageSize);

memset(bottomYEdge,-1,sizeof(long)*imageSize);

memset(topXEdge,-1,sizeof(long)*imageSize);

memset(topYEdge,-1,sizeof(long)*imageSize);

memset(zEdge,-1,sizeof(long)*imageSize);

memset(tempSlice,0,sizeof(short)*imageSize);

//計算某一層頂點和三角時所需要的一些變量

//一些循環(huán)變量

inti,j,k,w,r;

//cube類型

unsignedcharcubeType(0);

//計算法向量

floatdx[8],dy[8],dz[8],squaroot;

//記錄某個Cube生成

floatvertPoint[12][6];

intcellVerts[12];//whatuse

inttriIndex[5][3];//每個cube最多生成5條邊

//用于記錄已生成頂點索引的臨時變量

intoffset;

//當前cube8個頂點的灰度值

shortcubegrid[8];

long*edgeGroup;

//得到數(shù)據(jù)

pSliceD=m_volumeData;

pSliceB=tempSlice;

pSliceA=tempSlice;

inttt,tt1;

//掃描4層切片的順序

/*

-----------------------D|

-----------------------B|

-----------------------C|

-----------------------A|

V

*/

//marchingcubes算法開始實行？第一次循環(huán)時，只讀入一個切片？

for(i=0;i<=(sliceNumber);i++)

{

pSliceC=pSliceA;

pSliceA=pSliceB;

pSliceB=pSliceD;

if(i>=sliceNumber-2)

{

pSliceD=tempSlice;

}

else

{

pSliceD+=imageSize;

}

for(j=0;j<imageHeight-1;++j)

for(k=0;k<imageWidth-1;++k)

/*for(j=10;j<imageHeight-5;j++)//調(diào)試用

for(k=10;k<imageWidth-5;k++)*/

{

//得到八個頂點的灰度值step0

cubegrid[0]=pSliceA[j*imageWidth+k];

cubegrid[1]=pSliceA[j*imageWidth+k+1];

cubegrid[2]=pSliceA[(j+1)*imageWidth+k+1];

cubegrid[3]=pSliceA[(j+1)*imageWidth+k];

cubegrid[4]=pSliceB[j*imageWidth+k];

cubegrid[5]=pSliceB[j*imageWidth+k+1];

cubegrid[6]=pSliceB[(j+1)*imageWidth+k+1];

cubegrid[7]=pSliceB[(j+1)*imageWidth+k];

//計算cube的類型

cubeType=0;

for(w=0;w<8;w++)

{

if((cubegrid[w]>lowThreshold)&&(cubegrid[w]<highThreshold))//需要畫的點

{

cubeType|=(1<<w);

}

}//endfor計算cube的類型

if(cubeType==0||cubeType==255)

{

continue;

}

for(w=0;w<12;w++)//初始化cellVerts表到零

{

cellVerts[w]=-1;

}

//計算6個方向相鄰點的象素差值（用于計算法向量）

if(k==0)

{

dx[0]=pSliceA[j*imageWidth+1];

dx[3]=pSliceA[(j+1)*imageWidth+1];

dx[4]=pSliceB[j*imageWidth+1];

dx[7]=pSliceB[(j+1)*imageWidth+1];

}

else

{

dx[0]=pSliceA[j*imageWidth+k+1]

-pSliceA[j*imageWidth+k-1];

dx[3]=pSliceA[(j+1)*imageWidth+k+1]

-pSliceA[(j+1)*imageWidth+k-1];

dx[4]=pSliceB[j*imageWidth+k+1]

-pSliceB[j*imageWidth+k-1];

dx[7]=pSliceB[(j+1)*imageWidth+k+1]

-pSliceB[(j+1)*imageWidth+k-1];

}

if(k==imageWidth-2)

{

dx[1]=-pSliceA[j*imageWidth+imageWidth-2];

dx[2]=-pSliceA[(j+1)*imageWidth+imageWidth-2];

dx[5]=-pSliceB[j*imageWidth+imageWidth-2];

dx[6]=-pSliceB[(j+1)*imageWidth+imageWidth-2];

}

else

{

dx[1]=pSliceA[j*imageWidth+k+2]

-pSliceA[j*imageWidth+k];

dx[2]=pSliceA[(j+1)*imageWidth+k+2]

-pSliceA[(j+1)*imageWidth+k];

dx[5]=pSliceB[j*imageWidth+k+2]

-pSliceB[j*imageWidth+k];

dx[6]=pSliceB[(j+1)*imageWidth+k+2]

-pSliceB[(j+1)*imageWidth+k];

}

if(j==0)

{

dy[0]=pSliceA[imageWidth+k];

dy[1]=pSliceA[imageWidth+k+1];

dy[4]=pSliceB[imageWidth+k];

dy[5]=pSliceB[imageWidth+k+1];

}

else

{

dy[0]=pSliceA[(j+1)*imageWidth+k]

-pSliceA[(j-1)*imageWidth+k];

dy[1]=pSliceA[(j+1)*imageWidth+k+1]

-pSliceA[(j-1)*imageWidth+k+1];

dy[4]=pSliceB[(j+1)*imageWidth+k]

-pSliceB[(j-1)*imageWidth+k];

dy[5]=pSliceB[(j+1)*imageWidth+k+1]

-pSliceB[(j-1)*imageWidth+k+1];

}

if(j==imageHeight-2)

{

dy[2]=-pSliceA[(imageHeight-2)*imageWidth+k+1];

dy[3]=-pSliceA[(imageHeight-2)*imageWidth+k];

dy[6]=-pSliceB[(imageHeight-2)*imageWidth+k+1];

dy[7]=-pSliceB[(imageHeight-2)*imageWidth+k];

}

else

{

dy[2]=pSliceA[(j+2)*imageWidth+k+1]-pSliceA[j*imageWidth+k+1];

dy[3]=pSliceA[(j+2)*imageWidth+k]-pSliceA[j*imageWidth+k];

dy[6]=pSliceB[(j+2)*imageWidth+k+1]-pSliceB[j*imageWidth+k+1];

dy[7]=pSliceB[(j+2)*imageWidth+k]-pSliceB[j*imageWidth+k];

}

dz[0]=pSliceB[j*imageWidth+k]

-pSliceC[j*imageWidth+k];

dz[1]=pSliceB[j*imageWidth+k+1]

-pSliceC[j*imageWidth+k+1];

dz[2]=pSliceB[(j+1)*imageWidth+k+1]

-pSliceC[(j+1)*imageWidth+k+1];

dz[3]=pSliceB[(j+1)*imageWidth+k]

-pSliceC[(j+1)*imageWidth+k];

dz[4]=pSliceD[j*imageWidth+k]

-pSliceA[j*imageWidth+k];

dz[5]=pSliceD[j*imageWidth+k+1]

-pSliceA[j*imageWidth+k+1];

dz[6]=pSliceD[(j+1)*imageWidth+k+1]

-pSliceA[(j+1)*imageWidth+k+1];

dz[7]=pSliceD[(j+1)*imageWidth+k]

-pSliceA[(j+1)*imageWidth+k];

//計算三角形頂點的坐標和梯度

vertPos=0;

facePos=0;

for(w=0;w<12;w++)

{

if(g_EdgeTable[cubeType]&(1<<w))//what…..

{

//根據(jù)g_edgeTable[256]對應值判斷cube的那一條邊與等值面有交點

switch(w)

{

case0:

offset=j*imageWidth+k;

edgeGroup=bottomXEdge;

break;

case1:

offset=j*imageWidth+k+1;

edgeGroup=bottomYEdge;

break;

case2:

offset=(j+1)*imageWidth+k;

edgeGroup=bottomXEdge;

break;

case3:

offset=j*imageWidth+k;

edgeGroup=bottomYEdge;

break;

case4:

offset=j*imageWidth+k;

edgeGroup=topXEdge;

break;

case5:

offset=j*imageWidth+k+1;

edgeGroup=topYEdge;

break;

case6:

offset=(j+1)*imageWidth+k;

edgeGroup=topXEdge;

break;

case7:

offset=j*imageWidth+k;

edgeGroup=topYEdge;

break;

case8:

offset=j*imageWidth+k;

edgeGroup=zEdge;

break;

case9:

offset=j*imageWidth+k+1;

edgeGroup=zEdge;

break;

case10:

offset=(j+1)*imageWidth+k+1;

edgeGroup=zEdge;

break;

case11:

offset=(j+1)*imageWidth+k;

edgeGroup=zEdge;

break;

}//對應switch的{。。。endfor//根據(jù)g_EdgeTable對應值判斷cube的那一條邊與等值面有交點

//該邊上的頂點是否已經(jīng)在上一層中生成

if(edgeGroup[offset]==-1)

{

intindex1,index2;

shorts1,s2,s;

floatx1,y1,z1,nx1,ny1,nz1;

floatx2,y2,z2,nx2,ny2,nz2;

//得到該邊兩端點的索引進而得到兩點的灰度值

index1=g_CoordTable[w][3];

index2=g_CoordTable[w][4];

s1=cubegrid[index1];

s2=cubegrid[index2];

if(s1<highThreshold&&s1>lowThreshold)

{

if(s2>=highThreshold)

{

s=highThreshold;

}

elseif(s2<=lowThreshold)

{

s=lowThreshold;

}

}

elseif(s2<highThreshold&&s2>lowThreshold)

{

if(s1>=highThreshold)

{

s=highThreshold;

}

elseif(s1<=lowThreshold)

{

s=lowThreshold;

}

}

//計算兩端點實際坐標

x1=(k+g_CoordVertex[index1][0])*XSpace;

y1=(j+g_CoordVertex[index1][1])*YSpace;

z1=(i+g_CoordVertex[index1][2])*ZSpace;

x2=(k+g_CoordVertex[index2][0])*XSpace;

y2=(j+g_CoordVertex[index2][1])*YSpace;

z2=(i+g_CoordVertex[index2][2])*ZSpace;

//計算兩端點的法向量

nx1=dx[index1]/XSpace;

ny1=dy[index1]/YSpace;

nz1=dz[index1]/ZSpace;

nx2=dx[index2]/XSpace;

ny2=dy[index2]/YSpace;

nz2=dz[index2]/ZSpace;

floatfactor=((float)(s-s1))/((float)(s2-s1));

//插值計算交點坐標

vertPoint[vertPos][0]=factor*(x2-x1)+x1;

vertPoint[vertPos][1]=factor*(y2-y1)+y1;

vertPoint[vertPos][2]=factor*(z2-z1)+z1;

//計算法向量

vertPoint[vertPos][3]=factor*(nx1-nx2)-nx1;

vertPoint[vertPos][4]=factor*(ny1-ny2)-ny1;

vertPoint[vertPos][5]=factor*(nz1-nz2)-nz1;

//法向量歸一化

squaroot=sqrt(vertPoint[vertPos][3]*vertPoint[vertPos][3]+vertPoint[vertPos][4]*vertPoint[vertPos][4]

+vertPoint[vertPos][5]*vertPoint[vertPos][5]);

if(squaroot<=0)squaroot=1.0;

vertPoint[vertPos][3]/=squaroot;

vertPoint[vertPos][4]/=squaroot;

vertPoint[vertPos][5]/=squaroot;

//更新包圍盒數(shù)據(jù)

if(min[0]>vertPoint[vertPos][0])

{

min[0]=vertPoint[vertPos][0];

}

if(min[1]>vertPoint[vertPos][1])

{

min[1]=vertPoint[vertPos][1];

}

if(min[2]>vertPoint[vertPos][2])

{

min[2]=vertPoint[vertPos][2];

}

if(max[0]<vertPoint[vertPos][0])

{

max[0]=vertPoint[vertPos][0];

}

if(max[1]<vertPoint[vertPos][1])

{

max[1]=vertPoint[vertPos][1];

}

if(max[2]<vertPoint[vertPos][2])

{

max[2]=vertPoint[vertPos][2];

}

//記錄新生成的頂點索引

cellVerts[w]=m_vNumber;

edgeGroup[offset]=cellVerts[w];

m_vNumber++;

vertPos++;

}//endif(edgeGroup[offset]==-1)////

else

{

//若該點已經(jīng)在上一層生成，則直接得到其索引

cellVerts[w]=edgeGroup[offset];

}

}//end對應if(g_EdgeTable[cubeType]&(1<<w))//

}//對應for(w=0;w<12;w++)

//保存當前cubes頂點和法向量

tt1=m_vNumber-vertPos;

for(tt=0;tt<vertPos;tt++)

{

vPointNomal[tt1+tt][0]=vertPoint[tt][0];

vPointNomal[tt1+tt][1]=vertPoint[tt][1];

vPointNomal[tt1+tt][2]=vertPoint[tt][2];

vPointNomal[tt1+tt][3]=vertPoint[tt][3];

vPointNomal[tt1+tt][4]=vertPoint[tt][4];

vPointNomal[tt1+tt][5]=vertPoint[tt][5];

}

//memcpy(vPointNomal+6*(m_vNumber-vertPos),vertPoint,sizeof(float)*6*vertPos);

//記錄新生成的三角面片數(shù)據(jù)

w=0;

while(g_TriTable[cubeType][w]!=-1)

{

for(r=0;r<3;r++)

{

triIndex[facePos