CN102592267B

CN102592267B - 基于稀疏表示的医学超声图像滤波方法

Info

Publication number: CN102592267B
Application number: CN201210002260.3A
Authority: CN
Inventors: 郭翌; 汪源源
Original assignee: Fudan University
Current assignee: Fudan University
Priority date: 2012-01-06
Filing date: 2012-01-06
Publication date: 2014-09-03
Anticipated expiration: 2032-01-06
Also published as: CN102592267A

Abstract

本发明属于图像处理技术领域，具体为基于稀疏表示的医学超声图像滤波方法。本发明使用过完备字典将二维超声图像分解成最为稀疏的系数表示；采用基于子块相似度的全局滤波思想，对于每一个图像子块，利用全局冗余信息寻找整幅图像中与其最相似的子块；根据得到的所有最相似子块，对原子块内全部像素点一同滤波；对于各重叠区域内存在多个估值的像素点，对其进行均值处理，从而实现整幅图像的滤波。本发明非常适合于医学超声图像的滤波，可以有效地降低斑点噪声，并且较好地保持超声图像中的纹理信息和边缘特性。

Description

基于稀疏表示的医学超声图像滤波方法

技术领域

本发明属于图像处理技术领域，具体为一种基于稀疏表示的医学超声图像滤波方法。

背景技术

超声成像因具有无创、实时、可重复性强、诊疗费用低廉等优点，在医学临床诊断中有着重要的应用。但是超声图像中存在着大量的斑点，掩盖了某些细节信息，严重影响细微特征的分辨和图像质量，尤其是使边缘与细节变得模糊，增加了医学诊断的难度。

常用的斑点噪声滤波方法包括局部统计法[1]、中值滤波法[2]、各向异性滤波法[3]和小波滤波法[4]等。这些方法可以不同程度地降低噪声影响，但仍存在一些问题。对于基于局部特征的滤波算法而言，窗的尺寸和形状极大地影响滤波效果。此外，阈值的选取也十分关键，错误的阈值将导致边缘和细节模糊。

根据上述要求，本发明采用图像稀疏表示和基于子块相似度的全局滤波的思想，为医学超声图像滤波提供一种新的方法，可以有效地降低斑点噪声，且保持其细微结构、边缘和纹理特性。

发明内容

本发明的目的是提出一种可以有效降低医学超声图像中斑点噪声的医学超声图像滤波方法。

本发明提出一种的医学超声图像滤波方法，具体步骤为：

1、将二维超声图像分成子块，使用KSVD过完备字典学习算法[5]，对其进行最稀疏分解，求得每个子块的分解系数及灰度均值；

2、采用基于子块相似度的全局滤波思想，利用每个子块的稀疏分解系数，在整幅图像寻找灰度值和结构相似的候选块；

3、在所有候选块中，进一步利用灰度均值矩阵去除那些完全不相关的图像子块。通过最终得到的所有最相似子块，对原图像子块内全部像素点一同滤波。对于重叠区域多次滤波的像素点，对其进行均值处理，从而实现整幅图像的滤波。

下面就本发明方法的各个步骤涉及的相关技术细节作进一步的具体描述。

1、超声图像中的斑点噪声是一种乘性噪声，服从瑞利分布。设噪声图像Y∈R^M ^× ^N ，M×N为图像尺寸，f(•) 和η(•)分别代表无噪声图像和噪声，满足：

Y(i,j)=f(i,j)*η(i,j), i=1,…M, j=1,…,N (1)

η(•)独立于 f(•)，服从瑞利分布：

(2)

其中σ_η 为瑞利分布参数, η(•)的均值为1。

我们使用KSVD字典学习算法[5]，利用过完备的冗余字典对超声图像进行稀疏分解，使得分解后的系数中大部分分量为零，只有少数的非零大系数。这些非零大系数揭示了图像的内在结构与本质属性，用较少的系数捕获感兴趣目标重要信息，从而简化后续的处理任务。

如图1所示，当噪声图像Y∈R^M ^× ^N ，冗余字典D为L个K维原子的集合，经过稀疏分解后的系数为X，则图像Y的稀疏表示为：

(3)

其中T ₀为稀疏度。

首先，对图像Y进行分块。利用b×b的滑动窗口(b∈[5,15]的整数)，从上到下，从左到右遍历整幅图像，依次提取图像子块，并将其转开为b ²×1的列向量。这样，原图像Y分成了L个图像子块y_i ，L=(M-b+1)×(N-b+1)。因此L个图像子块的稀疏表示则为：

(4)

其中为初始化的随机字典。系数矩阵X∈，其第i列x_i 表示第i个图像子块y_i 的稀疏分解系数。

接着固定字典D，采用正交匹配跟踪算法(OMP)找到最优系数X的近似解。然后，利用奇异值分解(SVD)不断更新原子，求得最终的目标字典。设为的第j列，为经过第T次的迭代后X的第j行，假定更新第k个原子，则式(4)变为：

(5)

令记录{y_i }中使用原子的图像子块，即：

(6)

定义矩阵为L×||，在((i),i)处为1，其余全为0，令，，则式(5)改为：

(7)

对进行SVD分解：=U∆V^T ，利用矩阵U的第一列更新原子。经过多次迭代后，对所有原子进行更新，则可得到最终的冗余字典D。

同时，对于所有图像子块，计算其灰度均值，记为均值矩阵Y_mean∈R ¹ ^× ^L ，用于步骤2。

2、常用的超声滤波方法主要是基于局部特征的滤波算法，窗的尺寸和形状极大地影响滤波效果。若窗取得过大，滤波效果好，但边缘模糊；若窗取得小，则边缘保持良好，但滤波效果差。因此对于局部的滤波算法而言，滤波效果和边缘保持是一对矛盾的问题。

这里采用基于子块相似度的全局滤波思想。如图2所示，对于子块p ₁而言，在局部窗S中，没有与其结构相似的子块，若用S中的像素点对其滤波，则会模糊边缘。然而，观察整幅图像，子块p ₂，p ₃，p ₄和p ₅的灰度分布和结构信息与p ₁非常类似。若使用p ₂- p ₅对其滤波，可以在去除噪声的同时，很好地保持边缘。

判断子块相似性的方法有很多种，如均值、方差、灰度共生矩阵[6]和欧氏距离[7]等。这些方法都是基于图像灰度值本身。由于需要在整幅图像中寻找，对于一幅M×N的图像，计算复杂度为((M-b+1)×(N-b+1))²×b ²×(M×N)。计算量太大，不易实现。但是经过图像稀疏分解后，每个子块仅有少数不为零的大系数，利用这些大系数计算子块相似度，可以大大减少计算量，实现整幅图像的全局滤波。

这里根据每个图像子块的分解系数X，在整幅图像内寻找与其灰度值和结构最相似的其他子块。

经过稀疏分解后，原图像Y∈R^M ^× ^N 被分解为系数矩阵X∈。其中X的每一列大部分的系数分量为零，只有少数的非零大系数，这些大系数很好地表示了原图像的特征。

对X的每一列求和，得到每个图像子块的稀疏分量X_sum∈R ¹ ^× ^L 。

(8)

这里，我们将图像子块之间的相似性问题转化其对应稀疏分量之间的欧氏距离，即：

(9)

其中C_ij 为图像子块y_i 与y_j 的相似度。

对于需要滤波的图像子块y_i ，将其稀疏分量X_sum _i 减去其他图像子块的稀疏分量，得到差值矩阵∆(i)。

(10)

对矩阵∆(i)按升序排列，取前t个差值最小的列，将其对应的图像子块作为候选块。

3、进一步利用灰度均值矩阵Y_mean去除t个候选块(t∈[20,120])中那些完全不相关的图像子块。根据前t列，找到其所在图像子块的灰度均值Y_mean (k)，k=1,…t。若，保留该图像块；否则删除。故最终将保留n个最相似子块，即n列，n<t，记为Y_similar∈。

不同于点滤波，这里对整个图像子块内全部像素点一同滤波。考虑到Y_similar中每一列的灰度分布及结构均与y_i 十分相似，则可以直接对中Y_similar每一行取均值，作为滤波后的。这样可以提高算法速度且更好地保持边缘特性。图像子块y_i 滤滤波后的灰度值为：

(11)

将滤波后的列向量进行反列向量处理，即展开成b×b矩阵。由于图像子块在原图像中是全遍历的，对于重叠区域多次滤波的像素点，对其进行均值处理，从而实现整幅图像的滤波。

附图说明

图1 图像的稀疏表示。

图2 图像全局冗余信息。

图3仿真图像：(a)原始图像，(b)噪声图像。

图4 各种不同算法滤波后图像：(a)原始图像，(b)噪声图像，(c)KSVD算法滤波图像，(d)SRAD算法滤波图像，(e)小波算法滤波图像，(f)本发明算法滤波图像。

图5 不同信噪比的噪声图像滤波后的量化参数：(a)SNR值，(b)MSSIM值，(c)FOM值。

图6 对乳腺肿瘤超声图像滤波：(a)原始图像，(b)KSVD算法滤波图像，(c)SRAD算法滤波图像，(d)小波算法滤波图像，(e)本发明算法滤波图像。

图7对淋巴结超声图像滤波：(a)原始图像，(b)KSVD算法滤波图像，(c)SRAD算法滤波图像，(d)小波算法滤波图像，(e)本发明算法滤波图像。

具体实施方式

对本发明提出的滤波方法进行仿真图像和实际超声图像测试。分别采用本文方法、KSVD算法[5]、各向异性扩散滤波算法[3]和小波算法[4]滤波，并对结果进行比较。本文方法中，经过多次参数优化和比较，b=11，t=80时，滤波效果最好。

首先，比较仿真图像的滤波性能。对于一幅无噪声的原始图像，使用Yu等提出的超声模型[3]仿真，对其添加不同程度的噪声。图3(a)为原始图像，图3(b)为有噪声图像。

图4表示对于一幅噪声图像，各种不同算法的滤波效果。(a)为原始图像，(b)为噪声图像，(c)-(e)为分别经过KSVD、SRAD和小波算法滤波后的图像，(f)为本发明的滤波图像。可以得到，小波算法几乎不能去除斑点噪声，KSVD算法的滤波效果优于SRAD算法，但在背景部分仍存在部分噪声。而本发明取得了较为显著的噪声抑制效果，在去除噪声的同时，可以很好地保持边缘特性。

采用信噪比(signal to noise, SNR)、结构相似度(mean structure similarity，MSSIM)[8]和图像佳数(figure of merit，FOM)[3]评价算法滤波质量.

SNR反映滤波器对噪声的抑制作用，SNR越大，滤波效果越好：

(12)

MSSIM衡量两幅图像的结构相似度：

(13)

(14)

其中μ _Y , μ _f , σ _Y , σ _f , σ _Yf 表示Y(·)和f(·)内8*8窗的均值、方差和协方差。C ₁ =K ₁*L，C ₂ =K ₂*L，L 表示像素值的动态范围，对于医学超声图像，L取255，K ₁取0.01，K ₂取0.03。MSSIM为SSIM的窗平均。MSSIM∈[0,1]，值越大表示滤波效果越理想。

FOM表明滤波器边缘保持特性：

(15)

其中N和N _ideal为检测边缘和理想边缘的点数，d_i 表示第i个检测边缘与其最近的理想边缘之间的欧氏距离，e为常数，通常取1/9。FOM∈[0,1]，越接近1表明检测到的边缘越正确。

图5为不同信噪比的噪声图像滤波后的SNR、MSSIM和FOM值。这三个指标，本发明算法最好。随着SNR的提高，KSVD算法也能取得不错的效果。在图像相似性和边缘保持方面，本发明方法有明显优势，特别在SNR低的情况下，本发明仍然可以取得很好的效果，不仅可以去除同质区的噪声，还能较好地保持边缘和各类细节。

将上述四种滤波算法用于实际医学超声图像，如乳腺肿瘤图像和淋巴结图像。

从图6和图7可以看出，小波算法无法去除斑点噪声，KSVD和SRAD算法可以降低噪声的影响，但是图像对比度变差，边缘模糊，影响整体视觉。本发明能在有效抑制背景及感兴趣区域的斑点噪声，更好地保留边缘和纹理信息，改善图像质量。

综上所述，本发明非常适合于医学超声图像的滤波，可以有效地降低斑点噪声，并且较好地保持超声图像中的纹理信息和边缘特性。

参考文献

[1] LEE J S. Speckle analysis and smoothing of synthetic aperture radar images [J]. Computer Graphics and Image Processing, 1987, (17): 24–32.

[2] HUANG T, YANG G, TANG G. A fast two-dimensional median filtering algorithm [J]. IEEE Transaction on Acoustics Speech and Signal Processing, 1979, (27): 13-18.

[3] YU Y J, ACTON S T. Speckle Reducing Anisotropic Diffusion [J]. IEEE Transaction on Image Processing, 2002, (11): 1260-1270.

[4] Rabbani H, Nezafat R, Gazor S. Wavelet-Domain Medical Image Denoising Using Bivariate Laplacian Mixture Model [J]. I IEEE Transaction on Biomedical Engineering, 2009, (12): 2826-2837.

[5] Aharon M, Elad M, Bruckstein A. K-SVD: An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation [J]. IEEE Transaction on Signal Processing, 2006, (11): 4311-4322.

[6] Poonguzhali S, Ravindran G. A complete automatic region growing method for segmentation of masses on ultrasound images [C]. Proceedings of International Conference on Biomedical and Pharmaceutical Engineering, Singapore, 2006: 88-92.

[7] BUADES A, COLL B, MOREL M. A review of image denoising algorithms with a new one [J]. Multiscale Modeling and Simulation, 2005, (4): 490-530.

[8] WANG Z, BOVIK A C, SHEICKH H R et al. Image quality assessment: from error visibility to structural similarity [J]. IEEE Transaction on Image Processing, 2002, (13): 600-612.。

Claims

1.基于稀疏表示的医学超声图像滤波方法，其特征在于具体步骤为：

(1) 将二维超声图像分成子块，使用KSVD过完备字典学习算法，对其进行稀疏分解，求得每个子块的分解系数及灰度均值；

(2) 采用基于子块相似度的全局滤波思想，利用每个子块的稀疏分解系数，在整幅图像寻找灰度值和结构相似的候选块；

(3) 在所有候选块中，进一步利用灰度均值矩阵去除那些完全不相关的图像子块；通过最终得到的所有最相似子块，对原图像子块内全部像素点一同滤波；对于重叠区域多次滤波的像素点，对其进行均值处理，从而实现整幅图像的滤波；

步骤(1)中，所述使用KSVD字典学习算法，利用过完备的冗余字典对超声图像进行稀疏分解，使得分解后的系数中大部分分量为零，只有少数的非零大系数，具体步骤为：

设噪声图像Y∈R^M ^× ^N ，冗余字典D为L个K维原子的集合，经过稀疏分解后的系数为X，则图像Y的稀疏表示为：

(1)

其中T ₀为稀疏度；

首先，对图像Y进行分块，利用b×b的滑动窗口，b∈[5,15]的整数，从上到下，从左到右遍历整幅图像，依次提取图像子块，并将其转开为b ²×1的列向量；这样，原图像Y分成了L个图像子块y_i ，L=(M-b+1)×(N-b+1)，L个图像子块的稀疏表示则为：

i∈{1,2,…,L} (2)

其中为初始化的随机字典；系数矩阵X∈，其第i列x_i 表示第i个图像子块y_i 的稀疏分解系数；

接着固定字典D，采用正交匹配跟踪算法找到最优系数X的近似解；然后，利用奇异值分解不断更新原子，求得最终的目标字典；设为的第j列，为经过第T次的迭代后X的第j行，假定更新第k个原子，则式(2)变为：

(3)

令记录{y_i }中使用原子的图像子块，即：

(4)

定义矩阵为L×||，在((i),i)处为1，其余全为0，令，，则式(3)改为：

(5)

对进行SVD分解：=U∆V^T ，利用矩阵U的第一列更新原子；经过多次迭代后，对所有原子进行更新，则得到最终的冗余字典D；

同时，对于所有图像子块，计算其灰度均值，记为均值矩阵Y_mean∈R ¹ ^× ^L ；

步骤(2)中，所述采用基于子块相似度的全局滤波思想，利用每个子块的稀疏分解系数，在整幅图像寻找灰度值和结构相似的候选块的具体步骤为：

经过稀疏分解后，原图像Y∈R^M ^× ^N 被分解为系数矩阵X∈，其中X的每一列大部分的系数分量为零，只有少数的非零大系数，这些大系数很好地表示了原图像的特征；

对X的每一列求和，得到每个图像子块的稀疏分量X_sum∈R ¹ ^× ^L ：

(6)

这里，将图像子块之间的相似性问题转化其对应稀疏分量之间的欧氏距离，即：

(7)

其中C_ij 为图像子块y_i 与y_j 的相似度；

对于需要滤波的图像子块y_i ，将其稀疏分量X_sum _i 减去其他图像子块的稀疏分量，得到差值矩阵∆(i)：

(8)

对矩阵∆(i)按升序排列，取前t个差值最小的列，将其对应的图像子块作为候选块，t∈[20,120]；

步骤(3)中，所述进一步利用灰度均值矩阵去除候选块中那些完全不相关的图像子块，最终将保留n个最相似子块，n<t ，记为Y_similar∈；不同于点滤波，对整个图像子块内的全部像素点一同滤波；由于Y_similar中每一列的灰度分布及结构均与y_i 十分相似，所以直接对中Y_similar每一行取均值，作为滤波后的；

将滤波后的列向量进行反列向量处理，即展开成b×b矩阵；由于图像子块在原图像中是全遍历的，对于重叠区域多次滤波的像素点，对其进行均值处理，从而实现整幅图像的滤波。