CN113129235A - 一种医学图像噪声抑制算法 - Google Patents

Abstract

本发明提供了一种医学图像噪声抑制算法，对输入图像P₀进行3维高斯滤波得到P₁，利用差分图像D₀＝P₀‑P₁,对噪声标准差进行估计，得到σ，对图像P₀求梯度模图像T₁，再利用阈值函数H(T₁)进行变换得到图像H₀，利用图像H₀和噪声标准差σ，对图像P₀进行自适应NLM滤波得到滤波结果P₂。

Description

一种医学图像噪声抑制算法

技术领域

本发明涉及医学成像技术领域，尤其涉及一种医学数字化X线图像噪声抑制算法。

背景技术

近二十年来，X射线计算机断层成像技术迅猛发展，医疗领域也对图像的质量提高了要求。噪声是医疗图像不可避免的。现在的噪声抑制算法主要分为4类：第一类：空间域滤波，主要代表有均值滤波，中值滤波，高斯滤波，双边滤波，NLM，BM3D等，第二类：变换域滤波，主要有小波阈值去噪，傅里叶变换，第三类：全变分去噪以及偏微分方程去噪。第四类：深度学习。具体来说，由于第四类耗时较长，第三类及第二类去噪效果对参数选择敏感，因而我们选择第一类，鉴于其去噪效果以及运行时间，我们选择NLM去噪算法。

现有的图像去噪技术大多被各大公司视为商业机密，不轻易外泄。现有的公开技术多数只给出二维的去噪结果，对三维数据处理有借鉴意义。大部分现有技术对2维数据处理。在二维到三维的过程中，表面只增加了一维，但是，实际上，增加一维给传统的2维去噪带来了挑战，往往图像在2维看起来噪声抑制效果好，在3维空间上可能会出现新的伪影。

发明内容

针对背景技术中存在的问题，本发明致力于研究3维数据(但不仅限于3维，2维同样适用)的去噪算法，本发明的技术方案如下：

一种医学图像噪声抑制算法，包括如下步骤：

S1、对输入图像P₀进行3维高斯滤波得到P₁；

S2、利用差分图像D₀＝P₀-P₁,对噪声标准差进行估计，得到σ；

S3、对图像P₀求梯度模图像T₁，再利用阈值函数H(T₁)进行变换得到图像H₀；

S4、利用图像H₀和噪声标准差σ，对图像P₀进行自适应NLM滤波得到滤波结果P₂。

进一步地，所述S1具体包括：把CT图像加载至内存得到P₀，之后对P₀进行3维高斯滤波得到P₁，

其中，

为为波半径向量，需人为为置，

T表示转置；

滤波核定义如下：

进一步地，所述步骤S2具体包括：由差分图像D₀＝P₀-P₁，对噪声标准差进行估计得到σ，

D₀(i，j，k)＝P₀(i，j，k)-P₁(i，j，k)

N为图像的像素总数，

进一步地，所述步骤S3具体包括：对P₁求梯度模得到图像T₁,对T₁进行阈值变换H(T₁)得到H₀，

当下标越界时，利用镜像对称拓展。

进一步地，所述阈值转换函数H(x)还可为如下形式:

H(x)＝1+A log(1+B(1-e^-kt))+Ct

其中A,B,C,k均为认为设置的常数。

进一步地，所述阈值转换函数H(x)还可为如下形式:

参数k>0A控制最放大倍数，B控制增长速率变化。

参数选择原则如下：

为空气区域滤波半径放大倍数，进factor₁为边缘区域的最大滤波半径放大倍数，

所述空气区域为图像平坦区域，所述边缘区域为图像梯度模大的区域，

k控制阈值，图像梯度模大于多少时，需要缩小滤波半径，

B控制滤波半径放缩速率。

一步地，所述步骤S4具体包括：结合H₀和σ对图像P₀进行3维NLM自适应滤波得到P₂，

其中p,q∈N³,代表图像像素的三维位置向量，

p＝(p₁,p₂,p₃),q＝(q₁,q₂,q₃)，

U(p,r):{q∈N³|||p-q||＝max(|p₁-q₁|,|p₂-q₂|,|p₃-q₃|)≤r},

r为搜索半径，f为邻域半径，h为滤波半径，

r,f一般为向量形式，代表实际的图像在三个维度上相关性的不一致，此时，相应公式对应修改。

本发明的医学图像噪声抑制算法存在以下有益效果：

1.充分考虑医学3维图像的3维信息，充分考虑每一维度数据间的相关性及其差异。

2.参数设置依据医学图像的空间分辨率设置，可以最大限度保持图像的空间分辨率不下降。

3.自适应过程考虑图像的梯度信息，在保持图像边缘这方面较一般算法好。

附图说明

图1为本发明医学图像噪声抑制算法流程示意图；

图2为本发明医学图像噪声抑制算法中阈值变换函数的图像；

图3为本发明医学图像噪声抑制算法中第二种形式的阈值变换函数的图像；

图4为本发明医学图像噪声抑制算法中第三种形式的阈值变换函数的图像；

图5(1)、图5(2)为本发明医学图像噪声抑制算法在2维图像下去噪结果效果图；

图6(1)、图6(2)为本发明医学图像噪声抑制算法在3维图像下去噪结果效果图。

具体实施方式

下面结合附图1至6和实施方式对本发明作进一步说明，图1为本发明的医学图像噪声抑制算法流程示意图，图2至图4为本发明的医学图像噪声抑制算法中使用到的三种阈值变换函数的图像，下面例举实施例进行详细说明：

实施例一

一种医学图像噪声抑制算法，应用于三维医疗图像降噪技术(也适用于二维图像)，依次包括如下步骤：

S1、对输入图像P₀进行3维高斯滤波得到P₁；

S2、利用差分图像D₀＝P₀-P₁,对噪声标准差进行估计，得到σ；

S3、对图像P₀求梯度模图像T₁，再利用阈值函数H(T₁)进行变换得到图像H₀；

S4、利用图像H₀和噪声标准差σ，对图像P₀进行自适应NLM滤波得到滤波结果P₂。

步骤S1具体包括：把CT图像加载至内存得到P₀，之后对P₀进行3维高斯滤波得到P₁，

其中，

为为波半径向量，需人为为置，

T表示转置；

滤波核定义如下：

步骤S2具体包括：由差分图像D₀＝P₀-P₁，对噪声标准差进行估计得到σ，D₀(i，j，k)＝P₀(i，j，k)-P₁(i，j，k)

N为图像的像素总数，

步骤S3具体包括：对P₁求梯度模得到图像T₁,对T₁进行阈值变换H(T₁)得到H₀，

当下标越界时，利用镜像对称拓展。如P₁(-1，j)＝P₁(1，j)

H(x)＝A(1-e^-kt)+1

其中A,k均为控制参数，A控制滤波半径放大半径，

k控制梯度模平方阈值,均需要根据图像T1确定。

本实施例中，此时可使用第一种阈值变换函数的图像，参阅附图2医学图像噪声抑制算法中阈值变换函数的图像，前述步骤三中，公式中A控制滤波半径的放缩倍数(真实滤波半径放缩倍数为

k代表梯度模方的阈值倒数，若某个像素的梯度模大于该阈值

则表示该像素为边缘点或者角点，需要缩小滤波半径以保护该点的几何特性不被平滑掉。

在此步骤三中，阈值转换函数H(x)还可为如下形式:

H(x)＝1+A log(1+B(1-e^-kt))+Ct

其中A,B,C,k均为认为设置的常数。

该函数图像如附图3医学图像噪声抑制算法中第二种形式的阈值变换函数的图像所示，这些阈值转换函数H(x)都可以通过调节参数改变函数的形态。上述函数都具有如下性质：

第一，在x＝0附近，H(x)会把函数值快速抬升,从而可以保证在梯度不为0的地方尽可能把细节保留，

第二，在梯度值很大后，函数值趋于稳定(如第一个转换函数)或者增长缓慢(如第二个转换函数，可以调节C的值控制函数增长率，当C＝0时函数与第一个转换函数变化趋势相同)，保证在梯度大的像素点滤波半径不至于太小。

第三，可以通过控制参数调节阈值。两个转换函数中指数部分的k因子都是用于控制阈值的。

第四，可以通过调节参数控制函数增长趋势。第一个函数中A是用于控制最大放缩因子的，第二个函数中A,B也是用于控制最大放缩因子的。其中第二个转换函数中B也可以用于控制在x＝0附近的抬升速率。

阈值转换函数必须在x＝0处迅速增大，在增大到一定程度后趋于平缓。

本实施例中，步骤三还可以使用第三种形式的阈值变换函数的图像，且该阈值变换函数的图像如附图4所示，形式如下：

参数k>0A控制最放大倍数，B控制增长速率变化。

参数选择原则如下：

为空气区域(图像平坦区域)的滤波半径放大倍数

factor₁为边缘区域(图像梯度模大的区域)的最大滤波半径放大倍数

k控制阈值，图像梯度模大于多少时，需要缩小滤波半径，

B控制滤波半径放缩速率。

此步骤三中，三种形式的函数表达已示出，具体的参数数值在附图2至附图4中已定义清晰，横轴为X数值，纵轴为对应函数形式的H(x)的数值，参数选择已经在上面第三个函数定义里说明，其它两个函数类似。其中附图3中，曲线为函数曲线，直线则示出了函数的渐近线以清晰得出相关斜率。这些函数都可以作为转换函数使用。在这里我们使用了第一个转换函数。这里仅举例说明部分转换函数，其它具有类似(开关)性质的转换函数均可。

步骤S4具体包括：

结合H₀和σ对图像P₀进行3维NLM自适应滤波得到P₂，

其中p,q∈N³,代表图像像素的三维位置向量，p＝(p₁,p₂,p₃),q＝(q₁,q₂,q₃)，

r为搜索半径，f为邻域半径，h为滤波半径。

r，f一般为向量形式，代表实际的图像在三个维度上相关性的不一致。

此时，相应公式可对应修改。

本实施例中，结合医疗图像的特性，我们还可以把r设置为空间分辨率，f＝r,h＝factor1*σ,factor1一般取为0.1～10。

实验结果展示

对于2维和3维图像数据，本发明算法均可适用。下面展示在2维图像情况下，附图本文算法的去噪结果：(参阅附图5(1)、附图5(2)本发明医学图像噪声抑制算法在2维图像下去噪结果效果图)，供本领域技术人员参阅。

在应用于3维图像数据时，参阅附图6(1)、附图6(2)为本发明医学图像噪声抑制算法在3维图像下去噪结果效果图，供本领域技术人员参阅。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (7)

1.一种医学图像噪声抑制算法，其特征在于，包括如下步骤：

S1、对输入图像P₀进行3维高斯滤波得到P₁；

S2、利用差分图像D₀＝P₀-P₁,对噪声标准差进行估计，得到σ；

S3、对图像P₀求梯度模图像T₁，再利用阈值函数H(T₁)进行变换得到图像H₀；

S4、利用图像H₀和噪声标准差σ，对图像P₀进行自适应NLM滤波得到滤波结果P₂。

2.根据权利要求1所述的医学图像噪声抑制算法，其特征在于，所述步骤S1具体包括：把CT图像加载至内存得到P₀，之后对P₀进行3维高斯滤波得到P₁，

其中，

为为波半径向量，需人为为置，

T表示转置；

滤波核定义如下：

3.根据权利要求1所述的医学图像噪声抑制算法，其特征在于，所述步骤S2具体包括：由差分图像D₀＝P₀-P₁，对噪声标准差进行估计得到σ，D₀(i，j，k)＝P₀(i，j，k)-P₁(i，j，k)，

N为图像的像素总数，

4.根据权利要求1所述的医学图像噪声抑制算法，其特征在于，所述步骤S3具体包括：对P₁求梯度模得到图像T₁,对T₁进行阈值变换H(T₁)得到H₀，

当下标越界时，利用镜像对称拓展。

5.根据权利要求4所述的医学图像噪声抑制算法，其特征在于，所述阈值转换函数H(x)还可为如下形式:

H(x)＝1+Alog(1+B(1-e^-kt))+Ct，

其中A,B,C,k均为认为设置的常数。

6.根据权利要求5所述的医学图像噪声抑制算法，其特征在于，所述阈值转换函数H(x)还可为如下形式:

参数k>0A控制最放大倍数，B控制增长速率变化。

参数选择原则如下：

为空气区域滤波半径放大倍数，

factor₁为边缘区域的最大滤波半径放大倍数，

所述空气区域为图像平坦区域，所述边缘区域为图像梯度模大的区域，

k控制阈值，图像梯度模大于多少时，需要缩小滤波半径，

B控制滤波半径放缩速率。

7.根据权利要求1所述的医学图像噪声抑制算法，其特征在于，所述步骤S4具体包括：

结合H₀和σ对图像P₀进行3维NLM自适应滤波得到P₂，

其中p,q∈N³,代表图像像素的三维位置向量，

p＝(p₁,p₂,p₃),q＝(q₁,q₂,q₃)，

U(p,r):{q∈N³|||p-q||＝max(|p₁-q₁|,|p₂-q₂|,|p₃-q₃|)≤r},

r为搜索半径，f为邻域半径，h为滤波半径，

r,f一般为向量形式，代表实际的图像在三个维度上相关性的不一致，此时，相应公式对应修改。

Images