CN102651125A - 基于非局部全变差的核磁共振图像重构方法 - Google Patents

Abstract

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

Description

基于非局部全变差的核磁共振图像重构方法

技术领域

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

背景技术

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

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

发明内容

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

实现本发明目的的技术思路是，通过选择一些具有代表性的核磁共振MR图像作为原始图像，用傅里叶欠采样矩阵对原始核磁共振MR图像进行观测采样，对观测值逆傅里叶变换得到初始解，结合非局部全变差NLTV方法和小波域双变量阈值滤波方法对初始解进行恢复、凸集投影、交替迭代，最终实现图像重构。其具体实现方案包括如下步骤：

(1)对核磁共振MR图像x进行傅里叶变换，得到变换系数：x_f＝DFT(x)，用Monte-Carlo变密度欠采样观测矩阵Ф对图像x的傅里叶变换系数x_f进行观测采样，获得傅里叶域的观测值y＝Фx_f，对观测值y进行傅里叶反变换得到核磁共振MR图像x的初始值：x⁰＝DFT^-1(y)，其中DFT表示傅里叶变换，DFT^-1表示傅里叶反变换；

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

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

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

上述基于非局部全变差的核磁共振MR图像重构方法，其中步骤(2)所述的对核磁共振MR图像x的初始值x⁰用非局部全变差NLTV方法进行滤波，按如下步骤进行：

2a)令x表示一个Ω空间上的MR图像信号，在位置i处的非局部梯度的定义为：

$\left({\▿}_{\mathrm{NL}}x\right)(i,j):=(x\left(i\right)-x\left(j\right))\sqrt{w(i,j)},j\∈\mathrm{\Ω},$

其中 $w(i,j)=\exp \{-\frac{\left({G_{\mathrm{\α}}}^{*}\right|x(i+.)-x(j+.){|}^2)\left(0\right)}{h^2}\}$ 表示非局部权值

这里x(·)表示以某个像素为中心的一个邻域，G_α表示以α为标准差的高斯核函数，h为低通滤波算子，非局部权值w(i，j)是用来度量中心像素分别为i和j的两个图像块的相似度；

2b)基于非局部全变差NLTV方法进行滤波就是求解以下优化式：

$x^{\mathrm{NLTV}}=\underset{x}{\arg \min }{\∫}_{\mathrm{\Ω}}{\left|{\▿}_{\mathrm{NL}}x\right|}_1\mathrm{di}+\mathrm{\μ}{\left|\right|x-x_0\left|\right|}_2^2$

上式中μ为权重系数。由于非局部正则项具有不可微性，故引入一个变量d，将上式变形为一个约束优化式：

$\begin{array}{cc}{x}^{\mathrm{NLTV}}=\underset{x,d}{\min }{\∫}_{\mathrm{\Ω}}{\left|d\right|}_1\mathrm{di}+\frac{\mathrm{\μ}}{2}{\left|\right|x-x_0\left|\right|}_2^2,& s.t.d={\▿}_{\mathrm{NL}}x\end{array}$

2c)用Split-Bregman方法对上述问题进行求解，即可得到对初始解x⁰的滤波结果x^NLTV。

上述基于非局部全变差的核磁共振MR图像重构方法，其中步骤(2)所述的用小波域双变量阈值方法对非局部全变差NLTV滤波结果x^NLTV进行优化，按如下步骤进行：

3.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^NLTV是x^NLTV的小波变换系数，Ψ表示小波变换，i表示低i次迭代，f_p是f的父系数，如果

则取阈值为0，

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

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

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

附图说明

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

图2是本发明仿真实验使用的核磁共振测试图像；

图3是本发明中的Monte-Carlo变密度欠采样观测矩阵图；

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

具体实施方式

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

步骤1.对核磁共振图像x进行傅里叶变换，得到变换系数：x_f＝DFT(x)，用Monte-Carlo变密度欠采样观测矩阵Ф，对图像x的傅里叶变换系数x_f进行观测采样，获得傅里叶域的观测值y＝Фx_f，对观测值y进行傅里叶反变换得到核磁共振MR图像x的初始值：x⁰＝DFT^-1(y)，其中DFT表示傅里叶变换，DFT^-1表示傅里叶反变换。

所述的Monte-Carlo变密度欠采样观测矩阵Ф，是指采样点由中心到四周逐渐减少，如图3所示，即白点为采样点，黑点为非采样点。

步骤2.对得到的初始值x⁰进行滤波，得到滤波结果x^NLTV。

对得到的初始值x⁰进行滤波有多种方法，例如全变差TV滤波，非局部全变差NLTV滤波等，本发明采用非局部全变差NLTV滤波。

所述的非局部全变差NLTV方法，对初始值x⁰进行滤波，其实现步骤如下：

2a)令x表示一个Ω空间上的MR图像信号，在位置i处的非局部梯度的定义为：

$\left({\▿}_{\mathrm{NL}}x\right)(i,j):=(x\left(i\right)-x\left(j\right))\sqrt{w(i,j)},j\∈\mathrm{\Ω}$

其中非局部权值w(i，j)的定义如下：

$w(i,j)=\exp \{-\frac{\left({G_{\mathrm{\α}}}^{*}\right|x(i+.)-x(j+.){|}^2)\left(0\right)}{h^2}\}$

这里x(·)表示以某个像素为中心的一个邻域，G_α表示以α为标准差的高斯核函数，h为低通滤波算子。非局部权值w(i，j)是用来度量中心像素分别为i和j的两个图像块的相似度。

2b)基于非局部全变差NLTV方法进行滤波就是求解以下优化问题：

$x^{\mathrm{NLTV}}=\underset{x}{\arg \min }{\∫}_{\mathrm{\Ω}}{\left|{\▿}_{\mathrm{NL}}x\right|}_1\mathrm{di}+\mathrm{\μ}{\left|\right|x-x_0\left|\right|}_2^2$

上式中μ为权重系数，因为非局部正则项具有不可微性，所以我们引入一个变量d，将上式变形为一个约束优化公式：

$\begin{array}{cc}{x}^{\mathrm{NLTV}}=\underset{x,d}{\min }{\∫}_{\mathrm{\Ω}}{\left|d\right|}_1\mathrm{di}+\frac{\mathrm{\μ}}{2}{\left|\right|x-x_0\left|\right|}_2^2,& s.t.d={\▿}_{\mathrm{NL}}x\end{array}$

2c)用Split-Bregman方法对上述公式进行求解，即可得到非局部全变差NLTV方法滤波结果x^NLTV。

步骤3.对滤波结果x^NLTV进行优化，得到优化结果x′：

对得到的滤波结果x^NLTV进行优化有多种方法，例如小波域软阈值优化方法，小波域硬阈值优化方法，小波域双变量阈值优化方法等，本发明采用小波域双变量阈值优化方法。

所述的小波域双变量阈值优化方法，对滤波结果x^NLTV进行优化，其实现步骤如下：

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^NLTV是x^NLTV的小波变换系数，Ψ表示小波变换，i表示低i次迭代，f_p是f的父系数，f_p是f的父系数，如果

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

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

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

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

其中

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

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

1.仿真条件：

在CPU为Pentium 4 3.00GHz、内存2G、WINDOWS XP系统上使用MATLAB2010a进行了仿真。将图2所示三幅核磁共振MR图像作为测试图像，其中图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值

测试图像	采样率	TVCMRI(PSNR)	本发明(PSNR)
				brain	20％	25.82db	29.27db
chest	20％	27.30db	30.29db
				Renal Arteries	20％	30.51db	32.33db

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

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

测试图像	采样率	TVCMRI(PSNR)	本发明(PSNR)
				brain	25％	26.97db	30.87db
chest	25％	28.13db	31.07db
				Renal Arteries	25％	31.57db	33.56db

3.仿真结果：

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

Claims (3)

1.一种基于非局部全变差的核磁共振MR图像重构方法，包括如下步骤：

(1)对核磁共振MR图像x进行傅里叶变换，得到变换系数：x_f＝DFT(x)，用Monte-Carlo变密度欠采样观测矩阵Ф对图像x的傅里叶变换系数x_f进行观测采样，获得傅里叶域的观测值y＝Φx_f，对观测值y进行傅里叶反变换得到核磁共振MR图像x的初始值：x⁰＝DFT^-1(y)，其中DFT表示傅里叶变换，DFT^-1表示傅里叶反变换；

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

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

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

2.根据权利要求1所述的基于非局部全变差的核磁共振MR图像重构方法，其中步骤(2)所述的对核磁共振MR图像x的初始值x⁰用非局部全变差NLTV方法进行滤波，按如下步骤进行：

2a)令x表示一个Ω空间上的MR图像信号，在位置i∈Ω处的非局部梯度的定义为向量：

$\left({\▿}_{\mathrm{NL}}x\right)(i,j):=(x\left(i\right)-x\left(j\right))\sqrt{w(i,j)},j\∈\mathrm{\Ω}$

其中 $w(i,j)=\exp \{-\frac{\left({G_{\mathrm{\α}}}^{*}\right|x(i+.)-x(j+.){|}^2)\left(0\right)}{h^2}\}$ 表示非局部权值这里x(·)表示以某个像素为中心的一个邻域，G_α表示以α为标准差的高斯核函数，h为低通滤波算子，非局部权值w(i，j)是用来度量中心像素分别为i和j的两个图像块的相似度；

2b)基于非局部全变差NLTV方法进行滤波就是求解以下优化式：

$x^{\mathrm{NLTV}}=\underset{x}{\arg \min }{\∫}_{\mathrm{\Ω}}{\left|{\▿}_{\mathrm{NL}}x\right|}_1\mathrm{di}+\mathrm{\μ}{\left|\right|x-x_0\left|\right|}_2^2$

上式中μ为权重系数。由于非局部正则项具有不可微性，故引入一个变量d，将上式变形为一个约束优化式：

$\begin{array}{cc}{x}^{\mathrm{NLTV}}=\underset{x,d}{\min }{\∫}_{\mathrm{\Ω}}{\left|d\right|}_1\mathrm{di}+\frac{\mathrm{\μ}}{2}{\left|\right|x-x_0\left|\right|}_2^2,& s.t.d={\▿}_{\mathrm{NL}}x\end{array}$

2c)用Split-Bregman方法对上述问题进行求解，即可得到对初始解x⁰的滤波结果x^NLTV。

3.根据权利要求1所述的基于非局部全变差的核磁共振MR图像重构方法，其中步骤(2)所述的用小波域双变量阈值方法对非局部全变差NLTV滤波结果x^NLTV进行优化，按如下步骤进行：

3.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^NLTV是x^NLTV的小波变换系数，Ψ表示小波变换，i表示低i次迭代，f_p是f的父系数，如果

则取阈值为0，

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

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

Images