基于原始对偶字典学习的磁共振复数图像去噪

徐晓玲，刘沂玲，刘且根，张明辉

南昌大学信息工程学院，江西南昌330031



【电子与信息科学 / Electronics and Information Science】

基于原始对偶字典学习的磁共振复数图像去噪

徐晓玲，刘沂玲，刘且根，张明辉

南昌大学信息工程学院，江西南昌330031

针对磁共振(magnetic resonance, MR)幅度图像中带有不易去除的与信号相关的莱斯(Rician) 噪声问题，利用其复数图像中的实部与虚部所含噪声为不相关的加性高斯白噪声这一特性，代替对幅度图像直接去噪，提出将原始对偶字典学习(predual dictionary learning, PDL)算法用于对MR复数图像的实部与虚部分别进行去噪，然后组合得到幅度图像的方法.经仿真实验和在HT-MRSI50-50(50mm)1.2T小动物核磁共振系统中的实际应用，证明所提方法较直接对幅度图像去噪取得更好的效果，在有效去除MR图像噪声的同时能较好地保持图像中的细节.与经典的字典学习算法核奇异值分解(kernel singular value decomposition，K-SVD)相比，PDL算法去噪效果优于K-SVD算法，而运算速度提高约5倍.与经典的基于非局部相似块的三维块匹配滤波(block-matching and 3D filtering, BM3D)算法相比，在噪声水平较低时PDL算法略优于BM3D算法，噪声水平较高时BM3D算法略优于PDL算法，两者总体比较接近.

图像处理；字典学习；对偶字典学习算法；磁共振复数图像去噪；莱斯分布；核奇异值分解算法；三维块匹配滤波算法

磁共振成像(magnetic resonance imaging, MRI)因具有其无电离辐射、多参数多方位成像及高超的软组织对比度等优势，已成为临床诊断和科学研究的一种重要成像模式．然而，MRI属于低灵敏度成像，检测到的磁共振(magnetic resonance, MR)信号非常微弱，对噪声敏感，尤其是在快速和高分辨率成像时，信噪比大大降低，影响图像的视觉判读和后续的分析处理．图像去噪是图像处理的基本内容，寻求一种既能有效去除噪声又能很好地保留图像细节和边缘信息的去噪方法是人们一直追求的目标．目前，主流的图像去噪方法有小波阈值法[1]、各向异性扩散法[2]、非局部均值法(non-local means, NLM)[3]、三维块匹配滤波(block-matching and 3D filtering, BM3D)[4]、自适应字典学习法[5]和高阶奇异值分解法(high order singular value decomposition, HOSVD)[6]等．这些去噪方法基本上都是针对加性高斯白噪声而设计的，对高斯白噪声的去噪取得优良的效果．而在MRI中普遍采用幅度图像，所含噪声为与信号相关的莱斯(Rician)噪声[7]，加大了其去噪的难度．尽管一些研究对以上方法进行改进以适应MR幅度图像中的莱斯噪声[8-11]，但也增加了计算复杂度，且在低信噪比区域实际效果并不理想．

近年来，稀疏表示理论备受关注，其所应用的领域也愈加广泛，图像去噪就是其中一个方面．Elad等[12]将自适应字典学习用于图像去噪，通过核奇异值分解(kernel singular value decomposition，K-SVD)算法对图像信号进行自适应学习，针对高斯白噪声取得优良的去噪效果．本课题组提出原始对偶字典学习 (predual dictionary learning, PDL) 算法，将其用于各种自然图像的高斯噪声时取得较好效果[13]．该算法具有简洁的迭代格式，运算速度较快且具有很强的鲁棒性，其后相关改进算法先后用于MRI重建、椒盐噪声去噪、图像去模糊和图像填充等，取得良好效果[14-18]．

本研究基于课题组前期研究成果，研究MR图像的去噪问题和实际应用．磁共振成像首先通过机器扫描获取K空间数据，然后对得到的K空间数据进行傅里叶反变换重建复数图像,进一步取模得到幅度图像．K空间数据中的噪声主要源于人体和接收电路的热噪声及人体的生理噪声．K空间数据包括实部和虚部两部分，分别带有均值为0、方差相同且独立的高斯白噪声，经傅里叶线性变换所得复数图像的噪声依然为零均值加性高斯白噪声，又因求模运算为非线性运算，所得幅度图像的噪声不再是简单的加性高斯噪声，而是与信号相关的莱斯噪声．通过对磁共振成像技术从信号检测到图像形成过程中噪声特性的分析，提出以下MR图像去噪方法：代替对MR幅度图像直接进行去噪处理，采用对复数图像去噪的方法，即用PDL算法对MR复数图像的实部和虚部分别去噪，对去噪后的图像再进行取模运算，得到去噪后的幅度图像.将复杂的与信号相关的莱斯噪声去噪问题转化为两个简单的高斯噪声去噪问题．仿真结果和实际应用均表明，这种方法在用于MR图像的去噪时能获得较理想的结果，在有效去除图像噪声的同时，对MR图像中的细微组织结构依然能保持得很好．

1 理论基础

1.1 磁共振成像过程中的噪声分布

磁共振成像最常采用的成像方法为傅里叶变换法，依靠相位编码和频率编码获取K空间数据．K空间数据中的噪声主要源自被测人体和接收电路的热噪声及人体的生理噪声，经研究证实为加性、独立分布的零均值高斯白噪声.K空间数据为复数数据，即

S(kx, ky)=[SR(kx, ky)+nR(kx, ky)]+

i[SI(kx, ky)+nI(kx, ky)]

(1)

其中，x和y为空间域坐标；S(kx, ky)为含噪声复数K空间数据；SR(kx, ky)为不含噪声实部K 空间数据；nR(kx, ky)为K空间实部数据中的噪声；SI(kx, ky)为不含噪声K空间虚部数据；nI(kx, ky)为虚部K空间数据中的噪声；nR(kx, ky)和nI(kx, ky)均为加性和独立分布的零均值高斯白噪声．

对K空间数据进行傅里叶反变换得到重建的复数图像为

f(x, y)=[fR(x, y)+nR(x, y)]+

i[fI(x, y)+nI(x, y)]

(2)

其中，f(x, y)为含噪复数图像；fR(x, y)为不含噪声的实部图像；nR(x, y)为实部图像中的噪声；fI(x, y)为不含噪声的虚部图像；nI(x, y)为虚部图像中的噪声．由于傅里叶变换的线性和正交性质，不改变噪声的分布特性，nR(x, y)和nI(x, y)仍为加性、独立分布的零均值高斯白噪声．

鉴于复数图像不便观察和对相位变化导致的伪影敏感，磁共振成像中普遍采用幅度图像，通过取模运算得到幅度图像为

I(x, y)={[fR(x, y)+nR(x, y)]2+

(3)

由于求模运算的非线性，所得幅度图像的噪声不再是简单的加性高斯噪声，而是与信号相关的莱斯噪声．

1.2 PDL算法

PDL算法是在固定字典下的高效稀疏编码算法原始对偶近邻点算法(predual proximal point algorithm，PPPA)的延伸，表现为在PPPA算法的每一步小迭代后加上了一个字典更新步骤,目的是为了更好地解决式(4)在限定稀疏系数 α为非负情形下的问题．

(4)

其中，I2J为元素全为1的2J维向量；R2J+表示2J维正实数空间；α为稀疏系数；x为原图像；D为字典；τ为误差容限．

在PPPA算法的阐述中是使用原始-对偶框架求解式(4)的，具体方法可表述为

首先，证明(D)问题是(P)问题的对偶形式：

(5)

其中，u和α分别为原始和对偶变量．问题(P)的拉格朗日表达式为

(6)

然后，使用近似点算法(proximal point algorithm, PPA)求解式(5),得到迭代序列为

uTx+αT(DTu-I2J)

(7)

其中，k为迭代次数；L′为(P)问题的增广拉格朗日表达；βk为非负有界序列．

在此可注意到L′关于变量 u是βk模量强凸的，求解式(7)也就变得非常容易．因此可得解析解为

(8)

最后的困难是怎样确定变量α.Malgouyres[19]建议使用投影梯度算法进行求解来逼近原问题的解，即

max(0,αm+ρ(DTum-I2J))

(9)

其中，αm+1为第m次迭代后的稀疏系数, m为内层迭代次数；ρ为迭代步长．

可见，PPPA 是一个两循环的迭代算法.接下来将同时对字典进行更新,即PDL算法的改进之处.

由式(7)可见，若将式中的L′(u,α,ck,βk)拓展为包括字典变量D的新函数L′(u,α,ck,βk, D), 则更新字典会大大降低L′.最简单的方式是使用最速下降法更新字典D, 因此，对函数L′关于字典 D进行求导得到

Dk+μCk+1(Γk+1)T

(10)

其中，Ck+1=[um+11,um+12,…,um+1L]；Γk+1=[αk,m+11, αk,m+12,…,αk,m+1L]；μ为迭代步长．每次最速下降后，对字典的每1列进行归一化操作．

1.2 基于字典学习的图像稀疏去噪模型

(11)

其中，αij为稀疏表示系数；D为字典；x为列向量化的待恢复图像；b为列向量化的噪声图像；μij为正则化参数；Rij为图像块提取矩阵；第1项是对待恢复原始图像 x与带噪图像 b之间的总体相似程度的衡量，其中 b=x+n0，n0为列向量化的图像噪声；第2项对稀疏系数向量的稀疏度进行了约束，即稀疏约束项；在第3项中 Rij为提取子图像块的矩阵，Rijx为第ij张子图像块，Dαij是重构得到的近似子图像块．

式(11)中存在3个未知项，包括稀疏表示系数 αij、 字典 D和输出图像 x． 本研究采用分块坐标下降法求解式(11)，通过固定其中2个未知项，优化求取第3项．首先，预选一个初始字典 D(discrete cosine transform, DCT字典)，初始化 x=b，按式(11)求解以获取小图像块的最优稀疏表示系数 αij． 对所有图像小块稀疏表示，即稀疏编码．滑动窗口并对整个图像都进行处理，获得整幅图像的表示系数.然后，利用所获取的表示系数，更新字典 D，这是一个字典学习过程，字典的好坏决定了稀疏分解过程中图像信号能否获得最稀疏表示，继而影响重建图像的质量．最后，通过所获取的表示系数以及更新后的字典，求解式(11)优化问题从而更新 x，更新后的 x进入下一次迭代.此过程依次循环迭代直至满足终止条件，并输出去噪后的图像 x.

综上所述，图像去噪过程包括3部分：小图像块的稀疏分解(稀疏编码)、字典学习更新以及对去噪后的图像小块取平均处理从而得到最终去噪图像,这是一个迭代循环的过程．

2 仿真实验

为验证本研究提出的去噪方法的有效性，分别采用模拟数据和真实数据进行实验，并与经典的基于字典学习的K-SVD算法、基于非局部相似块的BM3D算法进行比较．对于PDL和K-SVD两种字典学习算法，设置子图像块大小为8×8像素，滑动步长为2，字典原子数目为128(即2倍冗余).对于BM3D采用其设定参数，子图像块大小为8×8像素，参考块间隔为3.

2.1 模拟数据实验

为模拟磁共振成像的实际过程，选取临床MR图像作为成像目标(基准参考图像)．首先对MR图像进行傅里叶变换，得到K空间数据后，分别向K空间数据的实部与虚部加入相同水平的高斯随机噪声，重建后得到复数图像．然后，分3种情况进行去噪处理：① 用PDL算法对复数图像的实部和虚部分别进行去噪处理，然后取模得到幅度图像；② 用K-SVD算法对复数图像的实部和虚部分别进行去噪处理，然后取模得到幅度图像；③ 先对复数图像进行取模得到幅度图像，再用PDL算法对幅度图像进行去噪处理．3种去噪处理结果的比较包括视觉主观比较、误差图、峰值信噪比(peak signal to noise ratio, PSNR)和结构相似度(structural similarity index, SSIM)，PDL和K-SVD算法性能的比较，还包括去噪时间的比较．

图1至图3分别为大脑、肝脏及膝盖的MR图像去噪效果比较图．其中，每组图中第1行从左到右依次为基准图像、在K空间数据中加入标准差σ=15的噪声后重建的幅度图像、用PDL算法进行复数图像去噪后的幅度图像、用K-SVD算法进行复数图像去噪后的幅度图像以及PDL算法对幅度图像进行去噪后的图像；第2行为对应的误差图.从图1至图3中不难看出，基于MR复数图像的去噪效果均十分理想，在去除噪声的同时边缘结构与图像的分辨率都保持得非常好，明显优于直接对幅度图像进行去噪，特别是在低信噪比区域．

为定量比较3种去噪情况，以肝脏图像为例，分别计算3种去噪方法在不同噪声水平下，对加噪声的肝脏MR图像去噪后所得的PSNR和SSIM值．表1为各噪声水平下，肝脏图像去噪后的PSNR值与SSIM值比较．由表1可见，无论是从PSNR值还是SSIM值来看，复数图像去噪结果均远优于幅度图像去噪的结果，而PDL算法略优于K-SVD算法．表2为加入不同噪声水平时PDL和K-SVD算法实现去噪所用时间的比较(处理平台为Intel(R)Core(TM)i7-6700CPU@3.40GHz, RAM为24.0Gbit的PC，使用Matlab R2013b软件) ．由表2可见，PDL算法的运算速度远快于K-SVD算法．

图1 当σ=15时，3种去噪方法对加噪声的大脑MR图像去噪效果对比Fig.1 The brain image’s results by three denoising methods in case of σ=15

图2 当σ=15时，3种去噪方法对加噪声的肝脏MR图像去噪效果对比Fig.2 The liver image’s results by three denoising methods in case of σ=15

图3 当σ=15时，3种去噪方法对加噪声的膝盖MR图像去噪效果对比Fig.3 The knee image’s results by three denoising methods in case of σ=15

表1 不同情况下PSNR值与SSIM值的比较Table 1 The PSNR and SSIM values under different cases

表2 不同噪声值下PDL与K-SVD去噪所用时间

Table2 The elapsed time for denoising of PDL and K-SVD under different noise levels

σPDL算法/sK-SVD算法/s528.63205.041014.9688.551312.5364.581511.3152.591810.5342.52209.7334.21258.8925.38

2.2 真实数据实验

对采集的颅脑K空间数据加入不同水平的噪声(σ= 3、5、7、10、15、20、25、30)，再分别用本研究提出的PDL算法和BM3D算法进行去噪处理．图4(a)到图4(d)依次为未加噪声时的重建图像、在K空间数据中加入σ=7的噪声后重建的图像、用PDL进行去噪处理后的图像以及用BM3D进行去噪处理后的图像．

表3为不同噪声水平下，颅脑图像去噪后的PSNR值与SSIM值比较．数据表明，在噪声水平较低时PDL略优于BM3D，噪声水平较高时BM3D略优于PDL，但总体比较接近．

图4 当σ=7时，PDL和BM3D对加噪声的颅脑MR图像去噪效果对比Fig.4 The brain image’s results of PDL and BM3D in case of σ=7

表3 不同噪声下PDL和BM3D的PSNR值与SSIM值比较

Table 3 The PSNR and SSIM values under different cases

σPSNR值/dBPDLBM3DSSIM值PDLBM3D341.830841.04780.96040.9097538.535938.28970.93010.9010736.465636.33590.90580.89201034.290134.38070.87380.88101531.984332.13420.82140.86632030.270730.65740.77660.85102529.112429.41540.74260.83433028.011528.18500.69270.7900

3 实例验证

基于所提出的去噪方法在仿真实验中取得良好去噪效果，本研究将其用于HT-MRSI50-50(50mm)1.2T小动物核磁共振系统，对实际磁共振成像过程中的图像进行去噪．为更好地观察所提出的去噪方法对细节的保持能力,对注射了血管内增强剂的小鼠进行高分辨率的3维数据采集，采集矩阵为512×512×32，扫描序列为自旋回波．去噪时先利用噪声估计函数分别估计重建所得的MR复数图像中实部与虚部的噪声水平；然后利用估计出来的噪声值确定PDL算法的参数；再分别对重建的MR复数图像的实部与虚部去噪；最后将去噪后的实部与虚部结合起来，得到MR幅度图像．图5为血管与其他信息均更丰富的小鼠第16层的未去噪和用本研究提出的去噪方法去噪后的图像比较．从图5可清楚地看到，未去噪小鼠图像受到了较大的噪声污染，而去噪图像在较好地抑制噪声的同时，各种细节与边缘均保持得很好，细小血管清晰可见．

图5 未去噪的小鼠MR图像与对其去噪后图像Fig.5 The initial rat image and the rat image denoised by PDL method

结 语

本研究提出将原始对偶字典学习算法PDL用于对MR复数图像的去噪．去噪时先利用噪声估计函数分别估计重建所得的MR复数图像中实部与虚部的噪声水平,然后利用估计出来的噪声值结合PDL算法分别对重建的MR复数图像的实部与虚部去噪，最后将去噪后的实部与虚部结合起来，得到MR幅度图像．通过仿真实验与实例验证，该方法优于直接对幅度图像去噪的效果，在去除噪声的同时能保持图像中的细节特征，达到较为理想的MR图像去噪效果．为在可便捷获得K空间数据或复数图像的情况下的MR图像去噪提供一种选择.

/ References：

[1] Chang S G, Yu Bin, Vetterli M.Adaptive wavelet thresholding for image denoising and compression[J].IEEE Transactions on Image Processing,2000, 9(9): 1532-1546.

[2] Bai Jian, Feng Xiangchu.Fractional-order anisotropic diffusion for image denoising[J].IEEE Transactions on Image Processing,2007, 16(10):2492-2502.

[3] Buades A, Coll B, Morel J M.A non-local algorithm for image denoising[C]// IEEE Conference on Computer Vision and Pattern Recognition.San Diego, USA：IEEE Computer Society Society,2005,2: 60-65.

[4] Lebrun M.An analysis and implementation of the BM3D image denoising method[J].American Society of Mechanical Engineers,2012,2(25): 175-213.

[5] Rubinstein R, Zibulevsky M, Elad M.Double sparsity: learning sparse dictionaries for sparse signal approximation[J].IEEE Transactions on Signal Processing,2010, 58(3): 1553-1564.

[6] Rajwade A, Rangarajan A, Banerjee A.Image denoising using the higher order singular value decomposition[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2013, 35(4): 849-862.

[7] Gudbjartsson H, Patz S.The Rician distribution of noisy MRI data[J].Magnetic Resonance in Medicine Official Journal, 1995, 34(6): 910-914.

[8] Liu Ryanwen, Shi Lin, Huang Wenhua, et al.Generalized total variation-based MRI Rician denoising model with spatially adaptive regularization parameters[J].Magnetic Resonance Imaging,2014, 32(6):702-720.

[9] Kang B, Choi O, Kim J D, et al.Noise reduction in magnetic resonance images using adaptive non-local means filtering[J].Electronics Letters,2013, 49(5): 324-326.

[10] Heydari M, Karami M R.A new adaptive diffusive function for magnetic resonance imaging denoising based on pixel similarity[J].Journal of Medical Signals and Sensors,2015, 5(4):201-209.

[11] Chen Geng, Zhang Pei, Wu Yafeng, et al.Denoising magnetic resonance images using collaborative non-local means[J].Neurocomputing,2016, 177:215-227.

[12] Elad M, Aharon M.Image denoising via sparse and redundant representations over learned dictionaries[J].IEEE Transactions on Image Processing,2006, 15(12): 3736-3745.

[13] Liu Qiegen, Wang Shanshan, Luo Jianhua.A novel predual dictionary learning algorithm[J].Journal of Visual Communication and Image Representation,2012,23(1): 182-193.

[14] Liu Qiegen, Wang Shanshan, Yang Kun, et al.Highly undersampled magnetic resonance image reconstruction using two-level Bregman method with dictionary updating[J].IEEE Transactions on Medical Imaging,2013, 32(7): 1290-1301.

[15] Wang Shanshan, Liu Qiegen, Xia Yong, et al.Dictionary learning based impulse noise removal via L1-L1minimization[J].Signal processing,2013, 93(9):2696-2708.

[16] Liu Qiegen, Liang Dong, Song Ying, et al.Augmented Lagrangian based sparse representation method with dictionary updating for image deblurring[J].SIAM Journal on Imaging Science,2013, 6(3): 1689-1718.

[17] Liu Qiegen, Zhang Minghui, Liang Dong.Two-level Bregmanized method for image interpolation with graph regularized sparse coding[J].Journal of Southeast University (English Edition),201329(4): 384-388.

[18] 张明辉，肖 凯，卢红阳，徐晓玲.基于加权双层Bregman及图结构正则化的磁共振成像[J].深圳大学学报理工版，2016，33(2)：119-126.Zhang Minghui, Xiao Kai, Lu Hongyang, et al.Weighted two-level Bregman method with graph regularized sparse coding for MRI reconstruction[J].Journal of Shenzhen University Science and Engineering,2016, 33(2):119-126.(in Chinese)

[19] Malgouyres F, Zeng T.A predual proximal point algorithm solving a non negative basis pursuit denoising model[J].International Journal of Computer Vision,2009, 83(3):294-311.

【中文责编：英 子；英文责编：子 兰】

A novel PDL denoising algorithm for magnetic resonance complex images

Xu Xiaoling, Liu Yiling, Liu Qiegen, and Zhang Minghui†

School of Information Engineering, Nanchang University, Nanchang 330031, Jiangxi Province, P.R.China

The noise in magnetic resonance (MR) magnitude images presents a signal-dependent Rician distribution, which is quite difficult to remove.We propose a novel denoising approach to MR images.Since the noise in MR complex images' real and imaginary parts is additive and uncorrelated zero-mean Gaussian noise, we first apply the predual dictionary learning algorithm (PDL) to respectively denoise the real and imaginary parts of the MR complex images.We, then, combine the two parts to be a denoised MR magnitude image.Extensive simulation experimental results and practical applications on small animal MRI system HT-MRSI50-50(50mm) 1.2T demonstrate that the proposed method is able to effectively remove noises while keeping the details of images.Compared with the classical dictionary learning algorithm of kernel singular value decomposition (K-SVD), the PDL algorithm not only attains better results but also requires less time, which is nearly five times faster than the K-SVD algorithm.Compared with block-matching and 3D filtering (BM3D), the proposed approach is slightly superior to BM3D in lower noise, and BM3D is slightly superior to ours in higher noise.

image processing; dictionary learning; predual dictionary learning; magnetic resonance complex images denoising; Rician distribution; kernel singular value decomposition (K-SVD); block-matching and 3D filtering

N 34

10.3724/SP.J.1249.2016.06578

国家自然科学基金资助项目(61362001)；江西省自然科学基金资助项目(20151BAB207007)

徐晓玲(1968—)，女，南昌大学副教授.研究方向：医学信号检测与处理，磁共振成像.E-mail： xuxiaoling98@163.com

引 文：徐晓玲，刘沂玲，刘且根，等．基于原始对偶字典学习的磁共振复数图像去噪[J].深圳大学学报理工版，2016，33(6)：578-585.