復剪切波變換域各向異性濾波DTI去噪

劉帥奇,李鵬飛,安彥玲,趙傳慶,耿 鵬

1(河北大學 電子信息工程學院,河北 保定 071000)2(河北省數字醫療工程重點實驗室,河北 保定 071000)3(河北省機器視覺工程技術研究中心,河北 保定 071000)4(石家莊鐵道大學 信息技術與工程學院,石家莊 050043) E-mail:pfli_hbu@163.com

1 引 言

擴散張量成像(DTI)是一種廣泛使用的腦成像方法,它作為測量連接大腦區域到分布式認知網絡的軸突纖維束的一個重要工具.DTI目前用于研究阿爾茨海默病、多發性硬化癥、亨廷頓病和其他神經退行性疾病中大腦纖維束的完整性[1,2].盡管DTI成像是無創的,并且具有為研究大腦神經的工作人員和臨床醫生提供重要信息的巨大潛力,但是眾所周知的是,其信噪比偏低,DTI圖像中的噪聲一方面會使得圖像中的信號產生不同類別的偏差,另一方面會導致DTI圖像出現失真,在一定程度上限制了DTI的發展與應用,因此開發新的DTI圖像去噪算法顯得尤為重要.

醫學圖像去噪是圖像處理任務中的一項核心問題,而DTI圖像去噪也是很重要的醫學圖像處理課題[3,4].隨著小波變換的迅猛發展,基于變換域的去噪算法也得到了發展[5],但在不同尺度的條件下,小波變換很難實現圖像中局部特征的精確逼近,因此在對噪聲進行預處理的過程中通常會造成原有的時域特征消失[6,7].為了更好地表示圖像的局部特征,眾多專家學者在基于小波變換的前提下提出了不同的多尺度幾何變換方法,其中應用最為廣泛有Contourlet[8]、Curvelet[9]和剪切波變換[10].其中,剪切波變換近年來已成為分析和處理多維數據的最有效的框架之一,該變換克服了Contourlet、Curvelet變換的缺點,其不僅有效地克服了其他多尺度幾何變換中的局限性并提供最優的稀疏表示,而且剪切系統由單個函數生成,該函數通過尺度矩陣和剪切矩陣進行膨脹,并在時域中平移,從而形成一個仿射系統,使得剪切波變換具有靈活的方向選擇性[11,12].但是,文獻[12]指出離散剪切波變換在實現多分辨率分析時運用了Laplace金字塔模型,使得變換不再具有移不變的性質,而非下采樣剪切波變換的冗余度較大,從而使得算法運行耗時巨大.因此,文獻[13]在引入雙樹復小波變換的基礎上提出了復剪切波變換,該變換只需在有限的冗余度的條件下便可輕松實現移不變性,并且隨著分解尺度的增加,其冗余度不會變大,從而在運行速度方面改進了現有的圖像處理算法.此外,復剪切波變換具有更加強大的方向選擇性,這也就意味著相比剪切波變換而言它對圖像的表示更為稀疏,顯然,這將對去噪結果產生積極的影響.基于復剪切波變換,文獻[14]提出了一種新的DTI去噪算法,該算法利用復擴散各向異性濾波對復Shearlet變換的高頻系數進行去噪,利用P-M各向異性濾波對復Shearlet變換的低頻系數進行去噪,不僅得到了良好的去噪結果,并且有效地保留了圖像的邊界信息.

各向異性濾波與傳統圖像濾波相比存在一些差異,根本原因在于其在濾波過程中的各向異性,在圖像中各個局部區域和方向上的濾波結果均與各向異性平滑的強度緊密聯系在一起[15].當圖像中局部區域或方向上的灰度變化較大時,擴散作用通常比較小,因此可以最優地保留圖像中的局部細節特征;當圖像區域中灰度變化較小或只有孤立噪聲點時,若應用較強的濾波過程便可以去除圖像中的噪聲,同時實現圖像邊緣細節有效保留.這一概念的發展改變了原有濾波方法的局限性,使得平滑技術在本質上顯著提升,一方面有效地降低了圖像中的噪聲,另一方面使得圖像結構邊界更加平滑.圖像中結構信息突出的區域和結構信息微弱的區域可以經由結構張量進行最優地劃分,即可得到圖像的邊緣輪廓區域與平坦光滑區域,為后期處理高標準的圖像邊界保留和去噪問題提供了便利.文獻[16]將結構張量應用于對圖像中邊緣、形狀角點特征等結構信息的提取與識別,從而實現了對圖像方向場的估計.文獻[17]提出了一種基于結構張量的線性濾波技術,將擴散方程轉變為偏微分方程的形式,使得原有的線性濾波過程成為了非線性過程,這樣改進后大大增強了邊緣檢測和角點檢測的準確性.文獻[18]通過結合結構張量和各向異性平滑技術,研究得到一種基于結構張量和各向異性平滑的DTI圖像去噪方案,該方法充分利用了結構張量區分圖像結構的特性,將原始圖像像素區分為均勻平坦區域和邊緣輪廓區域后分別進行各項同性濾波和各向異性濾波,很大程度上提升了算法的去噪能力及邊緣保留能力,與文獻[14]相比,無論是視覺效果還是客觀評價指標均顯示出更優的效果.鑒于復剪切波變換后得到的高頻分量主要包含了圖像中邊緣輪廓區域的信息,并且DTI噪聲主要存在于圖像的高頻分量中,因此本文將利用基于結構張量的各向異性濾波進行DTI圖像的高頻系數進行去噪.

基于歐氏度量的DTI圖像去噪過程將擴散張量圖像看作是嵌入到對稱矩陣構成的空間中的圖像,雖然在濾波處理中可以保證張量的正定性,但張量變形時,其會產生特征值膨脹效應.這些局限性產生的一個重要原因是歐氏度量并不適合于平滑擴散張量圖像,因此在DTI去噪領域產生了一個新的去噪框架——黎曼框架.黎曼框架將擴散張量圖像視為一個黎曼流形,其在仿射變換下具有一個固定度量,在文獻[19]中,作者在黎曼框架的相關理論基礎上總結得出了一種新穎的DTI圖像去噪算法,將原本為灰度圖像開發的幾種圖像處理技術,包括圖像插值、修復和非線性各向同性擴散濾波,推廣到擴散張量圖像.該框架是基于擴散張量的矩陣表示形式進行算法計算的,其中還大量的使用了計算量大的矩陣運算(如指數和對數運算).為了改進算法的有效性,文獻[20]提出了一種基于局部坐標的黎曼正則化方法,極大的降低了算法的復雜度.由于黎曼流形可以實現對DTI圖像的局部細節特征進行最優表示,而復剪切波變換后得到的低頻分量包含圖像的細節信息,其中也含有少量噪聲,因此本文根據文獻[19]和文獻[20]中提出的基于黎曼框架的去噪模型,利用基于黎曼幾何框架的各向異性濾波對DTI圖像的低頻系數進行去噪.

為了更加高效地抑制DTI圖像中的噪聲并且保留圖像的邊界特征信息,本文提出了基于復剪切波變換域各向異性濾波的DTI圖像去噪算法.算法首先對DTI張量圖像進行復剪切波分解,得到變換后的高頻系數和低頻系數;其次,利用基于結構張量和各向異性平滑的各向異性濾波對復剪切波變換的高頻系數進行去噪,利用基于黎曼幾何框架的各向異性濾波對復剪切波變換的低頻系數進行去噪;最后,利用復剪切波逆變換得到濾波后的DTI圖像.

2 基本理論

2.1 復剪切波變換

(1)

(2)

文獻[13]在雙樹復小波和剪切波變換的基礎上提出了一種新的多尺度幾何變換方法——復剪切波變換,下面簡述復剪切波變換的構造.標準的剪切波變換包含兩個過程,一是尺度分解,另外一個是方向分解.尺度分解通常采用拉普拉斯金字塔完成,而方向分解則采用剪切變換完成.為了克服拉普拉斯金字塔導致的剪切波變換缺乏移不變性的缺點,復剪切波被提了出來.在復剪切波變換中,在多尺度分解方面利用雙樹復小波變換取代拉普拉斯金字塔,從而使得復剪切波變換擁有平移不變性[13],其中雙樹復小波的高頻子帶系數是復數,這就意味著需要對高頻子帶系數的實部和虛部分別獨立地進行剪切變換.由于雙樹復小波具有移不變性,因此經過移不變的剪切方向濾波器后整個變換依舊具有移不變性.而且,圖像在進行雙樹復小波變換后會分解出6個高頻小波子帶,所獲得的圖像的表達系數更稀疏更有利于圖像去噪.

2.2 各向異性

Brox等人提出了一種各向異性擴散濾波技術[17],該技術在進行圖像去噪時可以有效的保留圖像的邊緣結構信息:

(3)

其中I表示圖像強度,而T提供了平滑過程所需的結構張量.div是指向量的散度,而表示向量的梯度,t是由公共梯度張量G構成迭代時間.通過圖像強度I的外積與高斯核Kρ的卷積可以得到公共梯度張量G:

G=Kρ*(I?I)

(4)

其中?為外積運算,*為卷積運算,Kρ表示標準差為ρ的高斯核.設λ1,λ2,λ3為張量G的特征值,其所對應的特征向量為v1,v2,v3.最大圖像強度梯度的方向由最大特征值表示,而最小圖像強度梯度的方向由最小特征值表示.在保留原始特征向量的條件下,通過調整張量G的特征值可構造出結構張量T,使得指向最大的圖像強度梯度方向的特征向量所對應的特征值最小.通過這種方式得到的結構張量可以明顯提升圖像中邊緣區域的平滑程度[18].

在含噪的DTI圖像中,由于圖像的平坦區域所含像素的數量高于邊界區域,在平滑過程中極易丟失圖像的邊界信息,因此文獻[18]提出利用結構張量的各向異性濾波算法,即在所有加權方向上使用一個公共梯度張量G,其可以通過公式(5)實現:

G=Kρ*∑(Ip?Ip)

(5)

其中ρ是高斯核Kρ的標準差,*為卷積運算,表示梯度運算,Ip是加權方向p的圖像強度,?表示外積運算.

為了對結構內部的平坦區域和結構邊界的邊緣輪廓區域分別實現各向同性平滑和各向異性平滑,通過求解公式(5)中公共梯度張量G的歸一化倒數便可以得到結構張量T.設λg1,λg2,λg3為G的三個特征值,vg1,vg2,vg3分別為G的三個特征值所對應的特征向量,因此結構張量T的特征向量和特征值可以表示為:

vti=vgi,i=1,2,3

(6)

(7)

在結構內部的平坦區域中,由于結構張量T與公共梯度張量G的特征值大小相同,從而使得在任意方向上對圖像進行各向同性濾波后所得到的結果基本保持一致;在結構邊界的邊緣區域中,由于最大圖像強度梯度的vt1方向對應的結構張量的特征值最小,而最小圖像強度梯度的vt3方向對應的結構張量的特征值最大,因此在垂直于和平行于邊界的方向上各向異性平滑的結果存在差異.為了解決在均勻區域和邊界區域上平滑結果不一致的問題,結構張量的特征值被歸一化為常數C,即λt1+λt2+λt3=C,此時對于平坦區域和邊界區域而言,各向同性平滑和各向異性平滑的總量在整個圖像區域上即可保持一致.

2.3 黎曼流形

在DTI中,張量可以看作是黎曼流形中的正定矩陣,距離對應于Fisher-Rao度量[19],這是基于統計學所產生的理論,因為正定矩陣表示模擬水分子擴散性的協方差矩陣.因此,與纖維結構相關的定向擴散性由相應張量的各向異性進行表示,通常纖維的方向由一個較大的特征值和其相應的特征向量來決定.我們用P(n)表示正定張量的空間,D表示P(n)中的正定矩陣,其黎曼度量表示為

(8)

其中對稱矩陣Λ1,Λ2表示在D∈P(n)處的切向量.利用這個度量即可將正定矩陣的空間轉換為一個Cartan-Hadamard流形[20].這意味著此流形是完整的,并且具有簡單的連接以及非正截面曲率.該類流形與一般意義上的流形相比具有特殊的優良特性,例如,測地線從以恒定速度移動時沒有固有加速度的曲線意義上來講代表局部最短路徑.此外,數學上非常重要的凸面概念可以以一種恰當的方式推廣到該空間.這種仿射不變的黎曼流形來源于信息幾何[21],它把張量分布的空間看作黎曼流形,用Fisher信息矩陣作為一個適當的度量.在具有固定均值的多元正態分布的情況下,Fisher信息矩陣歸結為仿射不變度量,其最初是作為度量固定均值正態分布之間的距離而開發的,現已廣泛地與張量值數據結合使用.

作為黎曼流形上的一般性質,測地線從切空間中任意點D到流形本身實現了一個局部微分同胚映射,稱為指數映射[22],這使得我們能夠利用切向量從局部上識別流形中的點.利用公式(9)可以將該微分同胚映射表示為矩陣指數的形式,如下所示

expD(A)=D1/2exp(D-1/2AD-1/2)D1/2

(9)

其中,切向量A∈TDP(n),正定矩陣D∈P(n).expD(A) 可以看作是起始于黎曼流形上的一點D在切向量A上的單位距離.同樣,通過公式(10)可以唯一地定義逆映射,即對數映射:

logD(E)=D1/2log(D-1/2ED-1/2)D1/2

(10)

實際上,我們想要通過各向異性濾波對DTI進行去噪,因此將去噪問題轉換為一個最小化問題.設f(∑)是一個目標函數,通常為梯度的風險函數.那么擴散張量圖像的去噪過程轉換為函數f(∑)最小化問題.令∑t表示對當前∑的估計值,Wt=?∑f=[?f/?σij]表示在∑處的導數,其中σij表示矩陣中第i行、第j列的元素.一階梯度下降的原理是朝向最快的下降方向,在與梯度相反的方向上前進一小段步長ε,并重復該過程.但是,標準算子∑t+1=∑t-εWt只對平坦的歐幾里得矩陣空間中的非常短的步長有效,而且在計算的過程中很容易地偏離正定張量所組成的錐形區域.而在黎曼流形中,可以通過沿著從∑開始的測地線反向地在Wt上前進步長為ε而進行計算.這種梯度下降確保迭代過程不離開流形.則黎曼流形映射為

∑t+1=exp∑t(-εWt)=∑1/2exp(-ε∑-1/2Wt∑-1/2)∑1/2

(11)

可以看到,該映射可以推廣到計算張量場上的偏微分方程,如?t∑(x,t)=-W(x,t).

一般情況下,在圖像處理領域采用梯度下降算法來解決偏微分方程的問題.為了得到更好的去噪效果,將目標函數f(∑)映射到黎曼空間中得到C(∑),從而將f(∑)最小化問題轉換為C(∑)最小化問題.黎曼流形下的各向異性濾波轉換為沿著測地線從∑t出發,在切向量-C方向上進行梯度下降.上述梯度下降方案可通過最小化C(∑)實現:

(12)

(13)

Λt+1=-εC(∑)

(14)

通過指數映射將Λt+1映射到流形上,即

∑t+1(x)=exp∑t(x)(Λt+1)=exp∑t(x)(-εC(∑t(x)))

(15)

并且通過迭代達到最大步長值εmax,最終即可得到去噪后的擴散張量圖像∑denoise=∑t+1.

3 基于復剪切波變換域各向異性濾波的DTI圖像去噪算法

本文在文獻[14]中復剪切波變換的基礎上,結合文獻[18]及黎曼流形的相關內容,提出了一種新的基于復剪切波變換域的DTI圖像去噪算法.該算法不僅可以利用復剪切波變換實現了對圖像的稀疏表示,使得算法在運行速度方面得到了進一步提升,去噪效果更優,并且可以更好的保留張量圖像的邊緣結構.

本文對DTI圖像噪聲去除的算法具體步驟為:

步驟1.對DTI圖像進行復Shearlet分解;

步驟2.對復Shearlet分解后得到的高頻分量進行基于結構張量和各向異性平滑的各向異性濾波;

步驟3.對復Shearlet分解后得到的低頻分量進行基于黎曼幾何框架各向異性濾波;

步驟4.將濾波后的復Shearlet系數進行復Shearlet逆變換得到最終的DTI去噪圖像.

基于復剪切波變換域各向異性濾波的DTI圖像去噪的流程圖如圖1所示.

圖1 本文去噪算法設計流程圖Fig.1 Flow chart of our denoising algorithm

4 實驗結果和分析

為了測試本文所提出的DTI圖像去噪算法的邊緣保留特性和去噪能力,我們分別對模擬和真實的圖像數據進行處理,將本文所提算法與文獻[14]提出的基于復Shearlet變換域和復擴散各向異性濾波的DTI圖像去噪算法(CST-CDAF)、文獻[18]提出的基于結構張量和各向異性平滑的DTI去噪算法(ST-AS)、文獻[19]提出的基于黎曼幾何的DTI正則化算法(RF)和文獻[23]提出的基于Shearlet變換域的復擴散DTI去噪算法(ST-CD)進行對比.

首先,我們檢驗了本文算法針對模擬張量值圖像的濾波效果.為了產生含噪擴散張量圖像,我們在無噪擴散張量圖像的切向量的每一個實體上將對稱矩陣與均值為零、標準差為σ的先驗高斯隨機變量相加.

圖2顯示了不含噪聲的模擬擴散張量圖像和三種不同噪聲水平(σ=0.1,0.3,0.5)的含噪的模擬擴散張量圖像,為了更加直觀地觀察實驗結果,在模擬實驗中,我們用橢球體對DTI圖像進行表示,將指向不同方向的張量用不同顏色的橢球體進行表示，顏色越深表示擴散張量受噪聲污染的程度越高.對圖2中的三幅含噪圖像應用ST-CD、RF、CST-CDAF、ST-AS和本文所提算法進行去噪處理,濾波結果如圖3所示.

圖2 模擬擴散張量圖像Fig.2 Simulated diffusion tensor image

由圖2可知,受到噪聲污染后,與原始DTI圖像相比,含噪圖像中張量的方向出現了明顯的偏離.由圖3可知,通過對含噪圖像進行去噪處理后,應用五種濾波方法后的實驗結果中張量方向的偏離程度在一定程度上均有所改善,與ST-CD、RF、CST-CDAF和ST-AS四種去噪算法相比,可以清楚地觀察到利用本文所提算法去噪后得到的DTI圖像更好地保留了張量圖像的邊緣結構,去噪效果更優.

下面用真實DTI數據測試本文去噪算法的濾波性能.本文使用的真實數據集采用GE公司的Signa Excite HD 1.5T磁共振系統對人腦進行DTI成像后的數據,其參數為:TR=2.4s,TE=65ms,b=1000s/mm2,FOV=22cm,像素尺寸4.0mm×1.7mm×1.7mm,矩陣為256×256.圖4給出了對真實DTI圖像數據應用ST-CD、RF、CST-CDAF、ST-AS和本文所提算法進行去噪后的實驗結果.從圖4可知,與原始張量場相比,經過本文算法處理過的張量場在圖像質量和視覺效果方面均優于其他四種算法,并且在邊界保留和噪聲去除方面都顯示出最佳性能.

圖3 模擬張量圖像去噪結果.(a-c)分別為對圖2中的3種含噪圖像應用ST-CD去噪算法濾波后的結果;而(d-f)為應用RF算法濾波后的結果;(g-i)為應用CST-CDAF算法濾波后的結果;(j-l)為應用ST-AS算法濾波后的結果;(m-o)為應用本文所提算法濾波后的結果Fig.3 Denoising results of simulated image.(a-c)are the results of applying the ST-CD algorithm to the three noised images in Fig.2;(d-f)are the results of applying the RF algorithm;(g-i)are the results of applying the CST-CDAF algorithm;(j-l)are the results of applying the ST-AS algorithm;(m-o)are the results of applying the proposed algorithm

(16)

圖4 5種算法去噪后的真實張量圖像Fig.4 Denoised real tensor images of five algorithems

表1 不同去噪算法中FA的絕對誤差實驗數據表
Table 1 Absolute error of FA for five different
denoising algorithms

去噪算法噪聲水平FA絕對誤差ST-CDσ=0.10.064σ=0.30.09σ=0.50.102RFσ=0.10.059σ=0.30.079σ=0.50.092CST-CDAFσ=0.10.051σ=0.30.073σ=0.50.084ST-ASσ=0.10.048σ=0.30.062σ=0.50.078本文算法σ=0.10.039σ=0.30.055σ=0.50.063

5 總 結

本文在結合多幾何尺度變換中剪切波變換與復剪切波變換的相關內容的基礎上提出了一種新的DTI去噪模型,首先利用復Shearlet變換將DTI圖像分解為高頻分量和低頻分量,然后將基于結構張量和各向異性平滑的DTI去噪算法與基于黎曼幾何框架各向異性濾波的DTI圖像去噪算法分別應用于去除高頻分量和低頻分量中的噪聲.模擬和真實的DTI實驗結果證明了本文提出的算法既能有效的去除噪聲,并且能完美地表示圖像,在醫學圖像處理領域提供了很大的實際應用潛力.