帶互異權(quán)值的B樣條曲線(xiàn)的最小二乘漸進(jìn)迭代逼近

王曾珍,劉華勇

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

1 引 言

數(shù)據(jù)擬合是幾何設(shè)計(jì)中廣泛使用的建模工具.曲線(xiàn)的擬合問(wèn)題通常出現(xiàn)在計(jì)算機(jī)圖形學(xué)、計(jì)算機(jī)輔助幾何設(shè)計(jì)和圖像處理等各個(gè)領(lǐng)域.漸進(jìn)迭代逼近法(Progressive Iterative Approximation,PIA)作為一種有效而直觀(guān)的數(shù)據(jù)擬合方法,在近幾年來(lái)受到廣泛關(guān)注.該方法起源于齊東旭等[1]和de Boor[2]針對(duì)均勻三次B-Spline曲線(xiàn)擬合的迭代性質(zhì)的研究,它不僅可以避免求解大型線(xiàn)性方程組的計(jì)算成本,還可以通過(guò)對(duì)控制頂點(diǎn)進(jìn)行迭代來(lái)調(diào)整曲線(xiàn),直到獲得滿(mǎn)足誤差要求的擬合曲線(xiàn).藺宏偉等[3]證明了非均勻三次B樣條曲線(xiàn)和曲面也具有PIA這一性質(zhì).此外,PIA方法還擴(kuò)展到了基于NTP基的混合曲線(xiàn)和曲面[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樣條曲線(xiàn)的漸進(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í),生成的曲線(xiàn)(曲面)序列的極限不是數(shù)據(jù)集的最小二乘擬合結(jié)果.因此,Deng等[10]提出了一種新的漸進(jìn)迭代最小二乘方法(progressive iterative approximation for least square fitting,LSPIA),LSPIA通過(guò)迭代調(diào)整控制點(diǎn)構(gòu)造一系列擬合曲線(xiàn)(曲面),極限曲線(xiàn)(曲面)是對(duì)給定數(shù)據(jù)點(diǎn)的最小二乘擬合結(jié)果.在每次迭代中,每個(gè)控制點(diǎn)的差向量是數(shù)據(jù)點(diǎn)和它們?cè)跀M合曲線(xiàn)(曲面)上對(duì)應(yīng)的點(diǎn)之間的差向量的加權(quán)和.Lin等[11]證明了最小二乘擬合方程組的系數(shù)矩陣為奇異矩陣時(shí),LSPIA方法是收斂的.為保證逼近生成的曲線(xiàn)(曲面)的極限是給定數(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)獲得一系列擬合曲線(xià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)整局部曲線(xiàn).為使LSPIA方法通過(guò)調(diào)整權(quán)值從而更加靈活地調(diào)整局部曲線(xiàn)形狀,本文將上述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)值μ滿(mǎn)足0<μ<2/λ0,λ0是矩陣BTB的最大特征值,由此可得新的控制頂點(diǎn)為：

(5)

新的曲線(xiàn)為：

(6)

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

(7)

(8)

(9)

可得第k+1次迭代后的曲線(xiàn)為：

(10)

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

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

(11)

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

(12)

第k+1次迭代后得到的曲線(xiàn)為：

(13)

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

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

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

定理2.當(dāng)μi滿(mǎn)足條件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算法是收斂的,且極限曲線(xiàn)收斂到初始數(shù)據(jù)點(diǎn)的最小二乘擬合結(jié)果.

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

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

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

例1.在雙紐線(xiàn)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樣條曲線(xiàn)逼近本例中所取11個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線(xiàn)圖.觀(guā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)整后的曲線(xiàn)圖.可見(jiàn),調(diào)整后的圖1(b)中迭代后的逼近曲線(xiàn)與位于圖中左下部分第7、8個(gè)初始控制頂點(diǎn)的距離更近,局部曲線(xià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樣條曲線(xiàn)逼近本例中所取11個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線(xiàn)圖.觀(guān)察圖2(a),第10個(gè)控制頂點(diǎn)附近所對(duì)應(yīng)的局部曲線(xiàn)逼近效果不佳,因此依據(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樣條曲線(xiàn)逼近本例中所取21個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線(xiàn)圖.根據(jù)圖3(a)結(jié)果,利用本文所提出方法調(diào)整第5個(gè)控制頂點(diǎn),令μ5=0.9,其余μ的取值不變,仍迭代5次后得到逼近曲線(xiàn)圖3(b),將五角星的右上角調(diào)整后,逼近曲線(xiàn)近似插值第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.在曲線(xiàn)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樣條曲線(xiàn)逼近本例中所取15個(gè)控制頂點(diǎn)連接的控制多邊形所得到的曲線(xiàn)圖.圖4(b)為依據(jù)本文提出方法調(diào)整第4個(gè)控制頂點(diǎn)對(duì)應(yīng)的權(quán)值μ4=0.001,保持其余μ值不變所得到的正弦曲線(xiàn)逼近圖.圖4(c)為4(a)的局部放大圖,圖4(d)為圖4(b)的調(diào)整局部放大圖.

圖4 μ取不同值時(shí)LSPIA的逼近正弦曲線(xiàn)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ì)曲線(xiàn)進(jìn)行局部調(diào)整.實(shí)例結(jié)果表明,調(diào)整后的局部曲線(xiàn)逼近誤差減小,達(dá)到調(diào)整局部曲線(xiàn)逼近效果的目的.本文依據(jù)LSPIA實(shí)驗(yàn)結(jié)果選取需要調(diào)整的權(quán)值,所給出權(quán)值為本文所提出方法迭代收斂范圍內(nèi)的實(shí)驗(yàn)值.因此,在權(quán)值選取和誤差控制范圍這方面仍待研究.