張桂梅 孫曉旭 劉建新 儲珺

基于分?jǐn)?shù)階微分的TV-L1光流模型的圖像配準(zhǔn)方法研究

張桂梅1孫曉旭1劉建新2儲珺1

圖像的非剛性配準(zhǔn)在計算機視覺和醫(yī)學(xué)圖像分析中有著重要的作用.TV-L1(全變分L1范數(shù)、Total variation-L1)光流模型是解決非剛性配準(zhǔn)問題的有效方法,但TV-L1光流模型的正則項是一階導(dǎo)數(shù),會導(dǎo)致紋理特征等具有弱導(dǎo)數(shù)性質(zhì)的信息模糊.針對該問題,將G-L(Grnwald-Letnikov)分?jǐn)?shù)階引入TV-L1光流模型,提出基于G-L分?jǐn)?shù)階微分的TV-L1光流模型,并應(yīng)用原始–對偶算法求解該模型.新的模型用G-L分?jǐn)?shù)階微分代替正則項中的一階導(dǎo)數(shù),由于分?jǐn)?shù)階微分比整數(shù)階微分具有更好的細(xì)節(jié)描述能力,并能有效地、非線性地保留具有弱導(dǎo)數(shù)性質(zhì)的紋理特征,從而提高圖像的配準(zhǔn)精度.另外,通過實驗給出了配準(zhǔn)精度與G-L分?jǐn)?shù)階模板參數(shù)之間的關(guān)系,從而為模板最佳參數(shù)的選取提供了依據(jù).盡管不同類型的圖像其最佳參數(shù)是不同的,但是其最佳配準(zhǔn)階次一般在1～2之間.理論分析和實驗結(jié)果均表明,提出的新模型能夠有效地提高圖像配準(zhǔn)的精度,適合于包含較多弱紋理和弱邊緣信息的醫(yī)學(xué)圖像配準(zhǔn),該模型是TV-L1光流模型的重要延伸和推廣.

分?jǐn)?shù)階微分,Grnwald-Letnikov,TV-L1模型,光流場,弱紋理,非剛性配準(zhǔn)

分?jǐn)?shù)階微積分理論在圖像處理中的應(yīng)用在最近幾年才引起學(xué)者的關(guān)注,越來越多的學(xué)者把分?jǐn)?shù)階微積分應(yīng)用到圖像處理中,并在圖像去噪、圖像增強、小波變換、圖像奇異性特征提取、圖像融合、圖像分割、圖像邊緣提取等方面取得了一些初步成果.文獻(xiàn)[1]認(rèn)為分?jǐn)?shù)階微分算子選擇合適的階次,在增強圖像的過程中可以大幅提升邊緣和紋理細(xì)節(jié),并且可以非線性保留圖像平滑區(qū)域的紋理信息.文獻(xiàn)[2]利用分?jǐn)?shù)階微分階次在(0,1)區(qū)間時,具有弱導(dǎo)數(shù)的性質(zhì),利用分?jǐn)?shù)階微分算子對含弱噪聲的圖像進(jìn)行邊緣檢測,成功解決了整數(shù)階梯度算子對噪聲敏感的問題,避免了噪聲的影響,準(zhǔn)確檢測噪聲圖像的邊緣.文獻(xiàn)[3]將分?jǐn)?shù)階積分理論引入圖像去噪,構(gòu)造分?jǐn)?shù)階積分模板對圖像進(jìn)行濾波去噪,該方法能夠在提高峰值信噪比的同時更好地保留圖像邊緣和紋理信息.文獻(xiàn)[4]將分?jǐn)?shù)階微積分引入奇異值分解用于人臉識別,該方法在人臉變化劇烈時相對于傳統(tǒng)方法具有更好的分類效果.文獻(xiàn)[5]將分?jǐn)?shù)階微分引入全變分模型用于圖像修復(fù)領(lǐng)域,該方法對空域和頻域損壞的圖像,特別是對圖像細(xì)節(jié)信息的修復(fù),都具有更好的修復(fù)效果,在提高視覺效果的同時還提高了峰值信噪比.文獻(xiàn)[6]將分?jǐn)?shù)階微分引入圖像分割領(lǐng)域,提出一種具有分?jǐn)?shù)階次擬合項的活動輪廓模型,分?jǐn)?shù)階擬合項能夠更加準(zhǔn)確地描述圖像,并且對噪聲具有魯棒性,該方法對噪聲圖像、醫(yī)學(xué)圖像和紅外圖像都有很好分割效果.受以上文獻(xiàn)啟發(fā),分?jǐn)?shù)階微積分能夠非線性地增強或者保留圖像的紋理細(xì)節(jié)信息,本文將分?jǐn)?shù)階微積分引入圖像配準(zhǔn)領(lǐng)域.

圖像配準(zhǔn)是圖像處理的一個基本問題,用于將不同時間、不同傳感器、不同視角及不同拍攝條件下獲取的兩幅或者多幅進(jìn)行匹配,目的是尋求一種空間變換,使一幅圖像與另一幅圖像上的對應(yīng)點達(dá)到空間上的一致,用以糾正圖像的形變.圖像配準(zhǔn)是近年來圖像處理領(lǐng)域的研究方向之一,在遙感圖像處理、醫(yī)學(xué)圖像分析,計算機視覺等領(lǐng)域有著廣泛的應(yīng)用.圖像配準(zhǔn)的方法根據(jù)待配準(zhǔn)目標(biāo)的類型可以將其分為剛性配準(zhǔn)和非剛性配準(zhǔn).剛性配準(zhǔn)只適用于不存在變形的配準(zhǔn),但現(xiàn)實生活中大多數(shù)形變是非剛性的,如在醫(yī)學(xué)圖像處理領(lǐng)域,大多數(shù)情況下,人體組織的形變是非剛性的,因而非剛性配準(zhǔn)更具有重要研究意義[7?9].目前最著名的非剛性配準(zhǔn)是基于光流場模型的Demons算法[7?9].其將圖像的配準(zhǔn)過程可以看作從源圖像流動到目標(biāo)圖像的過程,即配準(zhǔn)所求解的位移場可以看作光流場模型求解的速度場,因此可以通過求解光流場進(jìn)行圖像配準(zhǔn).如文獻(xiàn)[10]提出基于光流場理論的Demons算法,其思想是將配準(zhǔn)看作浮動圖像像素點在參考圖像灰度梯度驅(qū)動下向參考圖像移動的過程,文獻(xiàn)[11]在經(jīng)典Demons算法基礎(chǔ)上提出了主動Demons算法,該算法利用參考圖像和浮動圖像的灰度梯度作為共同驅(qū)動力驅(qū)動像素點向?qū)Ψ綄?yīng)像素點移動,配準(zhǔn)精度有所提高,但該類算法在圖像灰度均勻及弱紋理區(qū)域配準(zhǔn)精度較低,優(yōu)化易陷入局部最小從而導(dǎo)致配準(zhǔn)速度緩慢.文獻(xiàn)[12]將Horn-Schunck(H-S)光流場模型引入圖像配準(zhǔn),由于Horn-Schunck模型采用的光滑性約束無法保持圖像的不連續(xù)性,而在醫(yī)學(xué)圖像配準(zhǔn)中,允許位移場的不連續(xù)性是必要的,如呼吸運動,膈膜產(chǎn)生嚴(yán)重變形,而肋骨仍然是剛性形變,這時光滑位移場就不足以適用這種復(fù)雜的運動.另外,Horn-Schunck模型的二次數(shù)據(jù)項對亮度異常值不具有足夠的魯棒性.針對上述缺點,Pock等[13]提出了基于變分方法的TV-L1(Total variation-L1)光流模型,該模型能夠保持圖像的不連續(xù)性,即在圖像演化過程中可以有效保持圖像的邊緣等特征信息,而且數(shù)據(jù)項采用光流的L1范數(shù)增強了配準(zhǔn)的魯棒性.文獻(xiàn)[14]則具體描述了文獻(xiàn)[13]方法的實現(xiàn)過程.文獻(xiàn)[15]應(yīng)用可變形配準(zhǔn)方法得到比剛性配準(zhǔn)、快速Demons配準(zhǔn)和快速自由形狀配準(zhǔn)算法等更準(zhǔn)確的效果.但是文獻(xiàn)[13?15]對具有弱導(dǎo)數(shù)性質(zhì)的紋理特征等信息的保持仍然不夠理想.

雖然TV-L1模型相對H-S模型有了很大提升,能夠在擴散過程中保持邊緣信息,但對具有弱導(dǎo)數(shù)性質(zhì)的紋理細(xì)節(jié)等信息的保持仍不夠理想,并且在模型最小化過程中造成計算上的困難,主要原因是正則項和數(shù)據(jù)項在零點處都不可微,一個有效的方法是在零點處增加了一個極小值,然而新引入的項會造成收斂很慢并且會模糊位移場.本文的研究對象主要為存在弱邊緣和弱紋理等特征的圖像,如醫(yī)學(xué)圖像,在成像過程中由于受到電磁環(huán)境,成像設(shè)備,以及人體結(jié)構(gòu)的差異等客觀因素的影響,會造成圖像存在弱紋理、弱邊緣,所以要求配準(zhǔn)方法能較好保持圖像的邊緣和紋理特征.而分?jǐn)?shù)階微分算子具有弱導(dǎo)數(shù)的性質(zhì),通過合理選擇其階次,可以達(dá)到保持弱邊緣和弱紋理的特性,基于此,本文嘗試將G-L分?jǐn)?shù)階微分理論引入到TV-L1光流模型中,用于圖像的非剛性配準(zhǔn),以提高該模型的配準(zhǔn)精度.

1 基本理論

1.1 Horn-Schunck(H-S)光流場模型[16]

假設(shè)圖像上的像素點(x,y)在t時刻的灰度值為I(x,y,t),經(jīng)過較小的Δt時刻后該點位置為(x+Δx,y+Δy),該點灰度值為I(x+Δx,y+Δy,t+Δt),根據(jù)光流場的亮度恒定假設(shè)得到光流約束方程[16]

即Ixu+Iyv+It=0,其中

在圖像配準(zhǔn)中,只取Δt時刻前后的兩幀圖像,U=(u,v)T為求解的速度場矢量,即配準(zhǔn)要求解的位移場矢量,但是位移場矢量有兩個變量,而通過亮度恒定假設(shè)只能得到一個約束方程,不足以求解x和y方向的位移u和v,這就是光流場求解的病態(tài)問題.因此要求解光流場模型還需對位移場矢量附加另外的約束.

Horn等創(chuàng)造性地在光流約束的基礎(chǔ)上加上光滑性約束,該約束假設(shè)物體的運動速度在大多數(shù)情況下是局部光滑的,特別是目標(biāo)在無形變的剛性運動時,各鄰近像素點應(yīng)具有相同的運動速度,即相鄰點的速度變化率為0,因而將光流場的求解轉(zhuǎn)化為能量泛函極小值問題:

第一項為正則項,懲罰光流場中大的變化以獲得光滑的速度場,第二項為數(shù)據(jù)項,即基本光流約束,假設(shè)對應(yīng)點運動前后灰度值不變,λ為正則化參數(shù),是與正則項和數(shù)據(jù)項相關(guān)的一個權(quán)重參數(shù).

由于H-S模型正則項采用光滑性約束,導(dǎo)致無法保持位移場的不連續(xù)性,會在圖像的演化過程中產(chǎn)生嚴(yán)重模糊而丟失重要信息,數(shù)據(jù)項采用二次項懲罰偏差,對處理亮度異常值不具有魯棒性.

1.2 TV-L1光流場模型

為了克服Horn-Schunck光流模型的缺點,Pock[13]等提出了基于變分方法的TV-L1光流模型,并且由于在實際問題中,光流約束并不會一直成立,即圖像亮度恒定假設(shè)對于很多真實圖像并不適用,如在光照條件變化大,圖像亮度恒定假設(shè)存在較大誤差,Pock等引入變量w對光流約束加以改進(jìn),通過w對光照變化建模得到:

其中,ρ(u,v,w)=Ix(u?u0)+Iy(v?v0)+It+βw,參數(shù)β為光照變化項權(quán)重因子,u0和v0為已知的位移場,該模型用非二次正則項|?u|+|?v|替換了H-S模型的正則項|?u|2+|?v|2,并用L1范數(shù)的數(shù)據(jù)項|ρ(u,v,w)|替換了H-S模型的數(shù)據(jù)項(Ixu+Iyv+It)2.雖然與H-S模型相比改動較小,但是配準(zhǔn)精度有了很大的提升.首先全變分正則項保持了位移場的不連續(xù)性,在擴散過程中保護(hù)邊緣信息不被模糊,這個替換后的正則項與與著名的ROF(Rudin、Osher、Fatemi)模型[17]的正則項相同,而ROF模型在具有很好去噪效果的同時還能保持邊緣信息不被模糊,正是因為采用了非二次正則項.其次,數(shù)據(jù)項使用穩(wěn)健的L1范數(shù)對亮度的變化不及H-S模型敏感.

雖然TV-L1模型相對H-S模型有了很大提升,能夠在擴散過程中保持邊緣信息,但對具有弱導(dǎo)數(shù)性質(zhì)的紋理細(xì)節(jié)等信息的保持仍不夠理想,并且在模型最小化過程中導(dǎo)致了計算上的困難,主要原因是正則項和數(shù)據(jù)項在零點處都不可微,一個有效的方法是用可微的來近似替代|?u|,然而,極小值ε會導(dǎo)致收斂很慢并且會模糊位移場.本文嘗試將分?jǐn)?shù)階微分來代替TV-L1模型中的一階導(dǎo)數(shù),并采用文獻(xiàn)[18]的方法來求解基于分?jǐn)?shù)階微分的TV-L1光流場模型.

2 基于G-L分?jǐn)?shù)階微分的TV-L1光流場模型

分?jǐn)?shù)階微積分在17世紀(jì)提出,在19世紀(jì)得到了迅速發(fā)展,期間提出了許多關(guān)于分?jǐn)?shù)階微積分的定義,最著名的有R-L(Riemann-Liouville)定義、G-L(Grnwald-Letnikov)定義和Caputo定義.其中G-L和R-L定義在數(shù)值運算時可以轉(zhuǎn)化為卷積積分形式,因此適用于信號處理領(lǐng)域,而Caputo定義適用于分?jǐn)?shù)階微分方程的初邊值問題,多用于實際工程領(lǐng)域.而G-L定義在計算時較R-L定義更加準(zhǔn)確,因此我們從G-L定義出發(fā)推導(dǎo)出本文算法所需要的分?jǐn)?shù)階微分算子,對TV-L1算法進(jìn)行改進(jìn),用于圖像的非剛性配準(zhǔn).

2.1 離散分?jǐn)?shù)階微分算子

式中,f(x)是實函數(shù),為廣義二項式系數(shù),Γ(x) 為 Gamma 函數(shù),m為正整數(shù),m值越大,二項式系數(shù)值越趨于0.當(dāng)α=1,m≥2時,cαm=0,式(4)變?yōu)橐浑A差分,即為傳統(tǒng)梯度算子.

為簡化起見,所用圖像尺寸均為N×N 二維矩陣,用X 表示歐幾里得空間RN×N,對于X 中的任意變量定義X 中的內(nèi)積:

定義矢量空間Y=X×X,分?jǐn)?shù)階梯度算子?αu為矢量空間Y中的一個矢量:

其中

對于Y中任意的變量p=(p1,p2),分?jǐn)?shù)階散度算子的離散形式為

定義Y中的內(nèi)積:

容易得到式(5)和式(9)滿足如下性質(zhì):

為方便計算,本文根據(jù)G-L定義的分?jǐn)?shù)階導(dǎo)數(shù)在特定離散情況下與某種卷積積分的離散化形式相同這一個特點,將分?jǐn)?shù)階導(dǎo)數(shù)轉(zhuǎn)化為卷積積分的形式.若記卷積核函數(shù)為

則式(7)可看作為式(12)的一個離散化近似,

這樣算法中涉及分?jǐn)?shù)階微分的運算便轉(zhuǎn)化為卷積與相關(guān)運算.在具體計算時,采用4個方向的分?jǐn)?shù)階微分掩模算子W1、W2、W3、W4,從而可以控制圖像的上下左右4個方向的擴散程度,W1如圖1所示.其中W1 是沿x負(fù)方向的掩模算子,W2～W4依次順時針與x負(fù)方向夾角為 90°、180°、270°.

圖1 分?jǐn)?shù)階微分模板Fig.1 Diあerential template

2.2 分?jǐn)?shù)階TV-L1光流場模型

本文將G-L分?jǐn)?shù)階微分引入TV-L1光流模型,代替其中的一階微分,得到分?jǐn)?shù)階TV-L1光流場模型(Fractional TV-L1,記為FTV-L1):

將模型轉(zhuǎn)化為離散形式得到:

其中

參數(shù)α是個實數(shù),當(dāng)α=1時,該分?jǐn)?shù)階TV-L1光流場模型即可看作文獻(xiàn)[13]提出的一階TV-L1光流場模型,所以分?jǐn)?shù)階TV-L1光流場模型實際上是整數(shù)階TV-L1光流場模型的一個擴展和延伸.

2.3 分?jǐn)?shù)階TV-L1模型的求解

文獻(xiàn)[18]提出了有效且準(zhǔn)確的一階原始–對偶算法來求解變分問題,且證明了以O(shè)(1/N)的速率收斂.本文應(yīng)用文獻(xiàn)[18]的方法求解本文提出的FTV-L1光流場模型.文獻(xiàn)[18]的TV-L1光流模型迭代方法如下:對于一類問題的原始形式:

先變換為其的原始–對偶形式,再進(jìn)行迭代,得到下式:

其中,F?為F的凸共軛.

步驟1.初始化:τ,σ＞0,0=x0.

步驟2.迭代更新xn,yn,n:

分?jǐn)?shù)階TV-L1光流模型的離散形式為

引入對偶變量p,q,z∈Y,則分?jǐn)?shù)階TV-L1模型的原始–對偶形式如下:

根據(jù)式(12)的性質(zhì),式(25)可變換為

其中,P,Q,Z分別為變量p,q,z的凸集,定義如下:

式中,‖·‖∞表示無窮范數(shù):對‖q‖∞和‖z‖∞同理可得.

凸集P,Q,Z是L2球體內(nèi)的點集,式(26)中函數(shù)δP(p)、δQ(q)和δZ(z)分別為集合P,Q,Z的指示函數(shù),它們的表達(dá)式如下:

由式(21)和(22)知,該迭代方法的關(guān)鍵是求解(I+σ?F?)?1和 (I+τ?G)?1兩個算子. 對比式(26)和式(20)可知,F?(p,q,z)= δP(p)+δQ(q)+δZ(z),G(u,v,w)=λ‖ρ(u,v,w)‖1, 下面將應(yīng)用上述迭代方法來最小化式(15),求解分?jǐn)?shù)階TV-L1模型.具體求解步驟如下:

步驟1.

1)先求解式(21)中的(yn+σKˉxn):

2) 再求解式 (21) 中的 (I+σ?F?)?1,由于F?(p,q,z)為集合P,Q,Z 的指示函數(shù),把集合P,Q,Z內(nèi)的點按歐幾里得距離分別投影至L2球內(nèi):

步驟2.

1)先求解式(22)中的(xn?τK?yn+1):

2)再求解式 (22)中的 (I+τ?G)?1,因為G(u,v,w)=λ‖ρ(u,v,w)‖1是L1范數(shù), 而L1范數(shù)是連續(xù)不光滑的,為方便描述,定義ai,j=(?Ii,j,β),|a|2i,j=|?I|2i,j+β2:

步驟3.返回步驟1.

由于FTV-L1光流模型是個非凸優(yōu)化問題,光流約束僅對小位移有效,為了使其對大位移的情形也適用,本文對上述模型應(yīng)用從粗到精的配準(zhǔn)策略,避免模型在極小化過程中陷入局部極小值.具體方法為:取降采樣因子為2,建立高斯圖像金字塔,對每一層圖像極小化FTV-L1模型(式(15)),求其位移場.從最粗一層開始,直到最精一層.每一層求得的位移場作為下一層圖像求解的初值.第一層,初始化位移場(u,v)為0,光照變化變量w初始化為0,對偶變量(p,q,z)初始化為0.對于每一層金字塔圖像的位移場求解過程,如上述三個迭代步驟.由于本文使用4個方向的分?jǐn)?shù)階微分模板進(jìn)行求導(dǎo)操作,故計算過程中關(guān)于位移場(u,v,w)的離散梯度算子及關(guān)于對偶變量(p,q,z)的離散散度算子變化如下:

用式(33)替代原來的整數(shù)一階前后差分算法.

2.4 實驗參數(shù)的設(shè)置

本文實驗中迭代次數(shù)取50時,模型已經(jīng)達(dá)到收斂,此處迭代次數(shù)設(shè)置為50是因為經(jīng)過多次實驗發(fā)現(xiàn),在迭代次數(shù)為50時,灰度均方差(Mean square error,MSE)和峰值信噪比(Peak singnal to noise ratio,PSNR)均達(dá)到了收斂,故選擇了迭代次數(shù)為50.如圖2(a)和圖2(b)所示為Lena圖像的迭代次數(shù)與灰度均方差MSE和峰值信噪比PSNR的關(guān)系曲線,在迭代次數(shù)為42時,已達(dá)到了收斂.應(yīng)用其他圖像重復(fù)進(jìn)行多次實驗,發(fā)現(xiàn)迭代次數(shù)均在50次內(nèi)能達(dá)到收斂.

另外還有三個參數(shù)β、λ和α.β值控制著光照項對模型的影響,由于本文模型是基于光流強度保持恒定的假設(shè),但實際情況中,由于光照變化等原因,光流強度往往不會保持恒定,所以增加β項來稍微改進(jìn)光流約束項,β值應(yīng)取較小值,經(jīng)過多次實驗本文選取β為0.01.

正則化參數(shù)λ是影響著正則項和數(shù)據(jù)項的權(quán)重,由于光流模型的正則項起懲罰光流場中大的變化以獲得平滑光流場的作用,若λ值大,則對光流約束項的影響越大,即光流約束只要存在一點偏離,都會導(dǎo)致該泛函發(fā)散;反之,若λ越小,則對光流約束項的影響越小.總之,λ的大小決定正則化的強度,決定了模型解的傾向,較大的λ傾向于滿足正則項約束,較小的λ傾向于滿足光流約束,本文的λ值則在兩種極端中尋找某種平衡.一般λ值建議在10～100之間,本文經(jīng)過多次實驗,選取λ值為40.

圖2 迭代次數(shù)的選擇Fig.2 Choice of iteration number

下面分析分?jǐn)?shù)階微分階次α的選取,對于可積函數(shù)f(x),其傅里葉變換為根據(jù)傅里葉變換的性質(zhì),其階導(dǎo)數(shù)為把階次從整數(shù)階擴展到分?jǐn)?shù)階可得:

圖3 分?jǐn)?shù)階微分算子幅頻特性曲線Fig.3 The amplitude-frequency curve of fractional diあerentiator

從圖3可以看出,當(dāng)階次小于1時,與整數(shù)1階相比,分?jǐn)?shù)階微分算子對中低頻部分的幅值有所提升,而對高頻部分的幅值則有所壓縮,隨著階次的變大,對中低頻部分的提升幅度越來越小,對高頻部分的壓縮也減弱.當(dāng)階次大于1時,與整數(shù)1階相比,分?jǐn)?shù)階微分算子對中低頻部分的幅值由提升變成了壓縮,階次越大,壓縮越大,而對高頻部分的幅值則變成了提升,階次越大,提升越大.當(dāng)階次小于1時,可看作一個對中低頻增強放行,對高頻壓縮放行的濾波器,而當(dāng)階次大于1時,可看作一個高通濾波器.在圖像處理中,圖像的平滑區(qū)域?qū)?yīng)著信號的低頻部分,圖像的紋理信息對應(yīng)著信號的中低頻部分,圖像的輪廓和噪聲對應(yīng)著信息的高頻部分,總之根據(jù)實際圖像處理中目的不同選取不同的階次區(qū)域,而根據(jù)本文模型可知,我們希望阻止中低頻部分通過濾波器而非放行,從而達(dá)到保持紋理細(xì)節(jié)信息的目的,故階次的合理區(qū)間建議選擇在1～2之間.但針對不同的圖像,其最佳參數(shù)不相同.分?jǐn)?shù)階微分相對于整數(shù)階微分的優(yōu)勢是可以有目的地控制微分后的值大于或者小于整數(shù)階微分后的值,但這是通過某個像素點周圍的眾多像素點共同參與計算實現(xiàn)的,若模板寬度過小,則參與計算的像素點就少,體現(xiàn)不了分?jǐn)?shù)階微分的優(yōu)勢,但模板取得過大,又會起到光滑平均的作用,因此要在兩者之間平衡,需要經(jīng)過多次實驗取得最佳值.所以本文分?jǐn)?shù)階模板中的階次和模板寬度都是通過實驗確定,并且對于不同的圖像,其最佳的參數(shù)都是不同的,需要通過多次實驗確定.

3 實驗結(jié)果與分析

本文選用標(biāo)準(zhǔn)庫圖像Lena和真實大腦圖像Brain1、Brain2(兩種不同類型的)進(jìn)行測試,如圖4所示.圖4中第一行為參考圖像,第二行為與第一行對應(yīng)的浮動圖像,即變形圖像.為了證明本文方法的有效性,從定性和定量兩個方面來評價配準(zhǔn)精度.定性評價即觀察視覺效果,比較配準(zhǔn)后的圖像和原始圖像的差異以及通過二者差值圖像黑色區(qū)域所占大小來評價,定量評價則通過配準(zhǔn)后圖像和原始圖像的灰度均方誤差MSE(Mean square error)和峰值信噪比PSNR(Peak signal to noise ratio)來評價.

灰度均方誤差定義如下:

峰值信噪比定義如下:

其中,I1是t時刻的圖像,即參考圖像,I是經(jīng)過配準(zhǔn)后的圖像,m×n為像素個數(shù).理想情況下MSE值應(yīng)為0,表示兩幅圖像同一位置上點的灰度值相同.PSNR值越大越好,表示配準(zhǔn)后圖像越逼近參考圖像.

本文實驗所使用的計算機環(huán)境為:Intel Core i3-2130CPU,3.40GHZ,內(nèi)存4GB,操作系統(tǒng)為64位Windows 7.0,程序使用R2012b版Matlab實現(xiàn).

圖4 實驗圖像(其中第一行為參考圖像,第二行為與參考圖像對應(yīng)的浮動圖像)Fig.4 Experimental images(The fi rst line are reference images,the second line are fl oating images corresponding to reference images)

3.1 標(biāo)準(zhǔn)庫圖像實驗

本次實驗選用圖4(a)所示的Lena作為參考圖像進(jìn)行實驗,為了證明本文方法有效性,對Lena圖像進(jìn)行扭曲形變,如圖5(a)所示,Lena圖像頭發(fā)處發(fā)生了明顯扭曲變形.圖5展示了TV-L1算法和本文算法最佳階次配準(zhǔn)的視覺效果.圖5(b)為TV-L1算法配準(zhǔn)后的圖像,圖5(c)～(e)分別是本文算法取最佳階次三個不同的模板寬度所對應(yīng)的配準(zhǔn)結(jié)果圖,從圖中可得圖5(b)中的Lena圖像眼角處以及頭發(fā)處出現(xiàn)了信息丟失,而對比本文算法配準(zhǔn)后的配準(zhǔn)圖片圖5(c)、圖5(d)和圖5(e)在對應(yīng)處效果則要好很多,并未出現(xiàn)明顯信息丟失.圖5(f)是TV-L1算法配準(zhǔn)后的差值圖像(配準(zhǔn)后的圖像與輸入圖像的差值圖像),圖5(h)～(j)是本文算法取最佳分?jǐn)?shù)階階次和三個不同模板寬度時所對應(yīng)的差值圖像.如果兩幅圖像完全配準(zhǔn),理論上其差值圖像應(yīng)該全黑,即表示配準(zhǔn)后圖像與參考圖像相同,所以差值圖像黑色區(qū)域越多越好,通過對比可知,圖5(h)～(j)相對圖5(g)黑色區(qū)域明顯多很多,從而也證明了分?jǐn)?shù)階微分TV-L1配準(zhǔn)算法能夠較好地提高配準(zhǔn)精度.在分?jǐn)?shù)階階次為1.2時,模板寬度k為2(即為5×5)時,配準(zhǔn)效果達(dá)到最佳.但是,對于不同的圖片,最佳的實驗參數(shù)仍然需要進(jìn)行多次實驗人工選定.

為了得到最佳的模板參數(shù),我們進(jìn)一步進(jìn)行實驗,依據(jù)灰度均方誤差MSE和峰值信噪比PSNR來分析模板的參數(shù)α,k對配準(zhǔn)的影響.圖6(a)展示了均方誤差(MSE)和階次α與模板寬度k的關(guān)系,MSE值越小表示配準(zhǔn)精度越高,圖6(b)展示了峰值信噪比(PSNR)和階次α與模板寬度k的關(guān)系,PSNR值越大表示配準(zhǔn)精度越高.圖6(a)和圖6(b)中三條曲線分別是在模板大小參數(shù)k取1、2、3時的MSE值和PSNR值,當(dāng)階次α等于1時,即為TV-L1算法的實驗結(jié)果.為方便觀察,TVL1算法和本文算法最佳階次下不同模板的MSE值和PSNR值已列在表1和表2中.從圖6(a)中可以看出,MSE最小值并不是在階次α等于1的時候,階次從0.1開始,隨著階次的增大,MSE呈現(xiàn)減小的趨勢,當(dāng)α值在1.2處時,MSE值達(dá)到最小,這是因為當(dāng)α小于1時,紋理等中低頻信息丟失,圖像紋理細(xì)節(jié)模糊,導(dǎo)致圖片配準(zhǔn)精度降低,而當(dāng)α大于1時,紋理細(xì)節(jié)信息被阻止通過濾波器從而達(dá)到保持紋理的效果,但當(dāng)α過大時,會導(dǎo)致高頻信號如邊緣信息模糊,從而配準(zhǔn)精度降低,所以存在一個最佳階次,使中低頻信號和高頻信號達(dá)到一個平衡,而此次實驗的最佳階次為1.2.圖6(b)中PSNR值最高的階次也不在α為1處,原因同上.通過實驗得知,在分?jǐn)?shù)階階次α等于1.2,模板寬度k等于2(即為5×5)時,配準(zhǔn)效果達(dá)到最佳.不過,對于不同的圖片,最佳的實驗參數(shù)仍然需要進(jìn)行多次實驗人工選定.同時也能得到結(jié)論:一階微分并不能有效地刻畫位移場中的具有弱導(dǎo)數(shù)性質(zhì)的紋理細(xì)節(jié)等信息,而近年來的研究也表明分?jǐn)?shù)階微分比整數(shù)階微分具有更好的描述紋理信息能力,所以用分?jǐn)?shù)階微分來代替整數(shù)階微分在TV-L1算法中的應(yīng)用,能夠獲得更高的配準(zhǔn)精度.

圖5 Lena圖像實驗(第一行為浮動圖像和配準(zhǔn)后圖像,第二行為差值圖像.(a)為浮動圖像,(b)～(e)為配準(zhǔn)后的圖像;(b)TV-L1方法;(c)本文方法(α=1.2,k=1);(d)本文方法(α=1.2,k=2);(e)本文方法(α=1.2,k=3);(f)～(j)分別為第一行圖像與參考圖像(圖4(a))的差值圖像)Fig.5 Lena image(The fi rst line is fl oating image and registered image,the second line is diあerence image.(a)Floating image,(b)～(e)are registered images,(b)TV-L1,(c)Our method(α =1.2,k=1),(d)Our method(α =1.2,k=2),(e)Our method(α =1.2,k=3),(f)～(j)are diあerence images)

圖6 配準(zhǔn)精度與模板參數(shù)間的關(guān)系曲線Fig.6 Curve between registration accuracy with mask parameters

3.2 真實圖像實驗

圖7 Brain1圖像實驗(第一行為浮動圖像和配準(zhǔn)后圖像,第二行為差值圖像.(a)為浮動圖像,(b)～(e)為配準(zhǔn)后的圖像;(b)TV-L1方法;(c)本文方法(α=1.3,k=1);(d)本文方法(α=1.3,k=2);(e)本文方法(α=1.3,k=3);(f)～(j)分別為第一行圖像與參考圖像(圖4(b))的差值圖像)Fig.7 Brain1 image(The fi rst line is fl oating image and registered image,the second line is diあerence image.(a)Floating image,(b)～(e)are registered images,(b)TV-L1,(c)Our method(α =1.3,k=1),(d)Our method(α =1.3,k=2),(e)Our method(α =1.3,k=3),(f)～(j)are diあerence images)

本次做了兩組實驗,選用腦部圖像,如圖4中Brain1和Brain2,將參考圖像與浮動圖像配準(zhǔn).該兩組浮動圖像相對于參考圖像只有輕微變形,人眼辨識已稍顯困難,但通過參考圖像和浮動圖像的差值圖像還是可以看出大量白色區(qū)域存在,如圖7所示.圖7(g)是TV-L1算法配準(zhǔn)后圖像圖7(b)與參考圖像Brain1的差值圖像,圖7(h)、(i)和(j)分別是本文算法不同參數(shù)下的配準(zhǔn)圖像與參考圖像Brain1的差值圖像,通過比較可以發(fā)現(xiàn),本文算法的差值圖像白色區(qū)域有較大減少,并且當(dāng)α=1.3,k=2時,MSE值達(dá)到最小,PSNR值達(dá)到最大,也證明本文算法相比TV-L1算法具有更高的配準(zhǔn)精度,為了方便觀察,TV-L1算法和本文算法最佳階次下不同模板的MSE值和PSNR值已列在表1和表2中.同理,為了進(jìn)一步驗證本文方法的有效性,再次選用了Brain2圖像重復(fù)上述實驗,實驗結(jié)果如圖8所示和表1及表2所示.多組實驗結(jié)果表明本文的配準(zhǔn)方法較原始的TV-L1方法配準(zhǔn)精度得到較大的提高.

圖8 Brain2圖像實驗(第一行為浮動圖像和配準(zhǔn)后圖像,第二行為差分圖像.(a)為浮動圖像,(b)～(e)為配準(zhǔn)后的圖像;(b)TV-L1方法;(c)本文方法(α=1.3,k=1);(d)本文方法(α=1.3,k=2);(e)本文方法(α=1.3,k=3);(f)～(j)分別為第一行圖像與參考圖像(圖4(c))的差值圖像)Fig.8 Brain2 image.The fi rst line is fl oating image and registered image,the second line is diあerence image((a)lf oating image,(b)～ (e)are registered images,(b)TV-L1,(c)Our method(α =1.3,k=1),(d)Our method(α =1.3,k=2),(e)Our method(α =1.3,k=3),(f)～(j)are diあerence images)

表1 參考圖像和配準(zhǔn)圖像的均方誤差(MSE)Table 1 Mean square error(MSE)of reference image and registered image

表2 峰值信噪比(PSNR)Table 2 Peak signal to noise ratio(PSNR)

4 結(jié)論

將G-L分?jǐn)?shù)階微分引入到TV-L1光流模型中,提出了一種新的基于分?jǐn)?shù)階微分的全變分方法(FTV-L1)來解決圖像的非剛性配準(zhǔn)問題.結(jié)合全變分能量方程的對偶形式來極小化FTV-L1光流模型獲得位移場.融合G-L分?jǐn)?shù)階的TV-L1光流模型能夠解決圖像灰度均勻,弱紋理區(qū)域配準(zhǔn)結(jié)果中的信息模糊的問題,這是因為分?jǐn)?shù)階微分比整數(shù)階微分具有更好的細(xì)節(jié)描述能力,可以有針對性地選擇合適的階次對具有弱導(dǎo)數(shù)性質(zhì)的信息如紋理信息進(jìn)行抑制或者非線性保留,因此可以提高圖像的配準(zhǔn)精度.另外,通過實驗給出了配準(zhǔn)精度和階次、模板寬度的關(guān)系,從而為最佳模板參數(shù)的選取提供了依據(jù),盡管不同類型的圖像其最佳參數(shù)是不同的,但是其最佳配準(zhǔn)階次一般在1～2之間.理論分析和實驗結(jié)果均表明,本文的方法可用于弱紋理和弱邊緣圖像的非剛性配準(zhǔn),提高圖像配準(zhǔn)的精度,是TV-L1光流場模型的一個重要延伸和擴展.

但是不同圖像配準(zhǔn)的最佳階次需要不斷測試,因而是比較耗時和費力的.今后可以研究自適應(yīng)G-L分?jǐn)?shù)階TV-L1的非剛性圖像配準(zhǔn)算法.此外本文的分?jǐn)?shù)微分掩模是二維的,要完成三維圖像的配準(zhǔn),還需將其擴展到三維空間.

1 Pu Yi-Fei,Wang Wei-Xing.Fractional diあerential masks of digital image and their numerical implementation algorithms.Acta Automatica Sinica,2007,33(11):1128?1135(蒲亦非,王衛(wèi)星.數(shù)字圖像的分?jǐn)?shù)階微分掩模及其數(shù)值運算規(guī)則.自動化學(xué)報,2007,33(11):1128?1135)

2 Chen Qing,Liu Jin-Ping,Tang Zhao-Hui,Li Jian-Qi,Wu Min.Detection and extraction of image edge curves and detailed features using fractional diあerentiation.Acta Electronica Sinica,2013,41(10):1873?1880(陳青,劉金平,唐朝暉,李建奇,吳敏.基于分?jǐn)?shù)階微分的圖像邊緣細(xì)節(jié)檢測與提取.電子學(xué)報,2013,41(10):1873?1880)

3 Pu Y F,Siarry P,Zhou J L,Liu Y G,Zhang N,Huang G,Liu Y Z.Fractional partial diあerential equation denoising models for texture image.Science China Information Sciences,2014,57(7):1?19

4 Liu J,Chen S C,Tan X Y.Fractional order singular value decomposition representation for face recognition.Pattern Recognition,2008,41(1):378?395

5 Zhang Y,Pu Y F,Hu J R,Zhou J L.A class of fractionalorder variational image inpainting models.Applied Mathematics and Information Sciences,2012,6(2):299?306

6 Ren Z M.Adaptive active contour model driven by fractional order fi tting energy.Signal Processing,2015,117:138?150

7 Xue Peng,Yang Pei,Cao Zhu-Lou,Jia Da-Yu,Dong En-Qing.Active demons non-rigid registration algorithm based on balance coeきcient.Acta Automatica Sinica,2016,42(9):1389?1400(薛鵬,楊佩,曹祝樓,賈大宇,董恩清.基于平衡系數(shù)的 Active Demons非剛性配準(zhǔn)算法.自動化學(xué)報,2016,42(9):1389?1400)

8 Zhang Gui-Mei,Cao Hong-Yang,Chu Jun,Zeng Jie-Xian.Non-rigid image registration based on low-rank Nystr¨om approximation and spectral feature.Acta Automatica Sinica,2015,41(2):429?438(張桂梅,曹紅洋,儲珺,曾接賢.基于Nystr¨om 低階近似和譜特征的圖像非剛性配準(zhǔn).自動化學(xué)報,2015,41(2):429?438)

9 Yan De-Qin,Liu Cai-Feng,Liu Sheng-Lan,Liu De-Shan.A fast image registration algorithm for diあeomorphic image with large deformation.Acta Automatica Sinica,2015,41(8):1461?1470(閆德勤,劉彩鳳,劉勝藍(lán),劉德山.大形變微分同胚圖像配準(zhǔn)快速算法.自動化學(xué)報,2015,41(8):1461?1470)

10 Thirion J P.Image matching as a diあusion process:an analogy with Maxwell′s demons.Medical Image Analysis,1998,2(3):243?260

11 Wang H,Dong L,O′Daniel J,Mohan R,Garden A A,Ang K K,Kuban D A,Bonnen M,Chang J Y,Cheung R.Validation of an accelerated ‘demons’algorithm for deformable image registration in radiation therapy.Physics in Medicine and Biology,2005,50(12):2887?2905

12 Palos G,Betrouni N,Coulanges M,Vermandel M,Devlaminck V,Rousseau J.Multimodal matching by maximisation of mutual information and optical fl ow technique.In:Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society(IEMBS).San Francisco,CA,USA:IEEE,2004.1679?1682

13 Pock T,Urschler M,Zach C,Beichel R,Bischof H.A duality based algorithm for TV-L1-optical- fl ow image registration.Medical Image Computing and Computer-Assisted Intervention-MICCAI 2007.Berlin Heidelberg:Springer-Verlag,2007.511?518

14 P′erez J S,Meinhardt-Llopis E,Facciolo G.TV-L1 optical fl ow estimation.Image Processing on Line,2013,3:137?150

15 Yip S S F,Coroller T P,Sanford N N,Huynh E,Mamon H,Aerts H J W L,Berbeco R I.Use of registration-based contour propagation in texture analysis for esophageal cancer pathologic response prediction.Physics in Medicine and Biology,2016,61(2):906?922

16 Horn B K,Schunck B G.Determining optical fl ow.In:Technical Symposium East.Washington,D.C.,USA:International Society for Optics and Photonics,1981.319?331

17 Rudin L I,Osher S,Fatemi E.Nonlinear total variation based noise removal algorithms.Physica D:Nonlinear Phenomena,1992,60(1?4):259?268

18 Chambolle A,Pock T.A fi rst-order primal-dual algorithm for convex problems with applications to imaging.Journal of Mathematical Imaging and Vision,2011,40(1):120?145

19 Cafagna D.Fractional calculus:a mathematical tool from the past for present engineers.IEEE Industrial Electronics Magazine,2007,1(2):35?40

Research on TV-L1Optical Flow Model for Image Registration Based on Fractional-order Diあerentiation

ZHANG Gui-Mei1SUN Xiao-Xu1LIU Jian-Xin2CHU Jun1

In computer vision and medical image analysis,non-rigid image registration is a challenging task.TV-L1(Total variation-L1)optical fl ow model has been proved to be an eあective method in the fi eld of non-rigid image registration.It can solve the problem of fuzzy edge caused by smooth displacement fi elds of Horn-Schunck,but its fi rst-order derivative in regularization term leads to fuzzy texture information with weak derivative property.Aiming at the problem,this paper introduces G-L(Grnwald-Letnikov)fractional diあerentiation to TV-L1optical fl ow model,and proposes a new TV-L1optical fl ow model based on fractional diあerentiation,and then fi nds the solution of the model using primal-dual algorithm.In this paper we use Grnwald-Letnikov fractional order diあerential instead of the fi rst-order derivative in the regularization term for its better ability of detail description than fi rst-order′s.Then we purposefully control to retain or inhibit the texture information with weak derivative nature,thus improving the registration accuracy.Experimental results show that the proposed method has a better registration accuracy in registration of texture information with weak derivative,and that it can be considered an important extension and generalization of TV-L1optical fl ow modes.

Fractional diあerential,Grnwald-Letnikov,TV-L1(Total variation-L1)model,optical fl ow,weak texture,non-rigid registration

Zhang Gui-Mei,Sun Xiao-Xu,Liu Jian-Xin,Chu Jun.Research on TV-L1optical fl ow model for image registration based on fractional-order diあerentiation.Acta Automatica Sinica,2017,43(12):2213?2224

2016-04-29 錄用日期2016-10-05

April 29,2016;accepted October 5,2016國家自然科學(xué)基金(61462065,61661036),江西省自然科學(xué)基金 (20151BAB207036),江西省科技支撐計劃重點項目(20161BBF60091)資助

Supported by National Natural Science Foundation of China(61462065,61661036),Natural Science Foundation of Jiangxi Province(20151BAB207036),and the Key Science and Technology Support Program of Jiangxi Province(20161BBF60091)

本文責(zé)任編委張長水

Recommended by Associate Editor ZHANG Chang-Shui

1.南昌航空大學(xué)江西省圖像處理與模式識別重點實驗室 南昌330063 2.西華大學(xué)機械工程學(xué)院成都610039

1.Key Laboratory of Jiangxi Province for Image Processing and Pattern Recognition,Nanchang Hangkong University,Nanchang 330063 2.School of Mechanical Engineering,Xihua University,Chengdu 610039

張桂梅,孫曉旭,劉建新,儲珺.基于分?jǐn)?shù)階微分的TV-L1光流模型的圖像配準(zhǔn)方法研究.自動化學(xué)報,2017,43(12):2213?2224

DOI10.16383/j.aas.2017.c160367

張桂梅 南昌航空大學(xué)江西省圖像處理與模式識別重點實驗室教授.主要研究方向為計算機視覺,圖像處理與模式識別.

E-mail:guimei.zh@163.com

(ZHANGGui-Mei Professor at the Key Laboratory of Jiangxi Province forImage Processing and Pattern Recognition,Nanchang Hangkong University.Her research interest covers computer vision,image processing,and pattern recognition.)

孫曉旭 南昌航空大學(xué)江西省圖像處理與模式識別重點實驗室碩士研究生.主要研究方向為圖像處理與計算機視覺.

E-mail:sunxiaoxu@outlook.com

(SUN Xiao-Xu Masterstudent attheKey Laboratory ofJiangxi Province for Image Processing and Pattern Recognition,Nanchang Hangkong University.His research interest covers image processing and computer vision.)

劉建新 西華大學(xué)機械工程學(xué)院教授.主要研究方向為圖像處理與機器視覺.本文通信作者.

E-mail:jamsonliu@163.com

(LIU Jian-Xin Professor at the School of Mechanical Engineering,Xihua University. His research interest covers image processing and machine vision.Corresponding author of this paper.)

儲 珺 南昌航空大學(xué)江西省圖像處理與模式識別重點實驗室教授.主要研究方向為圖像處理與計算機視覺.

E-mail:chujun99602@163.com

(CHU Jun Professor at the Key Laboratory of Jiangxi Province for Image Processing and Pattern Recognition,Nanchang Hangkong University.Her research interest covers image processing and computer vision.)