基于C1自然鄰近插值的曲面擬合

聶志峰， 周慎杰， 王 凱， 孔勝利

（1. 山東大學(xué)機(jī)械工程學(xué)院，山東 濟(jì)南 250061；2. 山東科技大學(xué)機(jī)械電子工程學(xué)院，山東 青島 265510）

基于C1自然鄰近插值的曲面擬合

聶志峰1,2， 周慎杰1， 王 凱1， 孔勝利1

（1. 山東大學(xué)機(jī)械工程學(xué)院，山東 濟(jì)南 250061；2. 山東科技大學(xué)機(jī)械電子工程學(xué)院，山東 青島 265510）

以C0連續(xù)non-Sibsonian 插值作為三次單純形Bernstein-Bézier多項式的基坐標(biāo)，構(gòu)造C1連續(xù)自然鄰近插值函數(shù)。介紹了高階連續(xù)函數(shù)的構(gòu)建原理和性質(zhì)。將C1連續(xù)自然鄰近插值函數(shù)應(yīng)用于曲面擬合場合，由于Voronoi圖能夠自動調(diào)整數(shù)據(jù)點分布不規(guī)則和密度不均勻在空間上的差異，即使對于散亂數(shù)據(jù)點，也能獲得較好的擬合結(jié)果。

計算機(jī)應(yīng)用；曲面擬合；C1自然鄰近插值；Bernstein-Bézier多項式

對于規(guī)則網(wǎng)格數(shù)據(jù)點的插值問題，通常采用樣條曲線來構(gòu)建插值函數(shù)。Lai和Wenston采用連續(xù)三次樣條函數(shù)計算了橢圓型偏微分方程問題[1]。對于不規(guī)則數(shù)據(jù)點的曲線曲面擬合，通常采用最小二乘法[2]，通過使誤差的平方和最小，得到一個線性方程組，求解線性方程組可以得到擬合曲線或者曲面。但是，如果離散數(shù)據(jù)量比較大、形狀復(fù)雜，需要分段擬合和平滑處理，實際操作中存在一定的困難。為了克服以上困難，文獻(xiàn)[3]利用插值與逼近相結(jié)合的曲面擬合思路，較好地解決原始數(shù)據(jù)點分布不均勻的曲面造型問題。移動最小二乘法對最小二乘法作了改進(jìn)，能夠生成精度高﹑光滑性好的曲線曲面[4-5]。自然鄰近插值法是一種較新的插值理論，主要用于多變量不規(guī)則數(shù)據(jù)點的插值。由于Voronoi圖和對偶的 Delaunay三角化能夠自動調(diào)整數(shù)據(jù)點分布不規(guī)則和密度不均勻在空間上的差異，即使對于散亂數(shù)據(jù)點，也可以得到較好的曲面擬合結(jié)果。另外，自然鄰近插值法具有嚴(yán)格的插值特性，偏微分方程Garlerkin法可以直接施加本質(zhì)邊界條件，這是移動最小二乘法所無法相比的。Tily和Brace以汽車發(fā)動機(jī)半規(guī)則實驗數(shù)據(jù)為研究對象，對C0連續(xù)Sibson插值函數(shù)和MATLAB interp2 插值函數(shù)進(jìn)行了比較，發(fā)現(xiàn)當(dāng)實驗數(shù)據(jù)標(biāo)準(zhǔn)化處理后，Sibson 插值函數(shù)具有更高的插值精度[6]。高階偏微分方程的數(shù)值計算要求插值函數(shù)至少 C1連續(xù)，對 C1連續(xù)插值函數(shù)的插值性能進(jìn)行深入研究是有必要的。Farin將C0連續(xù)Sibson插值函數(shù)以自然鄰近坐標(biāo)的形式代入三次單純形Bernstein-Bézier曲面多項式，構(gòu)造出C1連續(xù)插值函數(shù)[7]。Sukumar和Moran提出了一種新的計算方法，將Bernstein基函數(shù)轉(zhuǎn)換為C1形函數(shù)，保證任意結(jié)點I的3個形函數(shù)分別與它的3個自由度、和直接相聯(lián)，從而適合于高階偏微分方程 Galerkin法[8]。與Sibson插值比較，non-Sibsonian插值具有計算簡單、計算效率高等優(yōu)點[9]，所以本文將non-Sibsonian 插值函數(shù)作為三次單純形Bernstein-Bézier曲面多項式的基坐標(biāo)，采用文獻(xiàn)[8]的計算方法，構(gòu)造出新的連續(xù)自然鄰近插值函數(shù)，該函數(shù)具有二次完備性、對結(jié)點函數(shù)值和梯度值的插值特性以及在分析域的邊界退化為三次多項式等性質(zhì)。為了考察新建連續(xù)插值函數(shù)的插值精度及其影響因素，作者將它應(yīng)用到曲面擬合場合。

1 C1自然鄰近插值

1.1 non-Sibsonian自然鄰近插值

Voronoi 圖及其對偶的Delaunay三角化源自計算幾何，是由一組不規(guī)則的點定義的最基本的幾何結(jié)構(gòu)?？紤]平面上由m個離散結(jié)點組成的集合，結(jié)點I的Voronoi圖是指與結(jié)點I對應(yīng)的區(qū)域，內(nèi)任意點到結(jié)點I的距離都小于該點到其它結(jié)點J的距離，用數(shù)學(xué)公式表示為

當(dāng)積分點x被引進(jìn)集合N的一階Voronoi 圖時，根據(jù)空圓準(zhǔn)則可以確定點x的自然鄰子。圖1是平面上7個結(jié)點的一階Voronoi 圖，點x共有4個自然鄰子，分別是結(jié)點1、2、3和4。

假設(shè)點x有n個自然鄰子，該點關(guān)于結(jié)點I的non-Sibsonian插值函數(shù)可以定義為

1.2 以non-Sibsonian插值作為自然鄰近坐標(biāo)的C1自然鄰近插值函數(shù)

將 non-Sibsonian插值φ代入三次單純形Bernstein-Bézier曲面多項式，沿用文獻(xiàn)[8]的計算方法，構(gòu)造出C1自然鄰近插值函數(shù)

對公式（3）作如下變換，得到

1.2.1 Bernstein-Bézier基函數(shù)

non-Sibsonian坐標(biāo)φ下由積分點x的n個自然鄰子組成的n邊形的三次基函數(shù)為

對于三次Bernstein-Bézier曲面，基函數(shù)的多指標(biāo)記法只會出現(xiàn)下面3種情況：和。表示數(shù)組中只有第α個元素為1，其余為0的多指標(biāo)記法。

1.2.2 變換矩陣的構(gòu)建

通常，可以把Bézier空間坐標(biāo)分為頂點、切面和中心Bézier空間坐標(biāo)三種。三種Bézier空間坐標(biāo)與結(jié)點函數(shù)值和梯度值具有不同的構(gòu)造關(guān)系：頂點Bézier空間坐標(biāo)等于結(jié)點函數(shù)值，即

切面Bézier空間坐標(biāo)為

中心Bézier空間坐標(biāo)的最佳方案為[7]

式中

由公式（4）可以看出：變換矩陣的作用是把結(jié)點函數(shù)值或者梯度值映射到Bézier空間坐標(biāo)上。利用公式（6）～（8），可以構(gòu)建出變換矩陣。一旦變換矩陣構(gòu)建完畢，Bernstein-Bézier基函數(shù)和變換矩陣相乘，可以得到形函數(shù)

圖3 結(jié)點網(wǎng)格及?形函數(shù)ψ3A-2(x)、ψ3A-1(x)

（1） 高階連續(xù)性

（2） 二次完備性

式中 a為切面Bézier空間坐標(biāo)的質(zhì)心，

c為頂點Bézier空間坐標(biāo)的質(zhì)心。

（3） 對結(jié)點函數(shù)值和結(jié)點梯度值的插值特性

（4） 分析域邊界積分點的插值函數(shù)是三次多項式

如果點x的兩個自然鄰子的 non-Sibsonian坐標(biāo)為，把它們代入方程（12），可以得到

顯然，分析域邊界上的插值函數(shù)是三次多項式。

2 曲面擬合

3個曲面函數(shù)分別定義為

圖4 單位正方形離散方案

為了衡量曲面擬合質(zhì)量，定義整體域?的相對誤差范Le

圖5 規(guī)則25結(jié)點的3個函數(shù)曲面與擬合曲面的對比（左為函數(shù)曲面，右為擬合曲面）

Le值越小，曲面擬合結(jié)果越好。3個函數(shù)在5種離散方案下的相對誤差范Le，見表1。

表1 曲面擬合相對誤差范Le

3 結(jié) 論

（2） 對于高次多項式函數(shù)（大于等于3次），曲面擬合質(zhì)量下降了。離散結(jié)點數(shù)量對曲面擬合質(zhì)量影響較大，加密結(jié)點可以提高連續(xù)插值函數(shù)的插值精度，得到質(zhì)量較高的擬合曲面。由于Voronoi圖能夠自動調(diào)整數(shù)據(jù)點分布不規(guī)則和密度不均勻在空間上的差異，所以對于不規(guī)則和密度不均勻數(shù)據(jù)點的情況，也得到了較好的曲面擬合結(jié)果。

[1] Lai M J, Wenston P. Report on numerical experi-ments with bivariate C1cubic splines for numerical solution of partial differential equations [R]. Technical Report, Department of Mathematics, University of Georgia, 1998.

[2] 彭芳瑜, 周 濟(jì), 周艷紅, 等. 基于最小二乘法的曲面生成算法研究[J]. 工程圖學(xué)學(xué)報, 1999, 20(3): 41-46.

[3] 彭芳瑜, 周云飛, 周 濟(jì). 基于插值與逼近的復(fù)雜曲面擬合[J]. 工程圖學(xué)學(xué)報, 2002, 23(4): 87-96.

[4] 曾清紅, 盧德唐. 基于移動最小二乘法的曲線曲面擬合[J]. 工程圖學(xué)學(xué)報, 2004, 25(1): 84-89.

[5] Lancaster P, Salkauskas K K. Surfaces generated by moving least squares methods [J]. Mathematics of Computation, 1981, (37): 141-158.

[6] Tily R, Brace C J. A study of natural neighbour interpolation and its application to automotive engine test data [J]. Proceeding of I mech E, Part D: Journal of Automobile Engineering, 2006, (220): 1003-1017.

[7] Farin G. Surfaces over Dirichlet tessellations [J]. Computer Aided Geometric Design, 1990, (7): 281-292.

[8] Sukumar N, Moran B.natural neighbor interpo-lant for partial differential equations [J]. Numerical Methods for Partial Differential Equations, 1999, 15(4): 417-447.

[9] Sukumar N, Moran B, Semenov A Y, et al. Natural neighbor Galerkin methods [J]. International Journal for Numerical Methods in Engineering, 2001, (50): 1-27.

Surface Fitting Based on C1Natural Neighbor Interpolation

NIE Zhi-feng1,2, ZHOU Shen-jie1, WANG Kai1, KONG Sheng-li1
( 1. School of Mechanical Engineering, Shandong University, Jinan Shandong 250061, China; 2. College of Mechanical and Electronic Engineering, Shandong University of Science and Technology, Qingdao Shandong 265510, China)

C1natural neighbor interpolation can be realized when C0non-Sibsonian interpolation is introduced in the Bernstein-Bézier surface representation of a cubic simplex. The principle and properties of C1natural neighbor interpolation are described and the surface fitting is given. The Voronoi diagram can adjust the spatial discrepancy caused by irregular data and data of varying density, so even for the scattered data, C1natural neighbor interpolation can get accurate fitting results.

computer application; surface fitting; C1natural neighbor interpolation; Bernstein-Bézier polynomial

TP 391

A

：1003-0158(2010)01-0110-06

2008-08-01

國家自然科學(xué)基金資助項目（10572077）；山東省自然科學(xué)基金資助項目（Y2007F20）；山東科技大學(xué)“春蕾計劃”資助項目（2009AZZ021）

聶志峰（1970-），男，山東日照人，博士研究生，主要研究方向為數(shù)值計算方法。