一種基于層選擇的彌散張量優化算法

劉 薇 周振宇 劉小征 嚴 序 楊 光 王遵亮 周永迪 Peterson Bradley S 徐冬溶，

(1華東師范大學腦功能基因組學教育部重點實驗室，上海 200062)

(2華東師范大學上海市磁共振重點實驗室，上海 200062)

(3哥倫比亞大學精神病學系，紐約 10032)

(4紐約州立精神疾病研究所，紐約 10032)

(5東南大學生物科學與醫學工程學院，南京 210096)

一種基于層選擇的彌散張量優化算法

劉 薇1，2周振宇3，4劉小征1，2嚴 序1，2楊 光2王遵亮5周永迪1Peterson Bradley S3，4徐冬溶1，2，3，4

(1華東師范大學腦功能基因組學教育部重點實驗室，上海 200062)

(2華東師范大學上海市磁共振重點實驗室，上海 200062)

(3哥倫比亞大學精神病學系，紐約 10032)

(4紐約州立精神疾病研究所，紐約 10032)

(5東南大學生物科學與醫學工程學院，南京 210096)

為了從帶偽影的彌散加權數據中重建更準確的張量，提出了一種基于層選擇的張量優化算法.首先對彌散加權圖像中常見的3種偽影(波狀、層間運動和對比度偽影)進行定性分析，分別提取3種對應的特征(小波標記、相似度和相關性)來識別這些偽影，從而區分出正常圖層和偽影圖層.然后，利用正常圖層數據進行張量重建.模擬實驗結果驗證了這3種特征對彌散加權圖像中相關偽影判斷的有效性，且對波狀偽影和層間運動偽影的判斷率均大于90%.真實數據實驗結果表明，與類似算法相比，所提算法可以更好地改善部分各向異性偽彩圖中的偏色現象，在白質結構分析中提供了更準確的方向信息.

彌散加權成像;彌散張量成像;張量優化;小波標記;互信息;灰度共生矩陣

彌散張量成像(diffusion tensor imaging，DTI)是一種測量不同腦組織區域內水分子局部彌散特性的新磁共振成像技術，可以用來顯示并定量分析腦白質結構.重建DTI圖像需要至少6組不同彌散梯度方向對應的彌散加權成像(diffusion weighted imaging，DWI)圖像［1］以及至少1組不施加彌散梯度的圖像數據.而在實際操作中，DWI數據的采集方向個數往往大于6，因此一般認為DWI數據具有冗余性.當DWI圖像中存在無法校正的偽影時，常用的解決方法是直接拋棄整個DWI數據，或者直接去除含有偽影層的梯度方向對應的DWI體數據，僅靠剩余的梯度方向對應的DWI體數據進行張量估計.后者在拋棄DWI偽影數據的同時也浪費了其中可能包含的很多正常圖層數據.為了能夠利用盡可能多的正常DWI數據進行準確高效的張量估計，研究者們提出了基于體素選擇或層選擇的張量優化算法.

張量數據通常采用最小二乘法(least square，LS)估計得到［1］.假設DWI不同體數據中同一位置處的信號對應于同一個腦解剖位置，由于偽影造成的異常信號會破壞該位置處的張量模型.為此，研究者們提出了基于LS的張量優化算法，如GMM算法［2］、RESTORE 算法［3］和 PATCH 算法［4］等.這些基于體素選擇的優化算法需要先計算原始的張量數據，通過多次迭代將異常點的權值降到最小，從而得到最佳的張量估計，因此計算時間較長.如果異常數據過多，迭代的初始值便會存在較大誤差，從而使得其后的張量估計不再穩定［5］.

張量數據的生成過程實際上是將二維圖像轉變到一維空間后再進行處理的，會造成二維圖像原來的空間關系以及相關統計特征的丟失，而空間性特征往往是有效識別圖像紋理模式的依據，為此研究者們提出了基于層選擇的張量優化算法［6－7］.但是這類算法所提特征無法對不同的偽影類型提供針對性信息，因此對偽影層的判斷不夠準確.

本文提出了一種基于層選擇的張量優化算法.首先對DWI圖像中的常見偽影進行定性分析;然后將小波標記、相似度和相關性分別作為不同偽影類型的識別判據，對偽影層進行較為精細的判斷;最后，利用正常圖層進行張量重建.

1 基于偽影特性的圖像特征提取

根據DWI中常見偽影的不同特性，分別采用3種不同的特征提取方法:① 利用小波標記來檢測具有方向特性的偽影(如渦流引起的波狀偽影);②利用圖像相似度來檢測非穩定偽影(如由不同層之間頭動引發的層間運動偽影);③利用相關性特征來檢測由于儀器缺陷所造成的對比度偽影.結合這些特征，選擇合適的閾值即可從原DWI圖像中辨別出偽影圖層.

1.1 小波標記特征

通常采用小波的近似子帶系數來反映圖像的灰度、輪廓等基本特性，采用小波的高頻子帶來反映圖像的紋理信息.在此基礎上，Pi等［8］將小波能量標記作為一種有效分析圖像紋理的特征.小波子帶Bk的能量定義如下:

式中，wk，n為k級分解后的第n個小波系數;Nk為wk，n的個數.

Bk中小波系數的平均絕對值MBk也可作為一種特征，通常與EBk同時使用，從而更好地表示圖像的紋理特性［9］.其表達式為

小波標記可以用特征向量{EH1，MH1，EV1，MV1，ED1，MD1，…，EHk，MHk，EVk，MVk，EDk，MDk}來表示，其中EHk，EVk，EDk分別表示水平、豎直和對角方向上第k級小波子帶能量，MHk，MVk，MDk分別表示水平、豎直和對角方向上第k級小波子帶系數的平均絕對值.

小波細節系數能很好地體現圖像的方向特征，故適合于檢驗DWI圖像中具有方向性的偽影(如波狀偽影).本文采用小波基db8對圖像進行一級小波分解，得到相應的小波標記特征(見圖1).如圖所示，當DWI圖像存在明顯的波狀偽影時，小波分解后3個方向上的細節圖像和正常DWI圖像小波分解后的結果有明顯差別，尤其是對角線方向.

1.2 相似度特征

互信息(mutual information，MI)可用于比較不同模態圖像之間的相似程度，故可作為衡量DWI圖像和相應參考圖像之間相似程度的特征指標.在計算MI時，為了使更能體現圖像特性的像素占有較大權重，以提高相似度計算的準確性，Luan等［10］將顯著性作為權值，提出了定量定性互信息(quantitative-qualitative measure of MI，QMI).顯著性測量是一種通過模擬生物體視覺選擇，對視覺數據進行高效篩選的算法，一般對圖像的噪聲不敏感，可以較好地反映圖像的局部特性.與傳統的MI相比，QMI對圖像位置的變化更加敏感，因此可以作為衡量圖像之間是否存在運動變化的一項重要特征指標.本文將QMI作為相似度判據來檢測DWI圖像中是否存在來自于層間頭動的偽影.在圖像域Ω中，假設Q為DWI圖像與對應參考圖像(即未加彌散梯度時的數據)之間的QMI，即

圖1 目標圖像的一級小波分解示意圖

式中，R，T分別為參考圖像和DWI圖像;IR，IT分別為R和T的灰度值;u(IR，IT)和P(IR，IT)分別為(IR，IT)的聯合效用和聯合概率密度;PIR，PIT分別為IR和IT的邊緣概率密度.

假設ARl為R在位置l的像素顯著性，ATl'為T在位置l'的像素顯著性［11］，IRl為R在位置l的灰度，ITl'為T在位置l'的灰度，則灰度對(i，j)的聯合效用可表示為

1.3 相關性特征

灰度共生矩陣(gray-level co-occurrence matrix，GLCM)不僅可用于描述圖像灰度的分布特性，也可以反映具有同樣亮度或者相似灰度的像素之間的位置分布特性［12］.2D加權灰度共生矩陣(2DWGLCM)C的定義如下:

式中，G(x，y)，G(x'，y')分別表示大小為S的圖像在(x，y)，(x'，y')位置處的灰度，(x'，y')表示以點(x，y)為中心3×3區域中的相鄰點;δ［·］表示邏輯判斷函數，根據其邏輯表達式的真假分別返回1和0;f(·)表示加權函數，且f(·)=1/［(xx')2+(y－y')2］1/2.

利用2DWGLCM計算灰度相關性統計特征［12］，即可判定DWI中的對比度偽影.

為了提高計算效率，去除大部分噪聲的影響，QMI和2DWGLCM的計算都是在小波分解后得到的近似圖像中進行的.實際應用中，可以簡單地利用上述3種特征的均值μ和均方差σ來確定離群值的閾值，如采用?F=aμ±bσ來區分偽影數據，其中a，b分別為預先設定的特征F的經驗修正系數.對于小波標記和相關性特征，如果其值不在上下閾值之間，則可認為該層是偽影數據;而對于相似度特征，如果其值小于下閾值，則可認為該層是偽影數據.本實驗中，采用小波標記特征時，a和b分別設置為1.5和2;采用相關性特征時，a和b分別設置為1和1.45;采用相似度特征時，a和b分別設置為1和1.5.

2 實驗

2.1 模擬實驗

根據DWI圖像中常見偽影的特性，模擬得到包含波狀偽影、層間運動偽影以及對比度偽影的模擬數據.然后，分別提取模擬圖像的小波標記、相似度和相關性特征，并對相應偽影類型進行判斷.最后，將判斷結果和真實情況進行比較，評估所提特征對其對應偽影類型判定的有效性.

2.1.1 波狀偽影模擬

為了簡化實驗，在正常DWI數據(彌散梯度方向個數為25)的前5個體數據中，對不同圖層施加大小和間距隨機的黑色平行線來模擬波狀偽影(亦作平行線偽影).此偽影的大小在圖像中心點周圍200×200的范圍內，間距為5～20.

2.1.2 層間運動偽影模擬

在同一組正常DWI數據的前5個體數據中，對不同圖層加入不同的剛體變化來模擬層間頭動引起的偽影.其中，旋轉變化是指沿層方向隨機旋轉體數據中的圖層，旋轉角度為－5°～5°;平移變化是指沿層平面任意方向平移體數據中的圖層，平移－5～5像素.

2.1.3 對比度偽影模擬

在同一組正常DWI數據的前5個體數據中，通過隨機改變不同圖層的對比度來模擬對比度偽影.假設當灰度拉升參數為s時，圖像直方圖平移大小為其半高寬度的s倍.圖像變化后超出原始范圍的像素點被設置為原始范圍內最相近的值.當s＞0時，圖像整體變亮;當s＜0時，圖像整體變暗.

2.2 真實數據實驗

實驗采集了320組成人和160組嬰兒的DWI數據(彌散梯度方向個數分別為25和11)，分別利用傳統的LS算法和本文算法對DWI數據進行張量估計，并將結果與文獻［7］進行比較.通過計算張量數據的部分各項異性(fractional anisotropy，FA)偽彩圖和張量的非正定性程度，分析優化后的張量是否可以得到更精確的FA值，以及是否可以重建出正確的張量方向.理論上正常張量矩陣都是對稱正定的，而DWI圖像中的偽影會使估計得到的張量矩陣不再正定，該張量矩陣即為病態張量.因此，可以將病態張量比例作為評價DTI圖像質量的關鍵指標.

3 實驗結果

模擬實驗結果顯示，所提3種特征對其相應的偽影類型均具有較好的選擇性.重復實驗100次，對波狀、層間運動和對比度偽影層的判斷率分別為92%，93%，63%.值得注意的是，在模擬對比度偽影的實驗中，當圖像的亮度拉升參數時，圖像的對比度變化并不明顯，而在人工檢查真實DWI圖像時，這種圖像通常也不會被作為偽影層，因此，無法區分此范圍內對比度偽影是可以理解的.去除這一區間的模擬數據，本文算法判斷對比度偽影的正確率為81%.

在真實實驗中，由一組成人數據和一組嬰兒數據重建得到的FA偽彩圖分別見圖2和圖3.由圖可知，傳統的重建結果中均出現了偏色情況，這是由于當某一方向上DWI體數據中的像素值發生異變時，在此方向上便會得到錯誤的彌散系數，從而使得最后的張量方向發生改變.由圖2可以看出，利用文獻［7］中算法得到的偽彩圖仍存在少量偏色，而本文算法則可完全消除該層的偏色現象.由圖3可以看出，利用文獻［7］中算法重建得到的FA圖中出現了局部噪點(如虛線箭頭所示)，這可能是由于去除過多偽影層時張量估計出現偏差所導致的.本文算法均可得到和已知腦解剖結構相比最為清晰準確的白質結構(如實線箭頭所示).

圖2 成人DTI數據的FA偽彩圖

圖3 嬰兒DTI數據的FA偽彩圖

對所有采集數據進行張量估計，人為檢測這些數據是否可以被挽回，從而大致估算出不同的張量優化算法成功挽救回來的DTI數據比例.實驗中人為檢測主要是通過多位專家檢測FA偽彩圖是否存在偏色及白質結構是否清晰準確來實現的.此外，專家們也對原始數據中的不同偽影類型進行了評定.采用不同張量估計算法對所有采集得到的成人和嬰兒DWI數據進行張量重建，統計結果見表1.由表可知，本文算法可以有效地降低病態張量的數量，大幅度提高DTI數據的利用率.本文算法對波狀、層間運動和對比度偽影的判斷率分別為84.6%，81.3%和 67.5%.

表1 不同算法的張量優化結果統計 %

4 結語

張量優化算法是DTI中的一個研究熱點，對于無法自身控制其行為的被測者(如嬰兒、多動癥患者等)來說尤為重要.本文提出了一種基于層選擇的張量優化算法.模擬實驗和真實數據實驗結果表明，該算法可以較好地剔除冗余DWI數據中的偽影圖層，有效地提高張量重構質量.張量重構算法的效果會直接影響DTI圖像中白質纖維追蹤的精確性，因此這種新的張量優化算法為纖維追蹤的準確性提供了有力的保障.

［1］Basser P J，Mattiello J，LeBihan D.Estimation of the effective self-diffusion tensor from the NMR spin echo［J］.Journal of Magnetic Resonance B，1994，103(3):247－254.

［2］Mangin J F，Poupon C，Clark C，et al.Distortion correction and robust tensor estimation for MR diffusion imaging［J］.Medical Image Analysis，2002，6(3):191－198.

［3］Chang L C，Jones D K，Pierpaoli C.RESTORE:robust estimation of tensors by outlier rejection［J］.Magnectic Resonance in Medicine，2005，53(5):1088－1095.

［4］Zwiers M P.Patching cardiac and head motion artefacts in diffusion-weighted images［J］.Neuroimage，2010，53(2):565－575.

［5］Chang L C.Improving RESTORE for robust diffusion tensor estimation:a simulation study［C］//Proceedings of the International Society for Optical Engineering.San Diego，California，USA，2010:762328.

［6］Liu Z，Wang Y，Gerig G，et al.Quality control of diffusion weighted images［C］//Proceedings of the International Society for Optical Engineering.San Diego，California，USA，2010:76280J.

［7］Zhou Z，Liu W，Cui J，et al.Automated artifact detection and removal for improved tensor estimation in motion-corrupted DTI data sets using the combination of local binary patterns and 2D partial least squares［J］.Magnetic Resonance Imaging，2011，29(2):230－242.

［8］Pi M H，Tong C S，Choy S K，et al.A fast and effective model for wavelet subband histograms and its application in texture image retrieval［J］.IEEE Transactions on Image Processing，2006，15(10):3078－3088.

［9］Hiremath P S，Shivashankar S.Wavelet based co-occurrence histogram features for texture classification with an application to script identification in a document image［J］.Pattern Recognition Letters，2008，29(9):1182－1189.

［10］Luan H X，Qi F H，Xue Z，et al.Multimodality image registration by maximization of quantitative-qualitative measure of mutual information［J］.Pattern Recognition，2008，41(1):285－298.

［11］Harel J，Koch C，Perona P.Graph-based visual saliency［C］//Proceedings of Neural Information Processing Systems.Vancouver，Canada，2006:545－552.

［12］Haralick R M，Shanmuga K，Dinstein I.Textural features for image classification［J］.IEEE Transactions on Systems Man and Cybernetics，1973，3(6):610－621.

Slice-wise optimization algorithm for diffusion tensor estimation

Liu Wei1，2Zhou Zhenyu3，4Liu Xiaozheng1，2Yan Xu1，2Yang Guang2Wang Zunliang5Zhou Yongdi1Peterson Bradley S3，4Xu Dongrong1，2，3，4

(1Key Laboratory of Brain Functional Genomics of Ministry of Education，East China Normal University，Shanghai 200062，China)
(2Shanghai Key Laboratory of Magnetic Resonance，East China Normal University，Shanghai 200062，China)

(3Department of Psychiatry，Columbia University，New York 10032，USA)
(4New York State Psychiatric Institute，New York 10032，USA)
(5School of Biological Science and Medical Engineering，Southeast University，Nanjing 210096，China)

A slice-wise optimization algorithm for diffusion tensor estimation is proposed to improve the accuracy of estimating tensors by using diffusion weighted imaging(DWI)data which contain artifacts.Firstly，three types of common artifacts(wavelike，motion-between-slice and contrast artifacts)are qualitatively analyzed in DWI data.Three types of features(wavelet signature，similarity and correlation)are extracted to identify these three artifacts，respectively.Thus，a slice with or without artifacts can be distinguished.Then，tensors can be reconstructed by using the slices tagged without any artifacts.The simulation results show that the three features are effective to identify related artifacts in DWI data.A high discrimination capability(＞90%)can be achieved for identifying wavelike and motion-between-slice artifacts.The experimental results using real datasets demonstrate that，compared with other similar algorithms，the proposed algorithm can improve the bias found in color-encoded fractional anisotropy map more effectively，and can provide more accurate directionality information to analyze white matter structure.

diffusion weighted imaging;diffusion tensor imaging;tensor optimization;wavelet signature;mutual information;gray-level co-occurrence matrix

TP391.41

A

1001－0505(2013)01-0030-05

10.3969/j.issn.1001－0505.2013.01.006

2012-05-08.

劉薇(1983—)，女，博士生;徐冬溶(聯系人)，男，博士，教授，博士生導師，dx2103@columbia.edu.

上海市科學技術委員會資助項目(10440710200).

劉薇，周振宇，劉小征，等.一種基于層選擇的彌散張量優化算法［J］.東南大學學報:自然科學版，2013，43(1):30－34.［doi:10.3969/j.issn.1001－0505.2013.01.006］