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,当满足
(qi−C)(qj−C)≤0 (公式 3-1 等值点判定条件式)则 Mi 和 M j 两点间庇
佑等值点 Mo ,另设等值点 Mo 的坐标为(x0,y0, z0),由 Mi 和 M j 两点根据
线性插值理论可得:
公式 3-2 等值点坐标公式
其中 k=(qi−C)(qj−C)≤0 根据等值点判定条件式 3-1,和等值点坐标公式
(3-2)可以按结构离散信息对网格棱边进行搜索判断,从而求出指定域中结构
体所有等值点。求出等值点以后,就可以将这些等值点连接成三角形或多边形形
成等值面的一部份。
1.1.6 等值面的法向量的计算
为了利用图形硬件显示等值面图像,必须给出三角面片个等值面的法向,选
择 适当的关照模型进行适当的光照计算,生成真实感图形。对于等值面上的每
一点,其沿面的切线方向的梯度分量应该是零,因此沿该点的梯度矢量方向也就
代表了等值面在该点的法向。而且等值面往往是具有不同密度物质的分界面,因
而其梯度矢量值不为零值即:
g(x,y,z)= ∇f(x,y,z)
公式 3-3 梯度矢量公式
直接计算三角面片的法向是费时的,而且,为了消除各三角面片之间的明暗
度的不连续变化,只要给出三角面片各顶点处的法向并采用哥罗德(Gouraud)模
型绘制各三角面片就行了。这里我们采用中心插分方法来计算各体素各角点的梯
度。在三角形的情况下,计算出每一个三角形面片的法向量,然后用三角面的法
向量求得每个顶点的法向量,最后用三角形三个顶点的三个法向量插值求出三角
形面上某一点的法向量。对于等值面来说有简单的方法计算顶点的法向量。考虑
等高线的情形。等高线的梯度方向与等高线的切线垂直,即可以用梯度代替等高
线的垂线。推广到三维情况,等值面的梯度方向就是等值面的法向方向。