Marching Cube 算法原理.doc

Marching Cube 算法原理1.1.1 Marching Cube 算法概述面绘制法则是根据设定的阈值，从体数据中提取出表面的三角面片集，再用光照模型对三角面片进行渲染，形成三维图像。面绘制法主要分为基于断层轮廓线的方法和基于体素的方法。基于断层轮廓线的方法是先在不同的断层上提取出感兴趣区的轮廓线，然后在相邻的断层的轮廓线间构造出三角面片的方法，这在同一断层上有多个轮廓线时会产生模糊性，上下两层的轮廓线不易对应。用户干预可以避免一定的模糊性，但是这样大大增加了操作的复杂性。因此不被广泛采纳使用。基于体素的方法以移动立方体法（Marching Cube，MC）为代表。Marching Cubes算法是面显示算法中的经典算法，它也被称为“等值面提取”（Iso-surface Extraction）。本质是将一系列两维的切片数据看作是一个三维的数据场，从中将具有某种域值的物质抽取出来，以某种拓扑形式连接成三角面片。算法的基本原理MC算法的基本思想是逐个处理体数据场中的各个体元,并根据体元各个顶点的值来决定该体元内部等值面的构造形式算法实现过程中,体元内等值面构造要经过以下两个主要计算:1、体元中三角面片逼近等值面的计算;2、三角面片各顶点法向量的计算。1.1.2预备知识介绍（体素模型和等值面介绍）1、体素模型的介绍体素一般有两种定义：一种与二维图像中像素定义相类似。直接把体数据中的采样点作为体素，另一种是把八个相邻的采样点包含的区域定义为体素。在三维空间某一个区域内进行采样，若采样点在x,y,z,三个方向上分布是均匀的。采样间距分别为x,y,z，则体数据可以用三维数字矩阵来表示。每八个相临的采样点相临的立方体区域就定义为一个体素。而这八个采样点称为该体素的角点。他们的坐标分别为：（i, j, k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), (i,j,k+1), (i,j,k+1), (i+1.j+k+1), (i,j+1,k+1) 和（i+1,j+1,k+1）如图-1所示图-1 移动立方体的体素对于体素内任一点P6(x, y,z)，其物理坐标可以转换为图像坐标 i6 , j6 ,k6， 其中i6=x/x, j6=y/y, k6=z/z.当把方向无关的三个线性插值作为体素模型时，其值可以表示为2、等值面（Iso-Surface）介绍在面重建算法中以重建等值面这一类算法最为经典。我们进行表面重建的目的就是用分割提取出的区域构建出对应组织或器官的三维几何模型。等值面的构造就是从体数据中恢复物体三维几何模型的常用方法之一。如果我们把体数据看成是某个空间区域内关于某种物理属性的采样集合，非采样点上的值用邻近采样点插值来估计，则该空间区域内所有具有某一个相同值的点的集合将定义一个或多个曲面，称之为等值面。因为不同的物质具有不同的物理属性，因此可以选定适当的值来定义等值面，该等值面表示不同物质的交界。也就是说，一个用适当值定义的等值面可以代表某种物质的表面。等值面是空间中所有具有某个相同值的点的集合，它可以表示成，(x, y,z), f (x, y,z) =c其中C为常数。并不是每个体素内都有等值面，当体素内角点都大于C或者都小于C时其内部不存在等值面只有那些即大于C又小于C的角点的体素才含有等值面，我们称这样的体素为边界体素。等值面在一个边界体素内的部分称为该体素的等值面片，等值面是一个三次曲面，它与边界体素面的交线是一条双曲线且这条双曲线仅由该面上四个角点决定。这些等值面片之间具有等值拓扑一致性，即它们可以构成连续的无孔的无悬浮面的曲面（除非在体数据的边界处）。因为对于任何两个边界共面的体素，如果等值面与他们的公共面有交线，则该交线就是两个边界体素中等值面片与公共面的交线，也就是说这两个等值面片完全吻合，所以可以认为等值面是由许多个等值面片组成的连续曲面。由于等值面是三次代数曲面，构造等值面的计算复杂，也不便于显示，而多边形的显示则非常方便，所以，等值面的三角面片拟合是常用的手段。我们本章论述的MC算法便是在边界体素中生成三角面片，以三角面片拟合成等值面。1.1.3 移动正方形法移动正方形法是一种二维算法，它是移动立方体法的依据。移动立方体法Marching Cube正是移动正方形法的三维引申发展起来的。移动正方形法也是找等值线的一种方法。首先找四个相邻的象素，编号为 1，2，3，4，如图-2所示。每个象素值有大于阈值和小于阈值两种情况，如果象素值大于阈值用代码 1 表示，用圆圈表示，如果小于阈值就用 0 表示。四个点就有 16 种组合形式，图-3列出了所有的可能组合形式。每一种形式就是等值线与正方形边之间的一种拓扑关系。图中的虚线就是等值线的路径。没有虚线的形式说明等值线不与正方形相交。以 0001 图为例，该图中左下角的象素值大于给定值，其它三个象素小于给定值，那么可以推断出等值线的一侧是圆圈代表的象素，另一侧是另外三个象素，那么等值线只能以图中虚线所示的这种方式与正方形相交。等值线与正方形边的交点坐标可以用线性插值来求得。这样当一幅图像中的所有正方形都求出了各自的一段等值线后，这些线段自然而然的就连成了一个闭合的等值线了。移动正方形算法如下： 在一幅图像中求出所有四个相邻象素点构成的正方形。 判断四个象素值与阈值的关系，生成 0101 的代码。 由上步生成的代码按照图-4所示的关系求出等值线与四个象素点间的拓扑关系。 由拓扑关系，用线性插值法求出等值线与正方形边的交点。 顺序连接等值线段就得到等值线了。图-2 移动正方形等值面的几种情况1.1.4移动立方体（Marching Cubes）算法MC算法的基本假设是沿着立方体的边的数据场是呈连续线形变化的，也就是说如果一条边的两个顶点分别大于小于等值面的值，在该边上庇佑且仅有一点是这条边与等值面的交点。确定立方体体素等值面的分布是该算法的基础。这里我们将理论与重建示例相结合使我们对MC算法进行更深一步的了解。首先我们将经过处理后的图片切片数据可以看做是一些网格点组成的，这些点代表了密度值。图-3 CT、MR图像的灰度单位网格每次读出两张切片，形成一层（Layer）图-4 连续读出两张图片所构成的Layer两张切片上下相对应的八个点构成一个Cube，也叫Cell,Voxel等。由相邻层上的各4 个像素组成立方体的8 个顶点，这 8个像素构成一个立方体。我们把这个立方体叫做体素。为了确定体元中等值面的剖分方式,因此所求等值面要的一个门限值,然后对体元的八个顶点进行分类,以判定顶点是位于等值面之内还是位于等值面之外;再根据顶点分类结果确定等值面的剖分模式。顶点分类规则为(1)如果顶点的数据值大于等值面的值,则定义该顶点位于等值面之内,记为“1”;顶点密度值域值，设为Outside（1）(2)如果顶点的数据值小于等值面的值,则定义该顶点位于等值面之外,记为“0”。顶点密度值域值，Inside （0）图-6 8个顶点内点外点在Cube中的表示那么这个等值面必定与三维图像的某些体素相交与其它体素不相交。对于某一个体素来说要么与等值面相交要么在等值面的某一侧，不与等值面相交。MC 方法的思想就是找出所有与等值面相交的体素，再分别找出每个体素与等值面相交的交面，这些交面连在一起就是所求的等值面。首先要确定等值面通过那些体素，然后在确定等值面如何与体素相交。当一个体素中一些象素的值大于阈值，而另一些象素小于阈值，那么等值面必然通过这个体素，一个体素的 8 个象素的值全都小于阈值或者全都大于阈值的话，那么该体素不与等值面相交，等值面不通过该体素。当一个体素与等值面相交的话，必然有一些象素值大于阈值，一些小于阈值。每个象素有两种状态，要么大于阈值，要么小于阈值确定包含等值面的体元。对于8个角点都为1或者都为0的体素，它属于“0”号结构没有等值面穿过该体素。当有1个角点标记为1时为1号结构我们用1个三角面片代表等值面它将该角点与其它七角点分成两部分。对于其余几种构型将产生多个三角面片。flag( i,j, k ) = 0 ( 1 )256种情况。因此共有256 种组合状态。每一种组合都对应一种等值面与体素相交的情况。因为 8 个点有对称关系，256 种组合经可简化为如图-6的15 种情况。每一种关系对应等值面如何与体素相交，知道了等值面如何与体素相交后就可以求得等值面与立方体边的交点，这些交点形成的面片就是等值面的一部分。当把所有与等值面相交的体素都找到，并求出相应的相交面后，等值面也就求出来了。黑点标记为（1）的角点。图-6 256中组合简化后的15种等值面与体素相交构成三角面片的情况应用上面的15种构型情况的具体方法是：对于每个体素根据它的索引在“构型-三角抛分形式然后再根据相应索引项中的旋转参数具体确定最终的三角剖分。从而根据上面的状态表当前体元素属于哪种情况以及等值面与哪一条边相交1.1.5 求等值面与体素边界的交点在确定立方体的三角剖分模式后，就要计算三角面片顶点位置。当三维离散数据场的密度较高时，即当体素很小时。可以假定函数沿体素边界呈线性变化，这就是MC算法的基本假设，因此根据这一基本假设可以用线性插值计算等值面与体素边界的交点。MC算法的基本假设是沿着体元的棱边数据场呈线性变化，即如果一条棱边的两个数据场值分布大于或小于等值面值，则这条边上由且仅有一点是等值面与该边界的交点。在已知结构空间网格结构点上某一物理量值的前提下。设对任意网格结构点M的直角坐标系下。以M（x,y,z,q）表示，其中x,y,z分别为M点的直角坐标值，q为该结点的物理量值（实际上等值面的等值点就是结构物体内所有具有相同量值得点。假定在结构体内等值点分布在离散网格的边棱上。空间等值点的判断;根据以上的假定，任取一离散网格边棱。设棱边上两结点分别为：Mi( xi, yi,zi,qi) 和Mj( xj, yj, zj,qj) ，取量值的等值为C，当满足(qiC)（qjC）0 （公式3-1 等值点判定条件式）则Mi 和M j 两点间庇佑等值点Mo ，另设等值点Mo 的坐标为（x0,y0, z0），由Mi 和M j 两点根据线性插值理论可得：公式 3-2 等值点坐标公式其中k=(qiC)（qjC）0 根据等值点判定条件式3-1，和等值点坐标公式（3-2）可以按结构离散信息对网格棱边进行搜索判断，从而求出指定域中结构体所有等值点。求出等值点以后，就可以将这些等值点连接成三角形或多边形形成等值面的一部份。1.1.6等值面的法向量的计算为了利用图形硬件显示等值面图像，必须给出三角面片个等值面的法向，选择 适当的关照模型进行适当的光照计算，生成真实感图形。对于等值面上的每一点，其沿面的切线方向的梯度分量应该是零，因此沿该点的梯度矢量方向也就代表了等值面在该点的法向。而且等值面往往是具有不同密度物质的分界面，因而其梯度矢量值不为零值即：g(x,y,z)= f(x,y,z)公式 3-3 梯度矢量公式直接计算三角面片的法向是费时的，而且，为了消除各三角面片之间的明暗度的不连续变化，只要给出三角面片各顶点处的法向并采用哥罗德(Gouraud)模型绘制各三角面片就行了。这里我们采用中心插分方法来计算各体素各角点的梯度。在三角形的情况下，计算出每一个三角形面片的法向量，然后用三角面的法向量求得每个顶点的法向量，最后用三角形三个顶点的三个法向量插值求出三角形面上某一点的法向量。对于等值面来说有简单的方法计算顶点的法向量。考虑等高线的情形。等高线的梯度方向与等高线的切线垂直，即可以用梯度代替等高线的垂线。推广到三维情况，等值面的梯度方向就是等值面的法向方向。公式3-4 插分梯度计算公式其中x,y,z都是体素长。1.1.7 MC 算法步骤 根据对称关系构建一个 256 种相交关系的索引表。该表指明等值面与体素的那条边相交。 提取相邻两层图片中相邻的 8 个象素，构成一个体素并把这 8 个象素编号。如图-7 所示。 根据每个象素与阈值的比较确定该象素是 1 还是 0。 把这 8 个象素构成的 01 串组成一个 8 位的索引值。如图 38 所示。 用索引值在上边的索引表里查找对应关系，并求出与立方体每条边的点。 用交点构成三角形面片或者是多边形面片。 遍历三维图像的所有体素，重复执行到。1.1.8MC算法分析根据上面对 MC 方法的介绍，可以认为 MC 方法是求等值面，进行三维表面重建的好方法，但 MC 方法还存在以下问题。 MC 方法构造的三角面片是对等值面的近似表示。首先，等值面与体素边界交点是用线性插值求得的。当体数据密度高时，体素很小，MC 方法中的三角面片与等值面比较接近，重建的效果好。当体数据比较稀疏时，这种重建方法将会产生较大的误差。 MC 方法中，在体素的一个面上，如果顶点值为 1 的象素和顶点值为 0的象素分别位于对角线的两端，那么有两种可能的连接方式，因此存在连接方式的二义性，如图-8 所示。这样的面称为二义面，包含 1 个以上的二义性面的体素就称为具有二义性的体素。在图-6 中的 15 中情况中，第 3，6，7，10，12，13 这几种情况都具有二义性。这种二义性如不解决，将造成等值面连接上的错误，从而形成空洞，如图-9 所示。等值面应该是连续的曲面，MC 中相邻体素重建出来的三角面片应该彼此相接，因为重建出来的三角面片有错误，相邻的三角面片不相接了，这样等值面上就出现了孔洞。图-8 MC算法中的二义性图-9 由于二义性引起所形成的空洞 在图 37 中，可以看出 7 号以后重建出来的三角面数多于 3 个，每个体素上都有几个三角形，重建出来的三角面片很多，也很零碎。如果 MC 方法的输入数据是灰度图，而不是分割后的二值图的话，那么所有体素都会重建，所有存在象素值小于和大于阈值的体素都重建出三角面片，就会重建出多余的等值面。而我们要的结果一般是一个物体的等值面。如果是用分割后的数据进行重建，那么所有的象素是二值的，要么是 0，要么是 255，MC 法求的等值面就没有意义了。1.1.9MC算法描述for(k=1;kNz;k+)读入k-1,k,k+1