基于結構張量和各向異性平滑的DTI去噪

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

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

1 引 言

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

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

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

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

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

2 基本理論

2.1 結構張量

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

初始結構張量ST(i,j)定義為

(1)

通常,我們在平滑處理時將結構張量定義為

(2)

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

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

2.2 各向同性平滑

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

(3)

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

(4)

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

2.3 各向異性平滑

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

(5)

(6)

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

3 基于結構張量和各向異性平滑的DTI去噪

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

3.1 結構張量對含噪DTI圖像的提取

我們將2.1節中的結構張量矩陣E記作

(7)

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

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

(8)

為了在結構內部各向同性地平滑圖像,在結構邊界各向異性地平滑圖像,結構張量T由上述梯度張量G構成,將T定義為G的歸一化倒數.vg1,vg2,vg3表示G的特征向量,λg1,λg2,λg3表示G的相應的特征值,T的特征向量和特征值是:

vti=vgi,i=1,2,3

(9)

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

(10)

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

圖1 算法設計流程圖Fig.1 Flow chart of our algorithm

為了總結本文所提出的算法的實現過程,對算法步驟給出了簡明地描述,算法設計流程圖如圖1所示.

基于結構張量和各向異性平滑的DTI去噪算法步驟:

步驟1.構造結構張量矩陣E,求出矩陣E的行列式K和跡H,然后根據行列式K和跡H的關系來區分圖像的均勻平坦區域和邊緣輪廓區域.

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

4 實驗結果和分析

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

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

圖2 結構張量提取結果Fig.2 Extraction results of structure tensor

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

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

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

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

圖4 數據集中兩個代表性區域及其PDD結構.左為掃描后獲得的人腦擴散加權圖像,右上方為人腦中含有白質區域的PDD結構,右下方為人腦中含有胼胝體區域的PDD結構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

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

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

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

(11)

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

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

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

5 總 結

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