基于扩散加权磁共振数据的脑纤维三维显示方法

摘要

本发明为一种基于扩散加权磁共振数据的脑纤维三维显示方法，可以快速精确重构脑纤维束，属于信息技术领域中的医学图像处理领域。提出了一种新的离散非负约束结合高斯函数拟合的球面反卷积算法来建立扩散加权磁共振数据的球面反卷积模型，同时我们提出了一种以纤维束为基本跟踪单位的群体纤维跟踪算法，应用该算法跟踪球面反卷积模型可以快速重构脑纤维。通过应用VTK可视化工具包将群体纤维跟踪算法生成的数据可视化，从而将脑纤维精确的显示出来。本方法简化了计算，提升了跟踪速度，同时有良好的抗干扰性能。

说明书

技术领域

本发明涉及信息技术领域下的医学图像处理领域。 

背景技术

利用扩散加权磁共振成像(DW-MRI)数据重构脑白质神经纤维是目前活体显示脑白质纤维走向的重要手段，它为研究脑白质的空间微结构开辟了崭新的途径，为脑的发育、精神分裂症、先天性与获得性脑白质病以及痴呆等的研究提供了新的方法。其基本思想是首先估计扩散张量磁共振图像每个体素的扩散张量，以此获得纤维方向分布函数,然后利用纤维跟踪算法形成具有解剖学意义的纤维空间微结构。由于该方法的直观性，并易于实现不同对象的纤维结构对比分析，从而引起了全世界研究者的关注。目前应用较广的是基于扩散张量模型（Diffusion Tensor Imaging,DTI）的Streamline纤维跟踪显示方法，但目前还存在一下问题： 

首先，现有扩散张量估计方法难以刻画复杂的脑纤维微结构。人脑结构的复杂性及纤维连接模式的多样性是目前面临困难的主要原因。目前常用的张量估计是利用扩散张量成像模型建立体素纤维扩散张量，并利用最小二乘法计算获得最优二阶张量矩阵。但是DTI成像模型的一个本质局限是其假设每个体素仅含一根神经纤维。实际上，DW-MRI数据的分辨率通常在3～15mm³之间，纤维跟踪算法中的纤维束直径为毫米级，而实际神经纤维直径为微米级。因此，应用DTI模型难以刻画诸如交叉、分叉、瓶颈型等复杂的纤维结构。 

其次，目前应用较广的纤维重构算法诸如Streamline算法、贝叶斯概率跟踪算法，其跟踪过程都是基于单根纤维进行的，而实际上一组普通DTI数据的脑纤维数量就达到1000-10000根，因此要跟踪出全体纤维计算耗时长，计算效率低下。最近，基于全局优化技术的纤维束跟踪算法得到极大的关注，成为快速重构神经纤维的新趋势。 

而在将脑白质纤维三维可视化方面，目前普遍采用的方法为应用Direct X显示控件，将脑纤维重构数据建模，通过将此模型导入虚拟三维空间来将其可视化。但Direct X对编程有较高要求，需要专业人员来建立纤维模型，不适合普通医务人员使用。 

发明内容

为了克服已有脑白质三维显示方法存在的跟踪精度低、跟踪速度慢、操作困难的问题，本发明提出了一种分辨率高、跟踪速度快、易于医务人员操作的脑纤维三维显示方法。 

本发明采用的技术方案为： 

（1）采用离散纤维方向密度函数建立球面反卷积（spherical deconvolution method,SD)模型： 

本文提出一种新的离散非负约束结合高斯函数拟合的球面反卷积算法：该算法首先采 用离散纤维方向点建立卷积模型，以此来隔离纤维方向密度函数对反卷积算法的影响，从而获得高角度分辨率识别；然后借助球面高斯函数拟合来补偿离散误差，从而消除由离散脉冲函数带来的离散误差，更精确地求出纤维方向密度函数。 

（2）进行基于球面反卷积模型的群体纤维重构： 

纤维跟踪可看作寻求在种子点与目标区域之间的一条最大后验概率路径的问题；在经典的跟踪算法如Streamline算法与贝叶斯概率跟踪算法中，纤维跟踪是基于单根纤维进行的，我们提出一种并行的基于纤维束进行的群体纤维跟踪算法；在我们的跟踪算法中，设定适当的粒子数目，粒子按照纤维概率密度方向随机得到纤维；将纤维分为不同的区段，通过应用局部函数模型为各区段内的各纤维的赋以不同的权重；选出各区段内权重最大的纤维路径为局部路径；遍历所有纤维区段，得到各段的局部路径；依次跟踪各区域内的局部路径得到纤维束的全局路径。 

我们定义给定的纤维束f的全局价值函数G为： 

$G\left(f(x_0\to x_n)\right)={\int }_{x_0}^{x_n}L(f\left(t\right),n,\alpha )\mathrm{dt}---\left(1\right)$

式中，L(f(t),n,α)为纤维束局部价值函数，f(t)为局部权重最大的某根纤维f的局部价值函数，n为纤维束中与纤维f的向量夹角小于参数α的纤维数；实际上，纤维束的路径可近似为纤维束中某根权重最大的纤维的路径，这样可极大地提升运算效率； 

一条脑白质纤维路径可看作是图像空间内的一些离散点的路径，即f={x₀,x₁,...,x_n}，假定所有向量的步长相同，即χ_i=χ,i=1,...,n；在离散时间内的路径可写成一个递归式： 

xⁱ(t+1)=xⁱ(t)+χvⁱ(t),i=1,...,n-1   （2） 

其中，vⁱ(t)是位置xⁱ上的纤维方向,假设A为起始区域，而B为目标区域，则从A到B的一条纤维路径的总值S可近似为： 

$S\left(f(x_0\to x_n)\right)=\chi \underset{i-1}{\overset{n}{\Sigma }}f\left(x_i\right)---\left(3\right)$

我们将纤维束中所有纤维方向都近似最优纤维方向，则可快速得出纤维束的走向。 

（3）数据的三维可视化显示，其主要技术特征为： 

我们通过调用VTK医学图像数据显示工具来可视化纤维的重构数据。VTK具有强大的三维体绘制、面绘制功能。采用模块化的面向对象的技术，将基本功能用类包装，系统结构简单，用户只需要调用对应功能的类，而不需从底层写算法就可以实现，从而使得用户可以很容易在很短时间内学习和掌握VTK。同时，较Direct X显示库，其曲面显示能力更强大，更适合纤维等结构的显示。 

具体纤维显示的实现过程如下： 

1．求取球面反卷积模型过程如下： 

1.1求取卷积模型。将观测信号表述为单条纤维信号响应函数R和纤维方向密度函P(v)在面上的卷积： 

$S=S\left(g\right)/S_0=R\otimes P\left(v\right)---\left(4\right)$

其中g为单位扩散脉冲梯度向量，S(g)为脉冲梯度方向g处的测量信号，S₀为无脉冲梯 度时的测量信号，v为单位方向向量，P(v)为球面上一系列基函数的加权组合。 

将S离散化即将单位球面离散为n个单位方向v_j(j=1,…,n)，该方向的密度函数值为P(v_j)；信号响应函数采用如下模型： 

$R(g,v_j)=e^{-l{\left(g^T{v}_j\right)}^2}---\left(5\right)$

其中v_j为球面上第j个采样方向向量，l是表征纤维各向异性程度和b值（扩散敏感系数）对信号衰减共同影响的参数。由(4)(5)式得离散化的卷积式为： 

$S\left(g\right)/S_0=\underset{j=1}{\overset{n}{\Sigma }}{e}^{-l{\left(g^T{v}_j\right)}^2}P\left(v_j\right)---\left(6\right)$

1.2基于公式（6）求解纤维方向密度函数。给定m个观测信号S(g_i)，i=1,…,m,定义能量函数： 

$E_n=\underset{i=1}{\overset{m}{\Sigma }}{(S\left(g_i\right)/S_0-\underset{j=i}{\overset{n}{\Sigma }}{e}^{-l{\left(g_i^T{v}_j\right)}^2}P\left(v_j\right))}^2---\left(7\right)$

其中，g_i表示第i个扩散梯度脉冲向量。因为增加非负性约束将显著提高反卷积的抗噪能力，本方法用非负最小二乘法求解(7)式。为了还原精度更高的连续的纤维方向密度函数，选取球面高斯函数： 

$f(p,\sigma )={\omega }_j{e}^{(-{\cos }^{-1}\left(\right|(p,v)\left|\right)/{\sigma }^2)}---\left(8\right)$

其中p_j为高斯函数的均值，σ为高斯函数的方差，v为三维单位向量，定义最小二乘曲面拟合函数： 

$E=\underset{i=1}{\overset{m}{\Sigma }}{({\lambda }_i-\underset{j=1}{\overset{n}{\Sigma }}{\omega }_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\sigma }^2)})}^2---\left(9\right)$

则问题转换为如下的方程组求解问题： 

Aω=λ     （10） 

其中ω为包含ω_j的n维未知数列向量，λ为包含λ_i的m维列向量，A为m×n的矩阵； 

$A_{i,j}=e^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\sigma }^2)}---\left(11\right)$

用非负最小二乘法求解(11)式得到ω，通过公式（12）我们可得到纤维的方向密度函（FOD）。 

$\mathrm{FOD}\left(v\right)=\underset{i=1}{\overset{n}{\Sigma }}{\omega }_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v)\left|\right)/{\sigma }^2)}---\left(12\right)$

2获得全局纤维束的后验概率路径的过程如下： 

2.1在第一步设定m个粒子：xⁱ(t),i=1,...,m； 

2.2利用最小二乘法计算扩散张量Dⁱ(t)； 

2.3将得到的纤维分为l个区间，在每个区间内，研究各条纤维的局部价值函数： 

假设我们研究第i个区间，定义如下纤维权重计算模型： 

W_i=W(f_j,θ,s)     （13） 

来统计对于某根纤维f_j，在其领域半径为s的柱体空间内，向量方向与f_j纤维向量方向的偏移角小于θ的纤维数目n_i，并将其作为纤维f_j的权重因子W_j； 

2.4重复步骤2.3直至区间i中的纤维全部被遍历，选出权重系数超过设定参数λ的d根纤维f_i,i=1,...,d并记录相应的系数n_i,i=1,...,d将其局部价值函数f(t)作为第i个区间内n_i根纤维组成的纤维束的局部价值函数L(f(t))_i； 

2.5重复步骤2.3、2.4直至l个区间的局部价值函数L(f(t))_i都被求出，我们规定这n_i个粒子沿此轨迹并列前进； 

2.6计算纤维路径的适当值。选择具有最合适值的β个粒子存为 并根据公式（14），分别计算出β个粒子在l个区段的纤维路径。 

xⁱ(t+1)=xⁱ(t)+χvⁱ(t),i=1,…,β   （14） 

3将重构数据数据三维显示化过程： 

3.1显示数据的获取。由群体纤维重构算法产生的纤维路径S是通过公式（15），由一系列的离散的数据点来拟合的，首先我们需要解析纤维束数据。为此建立一个在VTK库内用于存放数据的vtkData类的实例Data来保存由群体纤维重构算法生成的纤维数据，具体操作方法为：将纤维的各个离散数据点的坐标加一个偏移后存入Data类中，，并将其存储在内存中；我们可以通过修改Data中的数据来修改要显示的纤维； 

$S\left(f(x_0\to x_n)\right)=\chi \underset{i-1}{\overset{n}{\Sigma }}f\left(x_i\right)---\left(15\right)$

3.2建立三维空间显示实验数据。要将Data中的数据三维显示出来，我们需要为其建立一个虚拟空间，其中包括背景数据（VtkWindow）与作用物（VtkActor）两部分；我们建立一个名为Window的VtkWindow类实例来生成一个背景空间，通过设定参数可以调节背景颜色、是否有坐标显示、是否可旋转等内容；同时，我们建立一个名为Actor的VtkActor类，并将Data中的数据传递给Actor，通过配置Actor的参数，我们可以设定显示物是否被拉伸、显示物的颜色及位置的偏移量；然后调用VTK的显示控件来加载Window实例，为Window建立一个虚拟空间，同时将Actor中的数据信息显示到各数据点对应的坐标；通过此种方法，我们可以快速、精确地将脑纤维三维显示出来。 

本发明提出的基于球面卷积模型的群体纤维跟踪算法较传统的基于Q-ball成像模型的Streamline算法有显著优势。 

我们采用含有两束交叉纤维的Hough数据模型来验证本方法的优越性。其具体参数如下：体素量50×50×1，体素大小1×1×1mm，扩散系数b=1500s/mm²，76个扩散梯度方向。分别用SD与Q-ball两种方法在相同条件下进行数据模拟，如附图1中所示。 

我们选取了模型中的一块交叉区域作为研究样本，(a)为SD成像方法的区域模型图，(b)为Q-ball方法的区域模型图。可见，SD模型在非交叉处，保持单个特征方向，与Q-ball模型估计结果近似，而在纤维交叉处该模型走向清晰，较饼状的Q-ball模型更能准确地反映出纤维间的分布信息，符合实际的纤维情况。 

我们在起始位置m处设立50个种子点，产生满足正态分布μ=0，r=0、μ=0，r=5%、μ=0，r=10%的零均值高斯噪声（其中，r表示噪声水平）各50组，并将其施加于三组DW-MRI数据上，分别用基于DTI的Streamline算法、基于Q-ball的Streamline算法、基于SD群体跟踪算法对450组DW-MRI数据实施纤维跟踪，记录跟踪路径与到达目标位置n处的种子点数目，跟踪结果如附图2所示，并记录到达目标区域的种子点数目，计算其比例如表1所示。 

表1纤维跟踪结果比较表 

我们从实验数据可以看出，基于DTI模型的Streamline算法几乎不能如期到达目标区域，因此该模型不能解决复杂纤维结构；Q-ball模型的Streamline跟踪效果稍有改善，然而远不及SD模型下的群体跟踪算法。无噪声时，我们的方法纤维分布密集，种子点几乎都能到达目标区域，当噪声为5%时，种子区的纤维仍能如期到达目标区，而Q-ball模型则有部分纤维跑偏，当噪音干扰增大至10%左右时，由Q-ball模型得到的跟踪纤维效果不理想，而SD模型所得纤维受噪声干扰的影响不大，仍然保持较高的准确度。 

由此可见，我们提出的基于SD模型的群体纤维跟踪算法跟踪精确，同时由良好的抗干扰效果。 

附图说明

附图1是本发明的Hough数据的SD模型，（a）(b)分别为SD与Q-ball区域模型图 

附图2中的（a）（b）分别为基于DTI与Q-ball模型的Streamline纤维跟踪效果 

（c）为基于SD模型的群体纤维跟踪效果 

附图3是基于SD模型的群体纤维跟踪可视化方法效果图 

具体实施方式

参照附图1-3： 

（1）求取卷积模型。将观测信号表述为单条纤维信号响应函数R和纤维方向密度函P(v)在面上的卷积： 

$S=S\left(g\right)/S_0=R\otimes P\left(v\right)---\left(1\right)$

其中g为单位扩散脉冲梯度向量，S(g)为脉冲梯度方向g处的测量信号，S₀为无脉冲梯度时的测量信号，v为单位方向向量，P(v)为球面上一系列基函数的加权组合。 

将S离散化即将单位球面离散为n个单位方向v_j(j=1,…,n)，该方向的密度函数值P(v_j)；信号响应函数采用如下模型： 

$R(g,v_j)=e^{-l{\left(g^T{v}_j\right)}^2}---\left(2\right)$

其中v_j为球面上第j个采样方向向量，l是表征纤维各向异性程度和b值（扩散敏感系数）对信号衰减共同影响的参数。由(1)(2)式得离散化的卷积式为： 

$S\left(g\right)/S_0=\underset{j=1}{\overset{n}{\Sigma }}{e}^{-l{\left(g^T{v}_j\right)}^2}P\left(v_j\right)---\left(3\right)$

（2）基于公式（3）求解纤维方向密度函数。给定450个观测信号S(g_i)，i=1,…,450,定义能量函数： 

$E_n=\underset{i=1}{\overset{m}{\Sigma }}{(S\left(g_i\right)/S_0-\underset{j=i}{\overset{n}{\Sigma }}{e}^{-l{\left(g_i^T{v}_j\right)}^2}P\left(v_j\right))}^2---\left(4\right)$

其中，g_i表示第i个扩散梯度脉冲向量。因为增加非负性约束将显著提高反卷积的抗噪能力，本方法用非负最小二乘法求解(4)式。为了还原精度更高的连续的纤维方向密度函数，选取球面高斯函数： 

$f(p,\sigma )={\omega }_j{e}^{(-{\cos }^{-1}\left(\right|(p,v)\left|\right)/{\sigma }^2)}---\left(5\right)$

其中p_j为高斯函数的均值，σ为高斯函数的方差，v为三维单位向量，定义最小二乘曲面拟合函数： 

$E=\underset{i=1}{\overset{m}{\Sigma }}{({\lambda }_i-\underset{j=1}{\overset{n}{\Sigma }}{\omega }_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\sigma }^2)})}^2---\left(6\right)$

则问题转换为如下的方程组求解问题： 

Aω=λ     （7） 

其中ω为包含ω_j的n维未知数列向量，λ为包含λ_i的m维列向量，A为m×n的矩阵； 

$A_{i,j}=e^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\sigma }^2)}---\left(8\right)$

用非负最小二乘法求解(8)式得到ω，通过公式（9）我们可得到纤维的方向密度函（FOD）。 

$\mathrm{FOD}\left(v\right)=\underset{i=1}{\overset{n}{\Sigma }}{\omega }_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v)\left|\right)/{\sigma }^2)}---\left(9\right)$

其效果如附图1所示，我们建立了Hough数据的SD模型。 

（3）设定50个粒子：xⁱ(t),i=1,...,50； 

（4）利用最小二乘法计算扩散张量Dⁱ(t)； 

（5）将得到的纤维分为l个区间，在每个区间内，研究各条纤维的局部价值函数： 

假设我们研究第i个区间，定义如下纤维权重计算模型： 

W_i=W(f_j,θ,s)     （10） 

来统计对于某根纤维f_j，在其领域半径为s的柱体空间内，向量方向与f_j纤维向量方向的偏移角小于θ的纤维数目n_i，并将其作为纤维f_j的权重因子W_j； 

（6）重复步骤(5)直至区间i中的纤维全部被遍历，选出权重系数超过设定参数λ的d根纤维f_i,i=1,...,d并记录相应的系数n_i,i=1,...,d将其局部价值函数f(t)作为第i个区间内n_i根纤维组成的纤维束的局部价值函数L(f(t))_i； 

（7）重复步骤(5)、（6）直至l个区间的局部价值函数L(f(t))_i都被求出，我们规定这n_i个粒子沿此轨迹并列前进； 

（8）计算纤维路径的适当值。选择具有最合适值的β个粒子存为 并根据公式（11），分别计算出β个粒子在l个区段的纤维路径。 

xⁱ(t+1)=xⁱ(t)+χvⁱ(t),i=1,…,β    （11） 

效果如附图2（c）所示，设定的50个种子点并行排列跟踪，有良好的抗干扰能力。 

（9）显示数据的获取。由群体纤维重构算法产生的纤维路径S是通过公式（11），由一系列的离散的数据点来拟合的，首先我们需要解析纤维束数据。为此建立一个在VTK库内用于存放数据的vtkData类的实例Data来保存由群体纤维重构算法生成的纤维数据，具体操作方法为：将纤维的各个离散数据点的坐标加一个偏移后存入Data类中，，并将其存储在内存中；我们可以通过修改Data中的数据来修改要显示的纤维； 

$S\left(f(x_0\to x_n)\right)=\chi \underset{i-1}{\overset{n}{\Sigma }}f\left(x_i\right)---\left(12\right)$

（10）建立三维空间显示实验数据。要将Data中的数据三维显示出来，我们需要为其建立一个虚拟空间，其中包括背景数据（VtkWindow）与作用物（VtkActor）两部分；我们建立一个名为Window的VtkWindow类实例来生成一个背景空间，通过设定参数可以调节背景颜色、是否有坐标显示、是否可旋转等内容；同时，我们建立一个名为Actor的VtkActor类，并将Data中的数据传递给Actor，通过配置Actor的参数，我们可以设定显示物是否被拉伸、显示物的颜色及位置的偏移量；然后调用VTK的显示控件来加载Window实例，为Window建立一个虚拟空间，同时将Actor中的数据信息显示到各数据点对应的坐标；通过此种方法，我们可以快速、精确地将脑纤维三维出来，效果如图3所示。