王曾珍,劉華勇

(安徽建筑大學(xué) 數(shù)理學(xué)院,合肥 230601) E-mail：wangz_z1014@163.com

1 引 言

數(shù)據(jù)擬合是幾何設(shè)計(jì)中廣泛使用的建模工具.曲線的擬合問(wèn)題通常出現(xiàn)在計(jì)算機(jī)圖形學(xué)、計(jì)算機(jī)輔助幾何設(shè)計(jì)和圖像處理等各個(gè)領(lǐng)域.漸進(jìn)迭代逼近法(Progressive Iterative Approximation,PIA)作為一種有效而直觀的數(shù)據(jù)擬合方法,在近幾年來(lái)受到廣泛關(guān)注.該方法起源于齊東旭等[1]和de Boor[2]針對(duì)均勻三次B-Spline曲線擬合的迭代性質(zhì)的研究,它不僅可以避免求解大型線性方程組的計(jì)算成本,還可以通過(guò)對(duì)控制頂點(diǎn)進(jìn)行迭代來(lái)調(diào)整曲線,直到獲得滿足誤差要求的擬合曲線.藺宏偉等[3]證明了非均勻三次B樣條曲線和曲面也具有PIA這一性質(zhì).此外,PIA方法還擴(kuò)展到了基于NTP基的混合曲線和曲面[4].為加快PIA方法的收斂速度,陸利正[5]提出了一種加權(quán)PIA格式(Weighted Progressive Iteration Approximation,WPIA),給出了最優(yōu)權(quán)值取法.陳杰等[6]給出了WPIA的2種推廣方法：帶權(quán)值的局部PIA方法和均勻參數(shù)化下的三角Bézier曲面的帶權(quán)PIA方法.劉曉艷等[7]結(jié)合Jacobi迭代法,提出了非均勻三次B樣條曲線的漸進(jìn)迭代插值算法—Jacobi-PIA算法,證明了其收斂速度優(yōu)于經(jīng)典的漸進(jìn)迭代插值算法.張莉等[8]提出一種帶互異權(quán)值的漸進(jìn)迭代逼近算法,具有操作靈活且可以根據(jù)需要調(diào)整各控制頂點(diǎn)的優(yōu)勢(shì).當(dāng)參數(shù)取合適值時(shí),該算法的收斂速度比WPIA算法更快.

在經(jīng)典的PIA方法中,控制點(diǎn)的數(shù)量等于數(shù)據(jù)點(diǎn)的數(shù)量.然而,當(dāng)數(shù)據(jù)點(diǎn)的數(shù)量非常大時(shí),這是不可行的.因此,Lin等[9]給出了一種擴(kuò)展的PIA(Extended Progressive Iteration Approximation,EPIA)方法,此方法允許控制點(diǎn)的數(shù)量小于給定數(shù)據(jù)點(diǎn)的數(shù)量,但EPIA擬合給定的數(shù)據(jù)點(diǎn)集時(shí),生成的曲線(曲面)序列的極限不是數(shù)據(jù)集的最小二乘擬合結(jié)果.因此,Deng等[10]提出了一種新的漸進(jìn)迭代最小二乘方法(progressive iterative approximation for least square fitting,LSPIA),LSPIA通過(guò)迭代調(diào)整控制點(diǎn)構(gòu)造一系列擬合曲線(曲面),極限曲線(曲面)是對(duì)給定數(shù)據(jù)點(diǎn)的最小二乘擬合結(jié)果.在每次迭代中,每個(gè)控制點(diǎn)的差向量是數(shù)據(jù)點(diǎn)和它們?cè)跀M合曲線(曲面)上對(duì)應(yīng)的點(diǎn)之間的差向量的加權(quán)和.Lin等[11]證明了最小二乘擬合方程組的系數(shù)矩陣為奇異矩陣時(shí),LSPIA方法是收斂的.為保證逼近生成的曲線(曲面)的極限是給定數(shù)據(jù)點(diǎn)的最小二乘擬合結(jié)果,Yusuf Fatihu Hamza等[12]提出一種基于Gauss-Seidel迭代方法的快速PIA算法(GS-LSPIA),該算法比LSPIA算法需要更少的步驟和更短的運(yùn)算時(shí)間,因此收斂速度也得到提高.

此外,近幾年來(lái),PIA方法還被擴(kuò)展至細(xì)分格式,如Wang等[13]提出了用于細(xì)分曲面插值的Gauss-Seidel幾何迭代法(Gauss-Seidel Progressive Iterative Approximation,GS-PIA).最近,Zhang等[14]將PIA方法應(yīng)用于單幀圖像超分辨率重建.

LSPIA方法通過(guò)不斷調(diào)整控制頂點(diǎn)獲得一系列擬合曲線(曲面),具備高效、穩(wěn)定地?cái)M合大規(guī)模數(shù)據(jù)點(diǎn)及顯著減少計(jì)算量等優(yōu)點(diǎn).然而,在LSPIA的迭代過(guò)程中,計(jì)算差向量時(shí)的權(quán)值始終為確定范圍內(nèi)的固定值,無(wú)法根據(jù)需要單獨(dú)控制部分控制頂點(diǎn)從而調(diào)整局部曲線.為使LSPIA方法通過(guò)調(diào)整權(quán)值從而更加靈活地調(diào)整局部曲線形狀,本文將上述LSPIA方法進(jìn)行改進(jìn),在LSPIA方法的基礎(chǔ)上,控制部分?jǐn)?shù)據(jù)點(diǎn)的權(quán)值不變,賦予其余控制頂點(diǎn)不同的權(quán)值,對(duì)LSPIA方法和改進(jìn)后的帶互異權(quán)值的LSPIA方法進(jìn)行分析和比較.隨后證明本文所提出方法的收斂性.最后給出數(shù)值實(shí)例,計(jì)算出LSPIA和本文所提方法的誤差,驗(yàn)證了本文方法的可行性.

2 帶互異權(quán)值的LSPIA算法

2.1 LSPIA算法描述

(1)

其中,{Bi(t);i=0,1,…,n}為三次均勻B樣條基函數(shù),其配置矩陣為：

(2)

令差向量為：

(3)

且第i個(gè)控制頂點(diǎn)的第一個(gè)調(diào)整向量為：

(4)

其中,權(quán)值μ滿足0<μ<2/λ0,λ0是矩陣BTB的最大特征值,由此可得新的控制頂點(diǎn)為：

(5)

新的曲線為：

(6)

類(lèi)似地,假設(shè)第k次迭代后獲得曲線Pk(t),令：

(7)

(8)

(9)

可得第k+1次迭代后的曲線為：

(10)

2.2 帶互異權(quán)值的LSPIA算法

帶互異權(quán)值的LSPIA算法是在原LSPIA算法[10]上作出的改進(jìn).在該算法中,曲線和差向量的表達(dá)形式不變,改變調(diào)整向量.將1.1節(jié)中的LSPIA算法中的常數(shù)μ賦予不同的值,即將迭代過(guò)程中調(diào)整向量(式(8))改為：

(11)

假設(shè)已知第k次迭代后的曲線Pk(t),由式(7)和式(11)可得,第k+1次迭代后的控制頂點(diǎn)為：

(12)

第k+1次迭代后得到的曲線為：

(13)

由此可得到一系列迭代逼近曲線{Pk(t),k=0,1,…}(k為迭代次數(shù)).當(dāng)μi都取相同值且滿足0<μ<2/λ0時(shí),則改進(jìn)的算法退化為L(zhǎng)SPIA算法.

3 帶互異權(quán)值的LSPIA算法的收斂性

本節(jié)將證明,當(dāng)μi滿足一定條件時(shí),上述帶互異權(quán)值的LSPIA算法是收斂的,且極限曲線為初始數(shù)據(jù)點(diǎn)的最小二乘擬合曲線.

定理2.當(dāng)μi滿足條件0<μ<2/λ0,那么帶互異權(quán)值的LSPIA算法是收斂的.

(14)

根據(jù)式(9)和式(11)可得：

(15)

將式(15)寫(xiě)成矩陣形式：

Pk+1=Pk+AT(Q-BPk)

(16)

其中A=UB,B為三次均勻B樣條基函數(shù)的配置矩陣.

令D=I-ATB,I為n+1階單位矩陣,由式(16)可得：

Pk+1=Pk+AT(Q-BPk)=
(I-ATB)[Pk-(ATB)-1ATQ]+(ATB)-1ATQ

(17)

則：

Pk+1-(ATB)-1ATQ=(I-ATB)[Pk-(ATB)-1ATQ]=
(I-ATB)2[Pk+1-(ATB)-1ATQ]=…=Dk+1[P0-(ATB)-1ATQ]

(18)

設(shè){λi(D)(i=0,1,…,n)}是按非遞減順序排列的D的特征值.令λi(D)=1-λi,其中λi是按非遞減順序排列的ATB的特征值.

根據(jù)式(18)可得：

P∞=(ATB)-1ATQ+D∞[P0-(ATB)-1ATQ]=(ATB)-1ATQ

(19)

式(19)等價(jià)于(ATB)P∞=ATQ,意味著帶互異權(quán)值的LSPIA算法是收斂的,且極限曲線收斂到初始數(shù)據(jù)點(diǎn)的最小二乘擬合結(jié)果.

4 數(shù)值實(shí)例

在實(shí)例分析這一部分,本文選取一般的三次均勻B樣條曲線擬合給定點(diǎn)集.通過(guò)實(shí)驗(yàn)調(diào)整常數(shù)μ的取值來(lái)控制曲線形狀的變化,μ值改變后對(duì)應(yīng)的控制頂點(diǎn)隨之變化,進(jìn)而可以調(diào)整局部曲線的形狀,并觀察初始擬合曲線和本文提出算法迭代后的擬合曲線改變局部形狀后的變化.

圖1 μ取不同值時(shí)LSPIA的逼近雙紐線Fig.1 Approximation lemniscate of LSPIA with different values

例1.在雙紐線p(t)=(x(t),y(t))=(cost,sint·cost)上選取11個(gè)點(diǎn)p(ti)=(x(ti),y(ti))(ti=-π/2+i·π/5,i=0,1,…,10)作為初始控制頂點(diǎn),迭代次數(shù)為k=15.圖1(a)為μ取相同值0.03時(shí)依據(jù)LSPIA方法將控制頂點(diǎn)迭代15次后,三次均勻B樣條曲線逼近本例中所取11個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線圖.觀察圖1(a)迭代結(jié)果,將第7、第8和第10個(gè)控制頂點(diǎn)進(jìn)行調(diào)整,根據(jù)實(shí)驗(yàn)值,令μ7=μ8=0.08,μ10=0.1,圖1(b)為調(diào)整后的曲線圖.可見(jiàn),調(diào)整后的圖1(b)中迭代后的逼近曲線與位于圖中左下部分第7、8個(gè)初始控制頂點(diǎn)的距離更近,局部曲線得到調(diào)整.

表1 例1中調(diào)整第i個(gè)控制頂點(diǎn)前后的逼近誤差Table 1 Approximation errors before and after adjusting the ith control vertex in example 1

表1為μ0=μ1=…=μ11=0.03和μ7=μ8=0.08,μ10=0.1,μ0=μ1=…=μ6=μ9=μ11=0.03時(shí),第7、第8和第10個(gè)控制頂點(diǎn)調(diào)整前后的逼近誤差.根據(jù)圖表所示,誤差減小,說(shuō)明依據(jù)本文提出方法進(jìn)行局部調(diào)整后的逼近性更好.

例2.在圓弧p(t)=(x(t),y(t))=(sint,cost)上選取11個(gè)點(diǎn)p(ti)=(x(ti),y(ti))(ti=7π/6+i·π/6,i=0,1,…,10)作為初始控制頂點(diǎn),迭代次數(shù)k=20.圖2(a)為μ取相同值0.5時(shí)依據(jù)LSPIA方法將控制頂點(diǎn)迭代20次后,3次均勻B樣條曲線逼近本例中所取11個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線圖.觀察圖2(a),第10個(gè)控制頂點(diǎn)附近所對(duì)應(yīng)的局部曲線逼近效果不佳,因此依據(jù)本文所提出方法調(diào)整對(duì)應(yīng)控制頂點(diǎn)的位置,令μ10=0.01,得到圖2(b).

圖2 μ取不同值時(shí)LSPIA的逼近圓弧Fig.2 Approximation arc of LSPIA with different values

表2為μ0=μ1=…=μ11=0.5和μ10=0.01,μ0=μ1=…=μ9=μ11=0.5時(shí),第10個(gè)控制頂點(diǎn)調(diào)整前后的逼近誤差.

表2 例2中調(diào)整第i個(gè)控制頂點(diǎn)前后的逼近誤差Table 2 Approximation errors before and after adjusting the ith control vertex in example 2

圖3 μ取不同值時(shí)LSPIA的逼近五角星Fig.3 Approximate five-pointed star of LSPIA with different values

例3.在五角星上選取21個(gè)點(diǎn)作為初始控制頂點(diǎn)：(6.1232×10-16,10),(1.1226,6.5451),(2.2451,3.0902),(5.8779,3.0902),(9.5106,3.0902),(6.5716,0.9549),(3.6327,-1.1803),(4.7553,-4.6353),(5.8779,-8.0902),(2.9389,-5.9549),(-7.0166×10-16,-3.8197),(-2.9389,-5.9549),(-5.8779,-8.0902),(-4.7553,-4.6353),(-3.6327,-1.1803),(-6.5716,0.9549),(-9.5106,3.0902),(-5.8779,3.0902),(-2.2451,3.0902),(-1.1226,6.5451),(6.1232×10-16,10),迭代次數(shù)為k=5.圖3(a)為μ取相同值0.4時(shí)依據(jù)LSPIA方法將控制頂點(diǎn)迭代5次后,三次均勻B樣條曲線逼近本例中所取21個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線圖.根據(jù)圖3(a)結(jié)果,利用本文所提出方法調(diào)整第5個(gè)控制頂點(diǎn),令μ5=0.9,其余μ的取值不變,仍迭代5次后得到逼近曲線圖3(b),將五角星的右上角調(diào)整后,逼近曲線近似插值第5個(gè)控制頂點(diǎn).圖3(c)為圖3(a)的局部放大圖,圖3(d)為圖3(b)的局部放大圖.

表3為μ0=μ1=…=μ11=0.4和μ5=0.9,μ0=μ1=…=μ4=μ6=…=μ21=0.4時(shí),第5個(gè)控制頂點(diǎn)調(diào)整前后的逼近誤差.

表3 例3中調(diào)整第i個(gè)控制頂點(diǎn)前后的逼近誤差Table 3 Approximation errors before and after adjusting the ith control vertex in example 3

例4.在曲線p(t)=(x(t),y(t))=(t,sint)上選取15個(gè)點(diǎn)p(ti)=(x(ti),y(ti))(ti=i·π/7,i=0,1,…,14)作為初始控制頂點(diǎn),迭代次數(shù)為20次.圖4(a)為μ取相同值1.5時(shí)依據(jù)LSPIA方法將初始控制頂點(diǎn)迭代20次后,三次均勻B樣條曲線逼近本例中所取15個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線圖.圖4(b)為依據(jù)本文提出方法調(diào)整第4個(gè)控制頂點(diǎn)對(duì)應(yīng)的權(quán)值μ4=0.001,保持其余μ值不變所得到的正弦曲線逼近圖.圖4(c)為4(a)的局部放大圖,圖4(d)為圖4(b)的調(diào)整局部放大圖.

圖4 μ取不同值時(shí)LSPIA的逼近正弦曲線Fig.4 Approximation sinusoidal curve of LSPIA with different values

表4 例4中調(diào)整第i個(gè)控制頂點(diǎn)前后的逼近誤差Table 4 Approximation errors before and after adjusting the ith control vertex in example 4

表4為μ0=μ1=…=μ15=1.5和μ4=0.001,μ0=μ1=μ2=μ3=μ5…=μ15=1.5時(shí),第4個(gè)控制頂點(diǎn)調(diào)整前后的逼近誤差.

5 結(jié) 語(yǔ)

本文在權(quán)值相同的LSPIA算法的基礎(chǔ)上,提出帶互異權(quán)值的LSPIA算法.給出權(quán)值取值范圍,證明權(quán)值在該范圍內(nèi)取值時(shí)帶互異權(quán)值的LSPIA算法的收斂性.根據(jù)LSPIA算法結(jié)果調(diào)整不同控制頂點(diǎn)對(duì)應(yīng)的權(quán)值,對(duì)曲線進(jìn)行局部調(diào)整.實(shí)例結(jié)果表明,調(diào)整后的局部曲線逼近誤差減小,達(dá)到調(diào)整局部曲線逼近效果的目的.本文依據(jù)LSPIA實(shí)驗(yàn)結(jié)果選取需要調(diào)整的權(quán)值,所給出權(quán)值為本文所提出方法迭代收斂范圍內(nèi)的實(shí)驗(yàn)值.因此,在權(quán)值選取和誤差控制范圍這方面仍待研究.