韓 燮， 武敬民， 韓慧妍， 李定主， 孫福盛

(1.中北大學(xué)計(jì)算機(jī)與控制工程學(xué)院，山西 太原 030051；2.中國(guó)兵工集團(tuán)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)

在計(jì)算機(jī)仿真中，對(duì)散亂數(shù)據(jù)點(diǎn)進(jìn)行曲面擬合是一項(xiàng)基本任務(wù)。點(diǎn)云表面重建的常用方法有顯式重建和隱式重建兩種：顯式方式通常描述表面的位置，可表示為映射23S:R→R；隱式方式通常用標(biāo)量函數(shù)的等值面來(lái)描述，可表示為映射3F:R→R，其中顯式曲面一般又可分為參數(shù)曲面和三角化表面。三角化方法是通過(guò)點(diǎn)云的Delaunay三角剖分[1-3]或者Voronoi圖來(lái)構(gòu)造得到的拓?fù)潢P(guān)系，該方法的優(yōu)點(diǎn)是空間和時(shí)間復(fù)雜度僅取決于擬合的點(diǎn)數(shù)，能夠得到高質(zhì)量的三維網(wǎng)格曲面。三角化方法的缺點(diǎn)也很明顯，第一，點(diǎn)云數(shù)目增大時(shí)，網(wǎng)格化的效率會(huì)降低；第二，點(diǎn)云噪聲會(huì)對(duì)剖分后的效果造成嚴(yán)重的不保形；第三，對(duì)采樣密度很敏感等。參數(shù)曲面在繪制、細(xì)分等方面表現(xiàn)突出。隱式曲面的優(yōu)勢(shì)有以下幾個(gè)方面：第一、很容易辨別點(diǎn)在曲面的內(nèi)側(cè)或者外側(cè)，有利于解決碰撞檢測(cè)問(wèn)題；第二、擬合時(shí)對(duì)點(diǎn)云中的噪聲點(diǎn)不敏感，表面法向量易計(jì)算；第三、可以用簡(jiǎn)單的幾何體混合成表面復(fù)雜光滑的物體等。由于隱式曲面具有的優(yōu)點(diǎn)，尤其是簡(jiǎn)單的徑向基函數(shù)，使得隱式曲面在很多領(lǐng)域廣泛使用。

1982年Franke[4]第一次提出對(duì)散亂點(diǎn)云進(jìn)行徑向基函數(shù)插值擬合，該方法穩(wěn)定且精確，但因是全局支撐，對(duì)于大規(guī)模點(diǎn)云擬合時(shí)，效率低下且易變形。1995年Savchebko等[5]對(duì)Franke的方法加入了一次多項(xiàng)式約束，但效果并未有較大改善反而計(jì)算量加大。Carr等[6]先計(jì)算出表面法向量，依據(jù)法向量約束來(lái)完成曲面重建，但計(jì)算量大，且方向不易辨別。

本文綜合了徑向基函數(shù)理論和橢球理論[7-9]，利用帶有橢球約束的徑向基函數(shù)對(duì)點(diǎn)云進(jìn)行擬合。擬合點(diǎn)云時(shí)直接對(duì)散亂點(diǎn)云進(jìn)行擬合，不需要計(jì)算法向量等信息。不過(guò)本方法較適合擬合小規(guī)模來(lái)自封閉物體的點(diǎn)云，否則擬合后的模型會(huì)出現(xiàn)冗余部分，造成擬合結(jié)果不保特征且效率低下。此時(shí)需將點(diǎn)云進(jìn)行適當(dāng)分割并同時(shí)擬合分割的點(diǎn)云，然后對(duì)它們進(jìn)行光滑混合[10-12]。

1 徑向基函數(shù)

Blinn[13]提出元球模型，該技術(shù)采用具有等勢(shì)場(chǎng)值的點(diǎn)集來(lái)定義曲面。元球模型被定義為be-ar2，其中r是從空間特定點(diǎn)到等勢(shì)場(chǎng)的距離。這種模型有很多變形,一般情況下，影響場(chǎng)可用關(guān)于距離r的函數(shù)Φ(r)。對(duì)于Blinn模型中的Φ(r)=e-ar2是高斯徑向基函數(shù)，其他的基函數(shù)主要有以下幾種：

(1)Φ(r)=r(線性函數(shù))；

(2)Φ(r)=r3(三次函數(shù))；

(3)Φ(r)=r2log(r) (薄板樣條)；

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

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

通常徑向基函數(shù)隱式曲面的公式如下：

其中，P=(x,y,z)表示任意點(diǎn)，Pi表示采樣點(diǎn)，λi表示對(duì)應(yīng)每一個(gè)徑向基函數(shù)的權(quán)重，是歐幾里得距離，是一正定基函數(shù)，Ρ(Ρ)是一個(gè)多項(xiàng)式。多項(xiàng)式Ρ(Ρ)可看作一種約束，根據(jù)經(jīng)驗(yàn)選擇多項(xiàng)式時(shí)一般不超過(guò)四次，因四次以上，很難控制其擬合結(jié)果，本文選擇二次橢球。

2 基于橢球約束的RBF隱式曲面擬合

2.1 RBF隱式曲面擬合

可將式(1)分解為兩部分，并將這兩部分映射到兩個(gè)空間上：

令V1代表基函數(shù)空間：

令V2代表多項(xiàng)式部分的空間：

式中Ρj代表多項(xiàng)式函數(shù)，為使能量函數(shù)最小需讓其滿足正交條件即：

在空間V1和V2中，分別假設(shè)

則式(1)可寫(xiě)成下面的形式：

在隱式曲面上的點(diǎn)，均滿足F(Ρk)=0，對(duì)于k,k=1,2,…,n，有下式：

根據(jù)正交條件可將式(5)寫(xiě)成矩陣的形式：

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

假設(shè)矩陣M和ξ分別為：

利用式(1)對(duì)點(diǎn)云擬合現(xiàn)在等價(jià)于下面的線性系統(tǒng)：

為避免平凡解的存在，需對(duì)ξ稍作限制。一般限制形式如下：

式中D是半正定矩陣。通常情況下，可能沒(méi)有ξ同時(shí)滿足式(1)和式(10)。所以，用式(11)更能有效地解決這種問(wèn)題：

式(11)中的ξ必須滿足式(10)。

對(duì)于式(10)的解可以通過(guò)求特征值來(lái)解出：

當(dāng)D是一個(gè)單位矩陣時(shí)，隱式曲面的解是式(13)的特征值系統(tǒng)中最小特征值對(duì)應(yīng)的單位特征向量0ξ。

2.2 帶有橢球約束的RBF隱式曲面擬合

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

式中Ρ代表三維點(diǎn)數(shù)據(jù)。

令：

不管拋物面如何平移、旋轉(zhuǎn)，I、J、K的值是不變的，當(dāng)且僅當(dāng)J＞0且I、K＞0時(shí)，式(14)才表示為橢球。

使用J＞0且I、K＞0來(lái)保證式(14)為橢球時(shí)，會(huì)反復(fù)利用非線性約束來(lái)進(jìn)行迭代優(yōu)化等，使得計(jì)算復(fù)雜、效率低，需選用一種簡(jiǎn)單而直接的方法。文獻(xiàn)[14]指出只需滿足式(18)即可達(dá)到式(14)為橢球的目的。

對(duì)于式(18)從幾何意義來(lái)說(shuō)，它表示短軸的長(zhǎng)度至少要大于等于長(zhǎng)軸的一半，即24J-I必須是大于零的數(shù)。

對(duì)式(14)進(jìn)行調(diào)整：

由此可以得出式(6)中代表多項(xiàng)式部分的空間v2。

對(duì)應(yīng)過(guò)點(diǎn)Ρi(xi,yi,zi),i=1,2,…n的式(6)中對(duì)應(yīng)的矩陣B為：

式中

求帶橢球約束的隱式曲面方程實(shí)際上又變成求在24J-I＞0條件下的線性系統(tǒng)，如下式：

通過(guò)式(15)、(16)和(21)中的24J-I=1可以求出一個(gè)矩陣C0：

使得：

再定義一個(gè)矩陣C：

條件4J-I2=1可以轉(zhuǎn)化為：

式中的O11,O12,O21分別為(n-6)×(n-6),(n-6)×6,6×(n-6)大小的零矩陣?，F(xiàn)將式(7)、(8)式的矩陣重新分塊：

式中，M11,M12,M21分別是大小為(n-6)×(n-6),(n-6)×6,6×(n-6)，α,β的大小分別為(n-6)×1和6×1，式(21)現(xiàn)可簡(jiǎn)化為：

當(dāng)M11是非奇異時(shí)，式(26)的求解等價(jià)于求解下式中的β：

式中D=M22-MT12M11-1M12=-M12TMM11MM12,因M22是零矩陣，可以直接求出β，α。但當(dāng)D是非奇異矩陣時(shí)，式(27)無(wú)解，此時(shí)通常會(huì)選用最小化方法來(lái)求解。首先矩陣D是對(duì)稱的、正定的，利用下式求解：

帶橢球約束徑向基函數(shù)方程的求解變換為常見(jiàn)的特征值系統(tǒng)問(wèn)題：

對(duì)一般的徑向基函數(shù)，矩陣D不一定必須為正定的。不是正定時(shí)，需用至少是半正定的TDD來(lái)代替D，方程的解是最小非負(fù)特征值對(duì)應(yīng)的特征向量，但解不能保證唯一。

3 隱式曲面光滑拼接

本文選用文獻(xiàn)[10]提出的可自行控制參數(shù)的光滑絕對(duì)函數(shù)來(lái)對(duì)分塊后的重構(gòu)模型進(jìn)行光滑拼接。

定義1.表示當(dāng)x大于等于0時(shí)=x，當(dāng)x小于0時(shí)=-x，引進(jìn)如下絕對(duì)值函數(shù)：

對(duì)于光滑拼接我們只需保證光滑函數(shù)達(dá)到平方次光滑即可。式(30)的一階和二階光滑函數(shù)的顯式形式為：

式(30)可轉(zhuǎn)化為另一種等價(jià)形式：

式中的Gn(x)為：

依據(jù)式(33)可以計(jì)算出和：

定義2.對(duì)于δ＞0并且n＞0，我們定義：

簡(jiǎn)明地說(shuō)就是，當(dāng)＞δ時(shí),δ=,事實(shí)上，當(dāng)≥δ時(shí)有≥n。故：

由此可得出光滑拼接的函數(shù)為：

由以上定義得出由n,δ共同控制的光滑函數(shù)，要根據(jù)具體需要來(lái)選擇參數(shù)的值。一般情況選擇n=2,因?yàn)槎喂饣矢咔夜饣Ч诲e(cuò)。對(duì)于δ，對(duì)兩個(gè)半徑大小為0.5的圓球進(jìn)行實(shí)驗(yàn)，根據(jù)效果來(lái)確定，如圖1所示。

4 實(shí)驗(yàn)與分析

為驗(yàn)證本算法的效果和效率，與文獻(xiàn)[1-2]、[4-5]進(jìn)行了對(duì)比實(shí)驗(yàn)。硬件實(shí)驗(yàn)環(huán)境為：AMD X4 740，4 G內(nèi)存，軟件環(huán)境為Windows8 64位操作系統(tǒng)、VS2008 C++、OpenGL,Matlab2013等。

下面將對(duì)femur、head兩種點(diǎn)云進(jìn)行五種方法的重建工作，并對(duì)其重建效果進(jìn)行分析，詳見(jiàn)圖2、3中的(b)、(c)、(d)、(e)、(f)分別對(duì)應(yīng)文獻(xiàn)[1-2]、[4-5]和本文方法的效果圖。

圖2和表1表明各種算法對(duì)femur點(diǎn)云的重建效果，算法[1]、[4]、[5]三種在效果上、效率上都沒(méi)有優(yōu)勢(shì)，算法[2]雖能大體上保特征，但效率卻遠(yuǎn)低于本算法。

圖1 光滑函數(shù)操作

圖3和表1表明各種算法對(duì)head點(diǎn)云的重建效果，算法[1]、[4]、[5]明顯不保特征，算法[2]基本保特征，但五官信息、喉結(jié)特征沒(méi)有保留，本算法雖然出現(xiàn)腫脹部分，但細(xì)節(jié)特征依然保留。

圖2 femur點(diǎn)云擬合

總體來(lái)說(shuō)，算法[1]是對(duì)點(diǎn)云進(jìn)行三角剖分，重建模型外表比較粗糙。算法[2]對(duì)點(diǎn)云進(jìn)行Loop細(xì)分重建，重建模型表面雖然光滑了很多，但卻丟失了細(xì)節(jié)特征。算法[4]利用基于徑向基函數(shù)隱式曲面重建，該方法雖然重建模型表面光滑連續(xù)，但模型不飽和，有失真現(xiàn)象。算法[5]利用平面約束的徑向基函數(shù)隱式曲面重建，相對(duì)于算法[4]而言，重建模型稍有飽和，較精確地重建出模型，但有一個(gè)明顯的平面，使得整個(gè)模型重建效果不好。

對(duì)于本文提出的橢球約束徑向基函數(shù)隱式曲面重建方法?？梢钥醋鱿葘?duì)點(diǎn)云進(jìn)行基于徑向基函數(shù)方法重建，然后對(duì)整個(gè)模型曲面進(jìn)行橢球曲面約束。相當(dāng)于一個(gè)徑向基函數(shù)隱式曲面和一個(gè)橢球隱式曲面的線性光滑混合。利用橢球的幾何特性對(duì)徑向基函數(shù)隱式曲面進(jìn)行飽和處理，從而使得整個(gè)模型飽和精確。

實(shí)驗(yàn)表明基于橢球約束的徑向基函數(shù)方法重建后的模型比較精確，且圓潤(rùn)飽和。但對(duì)于有尖銳特征的點(diǎn)云重建時(shí)(如head)，因?yàn)闄E球的作用會(huì)形

表1 實(shí)驗(yàn)數(shù)據(jù)對(duì)比

圖3 head點(diǎn)云擬合

成腫脹部分。對(duì)于此缺陷，需對(duì)點(diǎn)云進(jìn)行簡(jiǎn)單重疊分割。首先觀察head點(diǎn)云，在耳朵處(y=1.5)出現(xiàn)腫脹，需在此進(jìn)行線性分割，分割為上下兩部分。這兩部分需要有重疊，以便拼接處理，所以將y=1.0以上部分分割為上半部分，將y=2.0以下部分分割為下部分。如圖4所示，圖中(a)表示上半部分，(b)表示下半部分。如圖5所示，將帶有尖銳特征點(diǎn)的云分割后，會(huì)使得尖銳特征性變小，然會(huì)對(duì)分割后的局部點(diǎn)云重建并進(jìn)行橢球約束，從而在不出現(xiàn)腫脹現(xiàn)象的情況下，局部點(diǎn)云生成的模型更具有飽和性，保特征性。最后通過(guò)對(duì)局部模型進(jìn)行光滑拼接。共消耗時(shí)間4.25 s，效率明顯提高。

圖4 head點(diǎn)云分割示意圖

圖5 光滑拼接效果

5 結(jié) 論

本文提出的基于橢球約束的徑向基函數(shù)隱式曲面重建算法，直接從散亂點(diǎn)獲取重建信息，重建結(jié)果精確。但對(duì)大規(guī)模點(diǎn)云集擬合時(shí)，需解大規(guī)模矩陣，效率低下甚至擬合模型會(huì)產(chǎn)生冗余部分。要通過(guò)對(duì)點(diǎn)云進(jìn)行簡(jiǎn)單的分割再進(jìn)行擬合光滑拼接，不僅效果好而且是并行擬合，這樣時(shí)間效率會(huì)飛速提高。實(shí)驗(yàn)效果表明：本算法相對(duì)于上面提到的效果保特征性好而且效率很高。

[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]陳慧群,黎景炎.基于變形網(wǎng)格法的高密點(diǎn)云曲面重建[J].工程圖學(xué)學(xué)報(bào),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 重疊區(qū)域作用效果混合[J].中國(guó)圖象圖形學(xué)報(bào),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.