全张量步长适应性脑白质纤维跟踪改进算法

212　北京师范大学学报(自然科学版)Journal of Beijing Normal University (Natural Science )2008204　44(2)全张量步长适应性脑白质纤维跟踪改进算法3刘晓东1)　蔡子凡2)　王婷婷2)　姚　力2) 卢　洁3)(1)北京师范大学认知神经科学与学习国家重点实验室,100875,北京;2)北京师范大学信息科学与技术学院,100875,北京;3)首都医科大学宣武医院医学影像学部放射科,100053,北京)摘要　由于成像设备及技术本身的限制,以及纤维跟踪算法的局限,准确刻画纤维结构遇到许多困难.本文在张量偏转(tensor deflection ,TEND )算法基础上,提出了一种基于全张量及历史跟踪信息的适应性步长纤维跟踪算法.该算法结合纤维跟踪历史,利用全张量信息与相应的适应性步长值来计算纤维跟踪的前进方向.关键词　磁共振成像;弥散张量成像;张量偏转;纤维跟踪;步长3国家自然科学基金资助项目(60628101);北京市自然科学基金资助项目(4061004) 通信作者收稿日期:2007207210 分子扩散是人体重要的生理活动,是体内重要的物质转运方式.分子的扩散是唯一能够使得我们探索物质微细结构的渠道,弥散对磁共振信号的影响也很早就提了出来[1],这是弥散磁共振成像的关键所在[2].生物组织中水分子的扩散能力在不同方向存有差异.在脑白质纤维中,沿着其长轴方向扩散的水分子由于阻力较小而扩散较快;而垂直于其长轴方向扩散的水分子由于阻力较大而扩散较慢,生物体内这样的现象称为各向异性.弥散张量成像(diff usion tensor imaging ,D TI )的生物物理学基础即是大脑白质中水分子的各向异性弥散[325].磁共振弥散张量可获得弥散的方向性信息,弥散张量纤维跟踪(diff usion tensor t ractograp hy ,D T T )即是利用该信息,假设弥散张量成像中的最大特征值对应的特征向量方向代表局部占优势的纤维轴索的走行方向,即纤维走向与最大特征弥散性相关的特征向量的方向共线.因而得以无创地重建人脑的白质纤维,实现活体测量脑白质纤维走形方向和显示纤维结构的完整性.然而,实施白质的纤维跟踪及结果的准确呈现受到诸多局限和实际困难,例如,磁共振成像设备与成像技术自身的限制;图像噪音影响和纤维跟踪算法的缺陷;弥散张量磁共振成像(diff usion tensor 2magnetic resonance imaging ,D T 2MRI )的数据是离散的,而所要追踪的纤维是连续的;相对较低的MRI 空间分辨率所引起的部分容积效应会影响D T 2MRI 的数据质量等等.鉴于此,在TEND 算法的基础上[6],本文提出一种更为可靠且具有较强灵活适应性的白质纤维跟踪算法———全张量步长适应性脑白质纤维跟踪算法(f ull tensor fiber tractograp hy wit h adaptive stepping ,F TF TAS ).1　纤维跟踪理论白质纤维是大脑功能区间联系的物质基础,精确描绘白质纤维的走向及关系,有助于加深对正常与非正常脑功能的理解,与功能性磁共振技术的结合,已经成为认知神经科学领域关于大脑结构与功能关系的最前沿的联合应用.神经解剖学采用侵入式技术研究活体动物[7],而人类相关研究数据则受到更多的限制.目前,白质纤维束跟踪是唯一能在活体非侵入式地显示白质纤维束走向及关系的成像技术,为研究脑和脊髓白质的空间结构开辟了崭新的途径,对于脑的发育、白质纤维的髓鞘化过程、先天性与获得性脑白质疾病的研究提供了新的应用前景,同时可以作为术前计划和术后评价的辅助工具.目前,有多种白质纤维跟踪的理论和计算方法,大致可分为以下2类:1)基于经验主义的线性延展法(line p ropagation techniques ).该方法以局部张量信息作为扩展的核心,将体素内最大特征向量连接起来,例如连续跟踪纤维分配(fiber assignment by continuous t racking ,FAC T )[8],张量线(tensor line )[9],EZ 2t racing [10]以及TEND [11],在此基础上衍生出其他改进的数值理论方法,包括欧拉方法(Euler )[12],龙格2库塔法(Runge 2Kutta )[12]以及弥散传递模型(diff usion 2convection　第2期刘晓东等:全张量步长适应性脑白质纤维跟踪改进算法213 model )[13]等.2)能量最小法(energy minimization techniques ).应用最小能量发现2个预先设定的像素间的最佳通道,如快速行进法(fast marching technique ,FM T )[14]、螺旋等级模型(spin class model )[15]、模拟退火法(simulated annealing approach ,SAA ),以及线性态空间模型(linear state space model )[16]等等.上述方法,各有其优势及局限,但有几个主要的因素影响白质纤维的跟踪和计算显示.首先,成像设备及成像数据采集时的各种干扰因素所造成的噪声会影响λ1、λ2、λ3计算的准确性,进而引起主特征向量方向的误差累积;其次,随着各向异性程度的降低,主特征向量的不确定性增加,这意味着在某些弥散张量未表示出强烈方向性信息的时候,跟踪会出现差误;同时,由于D TI 空间分辨率一般是mm 量级,远高于生物组织细胞的尺度(≤20μm ),在纤维束整齐规则排列的主干部分或是纤维束直径大于图像的体素大小时,体素弥散张量的主特征向量可近似估计局部的纤维方向,但在纤维交叉、汇合、分散处,一个体素仅能反映出其所包含的多方向或多类型纤维的平均弥散张量,跟踪就自然受到部分容积效应(partial volume effect s ,PV E )的影响.由于存在以上问题,跟踪结果会与真实的纤维走向有所偏差,除了外在因素的不断提升外,积极的改进方法也逐步提高了白质纤维跟踪的真实性和精确性.2　全张量步长适应性脑白质纤维跟踪算法211　TE ND 纤维跟踪方法　利用张量矩阵使相邻体素输入纤维方向向目标体素的主向量方向进行偏转,TEND 算法充分考虑了历史跟踪信息和当前体素弥散张量所有特征值的影响[6],表示如下:v out =D ・v in ,输入向量可用弥散张量3个特征向量表示:v in =α1e 1+α2e 2+α3e 3,因D e i =λi e i (λ1≥λ2≥λ3,i =1,2,3),等式可重写为:v out =λ1(α1e 1+λ2λ1α2e 2+λ3λ1α3e 3).采用TEND 纤维跟踪算法,纤维走向角度突变弱化,重建出的纤维更平滑且较长;穿过同一体素不同方向的纤维会根据其历史走向向该体元主向量方向偏转;追踪穿过各向异性程度较低的区域时,改变跟踪结束条件的域值,重建的纤维长度差距不大,相对比较稳定.但TEND 不足之处在于对所有的弥散张量都不加区别地采取偏转处理,对各向异性程度较高的区域,很可能使纤维跟踪方向出现偏差,并造成误差的累积;第二,该算法经过每一体素时单次偏转,相对白质纤维尺度而言对精细度的处理仍不足.212　全张量步长适应性纤维跟踪方法(FTFT AS )　对于TEND 算法,在每个体素内经过n 次偏转后输出向量的方向应为:v out =D n・v in .步长是每次偏转操作后纤维跟踪推进的长度,约为1/n 个体素宽度,输出向量与主特征向量方向靠近的程度与n 值成正比,即步长对偏转程度具有决定性影响.采用TEND 方法,对直线纤维束跟踪时,受噪音及各向异性程度的影响较小,但对于弯曲的纤维则会低估了其轨迹曲率.那么,减小纤维跟踪的步长显然也不能解决这个问题,这对于直线纤维会减小误差,但是增加了跟踪弯曲纤维时的系统误差和计算量.因此,基于以上步长和体素内偏转次数对跟踪的影响,提出基于弥散张量各向异性特征的全张量步长适应性跟踪方法.根据弥散张量特征值关系,可定义C L 、C P 及C S 如下:分别表示线型(linear case )、扁圆型(planar case )及球型(sp herical case )3种张量的形状特征,如图1(从左至右).图1　3种弥散张量形状特征同时可将弥散椭圆体分为2类,即线性各向异性(C L )和非线性各向异性(C P ,C S ),同时将特征值向量(λ1,λ2,λ3)标准化为(λ′1,λ′2,λ′3)=(1,λ2/λ1,λ3/λ1),定义S 如下:C L =λ1-λ2λ1(λ1µλ2,λ3),C P =λ2-λ3λ1(λ1≈λ2µλ3),C S =λ3λ1(λ1≈λ2≈λ3),S =2λ′21+λ′2=2(1-C L )2-C L . 体素内纤维跟踪的实际步长为体素宽度的S 倍.由上式可知,λ′2反映了弥散椭圆体的线型特征并直接决定纤维跟踪的步长大小.对于线性弥散张量,λ′2较小,各向异性程度较高,小步长使得体素内多次跟踪偏转以后,输出向量方向更靠近主特征向量方向,增强了对确定性纤维走向信息的呈现;而非线性弥散张量,其　214　北京师范大学学报(自然科学版)第44卷　各向异性程度较低,纤维跟踪易受噪声影响,尤其是在纤维交叉融合区域,此时采用大步长策略减少跟踪次数可有效地抑制不确定走向信息的呈现.体素内纤维跟踪采用上述步长适应性机制,如图2二维示例.输入向量v in =α1e 1+α2e 2+α3e 3,A 点为进入体素跟踪的初始点,O 点坐标为(0,0,0),那么在特征向量方向为(e 1,e 2,e 3)坐标系下A 点坐标为(α1,α2,α3),输入与主特征向量夹角arcco s (v in ・e 1),n 次偏转后,输出与主特征向量夹角arcco s (v out ・e 1),输出与输入向量方向夹角arcco s (v out ・v in ),分别如下:arcco s (v in ・e 1)=arccos (α1α21+α22+α23),arcco s (v out ・e 1)=arcco s (λn1α1α21λ2n 1+α22λ2n 2+α23λ2n3),arcco s (v out ・v in )= arcco s (α21λn 1+α22λn 2+α23λn3α21+α22+α23α21λ2n 1+α22λ2n 2+α23λ2n3). 虚线为体素边界.OA 为进入体素向量方向(O 为起始点,A 为进入体素时边界点),B 为偏转中间点,C 为体素内偏转终点.步长|A B |=|B C |=S ×L ,其中L 为体素边长.图2　改进TE ND 纤维跟踪的二维示例　考虑到体素内弥散张量的全部特征值信息,并采用上述步长适应性跟踪机制,与固定步长跟踪相比具有如下特点:1)如图3,输入、输出向量方向偏转角度的大小随输入向量方向与主特征值方向的夹角而变化,适应性步长相较固定步长机制而言,历史跟踪方向信息的作用体现的更加明显,输出向量方向更靠近输入向量方向,弱化了体素本身主特征值向量方向对跟踪偏转的影响.2)如图4,对于输入向量方向与主特征值方向不同的初始夹角,在高各向异性(fractional anisot ropy ,FA )值时(λ′2小,弥散张量线性程度高),输出与输入向量偏转角较大,更接近于主特征值向量方向;而在低FA 值时,输出与输入向量偏转角较小,更接近于输入向量方向,即历史跟踪方向.3)如图4,对于输入向量方向与主特征值方向不同的初始夹角,无论FA 值高低,输出与输入向量偏转角都小于固定步长跟踪时,即输出向量方向更贴近历史跟踪方向,而非一味向主特征向量方向偏转,增加了纤维跟踪的平滑性.固定步长S =012,适应性步长S =015,其中λ′1=1,λ′2=1/3,λ′3=012.图3　步长机制与偏转角变化的关系　 输入与主特征向量初始夹角方向分别为30°(令α1=3,α2=α3=1/2)和60°(令α1=1,α2=α3=3/2),λ′1=1,λ′3=012.图4　λ′2与向量夹角大小变化的关系(反映适应性步长对偏转角度的影响)　因此,适应性步长机制首先基于历史跟踪方向信息,其次充分体现弥散张量的线形特征,进而综合调整纤维走向的偏转,在实际跟踪中判断边界及步长关系以控制进入下一体素的纤维走向,从理论上达到更好的跟踪效果.3　讨论与总结从算法理论上看,F TF TAS 与传统算法相比具有　第2期刘晓东等:全张量步长适应性脑白质纤维跟踪改进算法215以下优势:1)F TF TAS充分考虑了历史跟踪信息及噪声对非线性扁长型弥散椭圆体的影响,每一次的迭代递进都包含有历史信息,体现了纤维的实际联系性;2)步长随2类弥散张量形状适应性变化,达到与实际的跟踪偏转角度和方向相近似的效果;3)F TF TAS 算法可增加纤维平滑度和精细度,减少纤维偏转角度突变.但是,F TF TAS算法在体素内计算时有n次偏转,并计算整个弥散张量,这在一定程度上增加了数据的计算量.对于白质纤维跟踪仍有许多要继续考证的问题.首先,无论使用何种跟踪技术,由于噪声的影响和部分容积效应,使3D重构结果总会包含一些数据误差.另外,活体人脑白质纤维成像没有金标准,这限制了D T T成像的发展.其次,精度问题的核心在于从一个像素到另一个像素推算出纤维走向的连续性.解决这一问题需要设计更为准确的神经纤维传导方向甄选算法以提高连续性.从物理或生物学的原理上看此问题,则需寻求神经纤维束方向场的几何性质,这将从根本上得以解决.第三,目前弥散张量方法的形式表达不能检测到几个肌束穿过同一像素的情况.第四,对于纤维交叉的情形,已有其他相关方法取得进展.Van Wedeen等人[17]提出弥散频谱影像(diff usion spect rum imaging,DSI),采用515个方向的弥散加权成像(b=0～10000s・mm-2)来侦测3维空间不同扩散距离的水分子,并将所有方向的扩散讯号填入q2space(描述三维空间水分子扩散讯号的空间),利用一些数学的运算便可以画出水分子的扩散方向函数(orientatio n distribution f unction,ODF),并正确地描绘出纤维交叉的情形.另外,QBI(Q2ball imaging)[18]采用252个方向的弥散加权成像(b= 4000s・mm-2),均匀分布在q2space中的一个球面上,同样可获得水分子的扩散方向函数,并描绘出纤维交叉的情形.利用DSI或QB I的数据进行大脑纤维跟踪,可免除纤维交叉的困境,也可大大提高纤维追踪的准确度,然而这2种方法所面临的问题是扫描时间很长,DSI 的扫描时间约1h,QB I的扫描时间则约015h,若能降低扫描时间达临床所能接受的范围,那么便可以提高其在临床上应用的可行性.4　参考文献[1]　Stejskal E O,Tanner J E.Spin diff usion measurements:spin echoes in the presence of a time2dependent fieldgradient[J].Chem Phys,1965,42:288[2]　Holmes A A,Scollan D F,Winslow R L.Directhistological Validation of diff usion tensor in formaldehyde2fixed myocardium[J].Magn Reson Med, 2000,44:157[3]　Basser P J,Mattiello J,LeBihan D.MR diff usion tensorspectroscopy and imaging[J].Biophysical Journal, 1994,66(1):259[4]　Basser P J,Mattieflo J,LeBihan D.Estimation of theeffective serf2diff usion tensor f rom the NMR spin echo [J].Journal of Magnetic Resonance,Series B,1994,103(3):247[5]　LeBihan D,Mangin J F,Poupon C.Diff usion tensorimaging:concepts and applications[J].Journal of Magnetic Resonance Imaging,2001,13(4):534[6]　Lazar M,Weinstein D M,Tsuruda J S,et al.Whitematter tractography using diff usion tensor deflection[J].Human Brain Mapping,2003,18:306[7]　Heimer L,Robards M J.Neuroanatomical tract2tracingmethods[M].New Y ork:Plenum Press,1981[8]　Mori S,Crain B J,Chacko V P,et al.Three2dimensionaltracking of axonal projections in the brain by magnetic resonance imaging[J].Annals of Neurology,1999,45(2):265[9]　Weinstein D M,K indlmann G L,L undberg E C.Tensorlines:advection2diff usion based propagation through diff usion tensor fields[C],IEEE Visualization Proceedings.San Francisco,USA,1999:2492253 [10]　Kenshi Terajima,Tsutomu Nakada.EZ2tracing:a newready2to2use algorithm for magnetic resonancetractography[J].Journal of Neuroscience Methods,2002(116):147[11]　Lazar M,Weinstein D M,Tsuruda J S,et al.Whitematter tractography using diff usion tensor deflection[J].Human Brain Mapping,2003,18(4):306[12]　Basser P J,Sinisa Pajevic,Carlo Pierpaoli,et al.Invivo fiber tractography using D T2MRI data[J].Magnetic Resonance in Medicine,2000,44:625 [13]　Batchelor P G,Hill D L,Atkinson D,et al.Fiber2tracking by solving the diffusion2convection equation[C]∥Proceedings of ISMRM,2002:1135[14]　Parker G J,Wheeler2K ingshott C A,Barker G J.Estimating distributed anatomical connectivity using fastmarching methods and diff usion tensor imaging[J].IEEE Trans Med Imaging,2002,21(5):505[15]　Mangin J F,Poupon C,Cointepas Y,et al.Af ramework based on spin glass models for the inferenceof anatomical connectivity f rom diff usion2weighted MRdata2a technical review[J].NMR Biomed,2002,15(7/　216　北京师范大学学报(自然科学版)第44卷　8):481[16]　G ssl C,Fahrmeir L,Putz B,et al.Fiber tracking f romD TI using linear state space models:detectability of thepyramidal tract[J].Neuroimage,2002,16(2):378 [17]　Van Wedeen J,Patric Hagmann,Wen2Y ih Isaac Tseng,et al,Mapping complex tissue architecture withdiff usion spectrum magnetic resonance imaging[J].Magnetic Resonance in Medicine,2005,54:1377 [18]　David S Tuch.Q2Ball imaging,[J].MagneticResonance in Medicine,2004,2:1358AN IMPROVED METH OD FOR FU LL TENSOR FIBERTRACTOGRAPH Y WITH ADAPTIVE STEPPINGL IU Xiaodong1)　CA I Zifan2)　WAN G Tingting2)　YAO Li2)　L U Jie3)(1)State Key Laboratory of Cognitive Neuroscience and Learning,Beijing Normal University,100875,Beijing,China;2)College of Information Science,Beijing Normal University,100875,Beijing,China;3)Depart ment of Radiology,Xuanwu Hospital,Capital University of Medical Sciences,100053,Beijing,China)Abstract　Due to limitations of imaging devices and techniques,and of fiber t ractograp hy algorit hms,it has been difficult to describe organization of white matter fiber accurately.An derivative tensor deflection algorit hm based on f ull tensor and historic tracking information was int roduced.This met hod could describe more accurately t he direction of fiber st ruct ure.K ey w ords　magnetic resonance imaging;diff usion tensor imaging;TEND;fiber t ractograp hy;stepping size我校成立全球变化与地球系统科学研究院日前,学校在主楼第三会议室召开全球变化与地球系统科学研究院成立暨全球变化与地球系统科学学科建设研讨会.会议宣布了学校成立“全球变化与地球系统科学研究院”的决定.刘川生书记、董奇副校长、史培军副校长、葛剑平副校长、科技部原部长徐冠华院士、刘昌明院士、张新时院士、符淙斌院士、李小文院士、中国气象局气候中心主任董文杰研究员、国家自然科学基金委员会地学部地学一处处长宋长青研究员以及20余位来自中国科学院、中国农业科学院、中国林业科学院、中国水利水电科学研究院、南京大学、同济大学、南京信息工程大学、香港中文大学的领导和专家出席了会议.我校学科规划与建设处处长陈丽教授、科学技术处处长高尚玉教授、地理学与遥感科学学院院长戴永久教授、全球变化与地球系统科学研究院常务副院长潘耀忠教授和地理学与遥感科学学院、民政部教育部减灾与应急管理研究院、物理系的专家参加了会议.刘川生书记致辞,对各位专家关心和支持我校全球变化与地球系统科学学科的建设表示感谢,希望各位专家为研究院的发展献言献策.徐冠华院士主持学科建设研讨会.会上,史培军副校长介绍了我校全球变化与地球系统科学学科的建设现状.与会专家针对该领域国际前沿及热点研究问题,面向我国可持续发展的重大需求,结合我校相关学科建设现状,为我校全球变化与地球系统科学学科建设的近期、中长期发展目标、重点研究内容提出了建议.(肖　希)。