CN102609905A - 基于迭代投影的mri图像重构方法 - Google Patents

Abstract

本发明公开一种基于迭代投影的核磁共振图像重构方法，主要克服已有技术重构图像速度慢，重构图像质量低，不利于硬件实现的问题。其实现步骤为：(1)对图像进行傅里叶变换，对变换系数用变密度欠采样矩阵进行观测采样，获得观测值，对观测值逆傅里叶变换得到初始解；(2)将得到的初始解用全变差方法进行滤波，再使用小波域双变量阈值方法对全变差方法滤波结果进行优化；(3)对优化结果进行投影，判断终止条件，最终得到最优解，输出重构图像。本发明具有重构图像时间短，重构图像边缘细节清晰，利于硬件实现的优点，可用于核磁共振仪器的成像系统中。

Description

基于迭代投影的MRI图像重构方法

技术领域

本发明属于图像处理领域，具体地说是一种在迭代投影的框架下，结合全变差TV滤波方法和小波域双变量阈值滤波方法来重构MRI图像的方法。该发明可用于解决核磁共振成像速度较慢，并且在采样率有限的情况下，提高重构图像质量和视觉效果的问题。

背景技术

在当前的医疗实践中，核磁共振成像MRI是继CT后医学影像学的又一重大进步。MRI成像是一种生物磁自旋技术，由于氢原子在人体内遍布全身，他在外加的强磁场中受到射频脉冲的激发，产生核磁共振现象。经过特殊的空间编码技术，用探测器检测到并接收以电磁形式放出的核磁共振信号。自80年代应用以来，它以极快的速度得到发展，技术日趋成熟，成为一项常规的医学检测手段，广泛应用于帕金森氏症、多发性硬化症等脑部与脊椎病以及癌症的治疗和诊断。MRI是一项非常重要的医疗成像工具，加快其成像速度一直是研究的热点问题，从技术和生理角度考虑，许多研究者正在寻求通过获得少量观测值来加快成像速度的方法，最近发展的压缩传感理论表明：若图像在某个变换域具有稀疏表示，则通过求解一个凸优化的L1最小化问题，就可由随机欠采样的傅立叶系数来进行重构。由于大部分核磁共振图像在某一变换领域都具有稀疏表示(如空间有限差分和小波变换域等)，满足了压缩传感图像重构的要求，因此在MRI中结合压缩感知理论来加快成像速度引起了人们的极大兴趣。

基于压缩感知的MRI重构算法利用MRI稀疏表示和局部光滑的先验知识，通过求解相应的优化问题来实现重构。目前已有多种算法解决此类优化问题。Lusting等提出了SparseMRI(Sparse MRI Reconstruction)算法，采用非线性共轭梯度和线性回搠的思想求解MRI重构问题，Shi等提出了TVCMRI(An Efficient Algorithm forCompressed MR Imaging Using Total Variation and Wavelets)算法，将MRI图像局部光滑特性和稀疏的先验知识相结合，利用凸函数和它共轭函数的关系特性将优化问题裂解，并采用固定点迭代算法来求解裂解后的优化问题实现MRI的图像重构。以上两种方法虽然都是将MRI图像局部光滑特性和稀疏的先验知识相结合，但仍然存在重构效果差，且求解过程复杂，不利于硬件实现的不足。

发明内容

本发明的目的在于针对上述已有技术的不足，在欠奈奎斯特采样的低速采样率下，将MRI图像局部光滑特性和稀疏的先验知识相结合，提出一种基于迭代投影的MRI图像重构方法，以缩短成像时间、提高图像质量、并利于硬件实现。

实现本发明目的的技术方案是先选择一些具有代表性的MRI图像作为原始图像，用傅里叶欠采样矩阵对原始MRI图像进行观测采样，对观测值逆傅里叶变换得到初始解，结合全变差TV滤波方法和小波域双变量阈值滤波方法对初始解进行优化，然后在凸集投影的原理下交替迭代，实现图像重构，具体实现步骤包括如下：

(1)对MRI图像x进行傅里叶变换，得到变换系数：x_f＝DFT(x)，用变密度欠采样观测矩阵Φ对图像x的傅里叶变换系数x_f进行观测采样，获得傅里叶域的观测值y＝Φx_f，对观测值y进行傅里叶反变换得到MRI图像x的初始值：x⁰＝DFT^-1(y)；

(2)结合MRI图像局部光滑特性和小波域稀疏的先验知识，将MRI图像x的初始值x⁰用全变差方法进行滤波，然后用小波域双变量阈值方法对该滤波结果x^TV进行优化，得到MRI图像x的优化结果x′；

(3)对优化结果x′进行投影，得到投影结果x^*＝x′+Φ^T(y-Φx′)，并判断x^*是否满足迭代终止条件：|D⁽ⁱ⁺¹⁾-D⁽ⁱ⁾|＜10^-5，其中

N是图像的像素个数，i表示第i次迭代，若满足条件，将投影结果x^*作为MRI图像x的重构结果并输出，否则返回步骤(2)继续执行。

所述的用小波域双变量阈值方法对全变差滤波结果x^TV进行优化，按如下步骤进行：

2.1)选取双变量阈值为： $\mathrm{Threshold}(f,\mathrm{\λ})=\frac{(\sqrt{f^2+f_p^2}-\mathrm{\λ}\frac{\sqrt{{3\mathrm{\σ}}^{\left(i\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}{\sqrt{f^2+f_p^2}}\·f,$ λ是控制算法收敛性的常数，取λ＝20，f＝Ψx^TV是x^TV的小波变换系数，Ψ表示小波变换，i表示低i次迭代，f_p是f的父系数，如果

则取阈值为0，

是系数ξ的3×3邻域边缘方差估计，σ⁽ⁱ⁾是中值估计函数，

2.2)根据确定的双变量阈值Threshold(f，λ)选择比阈值大的小波变换系数f_s保留，其余的系数置0，再对保留系数f_s进行小波反变换就得到对x^TV的优化结果x′＝Ψ^-1f_s。

本发明由于在对初始值进行全变差滤波的基础上，加入小波域双变量阈值方法对滤波结果进行优化，所以重构图像能更好的保持原图像的边缘和细节信息，从而能取得更好的图像质量和更好的视觉效果。

附图说明

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

图2是仿真实验使用的核磁共振测试图像Brain、Chest和Renal Arteries；

图3是变密度欠采样观测矩阵；

图4是用现有方法和本发明对图像Brain的重构效果比较。

具体实施方式

参照图1，本发明的具体实施步骤包括：

步骤1.对核磁共振图像x进行傅里叶变换，得到变换系数：x_f＝DFT(x)，用变密度欠采样观测矩阵Φ，对图像x的傅里叶变换系数x_f进行观测采样，获得傅里叶域的观测值y＝Φx_f，对观测值y进行傅里叶反变换得到MRI图像x的初始值：x⁰＝DFT^-1(y)。所述的变密度欠采样矩阵Φ是指采样点由中心到四周逐渐减少，和MRI图像的傅里叶变换系数分布相一致，从而能有效的采样，如图3所示，即白点为采样点，黑点为非采样点。

步骤2.将得到的初始值x⁰用全变差方法进行滤波，得到滤波结果x^TV，所述的全变差方法，参见Chambolle.“An Algorithm for Total Variation Minimization andApplications’IEEE Transactions on Signal Processing，vol.50，no.11，pp.2744-2756，November 2002。对初始值x⁰进行滤波，其实现步骤如下：

2a)定义dh(x⁰)为初始值x⁰与矩阵[1，-1，0的卷积，dv(x⁰)为初始值x⁰与矩阵[1，-1，0′]的卷积；

2b)初始化：令水平方向变差zh＝zeros(m，n)，垂直方向变差zv＝zeros(m，n)，(m，n)为初始值x⁰的大小，zeros为初始化函数；

2c)滤波过程：x¹＝dht(zh)+dvt(zv)-x⁰，计算权值w＝sqrt(dh(x¹).^2+dv(x¹).^2)，水平方向变差zh＝zh-τ*dh(x′)，垂直方向变差zv＝zv-τ*dv(x′)，

根据权值w重新计算水平方向变差zh＝zh./(1+2/λ*τ*w)和垂直方向变差zv＝zv./(1+2/λ*τ*w)，根据水平方向变差zh和垂直方向变差zv，计算残差dht(zh)为水平方向变差zh与[0，-1，1]的卷积，残差dvt(zv)为垂直方向变差zv与和[0，-1，1]′的卷积；其中参数τ＝0.249，λ＝1，*为数字相乘符号，/为数字相除符号；

2d)滤波结果为初始解减去残差，即：x^TV＝x⁰-dht(zh)-dvt(zv)。

步骤3.用小波域双变量阈值方法对x^TV进行优化，得到优化结果x′：

3a)选取双变量阈值为：选取双变量阈值为： $\mathrm{Threshold}(f,\mathrm{\λ})=\frac{\sqrt{f^2+f_p^2}-\mathrm{\λ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(i\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}}}{\sqrt{f^2+f_p^2}}\·f,$ λ是控制算法收敛性的常数，取λ＝20，f＝Ψx^TV是x^TV的小波变换系数，Ψ表示小波变换，i表示低i次迭代，f_p是f的父系数，如果

则取阈值为0，

是系数ξ3×3邻域的边缘方差估计，σ⁽ⁱ⁾是中值估计函数， ${\mathrm{\σ}}^{\left(i\right)}=\frac{\mathrm{median}\left(\right|x^{\mathrm{TV}}\left|\right)}{0.6745};$

3b)小波域双变量阈值优化：对全变差滤波结果x^TV进行小波变换，得到变换系数f＝Ψx^TV，Ψ表示小波变换；

3c)根据确定的双变量阈值Threshold(f，λ)选择比阈值大的变换系数f_s保留，其余的系数置0，再对f_s进行小波反变换就得到对x^TV的优化结果x′＝Ψ^-1f_s。

步骤4.对优化结果x′进行投影，得到投影结果x^*＝x′+Φ^T(y-Φx′)，并判断x^*是否满足迭代终止条件：|D⁽ⁱ⁾-D^(i-1))|＜10^-5，若满足该终止条件，输出MRI图像x的重构结果x^*，否则返回步骤2继续执行，

其中表示第i次迭代中优化结果和投影结果的均方误差，N是图像的像素个数，i表示第i次迭代，Φ是观测采样矩阵，T是转置符号，y是观测采样值。

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

1.仿真条件：

在CPU为Pentium 43.00GHz、内存2G、WINDOWS XP系统上使用MATLAB2010a进行了仿真。将图2所示三幅MRI图像作为测试图像，其中图2(a)为Brain图像、2(b)为Chest图像、2(c)为Renal Arteries图像。

2.仿真内容：

仿真1，对测试图像Brain，用本发明和TVCMR方法在采样率为20％和％25时进行重构仿真，结果如图4，其中图4(a)是用现有的TVCMRI方法在采样率为20％时的重构结果，图4(b)是用现有的TVCMRI方法在采样率为25％时的重构结果，图4(c)是本发明在采样率为20％时的重构结果，图4(d)是本发明在采样率为25％时的重构结果。

从图4可以看到，图4(a)，图4(b)有明显的伪影，图像边缘细节不清晰，而本发明重构结果图4(c)，图4(d)整体效果好，边缘细节清晰。

仿真2，对三幅测试图像Brain，Chest，Renal Arteries，用本发明和现有的TVCMR方法在采样率为20％时进行重构仿真，结果如表1。

表1本发明和TVCMR方法对三幅测试图像重构结果的峰值信噪比PSNR值(单位：db)：

测试图像	采样率	TVCMRI(PSNR)	本发明(PSNR)
				brain	20％	27.35	33.37
chest	20％	28.45	31.73
				Renal Arteries	20％	32.39	37.18

仿真3，对三幅测试图像Brain，Chest，Renal Arteries，用本发明和现有的TVCMR方法在采样率为25％时进行重构仿真，结果如表2。

表2本发明和TVCMR方法对三幅测试图像重构结果的峰值信噪比PSNR值(单位：db)：

测试图像	采样率	TVCMRI(PSNR)	本发明(PSNR)
				brain	25％	29.61	35.44
chest	25％	30.85	34.87
				Renal Arteries	25％	34.43	38.87

3.仿真结果：

评价图像处理效果最常用的技术指标是峰值信噪比PSNR值，峰值信噪比越大，表明重构图像和原图像的差异越小，从而图像重构效果就越好。由表1和表2可以看出，本发明在两种不同的采样率下，比TVCMRI方法的峰值信噪比都要高。说明本发明的重构图像相对原图像失真较小，从峰值信噪比和图像的视觉效果上来说，本发明比TVCMRI方法要更优越。

Claims (2)

1.一种基于迭代投影的核磁共振MRI图像重构方法，包括如下步骤：

(1)对MRI图像x进行傅里叶变换，得到变换系数：x_f＝DFT(x)，用变密度欠采样观测矩阵Φ对图像x的傅里叶变换系数x_f进行观测采样，获得傅里叶域的观测值y＝Φx_f，对观测值y进行傅里叶反变换得到MRI图像x的初始值：x⁰＝DFT^-1(y)；

(2)结合MRI图像局部光滑特性和小波域稀疏的先验知识，将MRI图像x的初始值x⁰用全变差方法进行滤波，然后用小波域双变量阈值方法对该滤波结果x^TV进行优化，得到MRI图像x的优化结果x′；

(3)对优化结果x′进行投影，得到投影结果x^*＝x′+Φ^T(y-Φx′)，并判断x^*是否满足迭代终止条件：|D⁽ⁱ⁺¹⁾-D⁽ⁱ⁾|＜10^-5，其中

N是图像的像素个数，i表示第i次迭代，若满足条件，将投影结果x^*作为MRI图像x的重构结果并输出，否则返回步骤(2)继续执行。

2.根据权利要求1所述的基于迭代投影的MRI图像重构方法，其中步骤(2)所述的用小波域双变量阈值方法对全变差滤波结果x^TV进行优化，按如下步骤进行：

2.1)选取双变量阈值为： $\mathrm{Threshold}(f,\mathrm{\λ})=\frac{(\sqrt{f^2+f_P^2}-\mathrm{\λ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(i\right)}}}{\mathrm{\σ\ξ}})}{\sqrt{f^2+f_P^2}}\·f,$ λ是控制算法收敛性的常数，取λ＝20，f＝Ψx^TV是x^TV的小波变换系数，Ψ表示小波变换，i表示低i次迭代，f_p是f的父系数，如果

则取阈值为0，

是系数ξ3×3邻域的边缘方差估计，σ⁽ⁱ⁾是中值估计函数，

2.2)根据确定的双变量阈值Threshold(f，λ)选择比阈值大的小波变换系数f_s保留，其余的系数置0，再对保留系数f_s进行小波反变换就得到对x^TV的优化结果x′＝Ψ^-1f_s。

Images