基于结构先验信息的人体心脏左心室扩散张量估计方法与流程

本发明涉及心脏磁共振扩散张量成像信号处理技术领域，尤其涉及一种 基于结构先验信息的人体心脏左心室扩散张量估计方法。

背景技术：

磁共振扩散成像是目前唯一能够测量活体组织内水分子扩散信息的无 创方法，在组织内部水分子的扩散运动受到组织微观结构的限制，因此可以 通过测量水分子的扩散信息进而重建组织的微观结构。该方法具有非侵入和 不需造影剂等优点，在理论研究和临床实践中均引起了极大关注。最初，磁 共振扩散成像主要是用于脑白质纤维结构的成像。近年来，随着相关技术的 发展，磁共振扩散成像已用于除脑之外的多种器官的成像，如心脏、肝脏、 舌头等。

扩散张量模型是磁共振扩散成像领域的经典模型，它是用于对体素内各 个方向水分子的扩散信息进行建模，从而描述体素内组织的微观结构特征。 虽然目前已有多种更为复杂的扩散模型，如多张量模型、高阶张量模型、扩 散方向分布函数模型等，但在心脏扩散张量成像领域，仍主要采用经典的张 量模型进行相关分析。

扩散张量估计的准确性直接影响后续的基于扩散张量提取的微观结构 信息的准确性，进而影响对相关组织或器官病变的诊断的可靠性。目前，常 用的扩散张量估计方法主要包括：线性最小二乘法、加权线性最小二乘法、 非线性最小二乘法等，这些方法计算简单、速度快，但均对噪声非常敏感， 而采用磁共振扩散成像方法得到心脏扩散加权图像通常噪声较强，采用这些 方法难以得到准确的扩散张量估计结果。

技术实现要素：

本发明提供了一种基于结构先验信息的人体心脏左心室扩散张量估 计方法，解决了现有扩散张量估计方法用于左心室扩散张量估计结果准 确性差且受噪声影响严重的问题。

为了达到上述目的，本发明包括以下技术方案：一种基于结构先验 信息的人体心脏左心室扩散张量估计方法，包括以下步骤：

步骤1：估计初始扩散张量：采用线性最小二乘法估计得到每个体素 的初始扩散张量DLLS，并采用特征值分解法得到每个体素的初始主扩散方 向dLLS；

步骤2：对于每层左心室短轴切片，计算每个体素的先验纤维走向Vf；

步骤3：构造基于邻域张量平滑限制和先验纤维走向限制的目标函数；

步骤4：求解带正定限制的目标函数最小化问题，获得扩散张量的估 计值。

本发明所提供的基于结构先验信息的人体心脏左心室扩散张量估计 方法，首先采用经典方法估计初始扩散张量，并根据初始扩散张量估计 结果拟合得到先验纤维走向，然后构造基于邻域张量平滑限制和先验纤 维走向限制的目标函数，最后求解带正定限制的目标函数最小化问题， 从而得到扩散张量的估计结果。

进一步地，所述步骤2具体为：

1)根据每个体素的初始主扩散方向，计算其对应的俯仰角(Elevation Angle)，选取俯仰角的绝对值小于10度的点拟合一个椭圆Ω；

2)计算每个体素p到椭圆Ω的距离dp，如果p在椭圆Ω外，则定义其 到椭圆Ω的距离为负值，即-dp；

3)计算每个体素的初始主扩散方向对应的螺旋角(Helix Angle)，对 所有体素到椭圆Ω的距离dp和对应螺旋角进行线性拟合，得到每个体素的 螺旋角的拟合值；

4)令所有体素对应的横向角(Transverse Angle)为0度，将螺旋角 的拟合值和横向角转换为向量表示的纤维走向Vf，即获得每个体素的先 验纤维走向Vf。

进一步地，所述步骤3中所构造的目标函数形式如下：

其中，S0为无扩散加权信号；Si为第i个扩散梯度方向的扩散加权信号；N 为扩散加权梯度方向数量；gi为第i个扩散梯度方向；D为扩散张量(3×3 对称矩阵)；d为扩散张量D的向量表示(若则 d＝[D1，D2，D3，D4，D5，D6]^T)；Vm为扩散张量D的主特征向量，即主扩散方向； α和β为两正则化参数，α的取值范围为8～120，β的取值范围为0.05～ 0.15；REG1和REG2分别为两正则化项，其定义如下：给定某体素p，则，

其中，q为p的3×3×3邻域Up内的体素；加权因子wpq＝||p-q||2为体素p和 q之间的欧氏距离；为加权因子之和；H为一N×6矩阵，其第i行 定义为：

其中，gix，giy和giz分别为第i个扩散梯度方向gi的三个分量，即 gi＝[gix，giy，giz]^T，b为扩散敏感因子；

其中Vf为步骤2所得到先验纤维走向。

进一步地，所述步骤4中求解的目标函数最小化问题如下：

arg min{f(D)} s.t. λ1＞0，λ2＞0，λ3＞0 (5)

其中，λ1，λ2和λ3分别为扩散张量D的三个特征值。

进一步地，所述目标函数最小化问题通过如下步骤求解：

1)对目标函数f(D)的形式进行变换：

根据特征值分解，将扩散张量D表示为如下形式：

其中，V1，V2，V3为扩散张量D的三个特征向量；λ1，λ2，λ3分别为扩散 张量D的三个特征值；在球坐标系统下，单位向量V1表示为如下形式：

其中，θ和分别为极角(polar angle)和方位角(azimuth angle)；

构造任意单位向量Vt满足即Vt⊥V1，则特征向量V2可以表示为：

V2＝T(V1，ω)Vt (8)

其中，T(V1，ω)表示绕V1轴旋转ω角度的旋转矩阵，其形式如下：

其中，vx，vy，vz为向量V1的三个分量，即V1＝[vx，vy，vz]^T；

根据三个特征向量V1，V2，V3的正交性可知：

V3＝±V1×V2 (10)

将式(6)、(7)、(8)、(10)代入(1)，可将目标函数f(D)转化 为自变量为θ，ω，λ1，λ2，λ3的新形式

2)求解新形式的目标函数最小化问题；新形式的目标函数最小化问 题表示为：

采用受限非线性优化问题求解方法对参数最小化问题进行求解，得到使 目标函数值最小的一组参数；

3)将上述步骤求解得到的参数转化为扩散张量，即根据式(6)、 (7)、(8)、(10)将参数θ，ω，λ1，λ2，λ3转化为扩散张量D。

采用上述技术方案，包括以下有益效果：本发明所提供的基于结构 先验信息的人体心脏左心室扩散张量估计方法，通过构造基于人体心脏 左心室结构先验信息(邻域张量平滑和先验纤维走向)的正则项，对目标 函数进行约束，同时对相应的优化问题进行正定限制，从而使得扩散张 量的估计结果更加符合实际，准确性更高；能够有效的提高人体心脏左 心室扩散张量估计的准确性，与经典的线性最小二乘方法相比，本发明 方法估计得到的扩散张量的主特征方向(主扩散方向)准确性提高了约 31％，扩散各向异性提高了约14％。

附图说明

图1为本发明实施例所提供仿真数据的实际扩散张量示意图；

图2为本发明实施例所提供仿真数据的实际主扩散方向示意图；

图3为本发明实施例所提供仿真数据的实际扩散部分各向异性示意图；

图4为本发明实施例采用线性最小二乘方法得到的扩散张量示意图；

图5为本发明实施例采用线性最小二乘方法得到的主扩散方向示意图；

图6为本发明实施例采用线性最小二乘方法得到的扩散部分各向异性示 意图；

图7为采用本发明所提供估计方法得到的扩散张量示意图；

图8为采用本发明所提供估计方法得到的主扩散方向示意图；

图9为采用本发明所提供估计方法得到的扩散部分各向异性示意图。

具体实施方式

下面通过具体的实施例并结合附图对本发明做进一步的详细描述。

实施例：

本实施例中采用扩散张量模型法模拟生成了一组分辨率为32×32的 心脏左心室短轴切片的磁共振扩散加权图像仿真数据，该组数据包含30 个不同扩散梯度方向的磁共振扩散加权图像，扩散敏感因子b值为1000， 张量的三个特征值为(1.2，0.3，0.3)×10^-3mm²/s，仿真的左心室短轴切 片从心外膜到心内膜的螺旋角为-60度到60度线性变化，横向角为0度。 其各个体素的实际扩散张量如图1所示，图中的每个椭球表示其所在体 素的扩散张量，椭球的长轴方向反应该体素内的纤维走向，椭球的三个 轴的长度反应了扩散张量的三个特征值的大小。对应的主扩散方向如图2 所示，扩散部分各向异性(Fractional Anisotropy)如图3所示，在生成 的仿真数据中加入一定强度的莱斯噪声，得到一组信噪比为5的仿真磁 共振扩散加权图像数据。

本实施例提供了一种基于结构先验信息的人体心脏左心室扩散张量 估计方法对该组含噪数据进行处理的步骤如下，

步骤1：估计初始扩散张量：采用线性最小二乘法估计得到每个体素 的初始扩散张量DLLS，并采用特征值分解法得到每个体素的初始主扩散方 向dLLS；其中，所获得的初始扩散张量如图4所示，相应的初始主扩散方 向如图5所示，扩散部分各向异性如图6所示。

步骤2：对于每层左心室短轴切片，计算每个体素的先验纤维走向Vf； 该步骤具体如下：

1)根据每个体素的初始主扩散方向，计算其对应的俯仰角(Elevation Angle)，选取俯仰角的绝对值小于10度的点拟合一个椭圆Ω；

2)计算每个体素p到椭圆Ω的距离dp，如果p在椭圆Ω外，则定义其 到椭圆Ω的距离为负值，即-dp；

3)计算每个体素的初始主扩散方向对应的螺旋角(Helix Angle)，对 所有体素到椭圆Ω的距离dp和对应螺旋角进行线性拟合，得到每个体素的 螺旋角的拟合值；

4)令所有体素对应的横向角(Transverse Angle)为0度，将螺旋角 的拟合值和横向角转换为向量表示的纤维走向Vf，即获得每个体素的先 验纤维走向Vf。

步骤3：构造基于邻域张量平滑限制和先验纤维走向限制的目标函数；

所构造的目标函数形式如下：

其中，S0为无扩散加权信号；Si为第i个扩散梯度方向的扩散加权信号；N 为扩散加权梯度方向数量；gi为第i个扩散梯度方向；D为扩散张量(3×3 对称矩阵)；d为扩散张量D的向量表示(若则 d＝[D1，D2，D3，D4，D5，D6]^T)；Vm为扩散张量D的主特征向量，即主扩散方向； α和β为两正则化参数，在本实施例中α的取值为100，β的取值为0.1； REG1和REG2分别为两正则化项，其定义如下：给定某体素p，则，

其中，q为p的3×3×3邻域Up内的体素；加权因子wpq＝||p-q||2为体素p和 q之间的欧氏距离；为加权因子之和；H为一N×6矩阵，其第i行 定义为：

其中，gix，giy和giz分别为第i个扩散梯度方向gi的三个分量，即 gi＝[gix，giy，giz]^T，b为扩散敏感因子，在本实施例中b＝1000；

其中Vf为步骤2所得到先验纤维走向。

步骤4：求解带正定限制的目标函数最小化问题，获得扩散张量的估 计值，具体如下：

求解的目标函数最小化问题为：

arg min{f(D)} s.t. λ1＞0，λ2＞0，λ3＞0 (5)

其中，λ1，λ2和λ3分别为扩散张量D的三个特征值。

所述目标函数最小化问题通过如下步骤求解：

1)对目标函数f(D)的形式进行变换：

根据特征值分解，将扩散张量D表示为如下形式：

其中，V1，V2，V3为扩散张量D的三个特征向量；λ1，λ2，λ3分别为扩散 张量D的三个特征值；在球坐标系统下，单位向量V1表示为如下形式：

其中，θ和分别为极角(polar angle)和方位角(azimuth angle)；

构造任意单位向量Vt满足即Vt⊥V1，则特征向量V2可以表示为：

V2＝T(V1，ω)Vt (8)

其中，T(V1，ω)表示绕V1轴旋转ω角度的旋转矩阵，其形式如下：

其中，vx，vy，vz为向量V1的三个分量，即V1＝[vx，vy，vz]^T；

根据三个特征向量V1，V2，V3的正交性可知：

V3＝±V1×V2 (10)

将式(6)、(7)、(8)、(10)代入(1)，可将目标函数f(D)转化 为自变量为θ，ω，λ1，λ2，λ3的新形式

2)求解新形式的目标函数最小化问题；新形式的目标函数最小化问 题表示为：

采用受限非线性优化问题求解方法对参数最小化问题进行求解，得到使 目标函数值最小的一组参数；

3)将上述步骤求解得到的参数转化为扩散张量，即根据式(6)、 (7)、(8)、(10)将参数θ，ω，λ1，λ2，λ3转化为扩散张量D。 本实施例估计得到的扩散张量如图7所示，相应的主扩散方向如图8所 示，扩散部分各向异性如图9所示。

本实施例所提供的基于结构先验信息的人体心脏左心室扩散张量估 计方法，首先采用经典方法估计初始扩散张量，并根据初始扩散张量估 计结果拟合得到先验纤维走向，然后构造基于邻域张量平滑限制和先验 纤维走向限制的目标函数，最后求解带正定限制的目标函数最小化问题， 从而得到扩散张量的估计结果。与线性最小二乘法的结果相比，可以看 出，采用本实施例所提供估计方法得到的扩散张量、主扩散方向和扩散 部分各向异性均与实际的扩散张量、主扩散方向和扩散部分各向异性更 接近。定量来讲，线性最小二乘法得到的主扩散方向的角度误差为13.2 度，部分各向异性的平均误差为0.14，而本实施例所提供估计方法得到 的主扩散方向的角度误差为9.1度，部分各向异性的平均误差为0.12。 因此，与线性最小二乘法相比，本实施例所提供估计方法得到的扩散张 量的主特征方向(主扩散方向)准确性提高了约31％，扩散各向异性提高 了约14％。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对 于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明 的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在 本发明的保护范围之内。