Marching Cube拓扑二义性解决方法改进

Marching Cube拓扑二义性解决方法改进1孙伟，张彩明,肖敏山东大学计算机科学与技术学院（250061）E-mail：sun_wei@摘要：在等值面的三角片逼近中，逼近等值面拓扑结构的正确性十分关键，直接关系到所得到的等值面在应用中的指导意义。

原始的Marching Cube方法是假设沿体元的边上数据由该边两顶点线性插值得到，渐近线消除面的二义性方法是基于在面上的双线性插值进行计算，后来Nielson在体元中通过计算三线性插值来解决三维空间中内部拓扑二义性的问题。

以上算法由于完全基于对采样数据线性变化的假设，所选择的采样点无法保证在非线性变化的数据场逼近的等值面具有正确的拓扑结构。

本文在原有算法的基础上，提出了新的解决面上二义性方法，假设原数据在二义性面上为一二次函数，通过引入周围立方体单元的信息来计算该二次函数，并用来解决面上的二义性。

文中用实例对新方法和已有方法的逼近结果作了比较。

关键词：Marching Cube，等值面逼近，拓扑二义性，临界点1 引言等值面技术在可视化中应用很广，许多标量场中的可视化问题都归纳为等值面的抽取和绘制，如医学图像中的三维重建；有限元计算中的标量场分析，分子化学中的分子表面显示；地质中矿藏分布的构造等等。

在散乱点的曲面重建[1，2]和关键帧动画中，也常常提取零等值面作为散乱点数据的重建曲面和中间帧。

近年来，等值面逼近问题得到了极为广泛的关注，并取得了许多重要的成果，其中，MC(Marching Cube)方法以其简单和高效的特点，成为应用最广泛的等值面逼近方法。

利用三角片对等值面进行逼近首先要保证有三角片构成的等值面与原曲面有相同的拓扑结构。

通常的逼近方法在等值面上提取有限个采样点，并将这些采样点进行三角化，以得到得三角片组合曲面作为逼近等值面。

因而采样点的选择是逼近等值面拓扑结构的正确性的关键，尤其是在医学图像可视化等处理中，逼近等值面拓扑结构的正确行更为重要，否则将影响对临床的指导意义。

现在的MC逼近方法，在处理二义性面时，由于基于在采样数据沿体元(voxel)边的变化是线性的假设，采用双曲线方法进行处理，这种方法导致了对许多二义性面的错误处理。

总结现在存在的MC方法的基础上，本文在处理二义性面时，通过计算二义性面四个顶点的梯度，引入周围体元的信息，较好的解决二义性，得到更可靠的结果。

本文的方法在增加计算量不大的情况下，保证了等值面具有正确的拓扑结构。

文中使用构造的二次函数实例对新的1本课题得到高等学校博士学科点专项科研基金（项目编号：20020422030）资助方法和已有的方法的逼近结果作了比较。

2 引入周围体元信息消除面的二义性文章中按照原始算法中的立方体体元顶点分类规则：对于大于等于等值面值的立方体顶点，定义该顶点位于等值面之外，称为正顶点，记为“+”；对于小于等值面值的立方体顶点，定义该顶点位于等值面之内，成为负顶点，记为“-”。

如果一个面上位于等值面内和在等值面外的顶点分别分布在对角线的两端，则称这个面为二义性面。

包含二义性面的立方体，即为二义性立方体。

本文中仍然假设等值面与立方体的边只有一个交点。

对于存在至少一个二义性面的体元，因为线性函数与体元面的交点为一条直线，不能将二义性面的正顶点和负顶点分开，所以面上的采样数据对应的原始函数至少为二次函数。

本文中设计如下两个二次函数对面上解决二义性的方法测试其判断结果正确性，函数如下：212222(,,)44()5(,,)4(1)2()2(3)1F x y z y x z F x y z y x z x z =+−−=−+−−+−+（1） （2）以上两个原函数当三维坐标均在实数域上取0到3时，函数的零等值面图形如下图：图1函数F1(左) 和F2(右)的零等值面对三个坐标分别在整数0到3上取样得到的点组成一个4*4*4的网格，其中中间立方体和其下面的立方体的各顶点的值在下图2中标出，（b ）（a ）图2 F 1和F 2中心以及中心下方的立方体各个顶点取值从图二中可以明显的看出，两个原始函数的中间立方体顶点取值是一样的，但从图一可以看出，函数F 1的等值面为一个连通的曲面，函数F 2的等值面为两个分离的双曲面片。

这就说明对于上面的这种非线性插值得到的原始函数，不可能使用一个立方体内的顶点信息得到正确的曲面拓扑结构，必须使用邻近体元的顶点信息，本文中通过计算立方体中二义性面的顶点梯度引入周围立方体的信息，称为梯度法。

由于梯度向量的方向为等值面的法方向，其大小为函数变化快慢的度量，所以估计二义性面的梯度，能够估计得出原始函数在二义性面上的变化情况，我们需要的仅仅是梯度在平面上的分量。

例如，如果我们计算梯度的二义性面平行于x-y 坐标轴，则我们只需要计算在x 和y 方向的梯度信息，用来估计原函数在这两个方向的变化情况。

由于在后面的光照明模型中仍然用到各个顶点处的梯度信息，即法向量，所以在此引入梯度来消除二义性，并不增加计算量。

本文假设立方体为单位立方体，各个边的长度为1，对于各个顶点处梯度值的计算，采用中心差分的方法，即对于图三中的中间的面，假设其平行于x －y 坐标轴我们采用下面的方法来求得顶点(i,j,k)在x,y 方向的梯度分量：x yg =(f(i+1,j,k)-f(i-1,j,k))/2g =(f(i,j+1,k)-(i,j-1,k))/2⎧⎪⎨⎪⎩（3）其中f(i,j,k)为点(i,j,k)处的函数值。

（i-1,j+1,k ） （i,j+1,k ） （i+1,j+1,k ） （i-1,j,k ） （i,j,k ） （i+1,j,k ） (i-1,j-1,k)(i,j-1,k)（i+1,j-1,k ）图3如上所述，由于两个变量的线性函数不能逼近一个二义性面，所以逼近二义性面的函数的最低次数为二次。

本文中就使用一个二次函数来逼近二义性面，通过对二义性面的四个顶点进行插值，对四个顶点处的梯度分量进行逼近来得到该函数。

由于四个顶点处的梯度分量是使用了图三中周围8个顶点计算出来，因此，二次函数也可以看作是对二义性面四个顶点进行插值，对周围8个点进行逼近的结果。

在对四个顶点处的梯度进行逼近时，选择最小平方误差作为标准。

假设所要处理的二义性面平行于x-y 轴，则不失一般性，可设二义性面在局部x-y 坐标系中的坐标为（0，0），（△x,0）,(0, △y),( △x,△y)，如图四所示，四个顶点分别标识为ll ，lr ，ul ，ur ，在顶点处x 方向的梯度分量表示为▽fx ，y 方向的梯度分量表示为▽fy 。

为了更方便地说明我们的梯度法，本文做从x-y 到t-u 的坐标变换，如图四所示， 其中1()21()2x t u x y t u y ⎧=−+∆⎪⎪⎨⎪=−+∆⎪⎩（4）图4 从x-y 坐标系到t-u 坐标系的变换在变换后的坐标空间t-u 中，梯度分量分别用▽ft ，▽fu 表示，可以通过下式（5）得到在t-u 空间的梯度分量▽ft ，▽fu 。

(5)t x u t f f x y f x y f ∇∇∆∆⎛⎞⎛⎞⎛⎞=⎜⎟⎜⎟⎜⎟∇−∆∆∇⎝⎠⎝⎠⎝⎠在下文中，我们用▽fx,ll, ▽fy,ll 表示在x-y 坐标系中顶点ll 处的梯度分量，相应的， 用▽ft,ll, ▽fu,ll 表示对应的t-u 坐标系中ll 处的梯度分量。

其它点处的梯度分量也可用相同的方法表示。

由于含有两个变量的二次函数有6个未知量，本文在t-u 坐标系中表示该二次函数如（6）所示。

2212111222(,)2(6)Q t u c b t b u a t a tu a u =+++++上式中存在6个未知量，由需要满足在四个顶点处的插值可知（7），ll lr u r u l Q (-1/2,0)=f Q (0,-1/2)=f (7) Q (1/2,0) =f Q(0,1/2) =f⎧⎪⎪⎨⎪⎪⎩由方程组（7），我们可以得出方程组（8），112212112214()2(8)1()2lllr ul ur ul ll ul lr ll lr ul urf f f f c a a f f b f f b f f f f a a ⎧+++=++⎪⎪−=⎪⎨−=⎪⎪−−+=−⎪⎩通过（8）式可以解出未知系数a 11,a 22,b 1,b 2。

其中a 11 =2(f ll +f ur )-4c, a 22 =2(f lr +f ur )-4c,二次函数Q 对于四个顶点的梯度采用的是最小平方误差方法进行逼近，本文使用函数E(t,u)定义每个顶点的平方误差，其表达式如（9）所示，221111222212(,)(22(,))(22(,))(9)t u E t u b a t a u f t u b a u a t f t u =++−∇+++−∇要使四个顶点处的平方误差达到最小，必须满足下面（10）所示的约束条件，11111111112222222222121212121111(,0)(0,(0,)(,0)2222()1111(,0)(0,(0,)(,0)2222()0(10)1111(,0)(0,)(0,)(,0)22220E E E E da a a a a dc E E E E da a a a a dc E E E E a a a a ⎧∂−∂−∂∂⎪+++⎪∂∂∂∂⎪⎪∂−∂−∂∂⎪⎪++++=⎨∂∂∂∂⎪⎪∂−∂−∂∂⎪⎪+++=∂∂∂∂⎪⎪⎩通过（10），解出系数c ，a 12 如下，,,,,12,,,,11()()416(11)1()4ll lr ul ur t ll u lr u ul u ur u ll t lr r ul u ur c f f f f f f f f a f f f f ⎧=++++∇+∇−∇−∇⎪⎪⎨⎪=−∇−∇+∇+∇⎪⎩这样即可得到二义性面的逼近函数Q （t,u ）的表达式。

当二次函数Q(t,u)=0表示的曲线为抛物线，椭圆，平行直线等凸集时，可以证明，函数在二义性面的中心的取值Q(0,0)=c 与在曲线中连通的对角顶点有相同的符号，即连通与Q （0，0）符号相同的两个节点。

当Q(t,u)=0表示的曲线为双曲线时，使用Nielson [7]的渐近线方法，此时仅使用二义性面四个顶点的就可以判断连通哪两个点，方法如下，如果Q(-1/2,0)Q(1/2,0）> Q(0，-1/2)Q(0，1/2）,即f ll f ur >f lr f ul ，则将ll,ur 节点连通；否则，连通lr,ul 节点。

3 新方法的体元内部曲面逼近过程MC 方法[3]将体元每个顶点的状态进行标记并根据不同的状态组合进行分类，其中标量值大于等值面阈值的顶点标为正，否则为负。