圖像中對(duì)數(shù)螺旋線的擬合

蒲平

（暨南大學(xué) 計(jì)算機(jī)科學(xué)系，廣東 廣州 510632）

對(duì)數(shù)螺旋線（等角螺旋線）是一種在自然界中經(jīng)常出現(xiàn)的曲線。由于它具有良好的幾何、數(shù)學(xué)、力學(xué)及美學(xué)性質(zhì),使得它在計(jì)算幾何、工業(yè)造型、天文及氣象研究、動(dòng)植物研究等現(xiàn)代學(xué)科中得到了廣泛的應(yīng)用[1]。

對(duì)數(shù)螺旋線是一條古老的曲線，很早以前人們就在幾何上對(duì)它進(jìn)行了一系列的理論研究。衛(wèi)星云圖中臺(tái)風(fēng)的形成路徑可以近似看作為一條螺旋線，進(jìn)而根據(jù)在圖中提取的部分曲線段所進(jìn)行的螺旋線擬合來確定臺(tái)風(fēng)中心[2]。銀河星系中的渦旋星系圖像中，旋臂結(jié)構(gòu)所構(gòu)成的形狀近似為對(duì)數(shù)螺旋線。由于星系演化和旋臂結(jié)構(gòu)之間的相互影響，通過研究旋臂結(jié)構(gòu)可以了解星系的演化[3]。

由于對(duì)數(shù)螺旋線是一條非閉合的曲線，不能簡單地按照傳統(tǒng)的閉合曲線擬合算法來對(duì)其進(jìn)行擬合。目前可以通過將對(duì)數(shù)螺旋線的極坐標(biāo)方程 r=aebθ變換為 lnr=bθ+lna的形式，也就是把局部的對(duì)數(shù)螺旋線上各點(diǎn)的坐標(biāo)通過在不同的極坐標(biāo)系下進(jìn)行變換使其成為能最佳地滿足直線方程的形式,可以通過最小二乘法[4]或Hough變換[5]來求得該直線方程的系數(shù)，進(jìn)而擬合出最佳逼近的對(duì)數(shù)螺旋線。

這些方法都是通過不斷掃描圖像中的每一個(gè)像素點(diǎn)來假設(shè)其為對(duì)數(shù)螺旋線極坐標(biāo)方程中的原點(diǎn)(x0,y0)，然后在以該點(diǎn)建立的極坐標(biāo)系中根據(jù)圖像中所獲取的對(duì)數(shù)螺旋線上的點(diǎn)來求出極坐標(biāo)系下的極半徑和旋轉(zhuǎn)角。由公式求出的旋轉(zhuǎn)角也只能是在0～180°的范圍。而實(shí)際給出的對(duì)數(shù)螺旋線的旋轉(zhuǎn)角度可以大于360°。

本文通過縮小極坐標(biāo)原點(diǎn)的遍歷區(qū)域來減少運(yùn)算次數(shù)，達(dá)到提高運(yùn)算效率的目的。同時(shí)考慮所獲得弧段的極角大小在任意范圍內(nèi)。

1 對(duì)數(shù)螺旋線的擬合

1.1 對(duì)數(shù)螺旋線的主要性質(zhì)

對(duì)數(shù)螺旋線的極坐標(biāo)方程為:

其中，θ為極角，r(θ)為極徑，α為極徑與切線間的夾角,如圖1所示。

根據(jù)式(1)可以看出,對(duì)每個(gè)θ值，都有一個(gè)對(duì)應(yīng)的r值,而且不同的 θ值所對(duì)應(yīng)的 r值也不同 (因?yàn)?cotα≠0)。這種現(xiàn)象表示：從對(duì)數(shù)螺旋線上某個(gè)點(diǎn)出發(fā)，隨著θ值的無限增大與無限減小，此曲線會(huì)環(huán)繞它的極點(diǎn)形成無數(shù)多圈，一面是越繞越遠(yuǎn)，一面是越繞越聚集在極點(diǎn)附近。若 cotα>0，則當(dāng) θ→-∞ 時(shí)，曲線聚集在極點(diǎn)附近；若 cotα<0，則當(dāng) θ→-∞ 時(shí)，曲線越繞越遠(yuǎn),如圖 2所示。

等角螺線上的一段弧經(jīng)伸縮若干倍后，必與該等角螺線上的另一弧全等。事實(shí)上，若對(duì)數(shù)螺旋線r=aeθcotα伸縮 成 r=ae(θ+φ)cotα，則 在 對(duì) 數(shù) 螺 旋 線 r=aeθcotα上 ，極 角 θ滿 足β≤θ≤γ的弧，經(jīng)伸縮后必與該等角螺線上極角θ滿足β+φ≤θ≤γ+φ的弧全等。對(duì)數(shù)螺旋線的這項(xiàng)特性使得自然界中許多物體都呈現(xiàn)對(duì)數(shù)螺旋線的形狀。例如許多貝殼都很接近對(duì)數(shù)螺旋線的形狀，因?yàn)樯钤跉?nèi)的動(dòng)物在成長過程中都是均勻地長大的,這就像相似地放大，所以，新生的部分所棲息的空間必與原有空間形狀相似。象鼻、動(dòng)物的角與毛等都呈對(duì)數(shù)螺旋線形。

1.2 坐標(biāo)變換

將參數(shù)方程 r=aebθ兩邊同時(shí)取對(duì)數(shù),得到：lnr=bθ+lna。可以將該式看作是關(guān)于θ和lnr的直線形式。由于從已知圖像中獲得的數(shù)據(jù)點(diǎn)的坐標(biāo)為直角坐標(biāo)系下的值，即需將直角坐標(biāo)系下的點(diǎn)(x,y)變換到(θ,lnr)坐標(biāo)系下。

選取圖像中任意一點(diǎn)(x0,y0)為極坐標(biāo)原點(diǎn)建立新的極坐標(biāo)系。根據(jù)在圖像中獲取的一系列數(shù)據(jù)點(diǎn)(xi,yi)(i=1,2,…,k),通 過 式(2)、 式(3)可 將 其 轉(zhuǎn) 換 為(θi,lnri)(i=1,2,…,k)形式的數(shù)據(jù)點(diǎn)。

1.3 對(duì)數(shù)螺旋線中心點(diǎn)的取值范圍

由于處理過程中需要得到r和θ的值，而它們需要在極坐標(biāo)系下獲取，圖像中的任意一點(diǎn)(xi,yi)在選擇不同的極坐標(biāo)原點(diǎn)所建立的極坐標(biāo)系下所獲得r和θ的值也不相同。由于給定圖像中并不知道該曲線的極坐標(biāo)原點(diǎn)，故常用的方法是假設(shè)極坐標(biāo)原點(diǎn)存在于該幅圖像中，進(jìn)而通過對(duì)整幅圖像進(jìn)行掃描來尋找該點(diǎn)。但這需要遍歷整個(gè)圖像，效率十分低下。

由對(duì)數(shù)螺旋線的每一個(gè)性質(zhì)可知，該曲線由內(nèi)到外逐漸發(fā)散或者是由外到內(nèi)逐漸聚集在極點(diǎn)附近。對(duì)于超過一圈的螺旋線而言,可以通過逐行與逐列地對(duì)整幅圖像進(jìn)行掃描，判斷掃描直線與圖像中螺旋線的交點(diǎn)個(gè)數(shù)，如圖3所示?？梢园l(fā)現(xiàn)越靠近螺旋線的中心其交點(diǎn)個(gè)數(shù)越多。分別在水平和豎直方向上統(tǒng)計(jì)其交點(diǎn)的個(gè)數(shù)，通過設(shè)定一定的誤差范圍可以獲得一個(gè)交點(diǎn)個(gè)數(shù)較大的區(qū)間。由于在水平和豎直方向上都能獲得一個(gè)區(qū)間,從而可以獲得一個(gè)包含中心點(diǎn)的矩形區(qū)域。則中心點(diǎn)的選取可以縮小到該區(qū)域來進(jìn)行遍歷。

1.4 旋轉(zhuǎn)角的獲取

對(duì)數(shù)螺旋線是一條非閉合的曲線，它能包含多圈。如圖4所示，點(diǎn)A和點(diǎn)B由式(2)可以得到不同的極徑值 r。由式(3)求得的極角的值相同，但是正確的應(yīng)該是如果點(diǎn)A的極角大小為 θa時(shí)，點(diǎn) B的極角大小應(yīng)該為(360°+θa)。

對(duì)于具有多圈并且曲線中只含有兩個(gè)端點(diǎn)的對(duì)數(shù)螺旋線，假設(shè)起始端點(diǎn)的極角為θ1，則沿著曲線進(jìn)行跟蹤,由式(3)得到的極角值的變化趨勢(shì)應(yīng)該為從 θ1遞增到360°，接著再從 0°遞增到 360°直到終點(diǎn)?？梢园l(fā)現(xiàn)每當(dāng)極角的值從 360°變?yōu)?0°時(shí)則剛好走完一圈，則剛好需要增加 360°。

通過上述分析可以得出，當(dāng)獲得的極角值前后兩個(gè)數(shù)據(jù)劇烈遞減且二者之差為一個(gè)較大的正數(shù)時(shí)，需要增加 360°,即 θi=(θ1+360k)(k=1,2,…)。

1.5 最小二乘擬合

由上面的分析可知,對(duì)于要搜索的每一個(gè)極坐標(biāo)原點(diǎn)(x0,y

0)都能由已知圖像中獲得的點(diǎn)經(jīng)變換而得到一系列數(shù)據(jù)點(diǎn)(θi,lnri)。通過這些點(diǎn)由最小二乘法[5]可以求出直線的參數(shù)b和lna。根據(jù)誤差公式可以求出對(duì)于每個(gè)極坐標(biāo)原點(diǎn)所對(duì)應(yīng)的誤差值大小，滿足誤差最小的直線參數(shù)即為所求。

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

為驗(yàn)證本文方法的有效性，分別采用仿真圖像及真實(shí)圖像對(duì)算法進(jìn)行測試。實(shí)驗(yàn)環(huán)境：CPU主頻為2.80 GHz，內(nèi)存為256 MB。實(shí)驗(yàn)結(jié)果如圖5所示。

圖5(a)為具有多圈的對(duì)數(shù)螺旋線人工仿真圖。圖5(b)是對(duì)原圖細(xì)化到單像素寬度[6]，分別找出起始點(diǎn)和終點(diǎn)。圖5(c)為分別從水平方向和豎直方向?qū)D像進(jìn)行掃描統(tǒng)計(jì)的交點(diǎn)個(gè)數(shù)。圖5(d)為在水平和豎直方向上獲得的交點(diǎn)個(gè)數(shù)較多的區(qū)域。圖5(e)為通過遍歷由圖5(d)獲得的區(qū)域內(nèi)每一像素點(diǎn)并以其為極坐標(biāo)原點(diǎn)建立極坐標(biāo)系，將原圖像中數(shù)據(jù)點(diǎn)變換為該坐標(biāo)系下，通過最小二乘擬合后求出的誤差最小的曲線。該方法耗時(shí)為2.111 s，如果在整個(gè)圖像中遍歷尋找極坐標(biāo)原點(diǎn)，則需要耗時(shí)58.467 s。

表1給出了對(duì)數(shù)螺旋線的原始數(shù)據(jù)和擬合后的參數(shù)，其中(x0,y0)為中心點(diǎn)坐標(biāo)，a和b分別為對(duì)數(shù)螺旋線中的參數(shù)值。由表1可知，擬合出的對(duì)數(shù)螺旋線與原始圖中心點(diǎn)相同，參數(shù)a和b值基本一致。

表1 人工仿真對(duì)數(shù)螺旋線擬合結(jié)果

圖6描述了對(duì)臺(tái)風(fēng)所形成的衛(wèi)星云圖進(jìn)行處理的結(jié)果。圖6(a)為真實(shí)衛(wèi)星云圖，圖中含有螺旋云帶。圖6(b)為提取原始圖像中的部分區(qū)域。圖6(c)為使用otsu方法得到的二值化圖像。圖6(d)為將二值化圖像通過形態(tài)學(xué)方法(膨脹和腐蝕等)將圖像中的螺旋云帶進(jìn)行獨(dú)立的分割所得到的圖像。圖6(e)為通過面積分割提取出來的螺旋云帶。圖6(f)為將螺旋云帶中的空洞去掉后的結(jié)果。圖6(g)為對(duì)圖像進(jìn)行細(xì)化后的結(jié)果。圖6(h)為通過找出曲線中的各個(gè)端點(diǎn)，提取出曲線中最長的一條曲線段,以此來對(duì)細(xì)化后的圖像去掉毛刺。圖6(i)為對(duì)提取出來的曲線段進(jìn)行樣條平滑后的結(jié)果。圖6(j)為分別在水平和豎直方向上進(jìn)行掃描統(tǒng)計(jì)出的交點(diǎn)個(gè)數(shù)。圖6(k)為找出統(tǒng)計(jì)的交點(diǎn)個(gè)數(shù)較多的相交區(qū)域。圖6(m)和圖6(n)為用最小二乘法對(duì)其進(jìn)行擬合后的結(jié)果。

本文通過約束對(duì)數(shù)螺旋線中心點(diǎn)的存在范圍，將對(duì)數(shù)螺旋線的擬合轉(zhuǎn)換為直線的擬合。實(shí)驗(yàn)結(jié)果表明,本文方法能夠快速、較為準(zhǔn)確地?cái)M合出圖像中存在的對(duì)數(shù)螺旋線。然而實(shí)際中，包含不同螺旋線的圖像需要具體地分析，尋找更為合理有效的方法。

[1]游世輝,李軍.對(duì)數(shù)螺線的應(yīng)用研究[J].九江學(xué)院學(xué)報(bào)(自然科學(xué)版),2005(3):6-9.

[2]謝俊元,艾早陽.臺(tái)風(fēng)中心定位中的螺旋線自動(dòng)識(shí)別算法[J].軟件學(xué)報(bào),1997,6(增刊):398-403.

[3]FRUHWIRTH R,STRANDLIE A.Helix fitting by an extended riemann fit[J].Nuclear Instruments and Methods in Physics Research(A),2002,490(1-2):366-378.

[4]劉凱,黃峰.無眼臺(tái)風(fēng)自動(dòng)定位方法研究[J].信息與控制，2001,30(6):543-546.

[5]孔秀梅.形成期臺(tái)風(fēng)螺旋云帶的提取、描述及中心定位的研究[D].天津:天津大學(xué),2003.

[6]STEFANELL R,ROSENFELD A.Some parallel thinning algorithms for digital pictures[J].Journal of ACM,1971,18(2):255-264. (收稿日期：2011-02-28)