基于椭球约束的径向基函数隐式曲面重建

韩 燮， 武敬民， 韩慧妍， 李定主， 孙福盛

(1.中北大学计算机与控制工程学院，山西 太原 030051；2.中国兵工集团207所，山西 太原 030051)

Han Xie1, Wu Jingmin1, Han Huiyan1, Li Dingzhu2, Sun Fusheng1

(1.School of Computer and Control Engineering,North University of China,Taiyuan Shanxi 030051,China;2.207,China Ordnance Group,Taiyuan Shanxi 030051,China)

在计算机仿真中，对散乱数据点进行曲面拟合是一项基本任务。点云表面重建的常用方法有显式重建和隐式重建两种：显式方式通常描述表面的位置，可表示为映射23S:R→R；隐式方式通常用标量函数的等值面来描述，可表示为映射3F:R→R，其中显式曲面一般又可分为参数曲面和三角化表面。三角化方法是通过点云的Delaunay三角剖分[1-3]或者Voronoi图来构造得到的拓扑关系，该方法的优点是空间和时间复杂度仅取决于拟合的点数，能够得到高质量的三维网格曲面。三角化方法的缺点也很明显，第一，点云数目增大时，网格化的效率会降低；第二，点云噪声会对剖分后的效果造成严重的不保形；第三，对采样密度很敏感等。参数曲面在绘制、细分等方面表现突出。隐式曲面的优势有以下几个方面：第一、很容易辨别点在曲面的内侧或者外侧，有利于解决碰撞检测问题；第二、拟合时对点云中的噪声点不敏感，表面法向量易计算；第三、可以用简单的几何体混合成表面复杂光滑的物体等。由于隐式曲面具有的优点，尤其是简单的径向基函数，使得隐式曲面在很多领域广泛使用。

1982年Franke[4]第一次提出对散乱点云进行径向基函数插值拟合，该方法稳定且精确，但因是全局支撑，对于大规模点云拟合时，效率低下且易变形。1995年Savchebko等[5]对Franke的方法加入了一次多项式约束，但效果并未有较大改善反而计算量加大。Carr等[6]先计算出表面法向量，依据法向量约束来完成曲面重建，但计算量大，且方向不易辨别。

本文综合了径向基函数理论和椭球理论[7-9]，利用带有椭球约束的径向基函数对点云进行拟合。拟合点云时直接对散乱点云进行拟合，不需要计算法向量等信息。不过本方法较适合拟合小规模来自封闭物体的点云，否则拟合后的模型会出现冗余部分，造成拟合结果不保特征且效率低下。此时需将点云进行适当分割并同时拟合分割的点云，然后对它们进行光滑混合[10-12]。

1 径向基函数

Blinn[13]提出元球模型，该技术采用具有等势场值的点集来定义曲面。元球模型被定义为be-ar2，其中r是从空间特定点到等势场的距离。这种模型有很多变形,一般情况下，影响场可用关于距离r的函数Φ(r)。对于Blinn模型中的Φ(r)=e-ar2是高斯径向基函数，其他的基函数主要有以下几种：

(1)Φ(r)=r(线性函数)；

(2)Φ(r)=r3(三次函数)；

(3)Φ(r)=r2log(r) (薄板样条)；

(4)Φ(r)=(1+r2)1/2(多二次函数)；

(5)Φ(r)=(1+r2)-1/2(逆多二次函数)。

通常径向基函数隐式曲面的公式如下：

其中，P=(x,y,z)表示任意点，Pi表示采样点，λi表示对应每一个径向基函数的权重，是欧几里得距离，是一正定基函数，Ρ(Ρ)是一个多项式。多项式Ρ(Ρ)可看作一种约束，根据经验选择多项式时一般不超过四次，因四次以上，很难控制其拟合结果，本文选择二次椭球。

2 基于椭球约束的RBF隐式曲面拟合

2.1 RBF隐式曲面拟合

可将式(1)分解为两部分，并将这两部分映射到两个空间上：

令V1代表基函数空间：

令V2代表多项式部分的空间：

式中Ρj代表多项式函数，为使能量函数最小需让其满足正交条件即：

在空间V1和V2中，分别假设

则式(1)可写成下面的形式：

在隐式曲面上的点，均满足F(Ρk)=0，对于k,k=1,2,…,n，有下式：

根据正交条件可将式(5)写成矩阵的形式：

其中矩阵A=(aij)n×n且aij=Φ，矩阵B=(bij)m×n且bij=Ρj(Ρi)。

假设矩阵M和ξ分别为：

利用式(1)对点云拟合现在等价于下面的线性系统：

为避免平凡解的存在，需对ξ稍作限制。一般限制形式如下：

式中D是半正定矩阵。通常情况下，可能没有ξ同时满足式(1)和式(10)。所以，用式(11)更能有效地解决这种问题：

式(11)中的ξ必须满足式(10)。

对于式(10)的解可以通过求特征值来解出：

当D是一个单位矩阵时，隐式曲面的解是式(13)的特征值系统中最小特征值对应的单位特征向量0ξ。

2.2 带有椭球约束的RBF隐式曲面拟合

通常抛物面的公式如式(14)所示：

式中Ρ代表三维点数据。

令：

不管抛物面如何平移、旋转，I、J、K的值是不变的，当且仅当J＞0且I、K＞0时，式(14)才表示为椭球。

使用J＞0且I、K＞0来保证式(14)为椭球时，会反复利用非线性约束来进行迭代优化等，使得计算复杂、效率低，需选用一种简单而直接的方法。文献[14]指出只需满足式(18)即可达到式(14)为椭球的目的。

对于式(18)从几何意义来说，它表示短轴的长度至少要大于等于长轴的一半，即24J-I必须是大于零的数。

对式(14)进行调整：

由此可以得出式(6)中代表多项式部分的空间v2。

对应过点Ρi(xi,yi,zi),i=1,2,…n的式(6)中对应的矩阵B为：

式中

求带椭球约束的隐式曲面方程实际上又变成求在24J-I＞0条件下的线性系统，如下式：

通过式(15)、(16)和(21)中的24J-I=1可以求出一个矩阵C0：

使得：

再定义一个矩阵C：

条件4J-I2=1可以转化为：

式中的O11,O12,O21分别为(n-6)×(n-6),(n-6)×6,6×(n-6)大小的零矩阵。现将式(7)、(8)式的矩阵重新分块：

式中，M11,M12,M21分别是大小为(n-6)×(n-6),(n-6)×6,6×(n-6)，α,β的大小分别为(n-6)×1和6×1，式(21)现可简化为：

当M11是非奇异时，式(26)的求解等价于求解下式中的β：

式中D=M22-MT12M11-1M12=-M12TMM11MM12,因M22是零矩阵，可以直接求出β，α。但当D是非奇异矩阵时，式(27)无解，此时通常会选用最小化方法来求解。首先矩阵D是对称的、正定的，利用下式求解：

带椭球约束径向基函数方程的求解变换为常见的特征值系统问题：

对一般的径向基函数，矩阵D不一定必须为正定的。不是正定时，需用至少是半正定的TDD来代替D，方程的解是最小非负特征值对应的特征向量，但解不能保证唯一。

3 隐式曲面光滑拼接

本文选用文献[10]提出的可自行控制参数的光滑绝对函数来对分块后的重构模型进行光滑拼接。

定义1.表示当x大于等于0时=x，当x小于0时=-x，引进如下绝对值函数：

对于光滑拼接我们只需保证光滑函数达到平方次光滑即可。式(30)的一阶和二阶光滑函数的显式形式为：

式(30)可转化为另一种等价形式：

式中的Gn(x)为：

依据式(33)可以计算出和：

定义2.对于δ＞0并且n＞0，我们定义：

简明地说就是，当＞δ时,δ=,事实上，当≥δ时有≥n。故：

由此可得出光滑拼接的函数为：

由以上定义得出由n,δ共同控制的光滑函数，要根据具体需要来选择参数的值。一般情况选择n=2,因为二次光滑效率高且光滑效果不错。对于δ，对两个半径大小为0.5的圆球进行实验，根据效果来确定，如图1所示。

4 实验与分析

为验证本算法的效果和效率，与文献[1-2]、[4-5]进行了对比实验。硬件实验环境为：AMD X4 740，4 G内存，软件环境为Windows8 64位操作系统、VS2008 C++、OpenGL,Matlab2013等。

下面将对femur、head两种点云进行五种方法的重建工作，并对其重建效果进行分析，详见图2、3中的(b)、(c)、(d)、(e)、(f)分别对应文献[1-2]、[4-5]和本文方法的效果图。

图2和表1表明各种算法对femur点云的重建效果，算法[1]、[4]、[5]三种在效果上、效率上都没有优势，算法[2]虽能大体上保特征，但效率却远低于本算法。

图1 光滑函数操作

图3和表1表明各种算法对head点云的重建效果，算法[1]、[4]、[5]明显不保特征，算法[2]基本保特征，但五官信息、喉结特征没有保留，本算法虽然出现肿胀部分，但细节特征依然保留。

图2 femur点云拟合

总体来说，算法[1]是对点云进行三角剖分，重建模型外表比较粗糙。算法[2]对点云进行Loop细分重建，重建模型表面虽然光滑了很多，但却丢失了细节特征。算法[4]利用基于径向基函数隐式曲面重建，该方法虽然重建模型表面光滑连续，但模型不饱和，有失真现象。算法[5]利用平面约束的径向基函数隐式曲面重建，相对于算法[4]而言，重建模型稍有饱和，较精确地重建出模型，但有一个明显的平面，使得整个模型重建效果不好。

对于本文提出的椭球约束径向基函数隐式曲面重建方法。可以看做先对点云进行基于径向基函数方法重建，然后对整个模型曲面进行椭球曲面约束。相当于一个径向基函数隐式曲面和一个椭球隐式曲面的线性光滑混合。利用椭球的几何特性对径向基函数隐式曲面进行饱和处理，从而使得整个模型饱和精确。

实验表明基于椭球约束的径向基函数方法重建后的模型比较精确，且圆润饱和。但对于有尖锐特征的点云重建时(如head)，因为椭球的作用会形

表1 实验数据对比

图3 head点云拟合

成肿胀部分。对于此缺陷，需对点云进行简单重叠分割。首先观察head点云，在耳朵处(y=1.5)出现肿胀，需在此进行线性分割，分割为上下两部分。这两部分需要有重叠，以便拼接处理，所以将y=1.0以上部分分割为上半部分，将y=2.0以下部分分割为下部分。如图4所示，图中(a)表示上半部分，(b)表示下半部分。如图5所示，将带有尖锐特征点的云分割后，会使得尖锐特征性变小，然会对分割后的局部点云重建并进行椭球约束，从而在不出现肿胀现象的情况下，局部点云生成的模型更具有饱和性，保特征性。最后通过对局部模型进行光滑拼接。共消耗时间4.25 s，效率明显提高。

图4 head点云分割示意图

图5 光滑拼接效果

5 结 论

本文提出的基于椭球约束的径向基函数隐式曲面重建算法，直接从散乱点获取重建信息，重建结果精确。但对大规模点云集拟合时，需解大规模矩阵，效率低下甚至拟合模型会产生冗余部分。要通过对点云进行简单的分割再进行拟合光滑拼接，不仅效果好而且是并行拟合，这样时间效率会飞速提高。实验效果表明：本算法相对于上面提到的效果保特征性好而且效率很高。

[1]Kuo C C,Yau H T.A Delaunay-based region-growing approach to surface reconstruction from unorganized points [J].Computer-Aided Design,2005,37(8):825-835.

[2]Cheng Fuhua,Fan Fengtao,Lai Shuhua,Huang Conglin,Wang Jiaxi,Yong Junhai.Loop subdivision surface based progressive interpolation [J].Journal of Computer Science and Technology,2009,24(1): 39-46.

[3]陈慧群,黎景炎.基于变形网格法的高密点云曲面重建[J].工程图学学报,2011,32(2): 64-67.

[4]Franke R.Scattered data interpolation: Tests of some methods [J].Mathematics of computation,1982,38(157):181-200.

[5]Savchenko V V,Alexander A P,Okunev OG,Kunii T L.Function representation of solids reconstructed from scattered surface points and contours [J].Computer Graphics Forum,1995,14(2): 181-188.

[6]Carr J C,Beatson R K,Cherrie J B,Mitchell T J,Fright W R,McCallum B C,Evans T R.Reconstruction and represe-ntation of 3D objects with radial basis functions [C]//Proceedings of the 28th annual conference on Computer graphics and interactive techniques.ACM,2001: 67-76.

[7]Klein P P.On the Ellipsoid and Plane Intersection Equation [J].Applied Mathematics,2012,3(11):1634-1640.

[8]Beck A,Sabach S.An improved ellipsoid method for solving convex differentiable optimization problems [J].Operations Research Letters,2012,40(6): 541-545.

[9]Reiner T,Mückl G,Dachsbacher C.Interactive modeling of implicit surfaces using a direct visualization approach with signed distance functions [J].Computers & Graphics,2011,35(3): 596-603.

[10]Li Qingde.Smooth piecewise polynomial blending operations for implicit shapes [C]// Computer Graphics Forum,2007,26(2): 157-171.

[11]Li Qingde,Tian Jie.Partial shape-preserving splines [J].Computer-Aided Design,2011,43(4): 394-409.

[12]李凌丰,谭建荣,赵海霞.Metaball 重叠区域作用效果混合[J].中国图象图形学报,2006,11(5): 695-699.

[13]Blinn J F.A generalization of algebraic surface drawing [J].ACM Transactions on Graphics (TOG),1982,1(3):235-256.

[14]Li Qingde.Registration techniques for computer assisted surgery [D].UK: University of Hull,2001.