[1] Lorensen W E, Cline H E .Marching cubes: a high-resoulution 3D suface construction algorithm [J], Computer Graphics,1987, 21(4):163~169
[2]集成化医学影像算法平台理论与实践田捷，赵明昌，何晖光 清华大学出版社2004年10月
[3]Polygonising a scalar field Also known as: "3D Contouring", "Marching Cubes", "Surface Reconstruction"
http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/Marching Cubes；
[4]www.3dmed.net
Marching Cubes算法工作原理
Marching Cubes算法是三维数据场等值面生成的经典算法，是体素单元内等值面抽取技术的代表。
等值面是空间中所有具有某个相同值的点的集合。它可以表示为， ，c是常数。则称F(f)为体数据f中的等值面。
在MC算法中，假定原始数据是离散的三维空间规则数据场。用于医疗诊断的断层扫描(CT)及核磁共振成像(MRI) 等产生的图像均属于这一类型。MC算法的基本思想是逐个处理数据场中的体素，分类出与等值面相交的体素，采用插值计算出等值面与体弯陪素棱边的交点。根据体素中每一顶点与等值面的相对位置，将等值面与立方体边的交点按一定方式连接生成等值面，作为等值面在该立方体内的一个逼近表示。在计算出关于体数据场内等值面的有关参数后山常用的图形软件包或硬件提供的面绘制功能绘制出等值面。
图2.1 离散的三维空间规则数据场中的一个体素
2.1.1 MC算法的主要步骤
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.2 MC算法流程
1、将三维离散规则数据场分层读入内存
2、扫描两层数据，逐个构造体素，每个体素中的8个角点取自相邻的两层
3、将体素每个角点的函数值与给定的等值面值c做比较，根据比较结果，构造
该体素的状态表
4、根据状态表，得出将与等值面有交点的边界体素
5、通过线性插值方法计算出体素棱边与等值面的交点
6、利用中心差分方法，求出体素各角点处的法向量，再通过线性插值方法，求出三角面片各顶点处的法向
7,根据各三角面片上各顶点的坐标及法向量绘制等值面图像。
========================================================
MC代码
MarchingCubes(float lowThreshold,float highThreshold,float XSpace,float YSpace,float ZSpace)
{
//记录生成的顶点数和面数初始时应该为0
m_vNumber=0
m_fNumber=0
//当前Cube中生成的顶点和面数
int vertPos,facePos
//包围盒的尺寸用于绘制程序计算整个场景的包围盒，用于调整观察位置，以使整个场景尽可能占满整个窗口。
float min[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
int imageWidth,imageHeight,imageSize,sliceNumber
imageWidth=imageHeight=512//我们是512×512的数据
imageSize=imageWidth*imageHeight
sliceNumber=m_FileCount-1
if((highThreshold*lowThreshold)==0)
{
return 0
}
pSliceD =new short [imageSize]
//因为等值面是每相邻两层切片为单位进行提取的，所以在处理后两层切片时难免生成前两层切片已经生成的顶点，这时候就用下面的数组记录哪些边上的顶点已经生成了，如果遇到已经生成的顶点就不再重复计算而是直接使用记录的索引，否则就生成新的顶点。
long *bottomXEdge=new long[imageSize]
long *bottomYEdge=new long[imageSize]
long *topXEdge=new long[imageSize]
long *topYEdge=new long[imageSize]
long *zEdge=new long[imageSize]
tempSlice=new short [imageSize]
if(bottomXEdge==NULL||bottomYEdge==NULL||
topXEdge==NULL||topYEdge==NULL||
zEdge==NULL||tempSlice==NULL)
{
return 0//错误
}
//初始化数据
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)
//计算某一层顶点和三角时所需要的一些变量
//一些循环变量
int i,j,k,w,r
//cube类型
unsigned char cubeType(0)
//计算法向量
float dx[8],dy[8],dz[8],squaroot
//记录某个Cube生成
float vertPoint[12][6]
int cellVerts[12]//what use
int triIndex[5][3]//每个cube最多生成5条边
//用于记录已生成顶点索引的临时变量
int offset
//当前cube8个顶点的灰度值
short cubegrid[8]
long *edgeGroup
//得到数据
pSliceD=m_volumeData
pSliceB=tempSlice
pSliceA=tempSlice
int tt,tt1
//扫描4层切片的顺序
/*
-----------------------D |
-----------------------B |
-----------------------C |
-----------------------A |
V
*/
//marching cubes 算法开始实行 ？第一次循环时，只读入一个切片？
for(i=0i<=(sliceNumber)i++)
{
pSliceC=pSliceA
pSliceA=pSliceB
pSliceB=pSliceD
if(i>=sliceNumber-2)
{
pSliceD=tempSlice
}
else
{
pSliceD+=imageSize
}
for(j=0j<imageHeight-1++j)
for(k=0k<imageWidth-1++k)
/* for(j=10j<imageHeight-5j++)//调试用
for(k=10k<imageWidth-5k++)*/
{
//得到八个顶点的灰度值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=0w<8w++)
{
if((cubegrid[w]>lowThreshold)&&(cubegrid[w]<highThreshold))//需要画的点
{
cubeType|=(1<<w)
}
}//end for计算cube的类型
if(cubeType==0||cubeType==255)
{
continue
}
for(w=0w<12w++) //初始化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=0w<12w++)
{
if(g_EdgeTable[cubeType]&(1<<w)) //what …..
{
//根据g_edgeTable[256]对应值判断cube的那一条边与等值面有交点
switch(w)
{
case 0:
offset=j*imageWidth+k
edgeGroup=bottomXEdge
break
case 1:
offset=j*imageWidth+k+1
edgeGroup=bottomYEdge
break
case 2:
offset=(j+1)*imageWidth+k
edgeGroup=bottomXEdge
break
case 3:
offset=j*imageWidth+k
edgeGroup=bottomYEdge
break
case 4:
offset=j*imageWidth+k
edgeGroup=topXEdge
break
case 5:
offset=j*imageWidth+k+1
edgeGroup=topYEdge
break
case 6:
offset=(j+1)*imageWidth+k
edgeGroup=topXEdge
break
case 7:
offset=j*imageWidth+k
edgeGroup=topYEdge
break
case 8:
offset=j*imageWidth+k
edgeGroup=zEdge
break
case 9:
offset=j*imageWidth+k+1
edgeGroup=zEdge
break
case 10:
offset=(j+1)*imageWidth+k+1
edgeGroup=zEdge
break
case 11:
offset=(j+1)*imageWidth+k
edgeGroup=zEdge
break
}//对应switch的{。。。end for//根据g_EdgeTable对应值判断cube的那一条边与等值面有交点
//该边上的顶点是否已经在上一层中生成
if(edgeGroup[offset]==-1)
{
int index1,index2
short s1,s2,s
float x1,y1,z1,nx1,ny1,nz1
float x2,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
}
else if(s2<=lowThreshold)
{
s=lowThreshold
}
}
else if(s2<highThreshold&&s2>lowThreshold)
{
if(s1>=highThreshold)
{
s=highThreshold
}
else if(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
float factor=((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++
} //end if(edgeGroup[offset]==-1) ////
else
{
//若该点已经在上一层生成，则直接得到其索引
cellVerts[w]=edgeGroup[offset]
}
} // end对应if(g_EdgeTable[cubeType]&(1<<w)) //
} //对应for(w=0w<12w++)
//保存当前cubes 顶点和法向量
tt1=m_vNumber-vertPos
for(tt=0tt<vertPostt++)
{
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=0r<3r++)
{
triIndex[facePos][r]=cellVerts[g_TriTable[cubeType][w++]]
if(triIndex[facePos][r]<0)
{
AfxMessageBox("有问题",MB_OK,0)
}
}
facePos++
m_fNumber++
} //end 对应while (g_TriTable[cubeType][w]!=-1)
//保存面数据
tt1=m_fNumber-facePos
for(tt=0tt<facePostt++)
{
pFace[tt1+tt][0]=triIndex[tt][0]
pFace[tt1+tt][1]=triIndex[tt][1]
pFace[tt1+tt][2]=triIndex[tt][2]
}
// memcpy(pFace+3*(m_fNumber-facePos)*sizeof(long),triIndex,sizeof(int)*3*facePos)
} memcpy(bottomXEdge,topXEdge,sizeof(short)*imageSize)
memcpy(bottomYEdge,topYEdge,sizeof(short)*imageSize)
memset(topXEdge,-1,sizeof(short)*imageSize)
memset(topYEdge,-1,sizeof(short)*imageSize)
memset(zEdge,-1,sizeof(short)*imageSize)
}
delete []tempSlice
delete []bottomXEdge
delete []bottomYEdge
delete []topXEdge
delete []topYEdge
delete []zEdge
return 1
}
在OnDraw
glBegin(GL_TRIANGLES)
for(i=0i<pDoc->m_fNumberi++)
{
glNormal3fv(&(pDoc->vPointNomal[pDoc->pFace[ i][0]][3]))
glVertex3fv(&(pDoc->vPointNomal[pDoc->pFace[ i][0]][0]))
glNormal3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][1]][3]))
glVertex3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][1]][0]))
glNormal3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][2]][3]))
glVertex3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][2]][0]))
}
glEnd()
以上代码只用于理解，未测试
包围盒层次的基本思想是通过建立对象的包围盒层次来逐渐逼近对象的几何模型,从而用体旦缓世积略大而形状简单的包围盒代替复杂的几何对象参加碰撞检测,通过包围盒间的相交测试快速地排除不相交的基本几何元素对,以减少相交测试的次数.对用于碰撞检测的包围盒有以下两方面的约束:(1)简单性:包围盒应该是简单的几何体,至少应该比被包围的几何对象简单.简单性不仅表现为几何形状简单、易于计算,而且包括相交测试算法的快速简单.(2)紧密性:包围盒应该尽可能地贴近被包围的几何对象.紧密性可以用包围盒B与被包围模肢对象G间的Hausdorff距离τ来衡(τ=maxb∈Bming∈Gdist(b,g)).τ越小,紧密性越好.紧密性哪握直接关系到需要进行相交测试的包围盒的数目.传统的包围盒类型有沿坐标轴的包围盒AABB(axis-aligned bounding boxes)[1,2]和球形包围盒(spheres).一个给定对象的AABB被定义为包含该对象且边平行于坐标轴的最小的正六面体,球形包围盒被定义为包含该对象的最小的球体.这两类包围盒的相交测试都是十分简单的,但它们的紧密性相对较差.近年来,比较著名的包围盒是带方向的包围盒OBB(oriented bounding box)[4].一个给定对象的OBB被定
义为包含该对象且相对于坐标轴方向任意的最小的正六面体.OBB最大的特点是其方向的任意
性,这使得它可以根据被包围对象的形状特点尽可能紧密地包围对象,但同时也使得它的相交测试
变得复杂.紧密性最好的包围盒应该是对象的凸包(convex hull).凸包的定义保证它是包围对象的
最小的凸多面体,但凸包自身的计算复杂性及其相交测试的困难使其难以用于碰撞检测.
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）