基于結(jié)構(gòu)張量和各向異性平滑的DTI去噪

劉帥奇,李鵬飛,安彥玲,扈 琪,趙 杰

1(河北大學 電子信息工程學院,河北 保定 071000)2(河北省數(shù)字醫(yī)療工程重點實驗室,河北 保定 071000)3(河北省機器視覺工程技術(shù)研究中心,河北 保定 071000)

1 引 言

擴散張量成像(DTI)已經(jīng)演變成為具有非侵入性的一種表征活體組織結(jié)構(gòu)的主要技術(shù)[1].這種表征技術(shù)是基于所測量的擴散張量的特征參數(shù),其中包括通過特征值的差異來評估組織結(jié)構(gòu)完整性,并且通過特征向量的方向來評估結(jié)構(gòu)特征.然而特征參數(shù)對噪聲極其敏感,因此隨機噪聲會影響結(jié)構(gòu)特征的系統(tǒng)偏差評估.由于DTI數(shù)據(jù)信噪比較低,且噪聲不僅會引起信號的隨機波動,還會導致信號與數(shù)據(jù)的偏差,特別是進行纖維結(jié)構(gòu)分布的分析時,被噪聲污染的張量通常表現(xiàn)為方向排列上雜亂、不規(guī)則,使得實驗得到的纖維結(jié)構(gòu)不夠平滑,甚至出現(xiàn)錯誤的結(jié)果,大大的限制了DTI的應(yīng)用.因此抑制擴散張量圖像中產(chǎn)生的噪聲有利于減少這些系統(tǒng)偏差,從而提高DTI組織表征的準確性.

對于所有矢量或張量圖像平滑來說,減少擴散張量圖像中的噪聲需要兩個基本要求:恢復包含在張量數(shù)據(jù)中的方向信息和保留結(jié)構(gòu)邊界.為了實現(xiàn)這兩個目的,Parker等人應(yīng)用了Perona和Malik提出的非線性擴散濾波技術(shù)[2]來減小擴散張量圖像中的噪聲.雖然可以有效地降低擴散各向異性的有偏估計,但是平滑對張量方向信息的影響只是停留在定性證明的階段.Pajevic等人應(yīng)用樣條擬合法從離散DTI數(shù)據(jù)中生成一個連續(xù)、平滑的張量場[3],但是這種方法不能有效地保留數(shù)據(jù)中的結(jié)構(gòu)邊界.Poupon等人提出的基于馬爾可夫模型[4]和隨機弛緩標記[5]的用于恢復主擴散方向(與最大特征值相對應(yīng)的方向)的正則化方法,也因其邊界保留能力而存在一定的局限性.

結(jié)構(gòu)張量技術(shù)可以很好地將結(jié)構(gòu)信息突出的部分和結(jié)構(gòu)信息微弱的部分進行區(qū)分,例如區(qū)分圖像中邊緣輪廓等細節(jié)與平坦光滑的部分,這樣可以大大地方便類似醫(yī)學圖像等對邊緣保留要求比較高的圖像去噪.Bigun和Granlund等人將結(jié)構(gòu)張量應(yīng)用于圖像方向場的估算[6]等方面,利用結(jié)構(gòu)張量快速提取圖像中的結(jié)構(gòu)信息,例如圖像中目標邊緣、目標形狀角點特征等.Weickert和Brox對結(jié)構(gòu)張量的線性高斯濾波技術(shù)做了很大的改進,利用基于偏微分方程形式的擴散方程的方法將線性濾波過程變成了非線性[7],這樣改進后很大程度上提高了邊緣檢測和角點檢測的精度.

迄今為止,國內(nèi)外科研人員已經(jīng)提出了大量的圖像平滑技術(shù)來減少圖像中的噪聲.其中包括B-樣條擬合[8]以及變分原理[9]、差分幾何[10]、變換域去噪算法[11]和經(jīng)驗分解去噪算法[12].然而,這些平滑技術(shù)的數(shù)據(jù)庫尚未實現(xiàn)其實用性,部分原因在于數(shù)值實現(xiàn)方面有些耗時,特別是計算復雜度會隨著加權(quán)方向數(shù)量的擴大而增加,或者是對其實用價值而言,缺乏嚴格的體內(nèi)DTI數(shù)據(jù)證明.各向異性平滑的概念最初由Weickert提出,用于增強標量圖像中的流狀結(jié)構(gòu)[7].各向異性平滑與傳統(tǒng)濾波方法的不同之處在于它的各向異性,即在不同局部區(qū)域、不同方向上,濾波與平滑的強度存在一定差異.它與圖像的局部結(jié)構(gòu)(灰度梯度)相關(guān),在灰度變化較大的區(qū)域或方向上,由于灰度梯度較大、擴散作用較小,從而最大限度地保留這些局部細節(jié)特征;若灰度變化較小或只有孤立噪聲點的區(qū)域,則采用較強的平滑處理以實現(xiàn)抑制噪聲、保留邊緣與細節(jié)的目的.這一概念的提出使得平滑技術(shù)有了質(zhì)的飛躍,其對減少圖像中的噪聲非常有效,并且使得結(jié)構(gòu)邊界更加平滑.

在本文中,我們結(jié)合結(jié)構(gòu)張量和各向異性平滑技術(shù)提出了一種新型DTI去噪方案,首先利用結(jié)構(gòu)張量區(qū)分DTI圖像中的邊緣結(jié)構(gòu)信息和平坦光滑信息,然后對圖像結(jié)構(gòu)內(nèi)部的平坦光滑區(qū)域進行各向同性平滑,對結(jié)構(gòu)邊界區(qū)域進行各向異性平滑,最終得到去噪后的DTI圖像.

2 基本理論

2.1 結(jié)構(gòu)張量

結(jié)構(gòu)張量是從圖像像素中提取結(jié)構(gòu)信息的張量,通常應(yīng)用于圖像和視頻處理[13,14].定義Ω?R2表示一幅二維圖像數(shù)據(jù)矩陣,該二維圖像的灰度I:Ω→R,那么這幅圖像中的每一個像素的初始結(jié)構(gòu)張量就是一個2×2的對稱矩陣場,通過該矩陣場可以求出包括邊緣、方向等有關(guān)灰度圖像的信息.

初始結(jié)構(gòu)張量ST(i,j)定義為

(1)

通常,我們在平滑處理時將結(jié)構(gòu)張量定義為

(2)

其中Gρ是標準差為ρ的高斯核函數(shù).因為和高斯核函數(shù)的卷積是一個線性過程,很多文獻都稱這種經(jīng)典的結(jié)構(gòu)張量為線性結(jié)構(gòu)張量[15-17].因為抑制了噪聲水平,經(jīng)過高斯核函數(shù)處理過的線性結(jié)構(gòu)張量能夠更好的體現(xiàn)出方向信息和邊緣信息.

對結(jié)構(gòu)張量的每個像素信息即對公式(2)進行特征分解,可以得到兩個特征值λ1≥λ2,其對應(yīng)的特征向量為v1和v2,v1近似平行于最小梯度的方向,而v2與該方向正交并且大致平行于最大梯度的方向,每個特征值的大小是其對應(yīng)的特征向量方向上的梯度量.

2.2 各向同性平滑

高斯濾波是各向同性平滑的經(jīng)典模型.其根據(jù)所選擇的高斯核函數(shù)來控制向周圍擴散的比重參數(shù)[10].高斯濾波對于去除高斯白噪聲具有良好的性能.均值為零、方差為σ的一維高斯函數(shù)可以表示為:

(3)

由于DTI圖像一般是二維或三維的數(shù)據(jù),所以通常采用二維高斯函數(shù)作為各向同性濾波器,其表達式為

(4)

其中A為歸一化參數(shù),σ為噪聲方差,其決定了濾波器的去噪效果,方差越大,去噪效果越好.但是這種濾波方法也存在一定的不足之處,那就是在去除圖像噪聲的同時產(chǎn)生圖像細節(jié)和邊緣模糊現(xiàn)象.

2.3 各向異性平滑

為了增強標量圖像中的流狀結(jié)構(gòu),Weickert提出了使用偏微分方程平滑圖像的各向異性擴散濾波技術(shù)[7]

(5)

(6)

其中?代表外積運算,參數(shù)ρ是高斯核的標準差,決定了梯度張量的空間尺度.G的特征值(λ1,λ2,λ3,以降序排列)和特征向量對應(yīng)于圖像強度梯度在空間尺度ρ上的大小和方向.與最大特征值相關(guān)聯(lián)的特征向量是最大強度梯度的方向,與最小特征值相關(guān)聯(lián)的特征向量是最小強度梯度的方向.通過將特征值不成比例地改變尺度(但保留特征向量),使得對應(yīng)于最小強度梯度方向的特征值是最大的,以這種方式構(gòu)造的結(jié)構(gòu)張量可以有效地提升圖像中沿邊緣切線方向區(qū)域的平滑程度.

3 基于結(jié)構(gòu)張量和各向異性平滑的DTI去噪

在本節(jié)中,詳細說明了本文提出的基于結(jié)構(gòu)張量和各向異性平滑的DTI去噪算法的實現(xiàn)過程.算法的主要步驟包括:(1)結(jié)構(gòu)張量對含噪DTI圖像的提取;(2)對含噪DTI圖像進行各向異性平滑.具體來說,首先利用結(jié)構(gòu)張量從DTI圖像像素中提取結(jié)構(gòu)信息,并得到圖像的均勻平坦區(qū)域和邊緣輪廓區(qū)域.然后在結(jié)構(gòu)內(nèi)部即平坦區(qū)域各向同性地平滑圖像,在結(jié)構(gòu)邊界即邊緣輪廓區(qū)域各向異性地平滑圖像.最終得到去噪后的DTI圖像.

3.1 結(jié)構(gòu)張量對含噪DTI圖像的提取

我們將2.1節(jié)中的結(jié)構(gòu)張量矩陣E記作

(7)

其中Ix,Iy就是對原圖像在x和y方向求得的偏導.然后我們求出矩陣E的行列式K=det(T)和跡H=tr(T),并且根據(jù)行列式K和跡H的關(guān)系來區(qū)分圖像的均勻平坦區(qū)域和邊緣輪廓區(qū)域.區(qū)分原則為:當H=0時,認為結(jié)構(gòu)張量提取的是圖像的均勻平坦區(qū)域:當H>0 &K=0時,則認為是邊緣輪廓區(qū)域.

3.2 對含噪DTI圖像進行各向異性平滑

(8)

為了在結(jié)構(gòu)內(nèi)部各向同性地平滑圖像,在結(jié)構(gòu)邊界各向異性地平滑圖像,結(jié)構(gòu)張量T由上述梯度張量G構(gòu)成,將T定義為G的歸一化倒數(shù).vg1,vg2,vg3表示G的特征向量,λg1,λg2,λg3表示G的相應(yīng)的特征值,T的特征向量和特征值是:

vti=vgi,i=1,2,3

(9)

λti=1/λgi,i=1,2,3

(10)

根據(jù)上述公式可以看出,在均勻區(qū)域(內(nèi)部束狀結(jié)構(gòu))中,結(jié)構(gòu)張量與梯度張量的三個特征值的大小是相同的,并且沿著所有方向的各向同性平滑結(jié)果是相似的.在結(jié)構(gòu)邊界,垂直于邊界方向(較大強度梯度的vt1方向)上的結(jié)構(gòu)張量的特征值是較小的,沿著邊界的方向(較小強度梯度的vt3方向)上的結(jié)構(gòu)張量的特征值是較大的,因此圖像結(jié)構(gòu)在沿著結(jié)構(gòu)邊界的切線方向上比垂直方向更平滑(各向異性平滑).為了實現(xiàn)在均勻區(qū)域和結(jié)構(gòu)邊界上產(chǎn)生同樣的平滑,結(jié)構(gòu)張量的特征值被歸一化為常數(shù)C,即λt1+λt2+λt3=C,因此無論是各向同性還是各向異性,平滑的總量在整個圖像上是相同的.

圖1 算法設(shè)計流程圖Fig.1 Flow chart of our algorithm

為了總結(jié)本文所提出的算法的實現(xiàn)過程,對算法步驟給出了簡明地描述,算法設(shè)計流程圖如圖1所示.

基于結(jié)構(gòu)張量和各向異性平滑的DTI去噪算法步驟:

步驟1.構(gòu)造結(jié)構(gòu)張量矩陣E,求出矩陣E的行列式K和跡H,然后根據(jù)行列式K和跡H的關(guān)系來區(qū)分圖像的均勻平坦區(qū)域和邊緣輪廓區(qū)域.

步驟3.使用公式(9)和公式(10)構(gòu)造結(jié)構(gòu)張量T.結(jié)構(gòu)張量T應(yīng)為一個3×3的對稱正定矩陣,然后分析并計算它的特征值和特征向量,從而實現(xiàn)在結(jié)構(gòu)內(nèi)部的均勻平坦區(qū)域各向同性地平滑圖像,在結(jié)構(gòu)邊界的邊緣輪廓區(qū)域各向異性地平滑圖像,最終得到去噪后的DTI圖像.

4 實驗結(jié)果和分析

本文算法及對比算法使用MATLAB編寫,實驗平臺為Intel Core i5 3.20GHz CPU和4.00 GB RAM的工作站.本文首先使用真實大腦的數(shù)據(jù)來驗證結(jié)構(gòu)張量是否能很好的提取圖像的平坦區(qū)域和邊緣區(qū)域.在實驗中,先求出第3節(jié)中提出的結(jié)構(gòu)張量矩陣T的行列式K和跡H,并根據(jù)行列式K和跡H的區(qū)分原則得出平坦區(qū)域和邊緣區(qū)域的分割結(jié)果.圖2給出了利用結(jié)構(gòu)張量區(qū)分出來的平坦區(qū)域和邊緣區(qū)域.這表明結(jié)構(gòu)張量可以有效提取圖像中像素并加以區(qū)分.

利用結(jié)構(gòu)張量將圖像區(qū)分并且得到平坦區(qū)域和邊緣區(qū)域之后,對含噪DTI圖像應(yīng)用本文所提出的各向異性平滑算法進行處理.為了驗證本文所提出的各向異性平滑算法在邊緣保留和去噪方面的良好性能,將本文所提的各向異性平滑算法與文獻[17]中基于復剪切波及復擴散各向異性濾波(CST-CDAF)的DTI去噪算法和文獻[18]中基于半隱式(S-I)方案的各向異性平滑算法進行比較.

圖2 結(jié)構(gòu)張量提取結(jié)果Fig.2 Extraction results of structure tensor

首先用模擬的DTI數(shù)據(jù)集來測試本文所提出的各向異性平滑算法的有效性.數(shù)據(jù)集由兩塊不同的“纖維”束組成,每個塊包含沿著不同取向的“纖維”片段.“纖維”的平均擴散系數(shù)(D)為0.7×10-5cm2/s,部分各向異性(FA)為0.9,與人腦體內(nèi)測量值相似.

上述各去噪算法的去噪結(jié)果如圖3所示.其中圖3(a)為模擬 “纖維”分布結(jié)構(gòu)的DTI圖像,圖3(b)為對模擬DTI數(shù)據(jù)添加零均值、SD=0.10的高斯噪聲的結(jié)果,圖3(c)為CST-CDAF算法的去噪結(jié)果,圖3(d)為S-I算法的去噪結(jié)果,圖3(e)為本文所提算法的去噪結(jié)果.可以看出在平滑后,由于噪聲導致的主擴散方向(PDD)的變化程度明顯降低,并且在任何“纖維”邊界附近沒有觀察到模糊現(xiàn)象,這也證明了我們所提出的算法在邊緣保留和去噪方面的有效性.

圖3 各向異性平滑算法對模擬DTI數(shù)據(jù)的影響Fig.3 Influence of anisotropic smoothing algorithm onsimulated DTI data

為了測試本文所提出的各向異性平滑算法對人腦DTI數(shù)據(jù)的去噪性能,使用磁共振掃描儀從健康的受試者大腦中獲得擴散加權(quán)圖像.DTI數(shù)據(jù)是基于自旋回波單次激發(fā)EPI成像獲得的,每次掃描中的成像參數(shù)b=1000sec/mm2,并且掃描體積為250×250×120mm2.掃描所獲取的數(shù)據(jù)矩陣的大小為128×128×34,然后運用插值法將其轉(zhuǎn)換為256×256×34.在掃描過程中,本文進行多次重復掃描并求其平均值,最終產(chǎn)生信噪比為SNR≥70%的數(shù)據(jù)集.選擇數(shù)據(jù)集中的兩個代表性區(qū)域分別應(yīng)用上述各去噪算法,并且與模擬數(shù)據(jù)中的參數(shù)保持一致.兩個代表性區(qū)域及其PDD結(jié)構(gòu)如圖4所示,第一個區(qū)域含有白質(zhì),其中存在具有不同方向的纖維束.第二個區(qū)域包含胼胝體的一部分,其中可以看見大部分平行的神經(jīng)纖維.

圖4 數(shù)據(jù)集中兩個代表性區(qū)域及其PDD結(jié)構(gòu).左為掃描后獲得的人腦擴散加權(quán)圖像,右上方為人腦中含有白質(zhì)區(qū)域的PDD結(jié)構(gòu),右下方為人腦中含有胼胝體區(qū)域的PDD結(jié)構(gòu)Fig.4 Two representative regions of the dataset and their PDD structure.The left image is a diffusion weighted image of the human brain obtained after scanning,the top right image is the PDD structure containing the white matter region in the human brain,and the bottom right image is the PDD structure containingthe corpus callosum region in the human brain

上述各去噪算法對人腦中兩個代表性區(qū)域的去噪結(jié)果如圖5和圖6所示.從圖中的實驗結(jié)果可以看出,本文所提出的各向異性平滑算法明顯降低了噪聲對兩個區(qū)域中主擴散方向(PDD)的影響,并且不會模糊結(jié)構(gòu)邊界.

圖5 各向異性平滑算法對人腦中含有白質(zhì)區(qū)域的DTI數(shù)據(jù)的影響Fig.5 Influence of anisotropic smoothing algorithm on DTI data containing white matter regions in human brain

為更加準確的說明本文所提算法的有效性,我們將上述各向異性平滑算法對PDD的影響進行定量評估.平滑后PDD的改善程度通過平均角度差改變的百分比來衡量,其定義為

(11)

其中θ0,θN分別是平滑之前和N次平滑之后DTI數(shù)據(jù)中PDD的角度差.上述各去噪算法對人腦中兩個代表性區(qū)域PDD的平均角度差變化的影響如圖7所示,圖中顯示了在大腦兩個代表性區(qū)域中PDD的平均角度差隨迭代次數(shù)變化的情況.從圖中可以很容易看出基于CST-CDAF去噪算法出現(xiàn)PDD的平均角度差隨著迭代次數(shù)的增加而增大的情況,這表明在平滑過程中出現(xiàn)了結(jié)構(gòu)模糊現(xiàn)象.基于S-I去噪算法和本文所提算法均呈現(xiàn)出PDD的平均角度差隨迭代次數(shù)的增加而減小的趨勢.對本文所提的各向異性平滑算法的定量分析表明,兩個代表性區(qū)域PDD的平均角度差均隨著迭代次數(shù)的增加而減小,其中迭代次數(shù)N= 1、2、4和 8對應(yīng)的兩個區(qū)域的平均角度差改善的百分比分別是27%、40%、60%和92%.綜上所述,定量評估的結(jié)果表明本文所提的各向異性平滑算法在去噪和保留結(jié)構(gòu)邊界等方面優(yōu)于CST-CDAF算法,并且可以得到和S-I算法不相上下的結(jié)果.

圖6 各向異性平滑算法對人腦中含有胼胝體區(qū)域的DTI數(shù)據(jù)的影響Fig.6 Influence of anisotropic smoothing algorithm on DTI data containing corpus callosum area in human brain

圖7 定量比較三種不同的平滑算法對人腦中兩個代表性區(qū)域的實驗結(jié)果Fig.7 Experiment results of the quantitative comparison between three different smoothing algorithms on the two representative regions of the human brain

5 總 結(jié)

在本文中,我們在結(jié)合結(jié)構(gòu)張量和各向異性平滑技術(shù)的基礎(chǔ)上提出了一種新的DTI去噪算法.該算法將結(jié)構(gòu)張量區(qū)分圖像結(jié)構(gòu)的特性和各向異性平滑的邊緣保留能力完美的結(jié)合在一起.基于模擬和體內(nèi)數(shù)據(jù)的實驗表明,本文所提出的去噪算法在邊緣保留和去除噪聲方面都顯示出了很好的性能,對醫(yī)學圖像后期處理和病理分析產(chǎn)生深遠的影響.