精確的快速極坐標(biāo)諧波變換*

沃焱 徐角

(華南理工大學(xué)計(jì)算機(jī)科學(xué)與工程學(xué)院，廣東廣州510006)

旋轉(zhuǎn)和縮放不變的模式表示是基于內(nèi)容的圖像檢索、模式識(shí)別等應(yīng)用領(lǐng)域的難點(diǎn)之一．在許多實(shí)際應(yīng)用中，圖像被縮放、旋轉(zhuǎn)后仍被視為內(nèi)容是相同的．研究者們已經(jīng)提出了許多旋轉(zhuǎn)不變特性的表示方法，其中包括比較廣泛使用的Zernike矩(ZMS)、偽 Zernike矩(PSMS)［1］．Zernike 矩、偽 Zernike 矩具有旋轉(zhuǎn)不變性，還具有對(duì)噪聲不敏感的特性，能較好地表達(dá)圖像特征，被廣泛地應(yīng)用于數(shù)字水印，人臉/字符識(shí)別，紋理分類、圖像檢索等領(lǐng)域［2-7］．盡管如此，由于計(jì)算Zernike矩、偽Zernike矩需多次階乘運(yùn)算，隨著階數(shù)增加，計(jì)算量和計(jì)算數(shù)值呈指數(shù)級(jí)增長(zhǎng)，不適合計(jì)算機(jī)存儲(chǔ)，且多次乘除運(yùn)算易產(chǎn)生精度損失，故對(duì)高階矩而言，不可避免地會(huì)存在迭代累計(jì)誤差，進(jìn)而影響圖像重建精度．針對(duì)Zernike矩在笛卡爾坐標(biāo)下存在的計(jì)算精度問(wèn)題，Xin等［8-9］通過(guò)像素重排的方式在極坐標(biāo)下計(jì)算Zernike矩，消除了幾何誤差和數(shù)值積分誤差，極大地提高了Zernike矩的計(jì)算精度．

Yap 等［10］提出了極坐標(biāo)諧波變換(PHT)，與Zernike、偽Zernike矩相比，PHT保留了正交性、旋轉(zhuǎn)不變性，且核函數(shù)形式簡(jiǎn)單，同時(shí)具有數(shù)值穩(wěn)定性．但PHT在極坐標(biāo)下定義，計(jì)算圖像的PHT矩時(shí)需將數(shù)字圖像由笛卡爾坐標(biāo)轉(zhuǎn)換到極坐標(biāo)，坐標(biāo)轉(zhuǎn)換離散化時(shí)會(huì)產(chǎn)生計(jì)算誤差．在計(jì)算PHT內(nèi)核系數(shù)時(shí)需要多次計(jì)算三角函數(shù)，從而使PHT計(jì)算速度較慢．Yang等［11］利用三角函數(shù)的對(duì)稱性減少PHT內(nèi)核系數(shù)的計(jì)算速度，但沒(méi)有解決離散PHT計(jì)算誤差問(wèn)題，且在計(jì)算PHT內(nèi)核系數(shù)時(shí)仍需多次三角函數(shù)計(jì)算．

針對(duì)上述問(wèn)題，文中將以極坐標(biāo)取代笛卡爾坐標(biāo)計(jì)算Zernike矩的方法推廣到PHT矩的計(jì)算中，提出在極坐標(biāo)下具有縮放旋轉(zhuǎn)不變的PHT矩計(jì)算方法，并對(duì)極坐標(biāo)進(jìn)行設(shè)計(jì)和規(guī)劃，以消除笛卡爾坐標(biāo)下PHT矩計(jì)算中存在的幾何誤差和積分近似誤差．為了提高運(yùn)算速度，在計(jì)算過(guò)程中利用三角函數(shù)的對(duì)稱性等相關(guān)性質(zhì)減少三角函數(shù)的計(jì)算次數(shù)，并利用查找表的方式進(jìn)一步提高計(jì)算速度，同時(shí)消除迭代累計(jì)誤差，提高圖像的重建精度．根據(jù)不同變換核，極坐標(biāo)諧波變換可分為極復(fù)指數(shù)變換(PCET)、極正弦變換(PST)和極余弦變換(PCT)3種變換．文中在算法推導(dǎo)的過(guò)程中選用PCET．

1 PHT極坐標(biāo)諧波變換

1．1 極坐標(biāo)諧波變換定義

對(duì)于一幅連續(xù)域上的圖像，將橫坐標(biāo)和縱坐標(biāo)投影到［-1，1］區(qū)間．PHT 在極坐標(biāo)下的定義式為［10］

式(1)中，［·］*表示復(fù)共軛，Mnl稱為圖像 f的 PHT(n，l)矩．式(2)表明基函數(shù) Hnl(r，θ)能分解成半徑分量Rn(r)和角度分量eilθ，n表示階數(shù)，l表示重復(fù)數(shù)，n、l為任意整數(shù)，==0，1，…，∞．極復(fù)指數(shù)變換PCET中

極余弦變換PCT中，Rn(r)=sin(nr2);極正弦變換PST中，Rn(r)=cos(nr2)．

PHT基函數(shù)具有正交性，且任何圖像都可用式(4)重建:

如果僅用n，l的一個(gè)子集重建，可以得到圖像的一個(gè)近似，隨著n、l增加可得到更加精確的重建圖像．

1．2 PCET 的離散計(jì)算

數(shù)字圖像都是離散的，因此需將式(1)轉(zhuǎn)換成離散形式．由于PCET是定義在極坐標(biāo)下的，在迪卡爾坐標(biāo)下定義PCET有如下形式:

式中，H'nl(x，y)=H'nl(rcosθ，rsinθ)≡Hnl(r，θ)，f'(x，y)=f(rcosθ，rsinθ)≡f(r，θ)．

假設(shè)圖f是一幅N×N的離散圖像，其定義域?yàn)?g［i，j］，其中 i=1，2，…，N，j=1，2，…，N．將圖像映射到區(qū)間(xi，yj)∈［-1，1］×［-1，1］．其中，xi=(2i-N-1)/N，yj=(2j-N-1)/N，則有

其中，xi，yj滿足+≤1，

其中Δ=2/N是像素的寬度．在笛卡爾坐標(biāo)系下無(wú)法得到式(7)的準(zhǔn)確值，只能由式(8)得到估計(jì)值，如此這就產(chǎn)生了數(shù)值積分誤差:

將估計(jì)值代入式(6)，得到PCET的離散形式［10］:

式(9)是笛卡爾坐標(biāo)下離散PCET的計(jì)算公式．利用式(9)計(jì)算PCET時(shí)，會(huì)產(chǎn)生數(shù)值積分誤差．另外，在計(jì)算PCET時(shí)如按數(shù)字圖像的像素離散時(shí)會(huì)產(chǎn)生幾何誤差，這與Zernike矩相似．幾何誤差產(chǎn)生的原因是正方形的單位像素組成的塊不可能恰好覆蓋一個(gè)單位圓．只要采用笛卡爾坐標(biāo)系，幾何誤差就是不可避免的，因此，要提高PCET的計(jì)算精度就需要減少幾何誤差和數(shù)值積分誤差．

2 快速極坐標(biāo)下的PCET離散計(jì)算

2．1 極坐標(biāo)劃分與像素插值

若圖像大小為N×N，則需計(jì)算PCET的圓型區(qū)域的半徑為N/2．先將圓等分成V個(gè)扇區(qū)，沿徑向?qū)A分成U層圓環(huán)區(qū)域，依次將第 u(u=1，2，3，…，U)層的V個(gè)扇區(qū)分別等分成為2u-1，個(gè)塊則第u個(gè)圓環(huán)被等分成V(2u-1)個(gè)像素塊．需計(jì)算PCET的圓型區(qū)域被劃分得到VU2個(gè)像素塊．文中取U=N/2，當(dāng)V=4時(shí)劃分結(jié)果如圖1(a)所示．每個(gè)扇區(qū)表示為 Ωuv，如圖 1(b)所示，u=1，2，…，U;v=1，2，…，V(2u-1)．扇區(qū) Ωuv表示極坐標(biāo)像素(ρuv，θuv)和分別表示 Ωuv的開(kāi)始和結(jié)束半徑和分別表示 Ωuv的開(kāi)始和結(jié)束角度，ρuv=(+)/2，=+)/2．

圖1 笛卡爾坐標(biāo)與極坐標(biāo)下計(jì)算PHT的像素分布Fig．1 Pixels distribution in Cartesian coordinates and polar coordinates for PHT computation

圖像單位圓被分成VU2個(gè)扇區(qū)，每個(gè)扇區(qū)的大小是/(VU2)，U和V的取值要適當(dāng)．VU2小則計(jì)算量小，對(duì)圖像信息的表達(dá)誤差大;反之則計(jì)算量大，對(duì)圖像信息的表達(dá)誤差小．

極坐標(biāo)中的像素位置(ρuv，θuv)與笛卡爾坐標(biāo)中的像素位置不吻合，需要通過(guò)插值的方式來(lái)解決這個(gè)問(wèn)題．文中采用雙三次插值方法，雙三次插值方法具有三階精度，可以滿足式(6)中定義的估計(jì)方法．對(duì)劃分后的像素塊Ωuv使用式(10)插值得到極坐標(biāo)下的圖像f'各像素的值:

式中:k=「ρuvcosθuv/Δ? +N/2;l=「ρuvsinθuv/Δ? +N/2;f為笛卡爾坐標(biāo)下的圖像;h(·)是雙三次插值的1-D核函數(shù)，

2．2 極坐標(biāo)下的PCET計(jì)算

由式(6)、(7)可知，插值得到的圖像f'的PCET可利用式(11)進(jìn)行計(jì)算:

式(13)是PCET基的半徑分量計(jì)算公式，式(14)是PCET基的角度分量計(jì)算公式．由復(fù)數(shù)與三角函數(shù)的關(guān)系，式(13)中，

式(14)中，

直接利用式(13)-(16)計(jì)算一個(gè)像素點(diǎn)的PCET基的半徑分量和角度分量時(shí)，分別需要4次三角函數(shù)計(jì)算、1次復(fù)數(shù)加法，這使得PCET的計(jì)算時(shí)間較長(zhǎng)．

2．3 PCET基的快速計(jì)算

借鑒快速傅里葉變換(FFT)的思想，根據(jù)三角函數(shù)的周期性和對(duì)稱性可提高PCET基的計(jì)算速度．cos(lθ)和 sin(lθ)是周期為 2/l的周期函數(shù)，由三角函數(shù)的周期性和對(duì)稱性可得:

同理可得表1．

表1 極坐標(biāo)像素對(duì)稱點(diǎn)的三角函數(shù)關(guān)系Table 1 Trigonometric relations of symmetric pixels in polar coordinates

表1描述了如圖2所示的8個(gè)對(duì)稱點(diǎn)a1，a2，…，a8的三角函數(shù)關(guān)系:

圖2 極坐標(biāo)像素對(duì)稱點(diǎn)Fig．2 Symmetric pixels in polar coordinates

由表1和式(14)、(15)可得，當(dāng)l mod 4=0時(shí)，

利用對(duì)稱性計(jì)算PCET基可以減少三角函數(shù)的計(jì)算次數(shù)，由式(19)可見(jiàn)，計(jì)算8個(gè)像素點(diǎn)的PCET基的角度分量時(shí)需要4次三角函數(shù)計(jì)算、1次復(fù)數(shù)加法．

2．3．1 PCET基中角度分量的快速計(jì)算在極坐標(biāo)劃分時(shí)，劃分扇形區(qū)域V取4，如圖1(a)所示，因而對(duì)于1≤u≤U，1≤v≤4(2u-1)，有

將式(23)代入(19)中，可得

計(jì)算 sin(vβu)、cos(vβu)(1≤u≤U，1≤v≤4(2u-1))，存儲(chǔ)在2個(gè)矩陣中，記為 A、B．A、B 中的元素 Au，v，Bu，v分別表示 sin(vβu)、cos(vβu)．當(dāng) l＞ 0時(shí)，sin(lvβu)、cos(lvβu)分別在矢量 Au和 Bu中找第lv個(gè)單元的值．第u個(gè)環(huán)被分成4(2u-1)個(gè)扇區(qū)，因此lv＞4(2u-1)時(shí)，lv需要對(duì) 4(2u-1)取模，則sin(lvβu)、cos(lvβu)分別在 Au、Bu中是第 lv mod 4(2u-1)個(gè)元素 Au，lvmod4(2u-1)，Bu，lvmod4(2u-1)．如果l＜0，因?yàn)?sin(- α)=-sinα，cos(- α)=cosα，則 sin(lvβu)為 -Au，lvmod4(2u-1)，cos(lvβu)為Bu，lvmod4(2u-1)．利用查找表，式(24)可表示為

其中，sgnl表示l的符號(hào)．

同理，當(dāng)l mod 4=1時(shí)，

當(dāng)l mod 4=2時(shí)，

當(dāng)l mod 4=3時(shí)，

當(dāng)v=u，0 ＜u＜U 時(shí)，對(duì)稱點(diǎn)為4個(gè)，分別是a1，a3，a5，a7．則利用式(25)-(28)時(shí)，將 a2，a4，a6，a8均視為0．

2．3．2 PCET基中徑向分量的快速計(jì)算

與角度分量的計(jì)算相同，半徑分量也可用查找表的方式快速計(jì)算．令γu=2u/U2，將式(22)代入式(15)中，則有

計(jì)算 sinγk、cosγk(1≤k≤U2)，存儲(chǔ)在 2 個(gè)一維向量中，記為 C、D．其中元素 Cu、Du分別表示sinγu、cosγu．當(dāng) n ＞0 時(shí)，sin(nuγu)、cos(nuγu)分別在矢量C、D中找第 nu2個(gè)單元的值．當(dāng) nu2＞U2時(shí)，nu2需要對(duì) U2取模，則 sin(nuγu)、cos(nuγu)分別在 C、D 中是元素 Cnu2modU2、Dnu2modU2．當(dāng) n ＜0 時(shí)，sin(nuγu)為-Cnu2modU2，cos(nuγu)為 Dnu2modU2．利用查找表，式(29)表示為

3 實(shí)驗(yàn)結(jié)果與分析

3．1 極坐標(biāo)下計(jì)算PCET的重建精度

對(duì)大小為256×256的Lena圖像分別在笛卡爾坐標(biāo)下和極坐標(biāo)下計(jì)算PCET后對(duì)圖像進(jìn)行重建，計(jì)算PCET時(shí)，≤≤K．K 分別為60、70、80 的基于笛卡爾坐標(biāo)的PCET重建圖像和基于極坐標(biāo)的PCET重建結(jié)果如圖3、4所示．從圖中可以看出，基于笛卡爾坐標(biāo)的PCET重建圖像沿單位圓邊界的一些錯(cuò)誤的像素是非常明顯的，隨著K的增大，錯(cuò)誤像素?cái)?shù)量迅速增加．然而，這種錯(cuò)誤的像素在基于極坐標(biāo)的PCET重建圖像中不存在．

圖3 基于笛卡爾坐標(biāo)的PCET重建結(jié)果Fig．3 Image reconstruction from PCETs based on Cartesian method

圖4 基于文中的極坐標(biāo)的PCET重建結(jié)果Fig．4 Image reconstruction from PCETs based on the proposed polar method

K為10～100時(shí)，重建圖像的峰值信噪比(PSNR)如圖5所示．由圖5可見(jiàn)，K較小時(shí)(K＜60)，基于笛卡爾坐標(biāo)的PCET重建圖像與基于極坐標(biāo)的PCET重建圖像的質(zhì)量相類似，重建圖像的PSNR隨著K增大而穩(wěn)步增長(zhǎng)．總體來(lái)說(shuō)，基于極坐標(biāo)的PCET重建圖像的PSNR要高約0．4dB．但當(dāng)K較大時(shí)，基于極坐標(biāo)的PCET重建圖像質(zhì)量明顯高于基于笛卡爾坐標(biāo)的PCET重建圖像．基于極坐標(biāo)的PCET重建圖像質(zhì)量隨K單調(diào)遞增．而在基于笛卡爾坐標(biāo)的PCET重建圖像中，當(dāng)K增大到某一點(diǎn)(大約為60時(shí))，圖像質(zhì)量達(dá)到最大值．當(dāng)K＞60時(shí)，隨著K值的增大，圖像質(zhì)量下降．這是因?yàn)橹亟ㄕ`差(包括幾何和數(shù)值積分誤差)隨著K的增加而增大，并在某些時(shí)候超過(guò)了因增加PCET矩而帶來(lái)的重建圖像質(zhì)量增益．

圖5 PCET重建結(jié)果的PSNRFig．5 PSNR of image reconstruction from PCETs

3．2 極坐標(biāo)下計(jì)算PCET的旋轉(zhuǎn)和縮放不變性

對(duì)256×256的Lena圖像進(jìn)行以下變換后分析不變性:(1)對(duì)圖像進(jìn)行旋轉(zhuǎn)，旋轉(zhuǎn)角度分別為5°、30°、60°、120°、135°;(2)對(duì)圖像進(jìn)行縮放，縮放因子分別為 300%、200%、150%、120%、90% 和 75%;(3)對(duì)圖像進(jìn)行旋轉(zhuǎn)和縮放，先對(duì)圖像旋轉(zhuǎn)30°，然后縮放120%．

對(duì)上述圖像及原圖分別在極坐標(biāo)和笛卡爾坐標(biāo)系下計(jì)算PCET基于極坐標(biāo)的原圖的PCET矩與不同變換類型的PCET矩平均如圖6所示．由圖6可見(jiàn)，變換后圖像的PCET矩與原圖的PCET矩基本重合，說(shuō)明基于極坐標(biāo)的PCET具有很好的旋轉(zhuǎn)和縮放不變性．

變換后圖像PCET矩與原始圖像PCET矩的誤差如圖7所示．由圖7可見(jiàn)，基于極坐標(biāo)的PCET矩比基于笛卡爾坐標(biāo)的PCET矩具有更好的縮放、旋轉(zhuǎn)不變性．

圖6 文中的極坐標(biāo)下PCET的旋轉(zhuǎn)和縮放不變性Fig．6 Rotation and scale invariance of PCET moments based on the proposed polar method

圖7 極坐標(biāo)與笛卡爾坐標(biāo)PCET的旋轉(zhuǎn)和縮放不變性Fig．7 Rotation and scale invariance of PCET moments based on Cartesian method and proposed polar method

3．3 速度和存儲(chǔ)空間分析

基于極坐標(biāo)的PCET快速計(jì)算方法，其速度提高的原因在于:(1)減少了三角函數(shù)和復(fù)數(shù)運(yùn)算的次數(shù);(2)利用查找表提高了速度．用Yap等［8］的方法直接在笛卡爾坐標(biāo)下計(jì)算PCET，為8個(gè)點(diǎn)生成內(nèi)核系數(shù)需要計(jì)算32次三角函數(shù)，16次復(fù)數(shù)加法，在極坐標(biāo)中用式(11)-(14)直接計(jì)算PCET，為8個(gè)點(diǎn)生成內(nèi)核系數(shù)需要計(jì)算64次三角函數(shù)、16次復(fù)數(shù)加法、8次復(fù)數(shù)乘法．在文中的快速算法中，為8個(gè)對(duì)稱點(diǎn)生成內(nèi)核系數(shù)需要計(jì)算2次復(fù)數(shù)加法、1次復(fù)數(shù)乘法，三角函數(shù)的計(jì)算轉(zhuǎn)換成了查找表的方式，因此計(jì)算速度比上述兩種方法直接計(jì)算PCET有很大提高．文中使用不同分辨率和內(nèi)容的圖像進(jìn)行實(shí)驗(yàn)，圖8給出了在PC環(huán)境(英特爾2．13 GHz的，2 G內(nèi)存)，K=10，20，…，100時(shí)3種方法的計(jì)算時(shí)間，算法實(shí)現(xiàn)采用Matlab 7．0．從圖中可見(jiàn):文中提出的快速算法的計(jì)算時(shí)間相比于直接的笛卡爾坐標(biāo)和直接的極坐標(biāo)下的計(jì)算時(shí)間都明顯減少;計(jì)算速度比笛卡爾坐標(biāo)的PCET計(jì)算方法快6倍以上;隨著K增大，計(jì)算速度提高越明顯．這主要是文中算法計(jì)算PCET內(nèi)核系數(shù)時(shí)利用查找表，無(wú)論K增加多大，計(jì)算時(shí)間都不會(huì)快速增加．

圖8 PCET的計(jì)算時(shí)間比較Fig．8 Comparison of computation time for PCET

4 結(jié)語(yǔ)

文中將以極坐標(biāo)取代笛卡爾坐標(biāo)計(jì)算Zernike矩的方法推廣到PHT矩的計(jì)算中，提出在極坐標(biāo)下具有縮放旋轉(zhuǎn)不變的PHT矩計(jì)算方法，并對(duì)極坐標(biāo)進(jìn)行設(shè)計(jì)和規(guī)劃，以消除笛卡爾坐標(biāo)下PHT矩計(jì)算中存在的幾何誤差和積分近似誤差．為了提高運(yùn)算速度，在計(jì)算過(guò)程中利用三角函數(shù)的對(duì)稱性等相關(guān)性質(zhì)減少三角函數(shù)的計(jì)算次數(shù)，并利用查找表的方式計(jì)算矩值，進(jìn)一步提高計(jì)算速度的同時(shí)消除迭代累計(jì)誤差，提高圖像的重建精度．文中還通過(guò)實(shí)驗(yàn)對(duì)提出的PCET方法進(jìn)行了驗(yàn)證，結(jié)果表明，該方法在重建精度和不變性上都優(yōu)于笛卡爾坐標(biāo)的PCET計(jì)算方法，計(jì)算速度比笛卡爾坐標(biāo)的PCET計(jì)算方法快6倍以上．提出的快速算法可以推廣到PCT和PST的計(jì)算中．

［1］Liao S，Pawlak M．On the accuracy of Zernike moments for image analysis［J］．IEEE Transactions on Pattern Analysis and Machine Intelligence，1998，20(12):1358-1364．

［2］Xin Y，Liao S，Pawlak M．Circularly orthogonal moments for geometrically robust image watermarking［J］．Pattern Recognition，2007，40(12):3740-3752．

［3］Gao Xinbo，Wang Qian，Li Xuelong，et al．Zernike-moment-based image super resolution［J］．IEEE Transactions on Image Processing，2011，20(10):2738-2747．

［4］Revaud J，Lavoue G，Baskurt A．Improving Zernike moments comparison for optimal similarity and rotation angle retrieval［J］．IEEE Transactions on Pattern Analysis and Machine Intelligence，2009，31(4):627-636．

［5］Li S，Lee M C，Pun C M．Complex Zernike moments features for shape-based image retrieval［J］．IEEE Transactions on Systems，Man and Cybernetics(Part A):Systems and Humans，2009，39(1):227-237．

［6］Chen Z，Sun S K．A Zernike moment phase-based descriptor for local image representation and matching［J］．IEEE Transactions on Image Processing，2010，9(1):205-219．

［7］Singh C，Walia E，Mittal N．Rotation invariant complex Zernike moments features and their applications to human face and character recognition［J］．IET Computer Vision，2011，5(5):255-265．

［8］Xin Y，Pawlak M，Liao S．Accurate computation of Zernike moments in polar coordinates［J］．IEEE Transactions on Image Processing，2007，16(2):581-587．

［9］Singh C，Walia E．Computation of Zernike moments in improved polar configuration ［J］．IET Image Processing，2009，3(4):217-227．

［10］Yap P，Jiang X，Kot A C．Two-dimensional polar harmonic transforms forinvariantimagerepresentation ［J］．IEEE Transactions on Pattern Analysis and Machine Intelligence 2010，32(7):1260-1270．

［11］Yang Zhuo，Kamata S．Fast polar harmonic transforms［C］∥2010 11th International Conference on Control Automation Robotics＆ Vision(ICARCV)．Singapore:IEEE Press，2010:673-677．