CN103445780A - 一种扩散加权磁共振成像多纤维重建方法 - Google Patents

Abstract

一种扩散加权磁共振成像多纤维重建方法，用于重建多纤维方向的。方法为计算扩散加权磁共振信号衰减函数的逆球面氡变换，并给出了一种保证扩散加权磁共振信号衰减函数的逆球面氡变换非负性的算法。附图(1)显示了在一定的扫描参数，纤维特性和计算参数下本算法重建出的多纤维结果。

Description

一种扩散加权磁共振成像多纤维重建方法

（一）技术领域

本发明涉及一种扩散加权磁共振成像多纤维重建方法

（二）背景技术

利用扩散加权磁共振成像（DWI,diffusion weighted magneticresonance imaging）数据重构脑白质神经纤维是目前活体显示脑白质神经纤维微结构的唯一方法，被越来越多的应用于脑认知功能的研究、脑外科手术的导航、精神类疾病的诊断等。

扩散加权磁共振成像能够把水分子的扩散运动反映到磁共振信号中来，施加扩散梯度脉冲时测得的信号和没有施加扩散梯度脉冲时测得的信号之间的关系为

S(g)=S₀e^-bA            （1）

其中g为单位三维向量，表示扩散梯度脉冲的方向；S(g)为施加扩散梯度脉冲时测得的信号，扩散梯度脉冲的方向为g代表的方向；S₀为不施加扩散梯度脉冲时测得的信号；b为扩散加权系数。

扩散张量成像（DTI,diffusion tensor imaging）由于其计算效率高，所需的扩散梯度脉冲方向少（最少只需要6个）是目前最常用的方法，但是其局限是不能解决体素内存在多于一条纤维的情况。而人脑的神经纤维往往存在交叉、分支或融合的复杂情况，这时DTI将还原不出正确的纤维方向。

为了解决DTI的缺点，多种多纤维重建方法被提出，这些多纤维重建方法可以分为模型相关和模型无关两种类型，模型相关的方法都对所求解的结果限制为某种固定的模型，模型无关的方法则不对所求解的结果进行任何限制。其中模型无关方法DSI计算水分子位移概率密度函数的径向积分获得多纤维方向，缺点是信号采集时间长；QBI方法通过计算扩散加权磁共振信号衰减函数的球面氡变换计算零阶贝塞尔曲线加权后的水分子位移概率密度函数的径向积分，解决了DSI信号采集时间长的缺点，但是角度分辨率低；EAP方法计算水分子一定位移处的概率密度函数，虽然角度分辨率可以较高，但是水分子位移量的选择困难。

（三）发明内容

本发明针对现有模型无关多纤维重建方法的缺点而提出一种多纤维重建方法。

一种扩散加权磁共振成像多纤维重建方法包括如下步骤：

1）获取数据，包括扩散加权系数固定的扩散加权磁共振成像数据和不施加扩散梯度脉冲时的磁共振成像数据；

2）计算扩散加权磁共振信号衰减函数的逆球面氡变换；

扩散加权磁共振信号衰减函数是指施加扩散梯度脉冲后的磁共振信号和不施加扩散梯度脉冲时的磁共振信号的比以单位扩散梯度脉冲方向g为自变量的函数，即

$h\left(g\right)=\frac{S\left(g\right)}{S_0}---\left(2\right)$

球面氡变换为

$\mathrm{\ξ}\left[L\left(x\right)\right]\left(y\right)={\∫}_{x\∈y^{\⊥}}L\left(x\right)\mathrm{dx}---\left(3\right)$

其中x、y为三维单位方向向量，L(x)为定义在球面上的函数。即L(x)球面氡变换后某一个方向上的函数值为L(x)在与该方向垂直的所有方向上值的和。

步骤2）中扩散加权磁共振信号衰减函数的逆球面氡变换的求解方法如下：

记扩散加权磁共振信号衰减函数的逆球面氡变换的结果记为f(x)，定义能量函数

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-{\∫}_{{x\∈g}_i^{\⊥}}f\left(x\right)\mathrm{dx})}^2---\left(4\right)$

其中g_i为单位三维向量，表示扩散梯度脉冲的方向；S(g_i)为扩散加权磁共振成像数据，施加的扩散梯度脉冲方向为g_i；S₀为不施加扩散梯度脉冲时的数据；m为扩散梯度脉冲方向的数目。

在单位半球面上均匀采样n个点，记为{c₁,c₂,…c_n}，用(x^Tic_j)^l,j=1…n的非负加权和估计f(x)以保证f(x)的非负性，即

$f\left(x\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{(x^T\·c_j)}^{l}s.t.{\mathrm{\ω}}_j\≥0---\left(5\right)$

其中l为偶数，ω_j(j=1…n)为权重。

将（5）式是带入（4）式

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l\mathrm{dx}\right))}^{2}s.t.{\mathrm{\ω}}_j\≥0---\left(6\right)$

整理得如下的等价形式

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l\mathrm{dx}\right))}^{2}s.t.{\mathrm{\ω}}_j\≥0---\left(7\right)$

其中的球面氡变换用均匀采样的函数值的和代替积分来近似

${\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l=\underset{v=1}{\overset{z}{\mathrm{\Σ}}}{(x_v^T\·c_j)}^l,x_v\∈g_i^{\⊥}---\left(8\right)$

其中z为采样的点数。

将（7）式是写成矩阵的形式为

$\mathrm{En}={\left|\right|\mathrm{A\ω}-y\left|\right|}_2^{2}s.t.\mathrm{\ω}\≥0---\left(9\right)$

其中ω为包含ω_j(j=1…n)的n维未知列向量，y为包含S(g_i)/S₀(i=1…m)的m维列向量，A为m×n维的矩阵，元素为

$A_{i,j}=\underset{v=1}{\overset{z}{\mathrm{\Σ}}}{(x_v^T,c_j)}^l,x_v\∈g_i^{\⊥}---\left(10\right)$

用非负最小二乘法最小化（9）式，得到ω_j(j=1…n)，将ω_j(j=1…n)代入（5）式。

3）寻找步骤2）结果的极值，极值方向作为还原的纤维方向。

本发明的优点是：不涉及与扩散特性有关的模型的选择问题，角度分辨率高，在纤维重建过程中信号采集时间短。

（四）附图说明

图1是本发明模拟数据结果图。

（五）具体实施方式

参照附图：

一种扩散加权磁共振成像多纤维重建方法包括如下步骤：

1）获取数据，包括扩散加权系数固定的扩散加权磁共振成像数据和不施加扩散梯度脉冲时的磁共振成像数据；

2）计算扩散加权磁共振信号衰减函数的逆球面氡变换；

扩散加权磁共振信号衰减函数是指施加扩散梯度脉冲后的磁共振信号和不施加扩散梯度脉冲时的磁共振信号的比以单位扩散梯度脉冲方向g为自变量的函数，即

$h\left(g\right)=\frac{S\left(g\right)}{S_0}---\left(11\right)$

球面氡变换为

$\mathrm{\ξ}\left[L\left(x\right)\right]\left(y\right)={\∫}_{x\∈y^{\⊥}}L\left(x\right)\mathrm{dx}---\left(12\right)$

其中x、y为三维单位方向向量，L(x)为定义在球面上的函数。即L(x)球面氡变换后某一个方向上的函数值为L(x)在与该方向垂直的所有方向上值的和。

步骤2）中扩散加权磁共振信号衰减函数的逆球面氡变换的求解方法如下：

记扩散加权磁共振信号衰减函数的逆球面氡变换的结果记为f(x)，定义能量函数

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-{\∫}_{{x\∈g}_i^{\⊥}}f\left(x\right)\mathrm{dx})}^2---\left(13\right)$

其中g_i为单位三维向量，表示扩散梯度脉冲的方向；S(g_i)为扩散加权磁共振成像数据，施加的扩散梯度脉冲方向为g_i；S₀为不施加扩散梯度脉冲时的数据；m为扩散梯度脉冲方向的数目。

在单位半球面上均匀采样n个点，记为{c₁,c₂,…c_n}，用(x^Tic_j)^l,j=1…n的非负加权和估计f(x)以保证f(x)的非负性，即

$f\left(x\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{(x^T\·c_j)}^{l}s.t.{\mathrm{\ω}}_j\≥0---\left(14\right)$

其中l为偶数，ω_j(j=1…n)为权重。

将（14）式是带入（13）式

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l\mathrm{dx}\right))}^{2}s.t.{\mathrm{\ω}}_j\≥0---\left(15\right)$

整理得如下的等价形式

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l\mathrm{dx}\right))}^{2}s.t.{\mathrm{\ω}}_j\≥0---\left(16\right)$

其中的球面氡变换用均匀采样的函数值的和代替积分来近似

${\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l=\underset{v=1}{\overset{z}{\mathrm{\Σ}}}{(x_v^T\·c_j)}^l,x_v\∈g_i^{\⊥}---\left(17\right)$

其中z为采样的点数。

将（16）式是写成矩阵的形式为

$\mathrm{En}={\left|\right|\mathrm{A\ω}-y\left|\right|}_2^{2}s.t.\mathrm{\ω}\≥0---\left(18\right)$

其中ω为包含ω_j(j=1…n)的n维未知列向量，y为包含S(g_i)/S₀(i=1…m)的m维列向量，A为m×n维的矩阵，元素为

$A_{i,j}=\underset{v=1}{\overset{z}{\mathrm{\Σ}}}{(x_v^T,c_j)}^l,x_v\∈g_i^{\⊥}---\left(19\right)$

用非负最小二乘法最小化（18）式，得到ω_j(j=1…n)，将ω_j(j=1…n)代入（14）式。

3）寻找步骤2）结果的极值，极值方向作为还原的纤维方向。

图1中，模拟数据由下式产生

$S\left(g\right)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{f}_i{S}_0{e}^{{-\mathrm{bg}}^T\mathrm{Dg}}---\left(20\right)$

其中f_i表示第i跟纤维所占的比例，f₁=0.5，f₂=0.5，S₀=1，b=2000s/mm²，扩散张量D的特征值为：λ₁=1.8×10^-3mm²/s,λ₂=0.3×10^-3mm²/s,λ₃=0.3×10^-3mm²/s。实验中l=8，z=60，n=321，81个在半球面内均匀分布的扩散梯度脉冲方向，图中第二行的黑线示意两条纤维的交叉情况。

Claims (1)

1.一种扩散加权磁共振成像多纤维重建方法，包括以下步骤：

1）获取获取数据，包括扩散加权系数固定的扩散加权磁共振成像数据和不施加扩散梯度脉冲时的磁共振成像数据；

2）计算扩散加权磁共振信号衰减函数的逆球面氡变换;

步骤2）中扩散加权磁共振信号衰减函数的逆球面氡变换的求解方法如下:

记扩散加权磁共振信号衰减函数的逆球面氡变换的结果记为f(x)，定义能量函数

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-{\∫}_{{x\∈g}_i^{\⊥}}f\left(x\right)\mathrm{dx})}^2---\left(2\right)$

其中g_i为单位三维向量，表示扩散梯度脉冲的方向；S(g_i)为扩散加权磁共振成像数据，施加的扩散梯度脉冲方向为g_i；S₀为不施加扩散梯度脉冲时的数据；m为扩散梯度脉冲方向的数目。

在单位半球面上均匀采样n个点，记为集合{c₁,c₂,…c_n}，用(x^Tic_j)^l的非负加权和估计f(x)即

$f\left(x\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{(x^T\·c_j)}^{l}s.t.{\mathrm{\ω}}_j\≥0---\left(3\right)$

其中l为偶数，ω_j(j=1…n)为权重。将（3）式带入（2）式得

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-{\∫}_{{x\∈g}_i^{\⊥}}\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{(x^T\·c_j)}^l\right)\mathrm{dx})}^{2}s.t.{\mathrm{\ω}}_j\≥0---\left(4\right)$

整理得如下的等价形式

$\mathrm{En}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_j{\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l\mathrm{dx}\right))}^{2}s.t.{\mathrm{\ω}}_j\≥0---\left(5\right)$

其中的球面氡变换用均匀采样的函数值的和代替积分来近似

${\∫}_{{x\∈g}_i^{\⊥}}{(x^T\·c_j)}^l=\underset{v=1}{\overset{z}{\mathrm{\Σ}}}{(x_v^T\·c_j)}^l,x_v\∈g_i^{\⊥}---\left(6\right)$

其中z为采样的点数。

用非负最小二乘法最小化（5）式，得到解ω_j(j=1…n),将ω_j(j=1…n)代入（3）式。

3）寻找步骤2）结果的极值，极值方向作为还原的纤维方向。