Marching Cubes算法工作原理

MarchingCubes算法工作原理

MarchingCubes算法是三维数据场等值面生成的经典算法，是体素单元内等值面抽取技术的代表。

等值面是空间中所有具有某个相同值的点的集合。它可以表示为，，c是常数。则称F(f)为体数据f中的等值面。

在MC算法中，假定原始数据是离散的三维空间规则数据场。用于医疗诊断的断层扫描(CT)及核磁共振成像(MRI)等产生的图像均属于这一类型。MC算法的基本思想是逐个处理数据场中的体素，分类出与等值面相交的体素，采用插值计算出等值面与体素棱边的交点。根据体素中每一顶点与等值面的相对位置，将等值面与立方体边的交点按一定方式连接生成等值面，作为等值面在该立方体内的一个逼近表示。在计算出关于体数据场内等值面的有关参数后山常用的图形软件包或硬件提供的面绘制功能绘制出等值面。

图2.1离散的三维空间规则数据场中的一个体素

2.1.1MC算法的主要步骤

1．确定包含等值面的体素

离散的三维空间规则数据场中的一个体素可以用图2.1表示。8个数据点分别位于该体素的8个角点上。MC算法的基本假设是:沿着体素的边其数据场呈局部连续线性变化，根据这个假设，可认为，如果两个相邻采样点一个为正点，一个为负点，则它们连成的边上一定存在且仅有一个等值点(设等值面值为c)。如果得到了体素各条边上的等值点，就可以以这些点为顶点，用一系列的三角形拟合出该体素中的等值面。因此确定立方体体素中等值面的分布是该算法的基础。

首先对体素的8个顶点进行分类，以判断其顶点是位于等值面之外，还是位于等值面之内。再根据8个顶点的状态，确定等值面的剖分模式。顶点分类规则为:

1.如体素顶点的数据值大于或等于等值面的值，则定义该顶点位于等值面之外，记为正点，即“1“

2.如体素顶点的数据值小于等值面的值，则定义该顶点位于等值面之内，记为负点，即“0"

由于每个体素共有8个顶点，且每个顶点有正负两种状态，所以等值面可能以=256种方式与一个体素相交。通过列举出这256种情况，就能创建一张表格，利用它可以查出任意体素中的等值面的三角面片表示。如果考虑互补对称性，将等值面的值和8个角点的函数值的大小关系颠倒过来，即将体素的顶点标记置反(0变为1,1变为0)，这样做不会影响该体素的8个角点和等值面之间的拓扑结构，可将256种方式简化成128种。其次，再利用旋转对称性，可将这128种构型进一步简化成15种。图3.2给出了这15种基本构型[131其中黑点标记为“1”的角点。

图2.2分布状态表

图2.3体素角点分布不同情况

基于上面的分析，MC算法中用一个字节的空间构造了一个体素状态表，如图2.2所示，该状态表中的每一位可表示出该体元中的一个角点的0或1的状态。根据这一状态表，就可知道当前体素属于图2.3中哪一种情况，以及等值面将与哪一条边相交。

2.求等值面与体元边界的交点

在确定体素内三角剖分模式后，就要计算三角片顶点位置。当三维离散数据场的密度较高时，即当体素很小时，可以假定函数值沿体素边界呈线性变化，这就是MC算法的基本假设。因此，根据这一基本假设，可以直接用线性插值计算等值面与体素边的交点。

对于当前被处理体素的某一条边，如果其两顶点，的标记值不同，那么等值面一定与此边相交，且仅有一个交点。交点为其中P代表等值点坐标,,代表两个端点的坐标，，代表两个端点的灰度值，v代表域值。求出等值面与体素棱边的交点以后，就可以将这些交点连接成三角形或多边形，形成等值面的一部分。

3.等值面的法向量计算

为了利用图形硬件显示等值面图象，必须给出形成等值面的各三角面片的法向分量，选择适当的局部面光照模型进行光照计算，以生成真实感图形。

对于等值面上的每一点，其沿面的切线方向的梯度分量应该是零，因此，该点的梯度矢量的方向也就代表了等值面在该点的法向量，当梯度值非零。所幸的是等值面往往是由两种具有不同密度的物质的分解面，因此其上的每点的梯度矢量均不为零，即

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光照模型公式，就可计算出小三角片表面的光强(灰度)。将其投影在某个特定的二维平面上进行显示，从而显示出物体富有光感的整个表面形态。其中我们在内存中保留四个切片来计算立方体中所有顶点梯度。

2.1.2MC算法流程

1、将三维离散规则数据场分层读入内存;

2、扫描两层数据，逐个构造体素，每个体素中的8个角点取自相邻的两层;

3、将体素每个角点的函数值与给定的等值面值c做比较，根据比较结果，构造

该体素的状态表;

4、根据状态表，得出将与等值面有交点的边界体素;

5、通过线性插值方法计算出体素棱边与等值面的交点;

6、利用中心差分方法，求出体素各角点处的法向量，再通过线性插值方法，求出三角面片各顶点处的法向;

7,根据各三角面片上各顶点的坐标及法向量绘制等值面图像。

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

MC代码

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

{

//记录生成的顶点数和面数初始时应该为0

m_vNumber=0;

m_fNumber=0;

//当前Cube中生成的顶点和面数

intvertPos,facePos;

//包围盒的尺寸用于绘制程序计算整个场景的包围盒，用于调整观察位置，以使整个场景尽可能占满整个窗口。

floatmin[3],max[3];

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

//当前扫描层的切片数据和一个临时的切片数据

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

pSliceA=pSliceB=pSliceC=tempSlice=NULL;

intimageWidth,imageHeight,imageSize,sliceNumber;

imageWidth=imageHeight=512;//我们是512×512的数据

imageSize=imageWidth*imageHeight;

sliceNumber=m_FileCount-1;

if((highThreshold*lowThreshold)==0)

{

return0;

}

pSliceD=newshort[imageSize];

//因为等值面是每相邻两层切片为单位进行提取的，所以在处理后两层切片时难免生成前两层切片已经生成的顶点，这时候就用下面的数组记录哪些边上的顶点已经生成了，如果遇到已经生成的顶点就不再重复计算而是直接使用记录的索引，否则就生成新的顶点。

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;//错误

}

//初始化数据

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);

//计算某一层顶点和三角时所需要的一些变量

//一些循环变量

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;

//得到数据

pSliceD=m_volumeData;

pSliceB=tempSlice;

pSliceA=tempSlice;

inttt,tt1;

//扫描4层切片的顺序

/*

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

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

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

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

V

*/

//marchingcubes算法开始实行？第一次循环时，只读入一个切片？

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++)//调试用

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…..

{

//根据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//根据g_EdgeTable对应值判断cube的那一条边与等值面有交点

//该边上的顶点是否已经在上一层中生成

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;

//更新包围盒数据

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

{

//若该点已经在上一层生成，则直接得到其索引

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);

//记录新生成的三角面片数据

w=0;

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

{

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

{

triIndex[facePos