本发明公开了一种基于双字典学习的4D‑MRI超分辨率重构方法，主要解决现有方法重构的4D‑MRI空间分辨率较低的问题。其主要步骤为：对多层sagittal2D动态MRI进行回顾性排序，导出4D‑MRI，在coronal方向上切出待超的低分辨率图像；从预先采集的多层coronal2D动态MRI中提取训练图像；再用KSVD算法对训练图像进行训练得到高、低分辨率字典；利用高、低分辨率字典之间的关系对待超的低分辨率图像进行超分辨率重构。本发明能够有效提高4D‑MRI的空间分辨率，可用于多个方向的MRI超分辨率重构。
1.一种基于双字典学习的4D-MRI超分辨率重构方法，其特征在于：包括如下步骤：(1)对一个slice的sagittal 2D动态MRI进行回顾性排序，得到该slice的M个呼吸相位的平均呼吸变化图，依次对p个slice的sagittal 2D动态MRI进行上述处理，可得到p个slice平均呼吸变化图，把这p个slice的平均呼吸变化图进行堆叠，导出4D-MRI，对于任一呼吸相位，假设sagittal 2D-MRI的大小为m×m，则可在coronal方向上切出m层，大小为m×p的低分辨率图像Lt,t＝1～m，即为待超的低分辨率图像；
(2)预先采集n，m>n层coronal 2D动态MRI图像Si,i＝1～n，在每个slice中抽取N幅不同呼吸相位的图像，得到n×N幅高分辨率图像Hi,i＝1～n×N，手动的在高分辨率图像Hi,i＝1～n×N中取一个窗口，该窗口要尽量包含待超的低分辨率图像中的信息，再在所取的窗口中等间隔抽取p列，得到低分辨率图像Li,i＝1～n×N，把这n×N对高分辨率图像Hi,i＝1～n×N和低分辨率图像Li,i＝1～n×N作为训练图像；
(3)分别输入高分辨率训练图像Hi和低分辨率训练图像Li，并采用不重叠的方式对每幅训练图像取4×4的小块，获得初始高分辨率字典H和初始低分辨率字典L；
(4)利用KSVD算法对初始高分辨率字典H和初始低分辨率字典L进行训练，得到新的高分辨率字典Dh和新的低分辨率字典Dl，以及高分辨率训练图像Hi的稀疏系数αhi和低分辨率训练图像Li的稀疏系数αli；
(5)输入待超的低分辨率图像Lt,t＝1～m，利用低分辨率字典Dl，求解待超的低分辨率图像Lt的稀疏系数αl；
(6)分别求待超的低分辨率图像Lt和n×N幅低分辨率训练图像Li的误差：得到待超的低分辨率图像Lt与n×N幅低分辨率训练图像中的第j幅训练图像Lj的最小误差：(7)求出待重构的高分辨率图像Ht的稀疏系数αh；
(8)利用高分辨率字典Dh和待重构的高分辨率图像Ht的稀疏系数αh，求得待重构的高分辨率图像：Ht＝Dh*αh；
所述步骤(1)所述的对一个slice的sagittal 2D动态MRI进行回顾性排序，得到该slice的M个呼吸相位的平均呼吸变化图，按如下步骤进行：输入一个slice的sagittal 2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周期中等间隔取M个呼吸相位的图像，作为导航图，然后求剩余图像和导航图的误差，判断它们属于哪一个呼吸相位，最后对排序后的每一个呼吸相位的图像进行平均，得到该slice的平均呼吸变化图；
所述步骤(2)所述的预先采集n，m>n层coronal 2D动态MRI图像Si,i＝1～n，在每个slice中抽取N幅不同呼吸相位的图像，按如下步骤进行：输入一个slice的coronal 2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周期中等间隔取N幅不同呼吸相位的图像；
所述步骤(4)所述的对初始高分辨率字典H和初始低分辨率字典L进行训练，按如下步骤进行：
4a)对KSVD算法的优化公式： Subject to 进行变形，即将其中的 表示为：
其中，Y为输入的初始字典，D为目标训练字典，X为稀疏分解矩阵， 为任意第l列，||Xl||0为Xl的0范数， 为求解Y-DX的2范数，T0为稀疏度控制系数；dm为D的第m列原子，为X的第m行，K为D的总列数，dk为目标训练字典D的第k列原子， 为X的第k行，Ek为不使用D的第k列原子dk进行信号稀疏分解所产生的误差矩阵；
4b)对变形后的公式 乘以矩阵Ωk，得到目标分解公式
其中 Ωk的大小为P*|ωk|，P为输入的初始字典Y的列数，
|ωk|为ωk的模值，且Ωk在(ωk(m),m)处为1，其它地方全为0，其中1≤m≤|ωk|，ωk(m)为ωk的第m个数；
4c)对目标分解公式 中的误差矩阵 进行奇异矩阵分解得到
其中U为左奇异矩阵，VT为右奇异矩阵，Φ为奇异值矩阵；
4d)依次取k＝1,2,…,K，用左奇异矩阵U的第一列更新目标训练字典D的第k列原子，求得更新后的字典D′，得到新的高分辨率字典Dh和新的低分辨率字典Dl；
4e)利用输入的初始字典Y和更新后的字典D′，求得稀疏分解矩阵X′，得到高分辨率训练图像Hi的稀疏系数αhi和低分辨率训练图像Li的稀疏系数αli；
所述步骤(5)所述的利用低分辨率字典Dl，求解待超的低分辨率图像Lt的稀疏系数αl；
所述步骤(7)所述的求出待重构的高分辨率图像Ht的稀疏系数αh，按如下步骤进行：
6a)求低分辨率差异矩阵：Δl＝αl-αlj，其中，αl是待超的低分辨率图像Lt的稀疏系数，αlj为低分辨率训练图像Lj的稀疏系数；
6b)由低分辨率差异矩阵Δl求出高分辨率差异矩阵Δh：
令Δh为一个元素全为零的矩阵，矩阵大小与Δl相等，求出Δl中所有不为零的元素的均值a；
找出高分辨率训练图像Hj的稀疏系数αhj中所有元素不为零的位置，令Δh中相同位置上的元素等于a；
6c)利用高分辨率差异矩阵Δh和高分辨率训练图像Hj的稀疏系数αhj，求得待重构的高分辨率图像Ht的稀疏系数：αh＝αhj+Δh。
一种基于双字典学习的4D-MRI超分辨率重构方法
技术领域
[0001] 本发明属于图像处理技术领域，涉及医学图像处理的方法，可用于4D-MRI超分辨率重构。
背景技术
[0002] 4D-MRI也就是立体动态MRI，将会为现在的放射治疗方案提供经得起检验的信息，现在，有两种的动态MRI技术，一种是3D扫描，另一种是2D多层扫描。3D扫描方法能够提供真实的时间信息，但是图像的质量不高。2D多层扫描技术在扫描过程中能够更好的适应运动，然而需要更多地后处理来导出3D信息。
[0003] 现在，研究者们大多利用多层2D动态MRI来产生4D-MRI。其中，Siebenthaletal重复采集了多层2DbSSFP序列数据，再评估这些序列数据与导航层的相似性来进行图像的排序，相似性评估是用肝脏中感兴趣区域的方向性转换的代价函数来进行的，最后进行层堆叠获得3D立体图。这种方法能够产生高连续性的腹部MRI立体图，但也存在许多缺点：第一，这种方法需要特定的MRI序列图，但这些MRI序列图在商业系统中无法应用；第二，它将被相当复杂的图像处理过程所阻碍；第三，怎样选择一个合适的肝脏区域来导出有代表性的平均呼吸周期还是一个不确定的问题。
[0004] 最近，Caietal.证明了一种2D-4D技术的可行性，他利用了多层bSSFP序列图像。这些动态MRI图像是在轴向上采集的，层的位置每个呼吸周期变化一次。图像的排序过程是利用感兴趣区域的面积进行的，也被叫做BA算法。通过减少整个BA轨迹的低频部分来导出代理呼吸轨迹，这些呼吸相位是相对可靠的。这种方法减化了动态MRI采集和后处理的过程，但所需成像时间较长。
[0005] ErikTryggestad提出了一种新的方法，通过回顾性的排序多层2D动态MRI来导出有代表性的4D-MRI，为放射治疗提供帮助。
[0006] 利用上述方法，我们能够获得有代表性的4D-MRI，但是，由于用来导出4D-MRI的多层2D动态MRI分辨率较低，每层的厚度大约为10mm，所以我们在其它方向上获得的连续的2D-MRI分辨率较低。
发明内容
[0007] 本发明的目的在于针对上述已有技术的不足，提出一种基于双字典学习的4D-MRI超分辨率重构方法，能够获得连续的高分辨率MRI图像。
[0008] 为实现上述目的，本发明提供了一种基于双字典学习的4D-MRI超分辨率重构方法，本发明包括如下步骤：
[0009] （1）对一个slice的sagittal2D动态MRI进行回顾性排序，得到该slice的M个呼吸相位的平均呼吸变化图，依次对p个slice的sagittal2D动态MRI进行上述处理，可得到p个slice平均呼吸变化图，把这p个slice的平均呼吸变化图进行堆叠，导出4D-MRI，对于任一呼吸相位，假设2D-MRI的大小为m×m，则可在coronal方向上切出m层，大小为m×p的低分辨率图像Lt(t＝1～m)，即为待超的低分辨率图像；
[0010] （2）预先采集n(m＞n)层coronal2D动态MRI图像Si(i＝1～n)，在每个slice中抽取N幅不同呼吸相位的图像，得到n×N幅高分辨率图像Hi(i＝1～n×N)，手动的在高分辨率图像Hi(i＝1～n×N)中取一个窗口，该窗口要尽量包含待超的低分辨率图像中的信息，再在所取的窗口中等间隔抽取p列，得到低分辨率图像Li(i＝1～n×N)，把这n×N对高分辨率图像Hi(i＝1～n×N)和低分辨率图像Li(i＝1～n×N)作为训练图像；
[0011] （3）分别输入高分辨率训练图像Hi和低分辨率训练图像Li，并采用不重叠的方式对每幅训练图像取4×4的小块，获得初始高分辨率字典H和初始低分辨率字典L；
[0012] （4）利用KSVD算法对初始高分辨率字典H和初始低分辨率字典L进行训练，得到新的高分辨率字典Dh和新的低分辨率字典Dl，以及高分辨率训练图像Hi的稀疏系数αhi和低分辨率训练图像Li的稀疏系数αli；
[0013] （5）输入待超的低分辨率图像Lt(t＝1～m)，利用低分辨率字典Dl，求解待超的低分辨率图像Lt的稀疏系数αl；
[0014] （6）分别求待超的低分辨率图像Lt和n×N幅低分辨率训练图像Li的误差：得到待超的低分辨率图像Lt与n×N幅低分辨率训练图像中的第j幅训练图像
Lj的最小误差：
[0015] （7）求出待重构的高分辨率图像Ht的稀疏系数αh；
[0016] （8）利用高分辨率字典Dh和待重构的高分辨率图像Ht的稀疏系数αh，求得待重构的高分辨率图像：Ht＝Dh*αh。
[0017] 所述步骤（1）所述的对一个slice的sagittal2D动态MRI进行回顾性排序，得到该slice的M个呼吸相位的平均呼吸变化图，按如下步骤进行：
[0018] 输入一个slice的sagittal2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周期中等间隔取M个呼吸相位的图像，作为导航图，然后求剩余图像和导航图的误差，判断它们属于哪一个呼吸相位，最后对排序后的每一个呼吸相位的图像进行平均，得到该slice的平均呼吸变化图。
[0019] 所述步骤（2）所述的预先采集n(m＞n)层coronal2D动态MRI图像Si(i＝1～n)，在每个slice中抽取N幅不同呼吸相位的图像，按如下步骤进行：
[0020] 输入一个slice的coronal2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周期中等间隔取N幅不同呼吸相位的图像。
[0021] 所述步骤（4）所述的对初始高分辨率字典H和初始低分辨率字典L进行训练，按如下步骤进行：
[0022] 4a)对KSVD算法的优化公式： 进行变形，即将其中的 表示为：
[0023]
[0024] 其中，Y为输入的初始字典，D为目标训练字典，X为稀疏分解矩阵， 为任意第l列，||Xl||0为Xl的0范数， 为求解Y-DX的2范数，T0为稀疏度控制系数；dm为D的第m列原子， 为X的第m行，K为D的总列数，dk为目标训练字典D的第k列原子， 为X的第k行，Ek为不使用D的第k列原子dk进行信号稀疏分解所产生的误差矩阵；
[0025] 4b)对变形后的公式 乘以矩阵Ωk，得到目标分解公式
[0026] 其中 Ωk的大小为P*|ωk|，P为输入的初始字典Y的列数，|ωk|为ωk的模值，且Ωk在(ωk(m),m)处为1，其它地方全为0，
其中1≤m≤|ωk|，ωk(m)为ωk的第m个数；
[0027] 4c)对目标分解公式 中的误差矩阵 进行奇异矩阵分解得到其中U为左奇异矩阵，VT为右奇异矩阵，Φ为奇异值矩阵；
[0028] 4d)依次取k=1,2,…,K，用左奇异矩阵U的第一列更新目标训练字典D的第k列原子，求得更新后的字典D′，得到新的高分辨率字典Dh和新的低分辨率字典Dl；
[0029] 4e)利用输入的初始字典Y和更新后的字典D′，求得稀疏分解矩阵X′，得到高分辨率图像Hi的稀疏系数αhi和低分辨率图像Li的稀疏系数αli。
[0030] 所述步骤（5）所述的利用低分辨率字典Dl，求解待超的低分辨率图像Lt的稀疏系数αl；其求解公式为：Lt=Dlαl；
[0031] 所述步骤（7）所述的求出待重构的高分辨率图像Ht的稀疏系数αh，按如下步骤进行：
[0032] 7a)求低分辨率差异矩阵：Δl=αl-αlj，其中，αl是待超的低分辨率图像Lt的稀疏系数，αlj为低分辨率训练图像Lj的稀疏系数；
[0033] 7b)由低分辨率差异矩阵Δl求出高分辨率差异矩阵Δh：
[0034] 令Δh为一个元素全为零的矩阵，矩阵大小与Δl相等，求出Δl中所有不为零的元素的均值a；
[0035] 找出高分辨率训练图像Hj的稀疏系数αhj中所有元素不为零的位置，令Δh中相同位置上的元素等于a；
[0036] 7c)利用高分辨率差异矩阵Δh和高分辨率训练图像Hj的稀疏系数αhj，求得待重构高分辨率图像Ht的稀疏系数：αh=αhj+Δh。
[0037] 本发明采取以上技术方案与现有的技术相比具有以下优点：
[0038] 1.本发明利用预先采集的几层2D动态MRI作为训练图像来获得高、低分辨率字典，训练样本图像与待超分辨率重构图像相关度高，重构图像的视觉效果较好，图像的细节信息都保持的比较完整。
[0039] 2.本发明在每层2D动态MRI都提取了一个完整呼吸周期的图像，使得重构图像能够保持呼吸变化规律，能够为医生的诊断和治疗提供较大的帮助。
[0040] 仿真结果表明，本发明能对4D-MRI图像进行高质量重构，提高了4D-MRI的空间分辨率。
附图说明
[0041] 图1是本发明的总流程图；
[0042] 图2是用本发明对第5个slice的sagittal图像进行回顾性排序的平均呼吸变化图；
[0043] 图3是用本发明对第132幅coronal测试图像的超分辨率重构效果图。
[0044] 具体实施方法
[0045] 如图1所示，本发明的具体步骤包括如下：
[0046] （1）导出4D-MRI
[0047] 1a)对一个slice的sagittal2D动态MRI进行回顾性排序，得到该slice的M个呼吸相位的平均呼吸变化图，其排序过程如下：
[0048] 输入一个slice的sagittal2D动态MRI，在图像中的横膈膜位置手动选取一个窗口，计算该窗口中像素值大于零的面积，得到该slice的呼吸曲线，从中提取一个完整的呼吸周期，在该呼吸周期中等间隔取M个呼吸相位的图像，作为导航图，然后求剩余图像和导航图的误差，判断它们属于哪一个呼吸相位，最后对排序后的每一个呼吸相位的图像进行平均，得到该slice的平均呼吸变化图；
[0049] 1b)依次对p个slice的sagittal2D动态MRI进行上述处理，可得到p个slice平均呼吸变化图，把这p个slice的平均呼吸变化图进行堆叠，导出4D-MRI，对于任一呼吸相位，假设2D-MRI的大小为m×m，则可在coronal方向上切出m层，大小为m×p的低分辨率图像Lt(t＝1～m)，即为待超的低分辨率图像；
[0050] （2）得到训练图像
[0051] 预先采集n(m＞n)层coronal2D动态MRI图像Si(i＝1～n)，在每个slice中抽取N幅不同呼吸相位的图像，得到n×N幅高分辨率图像Hi(i＝1～n×N)，手动的在高分辨率图像Hi(i＝1～n×N)中取一个窗口，该窗口要尽量包含待超的低分辨率图像中的信息，再在所取的窗口中等间隔抽取p列，得到低分辨率图像Li(i＝1～n×N)，把这n×N对高分辨率图像Hi(i＝1～n×N)和低分辨率图像Li(i＝1～n×N)作为训练图像；
[0052] （3）对训练图像进行预处理
[0053] 分别输入高分辨率训练图像Hi和低分辨率训练图像Li，并采用不重叠的方式对每幅训练图像取4×4的小块，获得初始高分辨率字典H和初始低分辨率字典L。
[0054] （4）训练高、低分辨率字典
[0055] 对初始高分辨率字典H和初始低分辨率字典L进行训练，得到新的高分辨率字典Dh和新的低分辨率字典Dl，以及高分辨率训练图像Hi的稀疏系数αhi和低分辨率训练图像Li的稀疏系数αli，训练字典的方法主要有两种，分别是主成分分析法PCA和K次奇异值分解法KSVD，本实例使用的是KSVD算法进行训练，其训练过程如下：
[0056] 4a)对KSVD算法的总优化公式： 进行变形，即将其中的优化公式 变形为：
[0057]
[0058] 其中，Y为输入的初始字典，D为目标训练字典，X为稀疏分解矩阵， 为任意第l列，||Xl||0为Xl的0范数， 为求解Y-DX的2范数，T0为稀疏度控制系数；dm为D的第m列原子， 为X的第m行，K为D的总列数，dk为目标训练字典D的第k列原子， 为X的第k行，Ek为不使用D的第k列原子dk进行信号稀疏分解所产生的误差矩阵；
[0059] 4b)给变形后的优化公式 乘以矩阵Ωk=P*|ωk|，得到目标分解公式
[0060] 其中 Ωk的大小为P*|ωk|，P为输入的初始字典Y的列数，|ωk|为ωk的模值，且Ωk在(ωk(m),m)处为1，其它地方全为0，
其中1≤m≤|ωk|，ωk(m)为ωk的第m个数；
[0061] 4c)对目标分解公式 中的误差矩阵 进行奇异矩阵分解得到T
其中U为左奇异矩阵，V为右奇异矩阵，Φ为奇异值矩阵；
[0062] 4d)依次取k=1,2,…,K，用左奇异矩阵U的第一列更新目标训练字典D的第k列原子，求得更新后的字典D′，得到新的高分辨率字典Dh和新的低分辨率字典Dl；
[0063] 4e)利用输入的初始字典Y和更新后的字典D′，求得稀疏分解矩阵X′，得到高分辨率图像Hi的稀疏系数αhi和低分辨率图像Li的稀疏系数αli。
[0064] （5）预处理待超的低分辨率图像
[0065] 5b)利用低分辨率字典Dl和待超的低分辨率图像Lt，求解待超的低分辨率图像Lt的稀疏系数αl，求解稀疏系数αl可使用匹配追踪算法，基追踪算法，正交匹配追踪算法等，本实例使用正交匹配追踪算法，求得初始重构图像Lt的稀疏系数αl，其求解公式为：Lt=Dlαl。
[0066] （6）求解高分辨率图像的稀疏系数
[0067] 6a)分别求待超的低分辨率图像Lt和n×N幅低分辨率训练图像Li的误差：
[0068]
[0069] 得到初始重构图像Lt与n×N幅低分辨率训练图像中的第j幅训练图像Lj的最小误差：
[0070] 6b1)求低分辨率差异矩阵：Δl=αl-αlj，其中，αl是待超的低分辨率图像Lt的稀疏系数，αlj为低分辨率训练图像Lj的稀疏系数；
[0071] 6b2）由低分辨率差异矩阵Δl求出高分辨率差异矩阵Δh：
[0072] 令Δh为一个元素全为零的矩阵，矩阵大小与Δl相等，求出Δl中所有不为零的元素的均值a；
[0073] 6b3）找出高分辨率训练图像Hj的稀疏系数αhj中所有元素不为零的位置，令Δh中相同位置上的元素等于a；
[0074] 6b4）利用高分辨率差异矩阵Δh和高分辨率训练图像Hj的稀疏系数αhj，求得待重构高分辨率图像Ht′的稀疏系数：αh=αhj+Δh。
[0075] （7）重构高分辨率图像
[0076] 利用高分辨率字典Dh和待重构的高分辨率图像Ht的稀疏系数αh，求得待重构的高分辨率图像：Ht＝Dh*αh；
[0077] 实施例1：
[0078] 本发明的效果可以通过以下实验进一步说明：
[0079] 1)实验条件
[0080] 本实验分别采用16个slice的coronal图像和20个slice的sagittal图像作为实验数据，每个slice有121幅图，每幅图像的大小为192×192。
[0081] 2)实验内容
[0082] 对20个slice的sagittal图像进行排序，其中每个slice提取10个不同呼吸相位的图像，从而导出4D-MRI，大小为192×192×20×10。然后，对于第一个呼吸相位，我们可以切出192幅192×20的coronal低分辨率图像，再对这192幅coronal低分辨率图像进行超分辨率重构：
[0083] 首先，对20个slice的sagittal图像进行回顾性排序，其中，图2为第5个slice的平均呼吸变化图；
[0084] 其次，对192幅coronal低分辨率图像进行超分辨率重构，其中，对第132幅coronal低分辨率图像的超分辨重构结果如图3所示，其中图3(a)为第132幅coronal低分辨率图像、图3(b)为真实的第132幅coronal高分辨率图像、图3(c)为本发明的重构结果；
[0085] 3)实验结果分析
[0086] 从图2中我们可以看出，本发明对sagittal图像进行回顾性排序时，能够提取出一个完成的呼吸周期，保持呼吸变化规律；从图3可以看出，本发明能够较好的超分辨率重构出图像的细节信息。
[0087] 以上的英文的中文解释：slice：层；sagittal：矢状；sagittal2D：二维矢状；coronal：冠状；coronal2D：二维冠状；MRI：磁共振图像；2D-MRI：二维磁共振图像；4D-MRI：
四维磁共振图像。
