（无线电物理专业论文）dti数据的重建与可视化.pdf

2 011s c h 0 0 1y e a rm a s t e r sd e g r e e1 1 1 e s i s i n s t i t u t i o nc o d e ：10 2 6 9 n u m b e r ：510 8 0 6 0 2 0 6 7 e a s tc h i n an o 肼a l u n i v e r s i t y d 印a n m e n t ：d e p t p h y s i c s m 旬o r ： r a d i op h y s i c s r e s e 2 u r c hd i r e c t i o n ：n u c l e a rm a g n e t i cr e s o n a n c es o 脚a r e t e c l u l i q u e s u p e r v i s o r ：1 y m gg u a n g g r a d u a t e ： w a n g l i l i m a y ，2 0 l 1 i i = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = ：= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = ：华东师范大学学位论文原创性声明 郑重声明：本人呈交的学位论文d t i 数据重建与可视化，是在华东师范 大学攻读硕士学位期间，在导师的指导下进行的研究工作及取得的研究成果。除 文中已经注明引用的内容外，本论文不包含其他个人已经发表或撰写过的研究成 果。对本文的研究做出重要贡献的个人和集体，均已在文中作了明确说明并表示 谢意。 作者签名：逊 日期：如1 1 年p 岁月够日 华东师范大学学位论文著作权使用声明 d t i 数据重建与可视化系本人在华东师范大学攻读学位期间在导师指导 下完成的硕士学位论文，本论文的研究成果归华东师范大学所有。本人同意华东 师范大学根据相关规定保留和使用此学位论文，并向主管部门和相关机构如国家 图书馆、中信所和“知网”送交学位论文的印刷版和电子版；允许学位论文进入 华东师范大学图书馆及数据库被查阅、借阅；同意学校将学位论文加入全国博士、 硕士学位论文共建单位数据库进行检索，将学位论文的标题和摘要汇编出版，采 用影印、缩印或者其它方式合理复制学位论文。 本学位论文属于( 请勾选) () 1 经华东师范大学相关部门审查核定的“内部 或“涉密”学位论 文木，于年月日解密，解密后适用上述授权。 a 2 不保密，适用上述授权。 导师签名 4 ， 爿、研 l i i i 本人签名至亟豳 伽f 1 年口夕月z 妇 王匦匦硕士学位论文答辩委员会成员名单 姓名职称单位备注 态馒甄亳毙吾眨孵纾争巷太芬抬l 箩曳物垤尔若韵耳甜 和萄 曷9 裱琨 。 释l 斥彬学瞧1 吃唬籀7 垤朱 癌盆目矮 蕃】敏髭 糯汗溅醒 龟氏构以朱 摘要 。 j d t i 作为一种新的m r 成像技术，能够利用水分子的扩散情况反映活体组织 的微观结构及病理组织的异常状态。通过对d t i 数据重建与可视化，得到组织 中神经纤维束的走行轨迹，能够非侵入式地对活体组织的微观结构进行观察与研 究，用于神经疾病的诊断及外科术前手术方案的制定。 通过运用c + + 语言和o p e i l g l 开发了一款d t i 数据重建与可视化的工具，工 具中包括多种可视化技术，如张量的计算与可视化，各向异性测度( f a v a ，m d 等) 的计算与可视化，图元( 椭球) 法可视化，多平面重建( m p r ) ，神经纤维 跟踪等。本文详细描述了这些可视化技术的实现过程，重点介绍了神经纤维跟踪 应用较多的几种算法_ f a c t 算法，t e n s o r l i n e 算法，b m t ef o r c e 算法。通过分 析这几种算法的利弊，将t e n s o r l i n e 算法与b m t ef o r c e 算法结合，得到了令人满 意的跟踪结果。 在以上工作的基础上，为了对d t i 可视化结果实时显示，运用了由n v i d 认 公司提供的g p u ( 图形处理器) 对跟踪算法进行加速计算。为了充分利用g p u 高效的并行计算的能力，文中对跟踪算法重新设计，提出了一种新的适合g p u 做并行处理的跟踪模型，将计算时间缩短到原来的十分之一。这种跟踪模型也适 用于其他跟踪算法的加速计算。 关键词：d t i 可视化，各向异性扩散，纤维跟踪，g p u 加速计算 a bs t r a c t d i f f i l s i o nt e n s o ri i i l a 酉n g ( d t i ) i saf a i r l yn e wm a g n “cr e s o n a n c ei n l 孵n g t e c h n i q u e i ts h o w sm i c r o s c o p i cs 饥l c t u r eo ft i s s u e a i l d1 e s i o nb ya n a l y z i n gt 1 1 e d i a h s i o no fw a t e rm 0 1 e c u l e si 1 1t h e m d t id a t ar e c o n 蛐m c t i o na n dv i s u a l i z a t i o nc a i l d i s p l a y 仃a c tn 匈e c t o d e si nt i s s u e ，w r h i c h c a nb eu s e di nn o n i n v a s i v ee x p l o r a t i o no f 也em i c r o s 旬m c t u r ei i l l i v i n gh m l a n ， d i a g n o s i so fn e u r o l o 百c a l d i s e a s e s柚d p r e o p e r a t i v es u 呼c a lp l a n 。av i s u a l i z a t i o nt 0 0 1f o rp r o c e s s i n ga n dv i s u a l i z a t i o nd t id a t ai sd e v e l o p e dw i t h c + + a n do p e n g l - i ti m p l e i i l e l l t sav 撕e t yo fv i s u a l i z a t i o nt e c ：h n 0 1 0 鼢f o ri n s t a l l c e ， v i s u a l i z a t i o nf o rt e n s o r ，a i l i s o 仃o p y ( f a ，v a ，m d ) m a p s ，p r i m i t i v er 印r e s e l l t a t i o n ， m u l t i _ p l a n a rr e c o n s t n l c t i o na i l d 舶e r 仃a c l ( i n g i n v 0 1 v e da l g o r i t l l l l l sa r ed e s 嘶b e di n d “l ，e s p e c i a l l y 丘b e r 仃a c k i n ga l g o 硼n ss u c ha sf a c t ，t h s o r l i n ea 1 1 db m t ef o r c e a l g o r i t h m s a d v a n t a g e sa i l dd i s a u a v a n t a g e so fd i 脑e 1 1 t 舶e r 仃a c k i n ga l g o d l m sa r e s u 础 1 1 撕z e d t 钿s o r l i l l ea n db r u t ef b r c ea l g o r i m m sa r ec o m b i n e dt op r e s e i l tt l l e 仃a c _ k i n gr e s u l t s f u r t h e 锄o r e ，t oa c h i e v er e a lt i m e 舶e r 仃a c b n g ，w eu s eg p u t os p e e du p6 b e r t r a c k i n gb ym e a l l so fc u d a an e wc o m p u t i n gm o d e li sp r o p o s e dt 0m a k e 如1 l u s eo f g p u sp a r r a l l e la b i l i 够o u r6 b e rt r a c k i n ga l g o r i m m si m p l e m e n t a t i o nh a v ea c h i e v e d a10 xp e r f - o r m a l l c eg a i n t h ep r o p o s e dm o d e lc a nb ee a s i l yu s e dt oa c c e l e r a t eo t h e r f i b e rt r a c k i n ga l g o r i t l h n s k e yw o r d s ：d t iv i s u a l i z a t i o n ，a n i s o 仃d p i cd i 髓s i o n ，丘b e rt r a c k j n g ，g p u a c c e l e r a t e dc o m p u t i n g 2 华东师范大学2 0 1 1 届硕士学位论文 ： 目录 摘要1 a b s t r a c t 2 目j 豪：3 第1 章绪论s 1 1 1 2 1 3 1 4 第2 章 2 1 2 2 第3 章 3 1 3 2 3 3 3 4 第4 章 4 1 4 2 4 3 4 4 4 5 4 6 第5 章 5 1 d t i ( d l f f u s i o nt e n s o ri m a g i n g ) 辛宇景5 d t l 数据可视化及意义5 本文的组织结构5 实验数据6 d t i 基本原理7 d t i 实验8 张量及模型：9 扩散张量及相关测度的可视化1 2 张量矩阵计算及显示1 2 各向异性测度计算及显示1 4 扩散张量重心空间的测度1 7 基于主特征向量颜色编码的椭球显示及交互1 8 神经纤维跟踪算法及显示2 l 相关工作2 1 多平面重建( m p r ) 2 1 f a c t 跟踪算法及结果2 3 t e n s o r l i n e 算法及结果2 7 b r u t ef o r c e 跟踪算法及结果3 0 程序结构3 3 c u d a 与加速算法实现。3 4 g p u 和c u d a 简介3 4 3 5 2 基于g p u 的纤维跟踪3 7 5 3数据结构3 9 5 4与c p u 计算结果比较4 1 第6 章总结及展望4 2 6 1工作总结4 2 6 2展望4 3 参考文献4 4 致谢4 7 4 第1 章绪论 1 1d ti ( dif f u sio nt e n s o rim a gin g ) 荦争景 m r i 技术的应用已经为医疗领域开启了新纪元。随着软硬件的飞速发展， m r i 技术也日趋成熟，各种新的成像技术也相继诞生，如扩散张量成像( d t i ) ， 血管成像( m r a ) 等。活体组织是不均匀介质，水的扩散运动会受到神经纤维 束的限制，d t i 正是以此作为物理基础，使用s t e j s k a l t a n n e r 脉冲序列对活体中 水分子的扩散运动进行编码，得到扩散加权磁共振图像( d w i ) 1 ，2 】。由于活 体组织中水分子的扩散运动是各向异性扩散，如果从多个方向获取d w i 图像， 即实验中在多个( 6 ) 不同方向施加扩散梯度，就可以通过计算得到每个像素 的扩散张量，这个过程叫d t - m r i 3 。作为非侵入式地探索活体组织细微结构的 重要工具，近十年来d t i 已经在临床诊断和科学认知领域得到重要的应用 4 。 1 2d tl 数据可视化及意义 d t i 数据中包含活体组织和病理组织的空间结构信息。但是从d t i 数据中很 难直接观察到这些信息，而是需要对d t i 数据进行处理，利用这些信息对组织 中神经纤维束重建并可视化，从而对神经纤维束的走向、连通和完整程度进行直 观的观察，并做病理分析。近几年d t i 数据重建与可视化已经成为研究的热点， 研究者和开发者提出多种可视化技术，这些技术可以大致分成颜色编码法，图元 显示法，多平面重建( m p r ) 和神经纤维跟踪技术 5 】，其中神经纤维的重建是 d t i 数据可视化技术中应用最多的。d t i 数据的重建与可视化不仅在神经科学和 一些疾病包括脑部肿瘤、多发性硬化症的研究和诊断中发挥了重要作用 6 ，而 且，对于活体组织细微结构和手术的规划等方面，也具有极大的潜质，因此，对 d t i 数据进行重建和可视化有重大的意义。 1 3 本文的组织结构 本文对d t i 可视化的各种方法做了深入的探索，分析了其优劣，并实现了 多种d t i 数据可视化算法，软件包含算法和显示两大部分，算法中由于数值计 算和矩阵计算的计算量较大，用c + + 实现，数据的3 d 显示用o p e n g l 完成，基 于g p u 的加速计算用c u d a 平台编程实现，软件有友好的交互和简洁现代的界 面。 第一章是绪论，简单介绍了d t i 的背景及可视化的意义； 第二章比较详细地介绍了d t i 的原理，介绍了各向同性扩散和各向异性扩 散，并详细讨论了扩散的张量表达形式和数学模型； 第三章实现了张量、扩散各向异性的几种测度以及扩散张量重心空间的测度 的计算流程与可视化的流程，并展示和分析了它们显示结果； 第四章实现了几种神经纤维束的跟踪算法，总结了其优劣，对不同算法的可 视化的结果进行了比较； 第五章在第四章的基础上，提出了一种适合并行计算的跟踪算法的实现模型， 并使用c u d a 进行加速计算，取得了一个1 0 倍的加速效果； 第六章对全文工作进行了总结，并对下一步的研究做了展望。 1 4 实验数据 为了避免计算结果受到噪声的影响，以及便于与其他d t i 工具包的计算结 果进行比较，本软件中使用的数据是一组d t i 模拟数据。数据分辨率是9 1 8 1 6 0 ，第三维表示图层数，每一层图对应6 0 个方向的加权梯度，= 7 0 0 如无特 殊说明，本文中的可视化结果都是利用这组数据进行计算的结果。数据下载地址： h t t p ：、) l ，w w d t i u a a c b e w i m o n 西n a l 2 0 a t l a s 2 0 a 1 1 d 2 0 9 r o u n d 2 0 仃1 l t h 6 第2 章d ti 基本原理 2 1分子扩散分子扩散即分子的布朗运动，是一种随机的、无规则的平移 运动，在不受限制的环境下，给定时间内分子扩散的距离满足高斯分布。想象一 滴黑墨水滴入纯净水中，f 时间后墨水会扩散开来，如图1 ( 出处不明) ，墨水分 子的扩散是无规则、完全随机的，最后扩散的形状会接近球状( 浅蓝色) 。 o 9、j 一、，；： ， ， ，| 一 “+ j 。? ：j。，叠 # ” ” ，一十 、。 一t“ ：。j ， ，。“j 。，。：i 。5 一o o 图1 分子的布朗运动，扩散形状是球状 f 时间后扩散的平均距离满足e i n s t e i n 方程 3 ： 仃：压瓦 ( 2 - 1 ) 其中仃是分子扩散的距离，与扩散时间f 的开方成正比，d 是扩散系数，用 来描述分子扩散的快慢。沿各个方向的扩散速率相同的扩散称为各向同性扩散。 但是，在复杂的生物组织中，存在大量的束状结构，如图2 7 中a 图是肌肉的微 观结构，肌肉由许多包含肌原纤维的纤维束组成，b 图是大脑中神经元的示意图。 在纤维束中水分子的扩散会受到细胞膜的限制，水分子沿着纤维束扩散比垂直纤 维束扩散要容易，这种情况下水分子的扩散不再是各向同性的，因此称为各向异 性扩散 8 】。此时的扩散结果也不再是球形，而是一个椭球 7 】，如图3 所示。 7 曼。 舀 ，0硌 k 0 。 一，、，一，9 口。：。，o 4 一誊 90，口位， ，留。气， 镑， 粕s c l e 移蛳b r n b 1 图2 ：a ) 肌肉组织结构：肌肉由大量神经纤维束组成；b ) 神经元的示意图，水分子在轴突中的扩 散是各向异性的 图3 图a ：水分子在神经纤维束中扩散的情形；图b ：水分子扩散的形状是一个椭球 2 2d tl 实验 使用s t e j s k a l t a n n e r 脉冲序列 1 ，9 】( 如图4 ) 对水分子的扩散运动效应进行 编码，能得到扩散加权磁共振图像( d w i ) 。s t e j s k a j t a i l e r 序列就是在常规自 选回波( s e ) 序列中加入一对大小和方向都相同的梯度脉冲，置于s e 序列的1 8 0 。 脉冲的两侧。第一个梯度脉冲引起质子自旋，使得磁化矢量相位发生变化，时 间后另一个梯度脉冲使磁化完全重聚，会导致信号衰减。信号( s ) 衰减的程度 可表示成 ；姜 坐丕嘘延太堂2 q ! ! 屋亟堂焦途塞 s = 驴x p 中双一害) g ，融】 ( 2 - 2 ) 其中而是非梯度加权图，2 万2 ( 一喜) h 2 是扩散加权因子，l e b i h a i l e ta l 将 其简写为b ，万是每个梯度脉冲的作用时间，是两个梯度脉冲的间隔时间，i 纠 是梯度脉冲的强度，d 是分子扩散系数，也叫做a d c ( a p p a r e n td i m s i o n c o e m c i e n t s ) ，其大小受梯度方向的影响 1 0 。 图4s t e j s l ( a l t a n n e r 脉冲序列 信号的衰减最终反映到d w i 图像中，体现为图像灰度的强弱，这就是d w i 的成 像原理。由前面的章节知道，各向异性扩散中，扩散系数d 是一个3 3 的对称 矩阵，有6 个独立的变量。要求出张量d ，至少需要6 个g 。也就是说，在实验 时，至少需要加6 个不同方向的梯度场，得到6 幅不同扩散方向的d w i 图像( s ) 和1 幅非梯度加权图像( 晶) ，才可以利用这些图像计算出每个提速的扩散张量， 这就是d t i 的基本原理。在实践中为了得到更精确的扩散张量数据，通常会使 用更多个不同方向的梯度场( 如1 2 ，2 5 ，6 0 个) ，得到更多d w i 图像3 1 。 2 3 张量及模型 描述各向异性的扩散需要使用椭球模型，由于椭球有六个自由度，对于各 向异性扩散的情况，扩散系数用一个标量d 来描述是不够的，而需要用一个二 阶张量来表示 3 ，9 ，1 1 ： 9 l 腑聊觑1 d = l 脚助助i j ( 2 3 ) d 勉d 弦d z z 张量模型是一个椭球 椭球方程 图s 张量模型，张量的特征值与椭球的主半径开根号成正比，特征 向量方向与主半径对应的方向一致 寿+ 专+ 专= 1( 2 4 ) 椭球的三个主轴长度分别与张量的三个特征值开根号即打、压、历成正比， a 、乃、忍的大小关系是a 乃 丑，三个主轴方向分别是张量的三个特征向量 写、乏、乏，其中最大特征值五对应的特征向量弓是主特征向量。椭球的形状 取决于扩散的情况，w e s t i ne t 龇1 2 将各向异性扩散分成下列三种基本情况，如 图6 7 ，其他的扩散情况介于这三者之间： 1 0 a )b )c ) 图6a ) 线性扩散，b ) 饼状扩散，c ) 球状扩散 1 ) 丑观五，椭球呈细长型，如图a ) ，扩散为各向异性扩散，也叫线性扩散， 白质中体现为这种扩散； 2 ) ，乃，椭球呈饼状，如图b ) ，扩散介于各向异性和各项同性之间，这 纠 5 。，椭球呈饼状，如图b ) ，扩散介于各向异性和各项同性之间，这 种扩散叫平面扩散，神经纤维交叉或分叉处表现为这种扩散； 3 ) 丑五乃，椭球呈球状，如图c ) ，扩散为各向同性扩散，也叫球状扩散， 脑灰质或脑脊液中体现为这种扩散。 第3 章扩散张量及相关测度的可视化 ： 3 1 张量矩阵计算及显示 前面提到过，扩散位移的布居数分布呈高斯型，而信号衰减满足： s = 而e x p ( - 6 d ) ( 3 1 ) 假设在d t i 实验中加6 个不同方向的梯度场，得到一组d w i 信号 s 。，毛，岛，s ，s 。，黾，s 。) ，其中是非加权图像，将公式( 3 1 ) 变换后得到 一组方程 1 】： 衄而) = 1 n ( ) _ 坞rd 9 1 警2 妖一6 9 2 h ( 3 - 2 ) 坂& ) = 妖而) 二喙r d 9 6 d 搿d 卿d 勉i 由6 个方程可以解出对称矩阵d = i 助功吵啦i 。如果实验条件允许，在 【脑助陇j d t i 实验中设置n ( n 6 ) 个不同方向的梯度场，d 的拟合公式变成 1 3 ： 一三1 1 1 ( 旦) 、s o 一三1 1 1 ( 兰) 6 0、7 一三1 i l ( 立) 6 0 、s o 7 式子( 3 3 ) 简化为 五2乃2毛22 五m2 y 1z 1 2 五z 1 x ? y ? z ? 2 x l y i2 y l z t2 x i z t 记乏2 x n y n2 y n z n2 x 。z ， 1 2 d x x d 列 d x z 功 d v z d z z ( 3 - 3 ) 垡丕! ! 亘蕉态堂2 q ! ! 昼亟堂焦途窒 a d c = g 车d ( 3 - 4 ) 为了保证矩阵可乘，需要在式子( 3 - 4 ) 等号左右同时乘以梯度矩阵g 的转置矩阵 g r ，即 再对式子( 3 5 ) 右边求逆 g ，宰么d c = g7 幸g 木d ( 3 - 5 ) ( g r 木回宰g 7 木么d c = ( g 7 木g ) 宰( g 7 车g ) 木d( 3 6 ) 便可以求出张量d 。 获得d t i 数据后，我们按照以上计算步骤计算出每个体素的张量d ，这里 使用的d t i 数据是由网站h t t p ：、 ，、 惭，搋u a a u c b e w i m 合成数据，数据的分辨率 是9 1 半8 1 枣6 0 ，第三维表示图层数目。图像中的每个体素对应一个3 掌3 的张量， 每个张量包含d 搬，珊，眈，助，助，陇六个独立分量。计算结果中， 每个独立分量都是一组9 1 木8 1 木6 0 大小的数据，下面分别对每个分量的第3 0 层图 像进行显示，如图7 所示，从图中我们能够看出，张量对角元素d 矗，d 。，见 的图像对比度更强，非对角元素的图像对比度较低。我们运用了c 十+ 语言对张量 d 进行计算和显示，并且，为了方便用户交互操作，每个分量的图像都可以通过 鼠标滚轮改变图层。 d捌dxz 聊d 弘 图7 张量d 的可视化结果，可以注意到张量的对角元素见，见的图像对比度更 强，非对角元素的图像对比度较低 另外，为了保证计算过程无误，在用c + + 实现算法前，首先用m a t l a b 编程 现实了以上计算过程，所有c + + 编程的计算结果都与m a t l a b 的计算结果进行了 比较，确保了结果的正确性。 3 2 各向异性测度计算及显示 扩散的各向异性 1 4 是d t i 的重要属性，对各向异性的测度的显示是可视化 的重要内容。d w i 图像中包含了分子各向异性扩散的信息，但这种信息并没有 直接反应在图像中，而是要对这些数据经过处理之后才显示出来。p j b a u s s e r 等 人提出了几个从数学的角度来描述扩散的各向异性的标量 1 5 ，这些标量的可视 化可以让我们直观地看到水分子扩散的各向异性，从而在不同程度地反映了生物 组织的细微结构，比分别可视化张量的6 个分量更有意义。这些标量的计算都基 于张量对应的特征值或特征向量。 。 1 4 f a ( 觚t i o n 甜a r i i s o 们p y ) 是最常用的各向异性的测度，取值范围是肛1 ，能 准确地描述扩散的各向异性程度，对于完全各向同性的组织剐= 0 ，而完全各向 异性的组织的剐= l 。剐也经常被用作纤维追踪时选择种子点的判据。计算公 杜 a ： 刚= 如型号特产 ( 3 - 7 ) 由于水分子在神经纤维中的扩散受限而常常表现出较大的各向异性，从f a 图像中能区分脑白质和灰质两种不同组织结构。 r a ( r e l a t i v ca i l i s o 仃o p y ) 刚= 压巫啄殍 l u 的表示组织中各向异性成分和各向同性成分的比值， 向异性扩散的分布情况，从而提供脑组织的微观结构信息。 ( 3 - 8 ) 能够较好地描述各 m d ( m e a i ld i h u s i o n ) 1 m d = ( 腑+ 助+ 陇) ( 3 9 ) j ， m d 表示平均扩散程度，值越大，表示组织内自由扩散的水分子越多，能很 好地区分大脑中的脑脊液成分。 在3 1 节中，我们已经得到了一组张量d 的矩阵，利用张量的特征值，能够 分别计算出以上各向异性测度。程序中，我们使用了rbd a v i e s 提供的n e w m a t l o ( n p ：r o b e r t n z n e 价吼1 0 h t n l ) 工具包来计算矩阵的特征值和特征向量，简化了 计算过程。图8 从左到右依次显示了以上各向异性测度的第3 0 层图像，跟张量 的显示一样，每个各向异性测度的图层可以通过滚动鼠标滚轮改变，跟其他通过 改变参数修改图层的软件相比，我们的软件操作起来更加便捷。考虑到篇幅，这 里只演示一层图像。 图8 从左到右依次为f a ，i 认，m d 第3 0 层的图像，从图中可以看出脑脊液部分的平均扩散系数较大， 但扩散的各向异性则较小。而扩散各向异性较大的部分基本上反映出了脑白质的结构 此外，为了更形象直观地描述扩散各向异性，我们还利用扩散的方向信息对 f a 图像进行颜色编码，得到二维的扩散各向异性彩色映射图【1 6 】，第2 5 层和第 3 0 层的图像： 图9f a 的彩色映射图，图左和图右分别是第2 5 层图像和第3 0 层图像 颜色由三个通道r g b 组成，分别对这三个通道进行颜色编码，主特征向量 不仅具有方向性，而且其方向能代表组织中扩散的方向，用它进行颜色编码能准 确地描绘组织中各向异性扩散的方向信息。具体的颜色编码方法：犬= 咖( 钆) ， g = 如心y ) ，b = 出( q ：) ，写是主特征向量，即r g b 对应每个体素的主特征向 1 6 华东师范大学2 0 1 1 届硕士学位论文 量巨的三个分量大小的绝对值。这样，红色表示左右方向，绿色表示上下方向， 蓝色表示前后方向。为了增强图像的对比度，每个通道的值需乘上对应图层的 f a 。图像中越明亮的区域，各向异性扩散程度越高。图层也能通过滚动鼠标滚 轮改变。 以上扩散各向异性测度的计算和显示用c + + 实现，另外用m a t l a b 编程验证 了计算过程的正确性。 3 3 扩散张量重心空间的测度 除了用朋、r 爿、m d 来描述脑组织中各向异性的分布， 另外三个测度 1 2 1 4 来描述张量的几何形状的情况： c ，：声鱼三l 。 2 + 如2 + 如2 旷耥 3 允 c = 1 。= = = = = 三兰= = = = o 5 2 + 如2 + 如2 w e s t i n 等提出用 ( 3 1 0 ) 当五五以时，椭球的形状呈细长，接近雪茄型，这种扩散叫线性( 1 i n e a r ) 扩散，用标量c ，来描述，q 值越大，线性扩散的程度越大，椭球越接近雪茄型； 当a 五忍时，主特征值约等于次特征值，椭球呈饼状，这种扩散叫平面( p l a i l a r ) 扩散，用标量0 来描述；当丑五乃时，椭球呈球状，为球状( s p h 谢c a l ) 扩 散，用g 表示。图l o 7 是由这三个标量组成的几何重心空间，三角形的三个顶 点分别是这三种情况的极端。三个标量的取值范围都是o 1 ，并且是归一化的， 满足c ，+ 0 + g = 1 。 图1 0 扩散张量重心空间，三角形的三个顶点对应线性、平面和球状扩散 运用上述公式，我们用c + + 编程对以上三个标量进行了计算和显示，帮助用 户做病理分析和诊断。图1 1 从左到右依次是q ，0 ，q 第3 0 层图像的显示结果， c ，0 都显示了较高的各向异性，能够清楚地反映线性扩散，g 反映的是各向同 性扩散的信息。 图1 1 q ，勺，g 第3 0 层图像，c ，0 都显示了较高的各向异性，能够清楚地反映线 性扩散，e 反映的是各向同性扩散的信息 3 4 基于主特征向量颜色编码的椭球显示及交互 图元显示法 1 4 以椭球作为基本模型，对每个体素对应的张量进行可视化， 方便用户查看脑局部的扩散信息。由于椭球是3 d 模型，我们进行可视化时使用 堡丕! ! 匝范太堂2 0 “届亟士堂僮论窒 了0 i p e n g l 软件包。o p e n g g l ( o p e n 研a p m c sl i b r a 巧) 是一个开放的3 d 图形软件 包，独立于操作系统，性能卓越，移植性好，是专业的图形处理、科学计算等高 端应用领域的标准图形库。关于o p e n g l 详细的介绍请参考网站 h 鲤；丛幽塑厶q p 鲤吐q 璐。o p e i l g l 中有超过7 0 0 个函数，这些函数可以用于指定2 d 或3 d 模型，对模型的颜色、材质、光照等进行设置，并可创建交互的3 d 应用 程序【9 】。 在我们的软件中，用椭球作为图元对每个体素对应的张量进行了可视化，图 1 2 是用图元法可视化张量的截图： 图1 2 张量的椭球可视化，左图是f a 彩色映射图，可拖动图中的自框到图层的任意空 间位置，右图是白框中体素对应的张量的椭球可视化结果，从图中能观察到脑局部的各向 异性扩散程度及神经纤维的大致走向 图中能够方便地查看脑组织局部各向异性的信息，椭球的形状和颜色形象地 描述了扩散的方向。左图是主特征向量的彩色映射图，右图是脑组织局部结构中 张量的椭球显示，每一个椭球对应一个体素的张量，椭球的半径分别是张量的三 个特征值，方向分别是椭球的三个特征向量，为了清楚地查看扩散的方向，椭球 的颜色采用左图中的颜色编码方法，用主特征向量来进行颜色编码，得到与左图 1 9 中体素颜色一致的结果。右图中椭球对应的体素就是左图中白色方框中的体素， 用户可以用鼠标拖动白色方框至图左中的任意区域，观察脑组织中局部各向异性 扩散信息，也可以将鼠标停放在左图用鼠标滚轮改变图层。 图1 2 中可以看出，我们为软件设计了简洁时尚的界面以及便捷的交互方式， 界面采用v s 2 0 l o 中m f c 提供的基于硒b b o n 的设计，交互过程用m f c 结合 o p e n g l 实现。 华东师范大学2 0 1 1 届硕士学位途文 一 第4 章神经纤维跟踪算法及显示 4 1相关工作 神经纤维束跟踪 1 7 1 8 1 9 是d t i 可视化的重要技术，利用像素之间的关 联在3 d 空间重建组织中的神经纤维束。这种直观的可视化技术已经在临床和神 经系统科学领域得到广泛的应用。神经纤维束跟踪中要面临的最大的挑战是神经 纤维结构的复杂性以及原始d t i 数据的低分辨率，为此，人们提出了多种不同 的扩散模型和跟踪算法。目前应用最多的跟踪算法是f a c t ( f i b e ra s s i 孕l i n e n tb y c o n t i n u o u st r a c b n g ) 算法 1 7 ，利用主特征向量进行跟踪，对于各向异性程度高 的组织能够获得准确稳定的跟踪结果；与f a c t 算法不同，t e n d ( t e n s o r d e f l e c t i o n ) 算法 2 0 】使用整个张量对向量进行偏转，重建出的纤维更长更光滑， 但对各向异性扩散程度高的组织使用偏转不合理；t e n s o r l i n e 算法 2 1 将f a c t 和t e n d 算法结合，各向异性扩散程度高的退化成f a c t 算法，各向异性程度 低的区域退化成t e n d 算法；交叉和分叉纤维是现有跟踪算法的一个瓶颈，b r u t e f o r c e 算法 1 9 基于f a c t 算法，以全脑体素作为种子点进行跟踪，能够较好地 跟踪出交叉纤维，但缺点是计算量大；随着d t i 数据的分辨率增高，多种基于 h a r d i ( h i 曲a n g u l a rr e s 0 1 u t i o nd i 自r l s i o ni l l l a g i n g ) 的跟踪算法 2 2 - 2 5 被提出来解 决分叉纤维的难题。 为了满足用户的不同需求，我们实现了多种跟踪算法供用户选择，包括f a c t ， t e n s o r l i n e ，b m t ef o r c e 几种准确性和稳定性较高的算法。其中，我们将1 h s o r l i n e 和b m t ef o r c e 算法结合，获得了令人满意的结果。这部分工作包括算法的实现 和图形显示及交互三部分，算法部分用c + + 编程实现，图形显示和交互用o p e n g l 结合m f c 实现。用户可以手动导入d i c o m 格式的d w i 数据，也可以直接导入 张量数据，得到m p r3 d 图像，用户可以手动定义感兴趣区域( r o i ) ，通过菜 单选择不同的跟踪算法，能够统计每种跟踪算法的计算时间，并且对跟踪结果进 行缩放，旋转等操作。下面将具体介绍这几种算法的原理、实现过程及可视化结 果。 4 2 多平面重建( m p r ) 2 1 为了方便用户选取感兴趣区域( r o i ) 及查看脑局都神经纤维束，我们实现了 多平面重建( m p r ) ，图1 3 演示了重建好的3 dm p r 。每一维的图层显示我们运 用o p e n g l 中的3 d 纹理贴图实现，白色的边框采用了o p e i l g l 中的去锯齿处理， r o i 的选取运用o p e n g l 的拾取原理实现。 图1 4 和1 5 显示的m p r 就是将冠状，矢状，横断面对应的特征向量映射图 正交组合的结果，方便用户手动选择感兴趣区域 2 6 】。用户可以通过鼠标进行缩 放、旋转操作，每个方向的图层都可以通过鼠标平移位置，不同空间位置的图层 包含的信息与大脑空间位置对应的特征向量映射图一致。图1 4 是用鼠标将图层 拖动平移后的结果，可以看到，图层包含的信息与大脑空间位置对应的扩散信息 。是一致的。用户能够通过拖动图层，改变视角等操作来定义r o i 。 m p r 中3 d 纹理的实现流程和主要函数： 启用纹理贴图功能 g l e n a b l e ( g l - t e t u r e 一3 d ) 绘制场景、提供纹理和几何坐标 g l d i s p i a y c o i o r m a p s 3 d ( ) 图1 3m p r 中3 d 纹理的实现流程和主要函数 一蔓互 一一一j 一一 耄帮 1 9 9 9 年m o r is 等人提出的f a c t ( l i n ep r o p a g a t i o na p p r o 犹h ) 是神经纤维束 跟踪算法中最简单直观的算法。该算法以种子点作为起始点，根据像素的主特征 向量向前追踪，如图1 6 。具体的跟踪步骤如下： 1 ) 图中的体素坐标系是连续的坐标系，选取中央像素点作为种子点，即跟踪的 起始点坐标为( 1 5 ，1 s ) ； 2 ) 从种子点出发，同时向两个相反的方向跟踪，得到与体素边界的两个交点分 别为( 1 8 ，2 o ) 和( 1 3 ，1 o ) ； 3 ) 分别从交点出发，根据邻近体素的主特征向量向前跟踪。 ( 1 8 ，2 0 ( 1 5 。1 5 ) ( 1 a1 o ) - 图1 6m 跟踪算法。以种子点的中央( 坐标是( 1 s ，1 5 ) ) 作为出发点，同时向两个相反的方 向跟踪，与相邻体素相交，交点坐标分别是( 1 3 ，1 0 ) 和( 1 8 ，2 o ) 跟踪过程停止的条件首先是不能超过数据边界。此外，图1 7 3 】a 图中，各 向异性测度( f a ) 大的像素椭球呈细长，趋近线性，很容易根据主特征向量确 定跟踪方向，但图中带星号的像素各向异性测度低，主特征向量和次特征向量大 小相近，很难根据主特征向量确定跟踪方向，此像素点跟踪方向最终会被噪声决 定，因而应停止跟踪。另外，两个相邻像素的主特征向量如果角度偏差太大 ( 4 5 。) ，同样的像素却会导致差异很大的结果，从而作为跟踪过程停止的第二 个判据。 o痧 po 。o0o b ，7 11_7 、 11 图1 7 左：每个体素用一个椭球来描述，不同的椭球形状表示各体素扩散的各向异性程度不 同；右：图中的箭头是每个体素的主特征向量，图中的星号表示该体素受到噪声影响 f a c t 跟踪算法的实现流程图： 图1 8f a c t 跟踪算法的实现流程图 图1 9 是我们用f a c t 算法进行跟踪的结果，用户可在m p r 中定义好l 的i ， 在上面的算法选框中选择f a c t 算法，点击“d r a wf i b e r 按钮，便可得到自己 关心区域的纤维束。图2 0 是对图1 9 中一束纤维放大后的结果。图1 9 中，种予 一一一一一7 7 d a t o o o 、 b 图2 1 同样的三个体素，由于各项异性程度不一致，用f a c t 算法跟踪得到不一样的结果， 各向异性程度越低，跟踪方向越不稳定 4 4t e n s o ri n e 算法及结果 图2 1 3 】中，从下往上看椭球的形状分别为线状、饼状、球状，椭球中的箭 头表示该像素点产生的跟踪方向，可以看出，细长型的椭球的方向由主特征向量 方向决定，但是饼状和球状椭球的方向比较随机，主特征向量和次特征向量相近， 导致跟踪过程极不稳定，三个同样的像素会产生完全不一样的跟踪结果。因而， 向异性测度低的像素的跟踪方向不能只依赖主特征向量，w e i n s t e i n 等人提出便 于处理饼状或球状各向异性的t e n s o r l i n e 算法。该算法中，每一步的跟踪方向由 公式( 4 1 ) 决定： 圪卸= c ，k + ( 1 一c ，) ( ( 1 一w ) + w 圪州) ( 4 - 1 ) 其中 k = d 圪 圪为前一步的跟踪方向，吃叫是当前点对应的张量d 对元偏转的结果 2 0 】。 ( 4 2 ) 由于d t i 数据天生分辨率很低，每一步的d 都需要用当前体素的八邻域张量 2 7 进行三线性插值 2 7 ，并计算出主特征向量露，并通过特征值计算当前点的线性 扩散测度q 。第三章中提到过，q 用于衡量当前点的扩散形状，取值在0 到1 之 间。该算法能够因地制宜，各向异性扩散程度不同，影响