一种基于层选择的弥散张量优化算法

刘 薇 周振宇 刘小征 严 序 杨 光 王遵亮 周永迪 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］