CN110570379A - 一种基于结构张量的非局部均值ct图像降噪方法 - Google Patents

Abstract

本发明涉及一种基于结构张量的非局部均值CT图像降噪方法，属于图像处理领域。该方法包括：S1：输入含噪CT图像，进行归一化处理，求解图像结构张量矩阵；S2：计算结构张量的迹；S3：根据结构张量的迹自适应地确定滤波系数；S4：计算图像块之间像素差值和差值积分图像，并利用积分图像进行加速；S5：根据邻域间相似度确定滤波权值，获得去噪图像。本发明解决了传统非局部均值算法在滤波过程中采用固定滤波系数从而导致的图像细节过度平滑问题，既提高图像降噪质量又提升运算效率。

Description

一种基于结构张量的非局部均值CT图像降噪方法

技术领域

本发明属于图像处理领域，涉及一种基于结构张量的非局部均值CT图像降噪方法。

背景技术

计算机层析成像(Computed Tomography,CT)利用X射线穿过不同物质的衰减信息，采用一定的重建算法得到被测物体内部密度分布，其图像清晰，分辨率高，是世界公认的先进的无损检测手段之一，已广泛应用在航天、航空、医学、生物、工业、农业、电子和考古等领域。快速、高分辨、高质量成像是CT领域需要攻克的重大难题，而噪声是影响成像质量主要因素之一，因此CT图像降噪对提高图像质量有重要现实意义。

2005年，Baudes提出了非局部均值滤波算法(Non-Local Means，NLM)，该算法利用图像中普遍存在的冗余信息来去除噪声，与双线性滤波、中值滤波等利用图像局部信息来滤波的方法不同，它是利用整幅图像信息进行去噪，以图像块为单位在图像中寻找相似区域，再对这些区域求平均来去除图像噪声。NLM算法主要不足在于：滤波系数无法自适应确定，整个降噪过程采用固定的滤波系数，而滤波系数控制滤波权值的衰减速度，是决定算法性能的关键参数，若滤波系数过大，该算法退化为均值滤波；算法计算复杂度高，运算耗时。

因此亟需一种能够自适应确定滤波系数的非均值滤波算法来解决图像降噪的问题。

发明内容

有鉴于此，本发明的目的在于提供一种基于结构张量(Structure Tensor,ST)的非局部均值CT图像降噪方法(以下简称ST-NLM)，用于克服现有技术应用NLM算法降噪时由于采用固定滤波系数而导致的图像细节过度平滑问题，同时引用积分图像对NLM方法进行加速，既提高图像去噪质量又提升运算效率。

为达到上述目的，本发明提供如下技术方案：

一种基于ST-NLM的CT图像降噪方法，先求取图像的结构张量及其结构张量的迹，然后根据迹的大小自适应确定滤波系数，并采用积分图像加速NLM运算效率。该方法具体包括以下步骤：

S1：输入含噪CT图像，进行归一化处理，求解图像结构张量矩阵；

S2：计算结构张量的迹；

S3：根据结构张量的迹自适应确定滤波系数；

S4：计算图像块之间像素差值和差值积分图像，并利用积分图像进行加速；

S5：根据邻域间相似度计算滤波权值，获得去噪图像。

进一步，所述步骤S1具体包括：定义一幅二维图像I，则图像中任一像素点x的初始结构张量定义为点x的梯度向量的外积空间域加权和，记作J(x)：

其中， 表示沿x和y方向的梯度算子；I_x(x)、I_y(x)分别表示沿x、y方向的梯度；G_σ是高斯卷积核，该过程是一个线性卷积过程。

进一步，所述步骤S2具体包括：对结构张量矩阵求特征值，特征值之和即为结构张量的迹。

进一步，所述步骤S3具体包括：建立一个控制滤波系数的权值函数c：

其中，H表示结构张量的迹，T为常数；

当结构张量的迹H较大时，采用较小滤波系数以保护图像的特征；反之，当H较小时，采用较大的滤波系数以去除噪声，由此得到自适应滤波系数函数：

其中，k₀为常数，σ为图像噪声方差。

进一步，所述步骤S4中，引入积分图像求解图像块之间像素差值，以避免逐点求像素差值计算效率低的问题。

进一步，所述步骤S5具体包括：利用欧式距离计算公式求解邻域相似度，欧式距离越小，相似度越大，反之，相似度越小；再根据邻域相似度大小，计算滤波权值，相似度越大，滤波权值越大，反之，滤波权值越小；最后根据滤波权值计算得到去噪图像。

本发明的有益效果在于：本发明针对NLM算法因采用固定滤波系数而导致图像细节过度平滑和计算耗时等问题，利用结构张量的迹表征图像几何结构信息，自适应地确定滤波系数，解决了NLM算法采用固定滤波系数所导致的细节过度平滑问题，同时引用积分图像对NLM算法进行加速，既提高图像降噪质量又提升运算效率。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：

图1为本发明所述图像去噪方法的流程图；

图2为本发明中搜索窗和相似窗几何关系图；

图3为本发明使用的结构张量能够表征图像几何结构示意图；

图4为本发明使用积分图像计算任意矩形的像素值示意图；

图5为本发明实验使用的昆虫纵切CT图像；

图6为分别使用NLM方法和本发明ST-NLM方法对图5进行降噪后的结果图；

图7为图5、图6的局部放大结果图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

请参阅图1～图7，图1为本发明所述基于结构张量的非局部均值CT图像降噪方法的一种优选实施例，具体包括以下步骤：

S1：输入含噪CT图像，进行归一化处理，求解图像结构张量矩阵。具体为：定义一幅二维图像I，则图像中任一像素点x的初始结构张量定义为点x的梯度向量的外积空间域加权和，记作J(x)：

其中， 表示沿x和y方向的梯度算子；I_x(x)、I_y(x)分别表示沿x、y方向的梯度；G_σ是高斯卷积核，该过程是一个线性卷积过程。参照图3，表示一幅图像结构张量分解示意图，图3(a)表示原图，图3(b)表示结构张量分解得到的图像，表明结构张量分解能够准确地表征图像的几何结构。

S2：计算结构张量的迹，具体为：利用步骤S1中所求的结构张量矩阵求其特征值，任意一点处的结构张量的特征值之和即为该点的迹，用H表示。

利用结构张量的迹表征图像几何结构信息，λ₁，λ₂分别表示像素点x的结构张量的两个特征值：

由上述分析可得，在图像的边缘及角点区域，特征值之和λ₁+λ₂值较大；在平坦区域，λ₁+λ₂则较小，而特征值之和即为结构张量的迹；当迹较大时，采用较小滤波系数以保护图像的特征；反之，采用较大的滤波系数以去除噪声。

S3：根据结构张量的迹自适应地确定滤波系数。具体为：建立一个控制滤波系数的权值函数c：

式中T为常数，实验中取0.1，用于控制c的变化趋势。

结构张量的迹H较大时，采用较小滤波系数以保护图像的特征；反之，当H较小时，采用较大的滤波系数以去除噪声，由此得到自适应滤波系数函数：

式中k₀为常数，σ为图像噪声方差。

S4：计算图像块之间像素差值和差值积分图像，利用积分图像进行加速。具体为：对含噪图像I边缘行列进行对称扩展。参照图2，设原始图像大小为N×N，相似窗大小设置为L×L，半窗宽表示为d_s，搜索窗大小设置为K×K，半窗宽表示为D_s，D_s≥d_s。对图像I的上下边缘区域均扩展D_s+d_s行、左右边缘区域均扩展D_s+d_s列，其扩展区域以镜像对称方式填充，记作I₀。在整个图像域I计算块的差值d_ifft(x)，V(x)表示搜索窗内以x为中心像素点的图像块，V(x)，V(x+t)表示平移前后对应图像块：

d_ifft(x)＝||V(x)-V(x+t)||²,x＝(x₁,x₂)∈I (5)

根据积分图像定义，积分图像中任意一点x＝(x₁,x₂)的值等于从图像的左上角到这个点所构成的矩形区域内所有点的灰度值之和，对d_ifft(x)分别做行列叠加构造积分图像S_t(x₁,x₂)：

参照图4，右下角所示区域的像素值可以由下式快速计算得到：

SDist＝S_t(x₁,x₂)+S_t(x₁-t,x₂-t)-S_t(x₁-t,x₂)-S_t(x₁,x₂-t) (8)

推广得到任意矩形区域内的像素值：

SDist＝S_t(x₁+d_s,x₂+d_s)+S_t(x₁-d_s,x₂-d_s)-S_t(x₁+d_s,x₂-d_s)-S_t(x₁-d_s,x₂+d_s) (9)

S5：根据邻域间相似度计算滤波权值，获得去噪图像。

计算相似度：

d(x,y)＝SDist/(L2) (10)

其中，y＝(y₁,y₂),y₁＝[x₁-D_s,x₁+D_s],y₂＝[x₂-D_s,x₂+D_s]。

根据相似度大小计算相似窗内每个像素点的权值：

w(x,y)＝exp(-d(x,y)/(h^2)) (11)

最后，根据滤波权值得到去噪图像。

下面结合附图对本发明ST-NLM方法加以详细说明，应指出的是，所描述的实施例仅旨在便于对本发明方法的理解，而对其不起任何限定作用。

如图5所示，原始的昆虫纵切CT图像噪声较大。

如图6所示，分别使用NLM方法和本发明ST-NLM方法对图5进行降噪，分别用图6(a)、图6(b)表示降噪结果，图6(a)整体上降噪效果明显，但局部细节出现不同程度模糊；而本发明有效降低噪声，同时边缘轮廓保持清晰完整，如图6(b)图所示。

如图7所示，取图5、图6的局部区域放大，从放大结果图7来看，经本发明所述ST-NLM方法降噪之后边缘保持性更好，而NLM方法降噪之后边缘较为模糊，如箭头处所示。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种基于结构张量的非局部均值CT图像降噪方法，其特征在于，该方法具体包括以下步骤：

S1：输入含噪CT图像，进行归一化处理，求解图像结构张量矩阵；

S2：计算结构张量的迹；

S3：根据结构张量的迹自适应确定滤波系数；

S4：计算图像块之间像素差值和差值积分图像，并利用积分图像进行加速；

S5：根据邻域间相似度计算滤波权值，获得去噪图像。

2.根据权利要求1所述的一种基于结构张量的非局部均值CT图像降噪方法，其特征在于，所述步骤S1具体包括：定义一幅二维图像I，则图像中任一像素点x的初始结构张量定义为点x的梯度向量的外积空间域加权和，记作J(x)：

其中， 表示沿x和y方向的梯度算子；I_x(x)、I_y(x)分别表示沿x、y方向的梯度；G_σ是高斯卷积核。

3.根据权利要求1所述的一种基于结构张量的非局部均值CT图像降噪方法，其特征在于，所述步骤S2具体包括：对结构张量矩阵求特征值，特征值之和即为结构张量的迹。

4.根据权利要求1所述的一种基于结构张量的非局部均值CT图像降噪方法，其特征在于，所述步骤S3具体包括：建立一个控制滤波系数的权值函数c：

其中，H表示结构张量的迹，T为常数；

当结构张量的迹H较大时，采用较小滤波系数以保护图像的特征；反之，当H较小时，采用较大的滤波系数以去除噪声，由此得到自适应滤波系数函数：

其中，k₀为常数，σ为图像噪声方差。

5.根据权利要求1所述的一种基于结构张量的非局部均值CT图像降噪方法，其特征在于，所述步骤S5具体包括：利用欧式距离计算公式求解邻域相似度，再根据邻域相似度大小，计算滤波权值，最后根据滤波权值计算得到去噪图像。