一种基于迭代p阈值投影算法的压缩感知磁共振成像方法与流程

本发明涉及磁共振成像领域，尤其涉及一种基于迭代p阈值投影算法的压缩感知磁共振成像方法。
背景技术：
：磁共振成像(magneticresonanceimaging,mri)在临床诊断中扮演着重要角色，尽管mri可以提供具有优异的软组织对比度信息的高质量图像，但是它的成像速度不尽人意，主要受到物理上的(如梯度脉冲的幅度和变化率)以及生理上的(神经刺激)的约束。冗余系统下的重建模型旨在利用先验信息从不完整的观测数据中恢复原始数据。目前，相关工作主要基于两个方面：分解型和综合型模型重建。对于综合型模型的求解，迭代软阈值算法(iterativesoft-thresholdingalgorithms,ista)被广泛应用。然而ista是在变换域进行迭代运算，因此变量的维度可能比图像本身的维度高很多。ista的加速算法命名为快速迭代软阈值算法(fastiterativesoft-thresholdingalgorithms,fista)。为了解决用fista来求解分解型模型的问题，用近邻映射的moreau包络来光滑化，然后再应用fista求解光滑后的近似模型，提出了光滑化的fista(smoothfista,sfista)。对于分解型模型的求解，交替方向乘子法(alternatingdirectionmultipliermethod，admm)常常被使用，但是admm本质上是一种原始-对偶算法，因此迭代过程中既需要存储图像也需要存储冗余系统下的系数。这些都会限制冗余系统在cs-mri重建中的应用。为了能够进一步解决磁共振图像重建的问题，专家学者们开始致力于研究基于平衡型重建模型的磁共振成像方法。刘运松等人提出了一个简单而有效的基于紧框架的压缩感知磁共振图像重建算法—快速迭代软阈值投影算法(pfista，projectedfastista)，它借鉴了fista，巧妙的编程可以使得pfista避免存储整个的冗余系统下的表示系数，因此pfista可以处理高冗余度系统下的大规模mri重建问题。同时pfista的重建模型和分解型模型近似。但快速迭代软阈值投影算法使用的阈值函数收敛速度慢、精度不高，无法对小系数进行更大的惩罚，从而导致最终的成像误差高；惩罚函数的近端映射可以解决正则化反问题，因此软阈值函数作为相应的惩罚函数的近端映射十分有效。技术实现要素：根据现有技术存在的问题，本发明公开了一种基于迭代p阈值投影算法的压缩感知磁共振成像方法，包括以下步骤：s1:将t2加权脑图进行傅立叶变换,将傅里叶变换后的t2加权脑图变换到k空间，得到t2加权脑图的k空间数据，将t2加权脑图的k空间数据按照非线性采样模板进行采样，得到t2加权脑图欠采样k空间数据；s2:将t2加权脑图欠采样k空间数据进行冗余变换，将t2加权脑图欠采样k空间数据进行稀疏表示；s3:将稀疏表示的t2加权脑图欠采样k空间数据，基于迭代p阈值投影算法或快速迭代p阈值投影算法进行图像重建。进一步地，所述迭代p阈值投影算法图像重建的过程如下：s3-1-1:确定最大迭代次数m，重构精度ε，参数阈值λ和步长γ；s3-1-2:初始化：初始化待重建欠采样mri图像x0，迭代次数i，且i＝0；s3-1-3：进行p阈值投影算法迭代，利用公式(1)、(2)进行重构；xk+1＝φtγλ(φh(xk+γfhut(y-ufxk)))(1)tγλ(t)＝sgn(ti)max{0,|ti|-λ|ti|p-1}(2)其中：y表示傅立叶变换，y表示采样算子，y表示欠采样得到的k空间数据，xk+1表示第k+1迭代的图像，φh表示一个分解矩阵，当φhφ＝i，称φ为标准紧标架，γ为步长，αk+1表示第k+1迭代的系数，fh为离散傅里叶反变换，ut为u的转置。tγλ(t)为迭代p阈值函数，λ表示阈值，γ为步长；sgn(αi)表示符号函数，即当αi>0时为1，当αi<0时为-1。s3-1-4:迭代次数i＝i+1；s3-1-5:迭代停止条件：若i>m或者结束迭代，输出图像x,否则，令i＝i+1返回s-1-3,继续迭代；s3-1-6:获得输出重构图像x。进一步地，基于快速迭代p阈值投影算法进行图像重建的过程如下：s3-2-1：确定最大迭代次数m，重构精度ε，参数阈值λ和步长γ；s3-2-2:初始化待重建欠采样mri图像x0，t0＝1迭代次数i，且i＝0；s3-2-3:进行p阈值投影算法迭代，利用公式(1)、(2)、(3)和(4)进行重构；xk+1＝φtγλ(φh(xk+γfhut(y-ufxk)))(1)tγλ(t)＝sgn(ti)max{0,|ti|-λ|ti|p-1}(2)其中：y表示傅立叶变换，y表示采样算子，y表示欠采样得到的k空间数据，xk+1表示第k+1迭代的图像，φh表示一个分解矩阵，当φhφ＝i，称φ为标准紧标架，γ为步长，αk+1表示第k+1迭代的系数，fh为离散傅里叶反变换，ut为u的转置。tγλ(t)为迭代p阈值函数，λ表示阈值，γ为步长；sgn(αi)表示符号函数，即当α>0时为1，当α<0时为-1。t0＝1,tk+1表示k+1次迭代的梯度下降的方向，第k+1次利用加速技术得到的图像；s3-2-4:迭代次数i＝i+1；s3-2-5:判断迭代停止条件，若i>m或者结束迭代，输出图像x，否则，令i＝i+1返回s-2-3继续迭代；s3-2-6:获得输出重构图像x。进一步地，所述非线性采样模板包括30％高斯采样模板和30％伪放射线模板。进一步地，所述冗余变换包括离散小波变换或contourlets变换。由于采用了上述技术方案，本发明提供的一种基于迭代p阈值投影算法的压缩感知磁共振成像方法，提出迭代p阈值投影算法(projectediterativep-thresholdingalgorithms,pipta)及其加速版快速迭代p阈值投影算法(projectedfastiterativep-thresholdingalgorithms,pfipta)，p阈值函数可以看作是具有稀疏约束更宽泛的惩罚函数的映射；针对软阈值函数收缩功能较差的问题，提出将收缩性能更好的p阈值函数替换pista中的软阈值，对图像边缘噪声进行了去噪，获得了更好的重建图像；通过灵活地改变p值来设计新的稀疏目标函数，获得更好的重建效果；利用p阈值对小系数的惩罚更大，提高重构mri图像的成像速度和成像精度；采用t2加权脑图作为数据集，sidwt作为实验中的紧标架，对比多种重建算法，实验结果表明：pfipta的重建误差比其他算法的重建误差要低，pfipta的rlne值比pfista的rlne值提高了大约0.01379个单位。附图说明为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。图1为本发明中所使用的p阈值函数图；图2(a)为本发明t2加权脑图第7层图；图2(b)为本发明t2加权脑图第10层图；图2(c)为本发明的30％高斯采样模板图；图2(d)为本发明的30％伪放射线模板图；图3为本发明在t2加权脑图第7层图像的重建结果图；图4为本发明在t2加权脑图第10层图像的重建结果图；图5为本发明在t2加权脑图第7层图像的收敛曲线图；图6为本发明在t2加权脑图第10层图像的收敛曲线图；图7为本发明在contourlets标架下的收敛曲线图；图8为验证本发明不同p值的重建误差图；图9为验证本发明sidwt和contourlets下，不同p值对应的重建误差图。具体实施方式为使本发明的技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：一种基于迭代p阈值投影算法的压缩感知磁共振成像方法，其特征在于：包括以下步骤：s1:将t2加权脑图进行傅立叶变换,将傅里叶变换后的t2加权脑图变换到k空间，得到t2加权脑图的k空间数据，将t2加权脑图的k空间数据按照非线性采样模板进行采用，得到t2加权脑图欠采样k空间数据；脑图是使用一台3t的西门子triotim磁共振成像仪扫描一个健康的志愿者头部得到的图像；s2:将t2加权脑图欠采样k空间数据进行冗余变换，将t2加权脑图欠采样k空间数据进行稀疏表示；s3:将稀疏表示的t2加权脑图欠采样k空间数据，基于迭代p阈值投影算法或快速迭代p阈值投影算法进行图像重建。进一步地，所述迭代p阈值投影算法(projectediterativepthresholdingalgorithms,pipta)图像重建的过程如下：s3-1-1:确定最大迭代次数m，重构精度ε，参数阈值λ和步长γ；s3-1-2:初始化：初始化待重建欠采样重建mri图像x0，迭代次数i，且i＝0；s3-1-3：进行p阈值投影算法迭代，利用公式(1)、(2)进行重构；其中：y表示傅立叶变换，y表示采样算子，y表示欠采样得到的k空间数据，xk+1表示第k+1迭代的图像，φh表示一个分解矩阵，当φhφ＝i，称φ为标准紧标架，γ为步长，αk+1表示第k+1迭代的系数，fh为离散傅里叶反变换，ut为u的转置。tγλ(t)为迭代p阈值函数，λ表示阈值，γ为步长；sgn(αi)表示符号函数，即当αi>0时为1，当αi<0时为-1；s3-1-4:迭代次数i＝i+1；s3-1-5:迭代停止条件：若i>m或者结束迭代，输出图像x,否则，令i＝i+1返回s-1-3,继续迭代；s3-1-6:获得输出重构图像x。表一为本发明中迭代p阈值投影算法算法一：pipta参数：λ，γ初始化：x0循环直到停止：xk+1＝φtγλ(φh(xk+γfhut(y-ufxk)))tγλ(α)＝sgn(αi)max{0，|αi|-λ|αi|p-1}输出：x进一步地，基于快速迭代p阈值投影算法(projectedfastiterativep-thresholdingalgorithms，pfipta)进行图像重建的过程如下：s3-2-1：确定最大迭代次数m，重构精度ε，参数阈值λ和步长γ；s3-2-2：初始化待重建欠采样mri图像t0＝1迭代次数i，且i＝0；s3-2-3：进行p阈值投影算法迭代，利用公式(1)、(2)、(3)和(4)进行重构；xk+1＝φtγλ(φh(xk+γfhut(y-ufxk)))(1)tγλ(t)＝sgn(ti)max{0，|ti|-λ|ti|p-1}(2)其中：y表示傅立叶变换，y表示采样算子，y表示欠采样得到的k空间数据，xk+1表示第k+1迭代的图像，φh表示一个分解矩阵，当φhφ＝i，称φ为标准紧标架，γ为步长，αk+1表示第k+1迭代的系数，fh为离散傅里叶反变换，ut为u的转置。tγλ(t)为迭代p阈值函数，λ表示阈值，γ为步长；sgn(αi)表示符号函数，即当αi＞0时为1，当αi＜0时为-1。t0＝1，tk+1表示k+1次迭代的梯度下降的方向，第k+1次利用加速技术得到的图像。s3-2-4：迭代次数i＝i+1；s3-2-5：判断迭代停止条件，若i＞m或者结束迭代，输出图像x，否则，令i＝i+1返回s-2-3继续迭代；s3-2-6:获得输出重构图像x。其中，利用下式定义的相对l2范数误差来作为图像重建质量的数字指标：其中，代表重建图像，x表示原始图像。数值越小，代表磁共振图像重建质量越高。针对紧标架下模型重建问题设计迭代p阈值算法，迭代p阈值投影算法求解紧标架下模型重建问题模型，迭代软阈值投影算法的迭代公式可以解决平衡型模型重建问题，将p阈值替换软阈值，并融入到pista的核心公式中，得到pipta的迭代公式：xk+1＝φtγλ(φh(xk+γfhut(y-ufxk)))(1)其中，tγλ(t)为迭代p阈值函数，λ表示阈值，γ为步长。进一步地，在紧标架为移不变小波变换和轮廓波下，验证不同p值对应的重建图像误差，从而得到了所提算法在p为0.7的情况下重建误差最小；同时，引入nesterov加速策略，nesterov方法采用从上一步迭代点处朝前走一步处的梯度，以极少的额外计算量大幅度提高了迭代p阈值投影算法的收敛速度，从而得到快速迭代p阈值投影算法pfipta。在紧标架为移不变小波变换和轮廓波下，验证不同p值对应的重建图像误差，从而得到了所提算法在p为0.7的情况下重建误差最小。同时，引入nesterov加速策略，nesterov方法采用从上一步迭代点处朝前走一步处的梯度，以极少的额外计算量大幅度提高了迭代p阈值投影算法的收敛速度，从而得到快速迭代p阈值投影算法pfipta。表1和表2为紧标架下mri重建的pipta和pfipta算法。进一步地，所述非线性采样模板包括30％高斯采样模板和30％伪放射线模板。进一步地，所述冗余变换包括离散小波变换或contourlets变换。图2(a)为本发明t2加权脑图第7层图；图2(b)为本发明t2加权脑图第10层图；图2(a)和图2(b)中的脑图均是使用一台3t的西门子triotim磁共振成像仪扫描一个健康的志愿者头部得到的，使用32线圈，脉冲序列为t2加权的turbo自旋回波序列，分别来自第7、10层，其中tr/te＝6100/99ms,fov＝220×220mm2，厚度为3mm。图2(c)为本发明的30％高斯采样模板图；图2(d)为本发明的30％伪放射线模板图；其中图2(c)模拟三维成像中的二维相位编码，图2(d)模拟二维成像中的伪放射性线采样，它的采样点都是在笛卡尔坐标的网格上最接近真实放射线采样轨迹的点。图3为t2加权脑图第7层图像重建结果图，第一行为原图和四种算法得到的重建图像；第二行第一幅为30％伪放射线采样模板，剩余四幅为四种算法的残差图。图中fipta、sfista、pfista、pfipta算法的rlne分别为0.132426、0.097078、0.097029和0.083288。图4为加权脑图第10层重建结果图，使用30％高斯采样模板进行采样，图中fipta、sfista、pfista、pfipta的rlne分别为0.104166、0.091703、0.091573和0.069123。从图3、4可以看出，fipta的重建图像含有明显的伪影，而sfista和pfista的重建图像中这些伪影被很好地压制了，pfipta中的伪影最少。pfipta的重建误差比其他三种算法的重建误差要低。如图5、6为对应的收敛曲线，所提出的pfipta比pfista、sfista提高了大约0.01379个单位。不同的稀疏化磁共振图像的紧标架对cs-mri的重建误差影响很大，为了全面地评估pfipta，我们在本节实验中选用轮廓波(contourlets)作为紧标架，contourlets应用图像的局部边缘方向性来进一步稀疏化图像，使得重建结果能够尽可能保留边缘信息。并和求解分解型模型的通用算法admm进行了对比，实验结果如下图7所示。图中结果显示，和使用sidwt作为紧标架效果一样。pfipta比其他的重建误差都要低，而且pfipta收敛更快。因此，pfipta具备pfista的优势不随使用的紧标架而变化，并且效果更好。因为pfipta求解的是平衡型模型，而非精确的分解型模型，这部分我们就通过对比求解分解型模型的通用算法admm，了解pfipta的重建结果与分解型的重建结果的效果。从图7中可知，sfista和pfista的重建误差的admm的都是基本一样的，但是我们所提pfipta算法收敛速度更快，精度更高。我们将从数值实验验证p值是怎么影响重建速度和重建误差的。图8中，我们选用p值分别为1、0.8、0.7、0.6时的收敛曲线，p＝1为pfista算法。图9为在紧标架为sidwt和contourlets时，取不同p值对应的rlne。从图8和图9我们可以看出，随着p值减小，收敛速度变快，但是rlne并不遵循这个规律。从图8可以看出在0.8、0.7出现拐点，0.8和0.7的收敛精度差不多，但是在0.7时，速度有所提升。图9中在contourlets下所代表的曲线，也表明在0.7时重建精度更好一点。综合考虑，我们建议将p的值设置为0.7。以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本
技术领域：
的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。当前第1页12