一种基于双字典学习的4D-MRI超分辨率重构方法

摘要

本发明公开了一种基于双字典学习的4D-MRI超分辨率重构方法，主要解决现有方法重构的4D-MRI空间分辨率较低的问题。其主要步骤为：对多层sagittal2D动态MRI进行回顾性排序，导出4D-MRI，在coronal方向上切出待超的低分辨率图像；从预先采集的多层coronal2D动态MRI中提取训练图像；再用KSVD算法对训练图像进行训练得到高、低分辨率字典；利用高、低分辨率字典之间的关系对待超的低分辨率图像进行超分辨率重构。本发明能够有效提高4D-MRI的空间分辨率，可用于多个方向的MRI超分辨率重构。

说明书

技术领域

本发明属于图像处理技术领域，涉及医学图像处理的方法，可用于4D-MRI超分辨率重构。

背景技术

4D-MRI也就是立体动态MRI，将会为现在的放射治疗方案提供经得起检验的信息，现在， 有两种的动态MRI技术，一种是3D扫描，另一种是2D多层扫描。3D扫描方法能够提供真实的 时间信息，但是图像的质量不高。2D多层扫描技术在扫描过程中能够更好的适应运动，然而需要 更多地后处理来导出3D信息。

现在，研究者们大多利用多层2D动态MRI来产生4D-MRI。其中，Siebenthaletal重复采集 了多层2DbSSFP序列数据，再评估这些序列数据与导航层的相似性来进行图像的排序，相似性 评估是用肝脏中感兴趣区域的方向性转换的代价函数来进行的，最后进行层堆叠获得3D立体图。 这种方法能够产生高连续性的腹部MRI立体图，但也存在许多缺点：第一，这种方法需要特定的 MRI序列图，但这些MRI序列图在商业系统中无法应用；第二，它将被相当复杂的图像处理过程 所阻碍；第三，怎样选择一个合适的肝脏区域来导出有代表性的平均呼吸周期还是一个不确定的 问题。

最近，Caietal.证明了一种2D-4D技术的可行性，他利用了多层bSSFP序列图像。这些动 态MRI图像是在轴向上采集的，层的位置每个呼吸周期变化一次。图像的排序过程是利用感兴趣 区域的面积进行的，也被叫做BA算法。通过减少整个BA轨迹的低频部分来导出代理呼吸轨迹， 这些呼吸相位是相对可靠的。这种方法减化了动态MRI采集和后处理的过程，但所需成像时间较 长。

ErikTryggestad提出了一种新的方法，通过回顾性的排序多层2D动态MRI来导出有代表性 的4D-MRI，为放射治疗提供帮助。

利用上述方法，我们能够获得有代表性的4D-MRI，但是，由于用来导出4D-MRI的多层2D 动态MRI分辨率较低，每层的厚度大约为10mm，所以我们在其它方向上获得的连续的2D-MRI 分辨率较低。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于双字典学习的4D-MRI超分辨率 重构方法，能够获得连续的高分辨率MRI图像。

为实现上述目的，本发明提供了一种基于双字典学习的4D-MRI超分辨率重构方法，本发明 包括如下步骤：

（1）对一个slice的sagittal2D动态MRI进行回顾性排序，得到该slice的M个呼吸相位的 平均呼吸变化图，依次对p个slice的sagittal2D动态MRI进行上述处理，可得到p个slice平均 呼吸变化图，把这p个slice的平均呼吸变化图进行堆叠，导出4D-MRI，对于任一呼吸相位，假 设2D-MRI的大小为m×m，则可在coronal方向上切出m层，大小为m×p的低分辨率图像 L_t(t＝1～m)，即为待超的低分辨率图像；

（2）预先采集n(m＞n)层coronal2D动态MRI图像S_i(i＝1～n)，在每个slice中抽取N幅不 同呼吸相位的图像，得到n×N幅高分辨率图像H_i(i＝1～n×N)，手动的在高分辨率图像 H_i(i＝1～n×N)中取一个窗口，该窗口要尽量包含待超的低分辨率图像中的信息，再在所取的窗 口中等间隔抽取p列，得到低分辨率图像L_i(i＝1～n×N)，把这n×N对高分辨率图像 H_i(i＝1～n×N)和低分辨率图像L_i(i＝1～n×N)作为训练图像；

（3）分别输入高分辨率训练图像H_i和低分辨率训练图像L_i，并采用不重叠的方式对每幅训 练图像取4×4的小块，获得初始高分辨率字典H和初始低分辨率字典L；

（4）利用KSVD算法对初始高分辨率字典H和初始低分辨率字典L进行训练，得到新的高 分辨率字典D_h和新的低分辨率字典D_l，以及高分辨率训练图像H_i的稀疏系数α_hi和低分辨率训练 图像L_i的稀疏系数α_li；

（5）输入待超的低分辨率图像L_t(t＝1～m)，利用低分辨率字典D_l，求解待超的低分辨率图 像L_t的稀疏系数α_l；

（6）分别求待超的低分辨率图像Lt和n×N幅低分辨率训练图像Li的误差：得 到待超的低分辨率图像L_t与n×N幅低分辨率训练图像中的第j幅训练图像L_j的最小误差： ${\mathrm{er}}_j=\underset{i=1}{\overset{n}{\min }}\left\{{\mathrm{er}}_i\right\};$

（7）求出待重构的高分辨率图像H_t的稀疏系数α_h；

（8）利用高分辨率字典D_h和待重构的高分辨率图像H_t的稀疏系数α_h，求得待重构的高分辨 率图像：H_t＝D_h*α_h。

所述步骤（1）所述的对一个slice的sagittal2D动态MRI进行回顾性排序，得到该slice的M 个呼吸相位的平均呼吸变化图，按如下步骤进行：

输入一个slice的sagittal2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗 口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周 期中等间隔取M个呼吸相位的图像，作为导航图，然后求剩余图像和导航图的误差，判断它们属 于哪一个呼吸相位，最后对排序后的每一个呼吸相位的图像进行平均，得到该slice的平均呼吸变 化图。

所述步骤（2）所述的预先采集n(m＞n)层coronal2D动态MRI图像S_i(i＝1～n)，在每个slice 中抽取N幅不同呼吸相位的图像，按如下步骤进行：

输入一个slice的coronal2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗 口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周 期中等间隔取N幅不同呼吸相位的图像。

所述步骤（4）所述的对初始高分辨率字典H和初始低分辨率字典L进行训练，按如下步骤进 行：

4a)对KSVD算法的优化公式：进行变形，即将其 中的表示为：

${\left|\right|Y-\mathrm{DX}\left|\right|}_F^2={\left|\right|Y-\underset{m=1}{\overset{K}{\Sigma }}{d}_m{x}_T^m\left|\right|}_F^2={\left|\right|(Y-\underset{m\ne k}{\Sigma }{d}_m{x}_T^m)-d_k{x}_T^k\left|\right|}_F^2={\left|\right|E_k-d_k{x}_T^k\left|\right|}_F^2,$

其中，Y为输入的初始字典，D为目标训练字典，X为稀疏分解矩阵，为任意第l列，||X_l||₀为X_l的0范数，为求解Y-DX的2范数，T₀为稀疏度控制系数；d_m为D的第m列原 子，为X的第m行，K为D的总列数，d_k为目标训练字典D的第k列原子，为X的第k行， E_k为不使用D的第k列原子d_k进行信号稀疏分解所产生的误差矩阵；

4b)对变形后的公式乘以矩阵Ω_k，得到目标分解公式 ${\left|\right|E_k{\Omega }_k-d_k{x}_T^k{\Omega }_k\left|\right|}_F^2={\left|\right|E_k^R-d_k{x}_R^k\left|\right|}_F^2,$

其中Ω_k的大小为P*|ω_k|，P为输入的初始字典Y的列数， |ω_k|为ω_k的模值，且Ω_k在(ω_k(m),m)处为1，其它地方全为0，其中 1≤m≤|ω_k|，ω_k(m)为ω_k的第m个数；

4c)对目标分解公式中的误差矩阵进行奇异矩阵分解得到其中U 为左奇异矩阵，V^T为右奇异矩阵，Φ为奇异值矩阵；

4d)依次取k=1,2,…,K，用左奇异矩阵U的第一列更新目标训练字典D的第k列原子，求得 更新后的字典D′，得到新的高分辨率字典D_h和新的低分辨率字典D_l；

4e)利用输入的初始字典Y和更新后的字典D′，求得稀疏分解矩阵X′，得到高分辨率图像H_i的稀疏系数α_hi和低分辨率图像L_i的稀疏系数α_li。

所述步骤（5）所述的利用低分辨率字典D_l，求解待超的低分辨率图像L_t的稀疏系数α_l；其 求解公式为：L_t=D_lα_l；

所述步骤（7）所述的求出待重构的高分辨率图像H_t的稀疏系数α_h，按如下步骤进行：

7a)求低分辨率差异矩阵：Δl=α_l-α_lj，其中，α_l是待超的低分辨率图像L_t的稀疏系数，α_lj为 低分辨率训练图像L_j的稀疏系数；

7b)由低分辨率差异矩阵Δl求出高分辨率差异矩阵Δh：

令Δh为一个元素全为零的矩阵，矩阵大小与Δl相等，求出Δl中所有不为零的元素的均值a；

找出高分辨率训练图像H_j的稀疏系数α_hj中所有元素不为零的位置，令Δh中相同位置上的元 素等于a；

7c)利用高分辨率差异矩阵Δh和高分辨率训练图像H_j的稀疏系数α_hj，求得待重构高分辨率 图像H_t的稀疏系数：α_h=α_hj+Δh。

本发明采取以上技术方案与现有的技术相比具有以下优点：

1.本发明利用预先采集的几层2D动态MRI作为训练图像来获得高、低分辨率字典，训练样 本图像与待超分辨率重构图像相关度高，重构图像的视觉效果较好，图像的细节信息都保持的比 较完整。

2.本发明在每层2D动态MRI都提取了一个完整呼吸周期的图像，使得重构图像能够保持呼 吸变化规律，能够为医生的诊断和治疗提供较大的帮助。

仿真结果表明，本发明能对4D-MRI图像进行高质量重构，提高了4D-MRI的空间分辨率。

附图说明

图1是本发明的总流程图；

图2是用本发明对第5个slice的sagittal图像进行回顾性排序的平均呼吸变化图；

图3是用本发明对第132幅coronal测试图像的超分辨率重构效果图。

具体实施方法

如图1所示，本发明的具体步骤包括如下：

（1）导出4D-MRI

1a)对一个slice的sagittal2D动态MRI进行回顾性排序，得到该slice的M个呼吸相位的平 均呼吸变化图，其排序过程如下：

输入一个slice的sagittal2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗 口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周 期中等间隔取M个呼吸相位的图像，作为导航图，然后求剩余图像和导航图的误差，判断它们属 于哪一个呼吸相位，最后对排序后的每一个呼吸相位的图像进行平均，得到该slice的平均呼吸变 化图；

1b)依次对p个slice的sagittal2D动态MRI进行上述处理，可得到p个slice平均呼吸变化 图，把这p个slice的平均呼吸变化图进行堆叠，导出4D-MRI，对于任一呼吸相位，假设2D-MRI 的大小为m×m，则可在coronal方向上切出m层，大小为m×p的低分辨率图像L_t(t＝1～m)，即 为待超的低分辨率图像；

（2）得到训练图像

预先采集n(m＞n)层coronal2D动态MRI图像S_i(i＝1～n)，在每个slice中抽取N幅不同呼 吸相位的图像，得到n×N幅高分辨率图像H_i(i＝1～n×N)，手动的在高分辨率图像H_i(i＝1～n×N) 中取一个窗口，该窗口要尽量包含待超的低分辨率图像中的信息，再在所取的窗口中等间隔抽取p 列，得到低分辨率图像L_i(i＝1～n×N)，把这n×N对高分辨率图像H_i(i＝1～n×N)和低分辨率图 像L_i(i＝1～n×N)作为训练图像；

（3）对训练图像进行预处理

分别输入高分辨率训练图像H_i和低分辨率训练图像L_i，并采用不重叠的方式对每幅训练图像 取4×4的小块，获得初始高分辨率字典H和初始低分辨率字典L。

（4）训练高、低分辨率字典

对初始高分辨率字典H和初始低分辨率字典L进行训练，得到新的高分辨率字典D_h和新的低 分辨率字典D_l，以及高分辨率训练图像H_i的稀疏系数α_hi和低分辨率训练图像L_i的稀疏系数α_li， 训练字典的方法主要有两种，分别是主成分分析法PCA和K次奇异值分解法KSVD，本实例使 用的是KSVD算法进行训练，其训练过程如下：

4a)对KSVD算法的总优化公式：进行变形，即将 其中的优化公式变形为：

${\left|\right|Y-\mathrm{DX}\left|\right|}_F^2={\left|\right|Y-\underset{m=1}{\overset{K}{\Sigma }}{d}_m{x}_T^m\left|\right|}_F^2={\left|\right|(Y-\underset{m\ne k}{\Sigma }{d}_m{x}_T^m)-d_k{x}_T^k\left|\right|}_F^2={\left|\right|E_k-d_k{x}_T^k\left|\right|}_F^2,$

其中，Y为输入的初始字典，D为目标训练字典，X为稀疏分解矩阵，为 任意第l列，||X_l||₀为X_l的0范数，为求解Y-DX的2范数，T₀为稀疏度控制系数；d_m为D的第m列原子，为X的第m行，K为D的总列数，d_k为目标训练字典D的第k列原子，为X的第k行，E_k为不使用D的第k列原子d_k进行信号稀疏分解所产生的误差矩阵；

4b)给变形后的优化公式乘以矩阵Ω_k=P*|ω_k|，得到目标分解公式 ${\left|\right|E_k{\Omega }_k-d_k{x}_T^k{\Omega }_k\left|\right|}_F^2={\left|\right|E_k^R-d_k{x}_R^k\left|\right|}_F^2,$

其中Ω_k的大小为P*|ω_k|，P为输入的初始字典Y的列数， |ω_k|为ω_k的模值，且Ω_k在(ω_k(m),m)处为1，其它地方全为0，其中 1≤m≤|ω_k|，ω_k(m)为ω_k的第m个数；

4c)对目标分解公式中的误差矩阵进行奇异矩阵分解得到其中U 为左奇异矩阵，V^T为右奇异矩阵，Φ为奇异值矩阵；

4d)依次取k=1,2,…,K，用左奇异矩阵U的第一列更新目标训练字典D的第k列原子，求得 更新后的字典D′，得到新的高分辨率字典D_h和新的低分辨率字典D_l；

4e)利用输入的初始字典Y和更新后的字典D′，求得稀疏分解矩阵X′，得到高分辨率图像H_i的稀疏系数α_hi和低分辨率图像L_i的稀疏系数α_li。

（5）预处理待超的低分辨率图像

5b)利用低分辨率字典D_l和待超的低分辨率图像L_t，求解待超的低分辨率图像L_t的稀疏系数 α_l，求解稀疏系数α_l可使用匹配追踪算法，基追踪算法，正交匹配追踪算法等，本实例使用正交 匹配追踪算法，求得初始重构图像L_t的稀疏系数α_l，其求解公式为：L_t=D_lα_l。

（6）求解高分辨率图像的稀疏系数

6a)分别求待超的低分辨率图像L_t和n×N幅低分辨率训练图像L_i的误差：

${\mathrm{er}}_i=\frac{{\left|\right|L_t-L_i\left|\right|}_2}{{\left|\right|L_t\left|\right|}_2}$

得到初始重构图像L_t与n×N幅低分辨率训练图像中的第j幅训练图像L_j的最小误差： ${\mathrm{er}}_j=\underset{i=1}{\overset{n}{\min }}\left\{{\mathrm{er}}_i\right\};$

6b1)求低分辨率差异矩阵：Δl=α_l-α_lj，其中，α_l是待超的低分辨率图像L_t的稀疏系数，α_lj为低分辨率训练图像L_j的稀疏系数；

6b2）由低分辨率差异矩阵Δl求出高分辨率差异矩阵Δh：

令Δh为一个元素全为零的矩阵，矩阵大小与Δl相等，求出Δl中所有不为零的元素的均值a；

6b3）找出高分辨率训练图像H_j的稀疏系数α_hj中所有元素不为零的位置，令Δh中相同位置 上的元素等于a；

6b4）利用高分辨率差异矩阵Δh和高分辨率训练图像H_j的稀疏系数α_hj，求得待重构高分辨 率图像H_t′的稀疏系数：α_h=α_hj+Δh。

（7）重构高分辨率图像

利用高分辨率字典D_h和待重构的高分辨率图像H_t的稀疏系数α_h，求得待重构的高分辨率图 像：H_t＝D_h*α_h；

实施例1：

本发明的效果可以通过以下实验进一步说明：

1)实验条件

本实验分别采用16个slice的coronal图像和20个slice的sagittal图像作为实验数据，每个slice 有121幅图，每幅图像的大小为192×192。

2)实验内容

对20个slice的sagittal图像进行排序，其中每个slice提取10个不同呼吸相位的图像，从而 导出4D-MRI，大小为192×192×20×10。然后，对于第一个呼吸相位，我们可以切出192幅192×20 的coronal低分辨率图像，再对这192幅coronal低分辨率图像进行超分辨率重构：

首先，对20个slice的sagittal图像进行回顾性排序，其中，图2为第5个slice的平均呼吸变 化图；

其次，对192幅coronal低分辨率图像进行超分辨率重构，其中，对第132幅coronal低分辨 率图像的超分辨重构结果如图3所示，其中图3(a)为第132幅coronal低分辨率图像、图3(b)为真 实的第132幅coronal高分辨率图像、图3(c)为本发明的重构结果；

3)实验结果分析

从图2中我们可以看出，本发明对sagittal图像进行回顾性排序时，能够提取出一个完成的呼 吸周期，保持呼吸变化规律；从图3可以看出，本发明能够较好的超分辨率重构出图像的细节信 息。

以上的英文的中文解释：slice：层；sagittal：矢状；sagittal2D：二维矢状；coronal：冠状； coronal2D：二维冠状；MRI：磁共振图像；2D-MRI：二维磁共振图像；4D-MRI：四维磁共振图 像。