一种基于外积有效和字典学习的改进灵敏度编码重建方法

1.本发明涉及一种基于外积有效和字典学习的改进灵敏度编码重建方法，属于磁共振成像技术领域。


背景技术：

2.磁共振成像技术(magnetic resonance imaging,mri)一直是图像处理领域的一个重要的研究热点，被广泛的应用于临床医疗和科学研究中。然而，mri的数据采集速度较慢，为了在保证磁共振成像质量的同时提高其成像速度，研究者们将压缩感知(compressed sensing,cs)理论应用到磁共振成像领域中，并提出了并行成像技术(parallel imaging,pi)。
3.cs理论利用信号的稀疏性，从欠采样数据中准确的恢复信号。2007年由lustig等提出的压缩感知磁共振成像理论(compressed sensing-magnetic resonance imaging,cs-mri)是其中的代表性著作之一，这是一种常用的减少数据采样量的方法，能有效加速扫描过程。
4.pi是一种利用多个频率线圈并行获取信号的空间灵敏度差异来定位信号的方法。灵敏度编码(sensitivity encoding,sense)通过对采集的欠采样k空间数据应用傅里叶逆变换获得多线圈混叠图像，为了与现有的mr临床条件相结合，研究人员提出了将压缩感知和并行成像技术相结合的方法。2008年，liang等提出了sparesense方法，将灵敏度编码(sense)与压缩感知相结合以实现快速mr成像。2017年，kim等提出了将全变分(total variation,tv)-loraks正则项(tv-loraks)和sense相结合的方法，可以在高加速因子下获得更好的图像质量。2020年，鲍中文等将lp伪范数全变分(lptv)正则项引入sense模型中，并使用算子分裂技术(operator splitting,os)来求解该问题，有效的提升图像信噪比。
5.正则项的选择对mri的重建质量至关重要，基于正则化器的sense模型算法主要有基于tv-loraks正则项的tv-loraks-sense算法和基于lp伪范数全变分(lptv)正则项的lptv-sense算法，虽然这些算法都能完成图像重建，但是重建图像存在阶梯伪影和模糊伪影，细节和边缘保留完整度较低，图像重建精度有待提高。


技术实现要素：

6.本发明的目的在于克服现有技术的不足，提供基于外积有效和字典学习的改进灵敏度编码重建方法，可进一步提高磁共振成像重建图像的质量。
7.本发明采用的技术方案是：一种基于外积有效和字典学习的改进灵敏度编码重建方法，包括以下几个步骤：
8.s0：初始化，令k＝0，x0＝(rfs)hy，y，
9.其中，k是循环变量，表示数据的第k次迭代，上标“0”表示初始值；表示由nv×
nh的待重建二维图像列向量化成的列向量，n＝nv×
nh为图像总像素个数，nv、nh分别为图
像垂直方向和水平方向的像素个数，x0表示x的初始值，表示欠就是采样k空间数据，m为单线圈k空间数据实际采样的点数，q表示接收线圈的总数，(
·
)h表示矩阵的共轭转置，m＜＜n，y的第q个线圈数据用yq表示，q＝1......q表示线圈索引，y1和yq分别表示列向量化的欠采样多线圈k空间数据y的第1个线圈数据和第q个线圈数据，y0表示y的初始值，为逐线圈的二维傅里叶变换算子，为二维傅里叶变换矩阵，uh和uv分别为nh、nv点傅里叶变换矩阵，表示克罗内克积；表示多线圈k空间欠采样算子，iq为q
×
q的单位矩阵，为单线圈欠采样矩阵，对采样k空间数据进行欠采样以减少扫描时间；是灵敏度算子，其中sq是一个n
×
n的对角矩阵，表示第q个线圈的灵敏度信息，其中，s1表示第1个线圈的灵敏度信息，sq表示第q个线圈的灵敏度信息，
10.表示中间变量，x为从单线圈图像x中抽取n
p
个大小为的图像块列向量化并水平拼接成的矩阵，p1x和分别表示从单线圈图像x中提取的第1个图像块列向量和第n
p
个图像块列向量，表示的初始值，x的第i列为列向量pix，表示第i个可重叠图像块提取矩阵，从单线圈图像x中提取大小为的图像块并列向量化成列向量，n
p
表示图像块的个数；
11.矩阵表示图像块的合成字典，j为变换系数的个数，若j＝n，则d为满秩；若j》n，则d为冗余字典，d1表示字典的第1列，dj表示字典的第j列，d的第j列用dj表示，表示dj的初始值，矩阵c1表示矩阵c的第1列，cj表示矩阵c的第j列，c的第j列用cj表示，表示cj的初始值，为fsx对应的辅助变量，为z对应的拉格朗日乘子，表示uz的初始值，
12.s1：初始化j＝1；
13.s2：计算第k+1次迭代的矩阵c
k+1
，计算公式如下：
[0014][0015]
当变量下标为“0”时，该数不存在；
[0016]
s3：计算第k+1次迭代的字典d
k+1
，计算公式如下：
[0017][0018]
s4：计算第k+1次迭代的中间变量b
k+1
，计算公式如下：
[0019][0020]
s5：计算第k+1次迭代的稀疏系数计算公式如下：
[0021]
[0022]
其中，硬阈值算子h
λ
(
·
)定义为下标“i”表示向量元素索引，bi表示向量b的第i个元素，l为一个标量，正数λ为控制整体稀疏性的权重，
⊙
表示逐元素乘法，z＝min(a,l)表示向量中每个元素与标量l的最小值，“∠”表示复数的相位；
[0023]
s6：计算第k+1次迭代的中间变量h
k+1
，计算公式如下：
[0024][0025]
s7：计算第k+1次迭代的字典原子系数计算公式如下：
[0026][0027]
其中，||
·
||2表示l2范数，v为n
×
n单位矩阵的第1列；
[0028]
s8：判断是否j＝j，当j＝j时进入步骤s10；否则j＝j+1，返回步骤s2；
[0029]
s9：计算第k+1次迭代的稀疏编码矩阵w
k+1
，计算公式如下：
[0030]wk+1
＝(c
k+1
)h；
[0031]wk+1
的第i列表示为wi即为矩阵w第i个图像块的稀疏系数向量，
[0032]
s10：计算第k+1次迭代的辅助变量z
k+1
，计算公式如下：
[0033][0034]
其中，r
t
r为对角矩阵，xk表示列向量化的待重建图像x的第k次迭代，上标“t”表示矩阵的转置，表示拉格朗日乘子uz第k次迭代的值，μ1＞0，(
·
)-1
表示求逆算子；
[0035]
s11：计算第k+1次迭代的待重建图像x
k+1
，计算公式如下：
[0036][0037]
其中，是对角矩阵，其对角元素等于对应的像素点出现在所有图像块中的次数，shs也为对角矩阵，μ＞0；
[0038]
s12：计算第k+1次迭代的中间变量x
k+1
，计算公式如下：
[0039][0040]
s13：更新第k+1次迭代的拉格朗日乘子
[0041]
s14：判断是否达到最大迭代次数，若达到则进入s15，否则令k＝k+1，返回s1；
[0042]
s15：输出x＝x
k+1
，得到重建图像x；
[0043]
本发明的有益效果是：sense是一种并行mri的重建技术，通过求解关于线圈灵敏度的固有编码信息来重建图像。本发明将包含l0范数的外积有效和字典学习(efficient sum of outer products dictionary learning with an l0 norm,soup-dillo)引入灵敏
度编码(sense)技术中，增强了灵敏度编码技术的重建性能，提高了磁共振成像精度。提出的算法主要基于块坐标下降法和admm(alternating direction method of multipliers,admm)技术求解，通过字典更新和图像更新两个步骤实现mri重建。实验结果表明，与其它基于灵敏度编码(sense)模型的算法(如tv-loraks-sense和lptv-sense)相比，本发明提出的新算法soup-dillo-sense能够有效的抑制噪声和去除图像的伪影，较完整地保留图像边缘部分。
附图说明
[0044]
图1为本发明方法流程图；
[0045]
图2为使用8通道的头部线圈进行完全采样来获取受试者的体内人脑切片数据即(dataset1)图；
[0046]
图3为具有5倍加速因子和24
×
24中心全采样自校准区域的二维泊松圆盘欠采样掩模图；
[0047]
图4～6分别为使用tv-loraks-sense算法、lptv-sense算法和soup-dillo-sense算法从含有5倍加速因子和24
×
24中心全采样自校准区域的二维泊松圆盘欠采样掩模进行欠采样的dataset1重建出的图像；
[0048]
图7～9分别为图4～6重建图像对应的误差图(白点表示误差，误差图越白表示误差越大)。
具体实施方式
[0049]
下面结合附图进一步详细描述本发明的技术方案：
[0050]
实施例1：本发明是基于sense框架提出的一种基于外积有效和字典学习的改进灵敏度编码重建方法。
[0051]
1)sense框架：
[0052]
假设表示由nv×
nh的待重建二维图像列向量化成的列向量，n＝nv×
nh为图像总像素个数，nv、nh分别为图像垂直方向和水平方向的像素个数；表示欠采样k空间数据，m为单线圈k空间数据实际采样的点数，q表示接收线圈的总数，(
·
)h表示矩阵的共轭转置，m＜＜n。y的第q个线圈数据用yq表示，q＝1......q表示线圈索引，y1和yq分别表示列向量化的欠采样多线圈k空间数据y的第1个线圈数据和第q个线圈数据。
[0053]
然后，得到sense重建模型为：
[0054]
y＝rfsx
ꢀꢀꢀ
(1)
[0055]
其中，为逐线圈的二维傅里叶变换算子，为二维傅里叶变换矩阵，uv和uh分别为nv、nh的点傅里叶变换矩阵，表示克罗内克积；表示多线圈k空间欠采样算子，iq为q
×
q的单位矩阵，为单线圈欠采样矩阵，可对k空间数据进行欠采样以减少扫描时间；是灵敏度算子，其中sq是一个n
×
n的对角矩阵，表示第q个线圈的灵敏度信息，s1表示第1个线圈的灵敏度图，sq表示第q个线圈的灵敏度图。
[0056]
2)算法推导：
[0057]
尽管结合了tv-loraks和lptv正则项方法的sense算法在一定程度上可以提高mri重建质量，但仍存在很大的不足之处。为了进一步提高mri重建质量，本发明将包含l0范数的外积有效和字典学习(soup-dillo)引入sense模型中，提出soup-dillo-sense算法，得到的优化问题表示如下：
[0058][0059]
其中，n
p
表示图像块的个数，表示第i个可重叠图像块提取矩阵，从单线圈图像x中提取大小为的图像块并列向量化成列向量。矩阵表示图像块的合成字典(j为系数个数，若j＝n，则d为满秩；若j》n，则d为冗余字典)，dj表示字典的第j列；矩阵w的第i列表示为wi，wi即为矩阵w第i个图像块的稀疏系数向量，正数λ为控制整体稀疏性的权重，l》0。最优化问题(2)的求解可分为字典更新和图像更新两个步骤，即：
[0060][0061][0062]
字典学习：令以及则问题(3)可以表示为：
[0063][0064]
令和cj是矩阵c的第j列，为了方便描述，下列推导中用λ2来表示则可以得出以下soup-dillo公式：
[0065][0066]
将式(6)中的矩阵dch表示为(稀疏)秩为一的矩阵或外积的和则问题(6)可以表示为：
[0067][0068]
问题(7)可以使用块坐标下降法来求解。可将问题(7)分解为稀疏编码和字典原子更新两个子问题：
[0069][0070]
[0071]
稀疏编码：在该过程中保持其它变量固定不变，求解问题(8)中的cj。
[0072]
假设是一个基于所有其它原子和系数的最新值的固定矩阵，则问题(8)可转化成：
[0073][0074]
当且仅当向量没有幅度为λ的元素时，(10)存在极小值且唯一。给定和稀疏编码问题(10)的全局极小值通过以下截断硬阈值操作获得：
[0075][0076]
其中，
⊙
表示逐元素乘法。假设l＞λ成立，对于向量z＝min(a,l)表示a中每个元素与标量l的最小值，“∠”表示复数的相位。硬阈值算子h
λ
(
·
)定义为：
[0077][0078]
其中，其下标“i”用来索引向量条目，使用bi来表示向量b的第i个元素。引入中间变量则优化问题(11)可以转换为：
[0079][0080][0081]
公式(13)和(14)中变量的上标“k+1”和“k”表示算法的第k+1次迭代与第k次迭代。b
k+1
表示第k+1次迭代的中间变量，xk表示第k次迭代的中间变量，表示第k次迭代的字典原子系数dj，表示第k次迭代的稀疏系数cj，表示第k+1次迭代的稀疏系数cj，表示第k+1次迭代的变量矩阵c，表示第k+1次迭代的变量字典d，当变量下标为“0”时，该数不存在。
[0082]
字典原子更新：保持其它变量固定，求解字典原子dj。借助前面的表达，问题(9)可转换为如下问题：
[0083][0084]
因此，问题(15)的一个全局极小值是：
[0085][0086]
此处将v设置为n
×
n单位矩阵的第一列，当且仅当cj＝0时，它存在唯一的解。引入中间变量h＝e
jcj
，则优化问题可以转换为：
[0087]
[0088][0089]
式(17)，(18)中h
k+1
表示第k+1次迭代的中间变量，表示第k+1次迭代的字典原子系数dj。
[0090]
则第k+1次迭代的稀疏编码矩阵w
k+1
由下式可得：
[0091]wk+1
＝(c
k+1
)hꢀꢀꢀ
(19)
[0092]wk+1
的第i列表示为
[0093]
图像更新：引入辅助变量以及相应的拉格朗日乘子则优化问题(4)可以转换为：
[0094][0095][0096]
式(21)中是第k+1次迭代的辅助变量z
k+1
对应的拉格朗日乘子，x
k+1
表示第k+1次迭代的待重建图像。
[0097]
使用admm方法求解，则问题(20)可以被分解成如下几个问题：
[0098][0099][0100][0101]
式(22)(23)(24)中z为fsx对应的辅助变量，z
k+1
表示第k+1次迭代的辅助变量，x
k+1
表示第k+1次迭代的中间变量，μ1＞0,μ＞0。xk表示第k次迭代的待重构图像，x
k+1
表示第k+1次迭代的待重构图像。
[0102]
保持其它变量固定，使得(22)的目标函数对z的导数为0，可得辅助变量z的方程：
[0103][0104]
由于对角矩阵r
t
r的逆易得，于是可得z的解析解：
[0105][0106]
控制其它变量固定，使得(23)的目标函数对x的导数为0，可得重建图像x的方程：
[0107][0108]
其中，f为二维傅里叶变换矩阵。由于shs与均为对角矩阵，则
也为对角矩阵，其逆易得，从而可得问题(23)的解析解为：
[0109][0110]
综上所述，获得了一种新的具有soup-dillo正则项的sense并行mri重建算法，称为souo-dillo-sense算法。
[0111]
具体流程图如图1所示，其中步骤如下：
[0112]
s0：初始化，令k＝0，x0＝(rfs)hy，y，
[0113]
其中，上标“0”表示初始值；是一个列向量，包含要重建的未知图像的样本，表示dj的初始值，表示cj的初始值，表示中间变量，x为从x中抽取n
p
个大小为的图像块列向量化并水平拼接成的矩阵，p1x和分别表示从单线圈图像x中提取的第1个图像块和列向量第n
p
个图像块列向量，表示的初始值；为fsx对应的辅助变量，为相应的拉格朗日乘子，表示uz的初始值。
[0114]
s1：初始化j＝1；
[0115]
s2：计算第k+1次迭代的矩阵
[0116]
s3：计算第k+1次迭代的字典
[0117]
s4：计算第k+1次迭代的中间变量b
k+1
，计算公式如(13)；
[0118]
s5：计算第k+1次迭代的稀疏系数计算公式如(14)；
[0119]
s6：计算第k+1次迭代的中间变量h
k+1
，计算公式如下(17)；
[0120]
s7：计算第k+1次迭代的字典原子系数计算公式如(18)；
[0121]
s8：判断是否j＝j，当j＝j时进入步骤s9；否则j＝j+1，返回步骤s2；
[0122]
s9：计算第k+1次迭代的稀疏编码矩阵w
k+1
，计算公式如(19)；
[0123]
s10：计算第k+1次迭代的辅助变量z
k+1
，计算公式如(26)；
[0124]
s11：计算第k+1次迭代的待重建图像x
k+1
，计算公式如(28)；
[0125]
s12：计算第k+1次迭代的中间变量x
k+1
，计算公式如(24)；
[0126]
s13：更新第k+1次迭代的拉格朗日乘子计算公式如(21)；
[0127]
s14：判断是否达到最大迭代次数，若达到则进入s15，否则令k＝k+1，返回s1；
[0128]
s15：输出x＝x
k+1
，得到重建图像x。
[0129]
实验结果：
[0130]
在以下实验中，为了测试soup-dillo-sense算法的重建性能，将提出的soup-dillo-sense算法分别与tv-loraks-sense和lptv-sense算法进行比较。所有算法均能对图像进行重建，所有实验均在matlab(mathworks，natick，ma)上实现。以下实验均在配置为intel core i7-6500u@3.1ghz的处理器、8gb内存和windows10操作系统(64位)的笔记本电
脑中进行的。
[0131]
为了比较各算法的重建性能，本发明使用1张人类脑部切片图进行仿真实验，命名为dataset1，如图2所示，dataset1为8通道头部线圈进行完全采样获取受试者的体内人脑切片图。使用加速因子为af＝5和24
×
24中心全采样自校准区域的二维泊松圆盘欠采样掩模，如图3所示，对全采样的数据集进行人工采样生成测试数据集。
[0132]
首先，对数据集dataset1下的三个重建算法进行视觉比较，本发明选择对加速因子为5倍(欠采样率为20％)时的重建图像进行比较，如图4～6所示。图4～6分别依次展示了tv-loraks-sense算法、lptv-sense算法和soup-dillo-sense算法的重建图像。图4是使用算法tv-loraks-sense进行重建得到的图像，可以看出，该图像存在模糊伪影、细节部分保留不太完整。图5是使用算法lptv-sense进行重建得到的图像，该算法的重建质量也不是很好不能很好的去除图像的伪影。图6是使用算法soup-dillo-sense进行重建得到的图像，该算法能有效的去除模糊伪影，保留了其它算法缺失的细节部分，边缘轮廓信息保留较完整，且重建图像更接近原始图像。
[0133]
为了进一步说明本发明提出的新算法soup-diilo-sense的重建性能优于算法tv-loraks-sense和lptv-sense，图7～9展示了dataset1在三种算法下的误差图(误差图白点越多表明误差越大)。由误差图可以看出，tv-loraks-sense和lptv-sense的误差图白点较多，图像重建质量较差，soup-dillo-sense误差图白点最少，重建误差较小，所提出的算法的重建性能较好。
[0134]
综上所述，对选用的数据集dataset1在不同的重建算法中进行实验，soup-dillo-sense算法在视觉上的表现效果良好，重建图像质量优于tv-loraks-sense算法和lptv-sense算法。
[0135]
以上结合附图对本发明的具体实施方式作了详细分析说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。