CN113313641A - 一种自适应中值滤波的ct图像去噪方法 - Google Patents

Abstract

本发明公开了一种自适应中值滤波的CT图像去噪方法，属于医学图像处理技术领域，特别适合于新冠肺炎的CT图像去噪。选取尺寸为n×n的方形滤波窗口，比较窗口内灰度值的自适应最大值

自适应最小值

与当前像素点灰度值f(i,j)，依据第一阈值T₀判断当前像素点是否为疑似噪声点，若是，再根据第二阈值T₁进一步精确判断是否为噪声点；若当前像素点不为疑似噪声点或噪声点，则遍历窗口内的下一像素点；通过中心加权中值滤波的方法处理噪声点；最后输出中值滤波去噪后的CT图像。在保持图像去噪的同时更好地保护了图像细节；通过改进熵权法修正传统熵权法传递信息存在偏差的问题，通过各项评价指标对去噪效果的贡献值确定中心加权滤波的最优权重，从而达到最优去噪效果。

Description

一种自适应中值滤波的CT图像去噪方法

技术领域

本发明涉及图像信号处理技术领域，具体涉及一种基于自适应中值滤波的CT图像去噪方法，适合于肺炎CT图像的滤波处理，特别适合于新冠肺炎的CT图像去噪，属于医学图像处理技术领域。

背景技术

新冠肺炎的病灶特点主要表现为各种形态的磨玻璃影，或者伴有实变影，CT图像在传输和获取过程中会受到脉冲噪声的干扰。新冠早期病灶变化不明显，病灶数目较少，病灶范围较小，密度较低，脉冲噪声干扰容易造成早期新冠患者漏诊的问题。

针对传统CT图像脉冲噪声去噪算法中值滤波器易受到窗口尺寸的影响，当窗口尺寸过小时，中值滤波容易受到周围像素点的干扰，从而使原来信号点的灰度值被替换成噪声点的灰度值，使得图像去噪能力受到限制，而无法正确的去除噪声。当窗口尺寸过大时，噪声点的干扰能力会大幅降低，能在很大程度上提高图像去噪能力，但会导致图像细节(例如边缘，线条，角等)遭到破坏。因此，传统中值滤波很难在图像去噪和保留图像细节中有着较好的表现。

发明内容

针对这些上述缺点，本发明提出多级阈值中心加权自适应中值滤波算法，通过采用多级阈值方法，更大程度上降低了对信号点的误判，在保持图像去噪的同时更好地保护了图像细节。通过对图像中心像素点赋予权重，权重的大小受窗口尺寸的影响，当该点为信号点时，通过权重的性质，极可能使该点排在中间位置，相反，当该点为噪声点时，会排在较后位置，通过修正像素权值后，中值运算在一定程度上加强了非噪声像素对噪声修复的主导作用。随着中值滤波迭代次数的不断增多，图像细节信息会逐渐减小，因此只有设置一个终止迭代条件，才能得到最佳滤波效果，拟提出变步长凸组合范数的去噪迭代终止准则，在减小图像去噪失真的情况下，加快原算法的收敛速度，减小中值滤波的计算量。

本发明根据CT图像在获取时易产生脉冲噪声的特点，以及为了能够适应新冠早期病灶变化不明显，病灶数目较少，密度较低的问题，提出了一种基于多级阈值中心加权自适应中值滤波的新冠肺炎CT图像去噪方法，包括以下步骤：

S10，选取尺寸为n×n的方形滤波窗口，比较窗口内灰度值的自适应最大值

自适应最小值

与当前像素点灰度值f(i,j)，依据第一阈值T₀，若

或

则初步判断该点为疑似噪声点；若不为疑似噪声点，则遍历窗口内的下一像素点；

S20，通过中值滤波方法处理所述疑似噪声点；

S30，输出中值滤波去噪后的CT图像。

其中，关于第一阈值T₀的选择，因第一阈值是噪声点检测的目的是为了找出疑似的噪声点。传统手段是基于极值判断噪声点，但一些图像存在亮度突变或颜色差异等，这种因素会导致像素灰度值的差异变大，若采用极值去判断噪声点，会将信号点误判为噪声点，故该方法不能成为噪声点的决定性判断依据。若噪声像素值趋于极值，可以先判断当前像素灰度值是否接近极值，需要去与该像素点周围的像素点灰度值去比较，计算每一个像素点灰度值与灰度自适应最大值和灰度自适应最小值、灰度中值、灰度平均值的绝对平均差。采用该平均差来作为阈值T₀，固定阈值无法适应图像的亮度突变或颜色差异等情况，采取使用随窗口多均值大小变化的可变动态阈值更好解决此类难点，因此，T₀通过下式来计算：

式中，n为滤波窗口大小相等的长和宽，

和

分别代表滤波窗口内像素点灰度值的自适应最大值和自适应最小值，f(i,j)为当前像素点灰度值。f_med为当前滤波窗口内的所有像素点的灰度中值，f_ave为当前滤波窗口内的所有像素点的灰度平均值。

传统的灰度最大最小值，用于检测噪声像素的像素值范围很大，会增大对噪声的错误检测概率。一般情况下，像素集合的均值通常与大部分像素值接近。所以可以通过分析一个像素与其所在集合的均值之间关系来判断该像素值是否被外界噪声干扰，通过9个均值来选取自适应最大最小值。

和

选取规则：为选取自适应最大最小灰度值，先对n×n滤波窗口中的所有像素值求取均值：

分别对待处理像素f(i,j)水平及垂直方向的像素值求取均值，得到2个均值：

分别在点(i，j)的上、下两个方向定义5×2的两个邻域矩形窗口，在点(i，j)的左、右两个方向定义2×5的两个邻域矩形窗口，示意如图6所示，同时定义这四个邻域所对应的灰度平均值分别为

其值计算公式为：

对n×n滤波窗口主副对角线上的像素值求取均值，得到2个均值：

以所求得的9个均值中的最大值作为滤波窗口的自适应最大值

以所求得的9个均值中的最小值来为滤波窗口的自适应最小值

即：

为了提高噪声点判断的精度，在步骤S10中初步判断该点为疑似噪声点之后，还包括通过第二阈值

对所述疑似噪声点的进一步精确判断，即当满足|f(i,j)-f_ave|-T₁＞0时，精确判断该点为噪声点；若|f(i,j)-f_ave|-T₁≤0，则继续针对下一疑似噪声点进行第二阈值精确判断；其中，

为阈值T₁的调节因子，

噪声密度

N_pn表示所述窗口内噪声点的个数。

进一步，步骤S20的中值滤波方法是滤波窗口内加权中值滤波，即用加权中值灰度值f(i,j)_CWM替换当前像素点灰度值，计算式为：

其中，median{}表示在一组数值中，返回居于中间数值的函数，k(r,s)×f(r,s)表示f(r,s)的权重是k(r,s)时灰度值，

是以像素点(i,j)为中心的n×n滤波窗口空心邻域内所有像素点的灰度值集合，f(r,s)表示

内所有像素点对应的灰度值，(r,s)表示

内任意一个像素点，即以像素点(i,j)为中心的n×n滤波窗口内除(i,j)外任意一个像素点。

f(r,s)的权重k(r,s)值计算：滤波窗口内，大部分的点应当与像素均值接近，与像素均值接近的灰度点的值，应当赋予像素加权较大的值；若像素灰度值远远大于像素均值，则可能为噪声点或者图像边缘，这时应当赋予像素加权较小的值。所以加权函数与灰度距离相关性公式呈负相关。设当前像素为(r,s)，当前滤波窗口的灰度均值为f_ave，根据当前像素f(r,s)与灰度均值f_ave的距离范数，赋予f(r,s)加权系数，如(14)所示，改进激活熵最优匹配混合L₁/L₂/L_∞范数灰度距离相关性公式为：

式中，

为像素点f(r,s)的L₁范数权重，

为像素点f(r,s)的L₂范数权重，

为像素点f(r,s)的L_∞范数权重，f_ave为当前滤波窗口内所有像素点的灰度平均值。

设

p＝1,2,∞；像素点(r,s)指标的p(p＝1,2,∞)范数信息熵为

式中：

n代表滤波窗口大小相等的长和宽。当

时，

计算激活熵：激活熵是信息熵通过激活函数计算得到的，激活函数主要用于人工神经网络，本文使用改进的双Sigmoid-Tan函数作为激活熵的算法，公式为：

式中：β是控制函数取值范围的常量，β∈[0,1]；α是控制函数形状的常量，α∈[-10,10]；h是控制步长变化速度的常量，h∈[0,5]。缩放系数T₂和位移参数θ决定了单个神经元的响应特性，T₂∈[0,3],θ∈[0,6]。在噪声比较大的环境中，评价指标变化率在算法收敛阶段仍然变化很大，较大的变化率在围绕最佳权系数值时会产生较大的波动。因此,当系统工作于低信噪比环境中，应用单Sigmoid算法及传统信息熵仍然会产生比较大的稳态误差。采用双Sigmoid-正切函数，相比Sigmoid函数及传统信息熵，误差相同时，该算法能够获取更快的最佳分解层数，使评价指标变化率在算法即将收敛阶段缓慢变化，从而降低稳态误差。

由于传统熵权法在计算

(5％范围内)的指标权重时，其微小的变化将引起熵权成倍数变化，因此由传统熵权法求出的指标权重存在较大误差。为改进这一缺陷，可利用改进后的熵权计算公式，即

式中

其中

为像素点(r,s)的L_p(p＝1,2,∞)范数权重，

ε＝0.05，

是所有位于[0 1-ε]的激活熵值的平均值，改进的权重公式减小权重误差。对于

均为0.005的f(r,s)像素点，根据实际空间距离(类比灰度距离)与人眼视觉感知(类比加权函数)是一种指数关系，据此关系得到如下加权函数：

式中，c为距离常数，其与滤波窗口的大小有关，

对于

非均为0.005的f(r,s)像素点，根据邻域像素之间的相关性随距离的增大而减小的特性，对邻域中的信号像素分别赋予不同的加权系数，取加权中值以滤除噪声。设当前像素点为(r,s)，离当前像素点最近的噪声点像素为(b₁,b₂)，根据待处理的噪声点像素(r,s)与(b₁,b₂)的距离，赋予f(r,s)加权系数，结合加权距离以及权值归一化，得到如下加权函数，保证负相关性，满足人眼视觉要求。设定加权函数为：

式中，当前像素点为(r,s)，离当前像素点最近的噪声点为(b₁,b₂)，c为距离常数，其与滤波窗口的大小有关，

为了进一步提升中值滤波的性能，在步骤S30输出中值滤波去噪后的CT图像之前，还包括：遍历下一像素点，采用相对信噪比(RSNR)变化率凸组合范数作为迭代终止准则，其定义如下：

式中：M是图像的行数，N是图像的列数，f(i,j)_z和f(i,j)_z+1分别是图像经过第z次和第z+1次迭代滤波后像素点(i,j)的灰度值大小。拟提出变步长凸组合范数的去噪迭代终止准则，公式如(19)所示：

当满足迭代终止准则时，即上式成立时，迭代过程终止，输出中值滤波去噪后的CT图像，式中的ε₂＝0.001是预定义的门限值。如果上式不满足，则继续寻找疑似脉冲噪声点和中值滤波处理。

公式(19)中：λ_z∈[0,1]，λ_z采用基于参数因子的改进反正切函数，即

其中，A＝2，

z＝0时，图像未进行滤波，a(z)＝0，此时图像为原始未经去噪的图像；z＝1，图像经过第一次滤波。迭代步长u(z)采用基于tanh融合Sigmoid函数的可变函数，即u(z)＝β₂(1+exp(-φtanh[E(z)]²))，取β₂＝0.6，φ＝0.08；

相比于现有技术，本方法有如下有益效果：

1、本方法使用多级阈值中心加权自适应中值滤波算法，通过采用多级阈值方法，更大程度上降低了对信号点的误判，在保持图像去噪的同时更好地保护了图像细节。通过对图像中心像素点赋予权重，权重的大小受窗口尺寸的影响，当该点为信号点时，通过权重的性质，极可能使该点排在中间位置，相反，当该点为噪声点时，会排在较后位置，通过修正像素权值后，中值运算在一定程度上加强了非噪声像素对噪声修复的主导作用。

2、本方法使用复合激活熵改进熵权法将像素点灰度值与平均值的灰度差L₁/L₂/L_∞范数做加权线性组合，使用改进的双Sigmoid-Tan函数激活熵替代经典熵权法中的信息熵，可以在不改变各指标信息嫡相对大小及各指标相对重要程度的情况下，降低部分异变样本数据对最终评价结果的影响。通过改进熵权法修正传统熵权法传递信息存在偏差的问题，可合理确定各项评价指标对去噪效果的贡献值，确定中心加权滤波的最优权重，从而达到最优去噪效果。

3、随着中值滤波迭代次数的不断增多，图像细节信息会逐渐减小，因此只有设置一个终止迭代条件，才能得到最佳滤波效果。利用相对信噪比变化率作为迭代终止条件，为了控制迭代次数，引入图像范数的概念，由于多次迭代后噪声已被去除，所以RSNR的变化率范数会很小，这时CT图像为最佳的滤波效果，减小中值滤波的计算量。变步长凸组合在保证系统的稳态性能、减小图像失真的情况下，加快原算法的收敛速度。其核心思想在于：将其中一个评价指标设点为快速动作，以保证算法快速收敛；而另一个则设定为慢速评价指标，保证收敛最终得到很小的稳态误差。

附图说明

图1为本发明的多级阈值中心加权自适应中值滤波流程图；

图2为新冠肺炎胸部CT原图像；

图3为加入脉冲噪声后的新冠肺炎胸部CT图像；

图4为本发明的多级阈值中心加权自适应中值滤波去噪后的胸部CT图像；

图5为传统中值滤波去噪后的胸部CT图像；

图6为像素点(i,j)上、下、左、右四个邻域矩形窗口示意图。

具体实施方式

下面根据附图和实例对本发明进行详细说明，但本发明的具体实施方式不仅于此。

本实施例阐述了本发明应用于新冠肺炎胸部CT图像脉冲噪声的去噪流程。下面结合附图对本发明的实施方式做具体说明。

图1为本发明算法的总体流程图，一种自适应中值滤波的CT图像去噪方法，具体包括：

中值滤波选取尺寸为5×5的方形滤波窗口，根据判别条件

和

来判断脉冲噪声的疑似噪声点。第一阈值T₀的计算原则为公式(1)，灰度自适应最大值最小值

和

的计算原则依据公式(2)～(12)，n＝5。遍历窗口内所有像素点，对像素点求和算平均，来计算滤波窗口内所有像素灰度值平均值f_ave。根据公式

计算第二阈值T₁，其中，

为第二阈值T₁的调节因子，

噪声密度

N_pn表示所述窗口内噪声点的个数。根据判断公式|f(i,j)-f_ave|-T₁＞0，来对该点进行中值替换。

对滤波窗口内所有像素点进行中心加权计算，加权计算原则为公式(13)。对于

非均为0.005的(r,s)像素点，滤波窗口像素点权值计算原则为公式(17)；对于

均为0.005的(r,s)像素点，加权函数为

对滤波窗口内进行加权计算后的所有像素点进行排序，排序后的中间值作为中值，来代替满足第二阈值判断条件的像素点。c为距离常数，其与滤波窗口的大小有关，可按式

来获得，其值为

进一步，中值滤波的加权函数中灰度距离相关性公式为公式(14)，采用改进激活熵加权混合L₁/L₂/L_∞最优范数匹配形式。其中范数权重选择计算利用改进后的熵权根据公式(16)来计算，公式(16)中的激活熵采用公式(15)的改进双Sigmoid-Tan函数作为激活熵，本方法中的β＝0.4，α＝4，h＝2，θ＝3，T₂＝0.8，n＝5。

遍历下一像素点，并检查是否满足终止条件，采用相对信噪比(RSNR)变化率凸组合范数作为迭代终止准则，计算公式采用(18)、(19)所示；若满足中值滤波迭代终止准则公式(19)所示，则输出改进中值滤波去噪后的新冠肺炎胸部CT图像。如果公式(19)不满足，则继续寻找疑似噪声点。

具体到本实施例，本发明采用Matlab R2018a软件平台，对算法进行编程。图2为新冠肺炎胸部CT图像，箭头指向的红色圆圈部位是新冠肺炎的磨玻璃影病灶；图3为加入噪声(噪声强度为0.3的脉冲噪声)后的新冠肺炎胸部CT图像，图片大小均为512*512，滤波窗选用为5*5的方形结构元。使用多级阈值中心加权自适应中值滤波对带有脉冲噪声的新冠肺炎胸部CT图像进行去噪，得到去噪后的的新冠肺炎胸部CT图像，如图4所示。

为了进一步验证本发明的新冠肺炎CT图像去噪方法去噪效果，从客观数据比较分析所提出的方法性能，使用均方误差(mean squared error，MSE)、峰值信噪比(peaksignal-to-noise ratio，PSNR)和信噪比(SNR)对去噪后的图像进行数值计算，均方误差数值越小说明去噪图像质量越好；峰值信噪比、信噪比数值越大，说明去噪效果越好。均方误差的表达式为：

峰值信噪比的表达式为：

信噪比的表达式为：

其中f(j,l)表示去噪后的信号，

代表带脉冲噪声的输入信号，M、N表示输入图像的长度和宽度，图片大小均为512*512。

采用传统中值滤波做对比去噪仿真实验，进一步比较分析不同去噪方法下的MSE、PSNR和SNR，图5为传统中值滤波去噪后新冠肺炎CT图像，图4为本方法去噪后新冠肺炎CT图像，评定指标数值如表1所示。

表1本方法与传统方法去噪后新冠肺炎CT图像性能对比

去噪方法	MSE	PSNR/dB	SNR/dB
				含噪声图像	580.5461	13.0473	9.4051
传统中值滤波方法	166.1739	31.3770	31.5043
				本专利方法	67.4723	39.7349	43.6808

从表中数据的变化趋势可知，与对比实验去噪方法相比，当前已改进的去噪方法下加噪图像的PSNR增加了约8dB，MSE值大大降低，SNR增加了约12dB。本专利改进的CT图像去噪方法下的峰值信噪比值最大值为39.7349，与传统去噪方法下的PSNR相比，增加了约8dB，由此可知，采用本专利方法下的去噪效果是最佳的。

至此，完成了本实施例自适应中值滤波的新冠肺炎CT图像去噪过程。

Claims (10)

1.一种自适应中值滤波的CT图像去噪方法，其特征在于，对所述CT图像，选取尺寸为n×n的方形滤波窗口，比较窗口内灰度值的自适应最大值

自适应最小值

与当前像素点灰度值f(i,j)，依据第一阈值T₀，若

或

则初步判断当前像素点为疑似噪声点；若当前像素点不为疑似噪声点，遍历窗口内的下一像素点；通过中值滤波方法处理所述疑似噪声点；输出中值滤波去噪后的CT图像；所述第一阈值T₀：

式中，f_med为当前滤波窗口内的所有像素点的灰度中值，f_ave为当前滤波窗口内的所有像素点的灰度平均值；

选取自9个均值f_ave,

中的最大值，

选取自9个均值f_ave,

中的最小值，其中，当前滤波窗口内的所有像素点的灰度平均值

当前像素点(i,j)在水平方向像素值的均值

当前像素点(i,j)在垂直方向像素值的均值

当前滤波窗口主对角线上像素值的均值

当前滤波窗口副对角线上像素值的均值

在当前像素点(i,j)上方5×2的邻域矩形窗口内像素灰度平均值

在当前像素点(i,j)下方5×2的邻域矩形窗口内像素灰度平均值

在当前像素点(i,j)左方2×5的邻域矩形窗口内像素灰度平均值

在当前像素点(i,j)右方2×5的邻域矩形窗口内的像素灰度平均值

2.如权利要求1所述的CT图像去噪方法，其特征在于，在所述初步判断该点为疑似噪声点步骤之后，还包括通过第二阈值

对所述疑似噪声点的进一步精确判断，即当满足|f(i,j)-f_ave|-T₁＞0时，精确判断该点为噪声点；若|f(i,j)-f_ave|-T₁≤0，则继续针对下一疑似噪声点进行第二阈值精确判断；其中，

为阈值T₁的调节因子，

噪声密度

N_pn表示所述窗口内噪声点的个数。

3.如权利要求1所述的CT图像去噪方法，其特征在于，所述中值滤波方法为滤波窗口内加权中值滤波，即用加权中值灰度值

替换当前像素点灰度值；其中，median{}表示在一组数值中，返回居于中间数值的函数，k(r,s)×f(r,s)表示f(r,s)的权重是k(r,s)时灰度值，

是以像素点(i,j)为中心的n×n滤波窗口空心邻域内所有像素点的灰度值集合，f(r,s)是

内所有像素点对应的灰度值，(r,s)是以像素点(i,j)为中心的n×n滤波窗口内除(i,j)外任意一个像素点。

4.如权利要求3所述的CT图像去噪方法，其特征在于，所述权重

式中，距离常数

Q(r,s)依据不同的范数情况取不同的距离范数，即：当

且

时，Q(r,s)＝||f(r,s)-f_ave||₁，|| ||₁为L₁范数；当

且

时，Q(r,s)＝||f(r,s)-f_ave||₂，|| ||₂为L₂范数；当

且

时，Q(r,s)＝||f(r,s)-f_ave||_∞，|| ||_∞为L_∞范数；其中，

为像素点f(r,s)的L₁范数权重，

为像素点f(r,s)的L₂范数权重，

为像素点f(r,s)的L_∞范数权重，

5.如权利要求3所述的CT图像去噪方法，其特征在于，所述权重

式中，当前像素点为(r,s)，离当前像素点最近的噪声点像素为(b₁,b₂)，距离常数

Q(r,s)依据不同的范数情况取不同的距离范数，即：当

且

时，Q(r,s)＝||f(r,s)-f_ave||₁；当

且

时，Q(r,s)＝||f(r,s)-f_ave||₂；当

且

时，Q(r,s)＝||f(r,s)-f_ave||_∞；其中，

为像素点(r,s)的L₁范数权重，

为像素点(r,s)的L₂范数权重，

为像素点(r,s)的L_∞范数权重，

和

中至少一个不等于0.005。

6.如权利要求4或5所述的CT图像去噪方法，其特征在于，所述范数权重

的计算式为

其中，误差因子ε＝0.05，参数p＝1,2,∞，

是所有位于[0 1-ε]的激活熵值的平均值，激活熵

α是控制函数形状的常量，β是控制函数取值范围的常量，h是控制步长变化速度的常量，T₂是缩放系数，θ是位移参数；范数信息熵为

式中：

k₁＝1/ln(n×n×3)，

当

时，

7.如权利要求1所述的CT图像去噪方法，其特征在于，在所述输出中值滤波去噪后的CT图像之前，还包括以下步骤：若满足

则输出中值滤波去噪后的CT图像，否则，重新进行疑似噪声点的判断和中值滤波处理，其中，

ε₂是预定义的门限值，

M是图像的行数，N是图像的列数，f(i,j)_z和f(i,j)_z+1分别是图像经过第z次和第z+1次迭代滤波后像素点(i,j)的灰度值大小，A＝2，

z＝0时，图像未进行滤波，a(z)＝0，此时图像为原始未经去噪的图像；z＝1时，图像经过第一次滤波；迭代步长u(z)＝β₂(1+exp(-φtanh[E(z)]²))，

β₂、φ和ε₂为调节参数。

8.如权利要求1所述的CT图像去噪方法，其特征在于，n＝5。

9.如权利要求6所述的CT图像去噪方法，其特征在于，α∈[-10,10]，β∈[0,1]，h∈[0,5]。

10.如权利要求7所述的CT图像去噪方法，其特征在于，β₂＝0.6，φ＝0.08，ε₂＝0.001。

Images