本发明属于图像处理领域，尤其涉及一种基于自适应快速迭代收缩阈值算法的图像处理方法。
背景技术：
磁共振(MR)成像是一种安全、快速和准确的图像获取技术。它具有多个方向、参数和模式的优点，对人体亦无害。它能显示人体组织的解剖和功能信息。MR成像具有广泛的应用。然而，MR成像的扫描时间长，扫描速度慢，由于器官运动可能导致图像模糊，无法提供动态的实时图像和导航。所以，MR成像的缺点限制了功能成像的推广并给使用者带来额外的痛苦。
为此，该领域技术人员做了深入研究，并提出了压缩感知理论：在非常低于奈奎斯特采样率时，通过随机采样获得离散信号。在一些已知变换域中，根据信号的稀疏性通过非线性重构算法重建原始信号。
目前，构建一种基于压缩感知的稳定高效的重建算法是非常重要的。重构算法主要包括贪心算法和凸松弛算法。对于低维小规模信号，贪婪跟踪算法快速且质量好，如匹配追踪，正交匹配追踪和正则化正交匹配追踪。但是，对于这种高维度的大规模信号，很难满足重构精度的要求。凸松弛算法在重建中花费的时间更少。经典的凸函数优化算法主要包括共轭梯度，bre gman迭代和迭代重加权最小二乘。
为了提高磁共振图像重建的速度，提出了一种迭代收缩阈值算法和一系列改进的算法。这些方法直接解决了L1最小化问题。Daubechies等人提出了迭代算法(ISTA)，该迭代算法在每次迭代步骤中应用阈值(或非线性收缩)的Landweber迭代来解决线性反问题。BeckA等人提出了一种新颖的迭代收缩阈值算法(FISTA)，该算法保留了ISTA的计算简单性并提高了全局收敛速度。Xiaobo Qu等使用contourlet变换来稀疏表示曲线和边以解决CS-MRI的L1范数优化问题。Bayram I等研究了流行的迭代收缩/阈值算法的子带自适应版本，其针对每个子带采用不同的更新步骤和阈值。他们还给出了一种算法，以便为失真算子的线性和时间不变性选择适当的更新步骤和阈值。Wang提出了一种新颖的稀疏配音指数小波变换(EWT)，它提供了比传统小波变换更稀疏的系数。Elahi等人。介绍了一种改进的基于P阈值的CS-MRI图像重建迭代算法(GTIA)，使用P阈值促进了图像中的稀疏性，这是基于CS的图像重建的关键因素。
然而，上述算法是基于收缩因子在迭代中已经减小了固定步长的事实，这不适用于整个迭代过程。事实上，我们希望在早期的迭代中快速收敛，并且降低收缩速度以在后面的迭代中保持重建精度。
因此，本申请提出了一种自适应快速迭代收缩阈值算法(SAFISTA)来自适应地调整迭代收缩算子。它使用反馈来动态调整迭代步长以提高收敛速度。
技术实现要素：
：
针对现有技术的不足，本发明提供了一种基于自适应快速迭代收缩阈值算法的图像处理方法。
本发明的的技术方案为：一种基于自适应快速迭代收缩阈值算法的图像处理方法，包括如下步骤。
步骤一：采集K空间欠采样的MR图像数据，重构模型如下；
Fu＝kF (2)
其中，m∈CN是重构图像，y∈CM是MR的K空间欠采样信号，Fu是观察矩阵，k i是采样模式，F代表二维傅里叶变换。
步骤二：引入原图像在稀疏域的投影系数θ，对m进行稀疏表示m＝Ψθ，Ψ＝[ψ1,ψ2,...,ψN]∈RN×N，将上述模型的L0范数最小化问题转化为范数L1最小化问题；
m＝Ψθ (4)
步骤三：结合成像过程中的噪音，将范数L1最小化问题转化为求解投影系数θ的问题：
其中，A＝FuΨT，ε是允许误差
步骤四：结合稀疏性和K空间数据一致性，转化为拉格朗日约束条件
其中，A＝FuΨT，第一项是数据保真度项，第二项是正则化项；λ是平衡两个项的比例的正则化参数；
步骤五：令g(θ)＝||θ||1，将式(6)转变为
步骤六：对θ采用梯度下降法进行修改，使得式(7)和式(6)进行转变分别得到
和
考虑θ的连续性，得到
考虑L1范数和L2范数的平方是可分离的，将问题(8)转化为多个最小化问题，且通过阈值收缩法求解
其中，shrink(x,β)＝sign(x)·max{|x|-β,0})，λ＝max(ATy)，λk+1＝ρλk；ρ是收缩参数，为常数；
步骤七：定义自适应收缩算子ρA＝Rρ，进行迭代直至到达预设值时停止迭代，其中
步骤八：利用步骤五-八得到的θ，利用步骤一中的式(1)获得重构图。
有益效果：本方法获取的图像相比现有的三种算法而言，可以显示出更多的局部细节和更清晰的轮廓；产生较少的误差点，提供更准确的重建；收敛速度快，具有相较其他方法更高的迭代效率；对于不同部位的图像重建具有稳定性。
附图说明：
图1为原始图像。
图2为头部重建后的图像、细节图像和误差图像，采样率为20％。
图3为血管重建图像、细节图像和误差图像，采样率为20％。
图4为膝盖重建图像，细节图像和误差图像，采样率为20％。
图5为人膝盖头部不同方法的性能比较
图6为血管不同方法在的性能比较
图7为膝关节不同方法的性能比较。
图8为头部的重建后的图像、细节图像和误差图像，采样率为50％。
图9为不同采样率下的性能比较。
具体实施方式：
下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
MR图像数据采集不完整的常见模型用于解决最稀疏的解决方案是根据以下公式：
这里m∈CN是重构图像.y∈CM是MR的K空间前采样信号.Fu是观察矩阵由下式表达：
Fu＝kF (2)
这里k是采样模式and F代表二维傅里叶变换。
原始信号的重构可以看作是非确定性多项式时间(NP)的一个L0范数最小化问题。为了解决这个问题，列出稀疏系数中所有可能的非零值。然而，如果观察矩阵满足约束等距条件，则可以稀疏表示MR图像。根据下式(3)，方程式将L0范数最小化问题转化为L1范数最小化问题：
这里m由稀疏变换Ψ进行稀疏表示：
m＝Ψθ (4)
这里Ψ＝[ψ1,ψ2,...,ψN]∈RN×N，本文可采用haar小波稀疏变换。θ是原图像在稀疏域的投影系数.
等式(3)改写为：
这里FuΨT为编码矩阵，
A＝FuΨT (6)
结合成像过程造成的噪音影响，等式(5)和(6)等价于：
这里ε为允许误差，MRI重建问题包括稀疏性和k空间数据一致性，由以下拉格朗日约束条件给出：
等式(8)中,第一项是数据保真度项，第二项是正则化项。λ是平衡两个项的比例的正则化参数.
用最优化算法通过等式(8)求解θ。然后，重建图像可以通过等式(1)获得。.下面运用自适应收缩阈值算子的阈值收缩算法(SAFISTA)，用于CS-MRI重建。通过SAFISTA求解等式(8)的方法如下：
令g(θ)＝||θ||1,然后式(8)可以写为：
首先考虑解下面的等式：
为了求解等式(9),在每步迭代过程中，采用梯度下降法用来修改θ：
这里a＞0,设f(θ)满足Lipschitz-continuous条件，
采用同样步骤解式(8)，迭代公式如下：
因为θ是连续的,忽略常数项，式(13)等价于：
由于L1范数和L2范数的平方是可分离的，所以式13的解也被分离成每个元素的最小化问题，这可以通过阈值收缩法来解决：
这里shrink软阈值收缩算子，这里定义为：
shrink(x,β)＝sign(x)·max{|x|-β,0}) (16)
这里sign是符号函数。为了提高收敛速度，引入参数t和z。具体的迭代步骤如下：
这里θk-1和θk是θ的最近两次迭代值,，等式(15)改写为：
λ是用来平衡保真项f(θ)和正则项g(θ)的参数。这里λ值越大，g(θ)所占的比重越小，and f(θ)所占比用越大。由于θ是稀疏的，取初值为0。也就是说,θ在初期比较小，λ应该选择一个较大的参数。通常λ按照以下方法选择：
λ＝max(ATy) (20)
这里A由等式(6)定义。Andλ由下式更新：
λk+1＝ρλk (21)
这里ρ是收缩参数，通常取一个常数，and 0＜ρ＜1。更新后的lambda决定了式(19)中阈值收缩函数的每个迭代步长。为了进一步加快收敛速度，引入自适应收缩算子来加快初始收敛速度并保持最终的收敛精度。自适应收缩算子定义如下：
ρA＝Rρ (22)
这里R是一个参数，我们希望它在初期迭代中取值较小，在后期迭代中取值较大。根据这个规则，选择迭代最后两次图像总变分的值的比来定义R：
这里图像的TV总变分定义为：
在MRI图像重建的迭代过程中，R在初期小于1，通常在0.6左右，在迭代过程中逐渐收敛到1。λ由式(25)迭代：
迭代的停止标准定义如下：
当ε减小到预设值，停止迭代。
实施例1：
使用人体各部分的不同MR图像，包括人体头部，血管和膝关节MR图像作为实验图像。用于测试的MR图像在江苏省人民医院的1.5T Philips Achieva磁共振上使用8通道接收器线圈获得，梯度回波序列具有表1中的以下参数。
本实施例中，式(25)中q＝2，在本申请中，使用三个参数用来评估重建质量：均方误差(MSE)，峰值信噪比(PSNR)结构化相似度(SSIM)。MSE反映了估计值与原始值之间的差异程度。PSNR表示信号的最大可能功率与噪声功率的比值。SSIM是从亮度、对比度和结构信息两个图像的相似度的测量。与PSNR相比，它更符合人眼视觉特性。
图1显示了从完全采样的K空间中重建的MR图像，以人头、血管和膝关节数据作为参考。
图2-4用不同的四种方法对人体三部分进行了重建。图2-4(a)-(d)示出了使用迭代收缩算法与ISTA、FISTA、GTIA和SAFISTA的20％采样率的重建图像。图2-4(E)-(H)通过不同的算法显示了重构图像的相同部分细节。图2-4(I)-(L)为用不同的方法产生的重建误差。因为重建误差较小，为了便于观察，增加了每个误差图像的对比度。
图5-7在30次强制迭代后，用20％的采样率，用不同的方法显示了每次迭代的三个评价参数。图5-7(a)-(c)分别示出了MSE、PSNR和SSIM的性能。表2列出了使用ISTA、FISTA、GTIA和SAFISTA的实际MRI数据的MSE、PSNR和SSIM，其采样率为20％。
对人头图像在其他采样率下的重建性能进行试验。图8在50％的采样率下用不同的方法显示了重建的人头图像。此外，图9示出了MSE、PSNR和SSIM在采样率分别为30％、40％和50％的重建性能。
实施例针对三个评价参数对本发明进行了定量和定性分析。为了验证本发明有效性，比较了迭代收缩阈值算法(ISTA)、快速迭代收缩阈值算法(SUTA)〔和广义阈值迭代算法(GTIA)。
附图2-4(a)-(d)。用四种算法对人体三个部分进行重建。图2-4(e)-(h)示出了重建图像的细节。说明本发明比其他三种算法显示出更多的局部细节和更清晰的轮廓。图2-4(i)-(l)示出了重建误差的比较。从图中可以对比出，本发明产生较少的误差点，可以提供更准确的重建。图8示出了较高的采样率，本发明提供了更丰富的细节的重建图像。
在实施例中，针对迭代性能，对本发明与其他三个算法进行了比较。图5-7示出了每次迭代的MSE、PSNR和SSIM比较图。可以看出，随着迭代次数的增加，MSE逐渐减少，PSNR和SSIM逐渐增加。本发明的收敛速度明显优于其他三种方法，收敛速度明显不同。这些结果表明，本发明具有较高的迭代效率。
表2表明，该方法在人体不同部位的MSE、PSNR和SSIM方面具有较好的重建性能。这表明所提出的方法对于不同类型的MRI图像具有良好的稳定性。
表1磁共振成像参数
简写说明:
TE:echo time回声时间
TR:repetition time重复时间
FOV:field of view视矩
SL:slice切片
表2重建图像评价参数表
简写说明：
MSE:均方误差
PSNR:峰值信噪比
SSIM:结构化相似度
ISTA:迭代收缩阈值算法
FISTA:快速迭代收缩阈值算法
GTIA:广义阈值迭代算法
SAFISTA:自适应快速迭代收缩阈值算法