CN103295190A - 数字x射线图像降噪方法 - Google Patents

Abstract

一种数字X射线帧降噪方法，包括获取原始图像、估计信号依赖噪声方差、组合双向非局部过滤、降低残留LF噪声、过平滑校正。估计信号依赖噪声方差包括：形态学提取噪声图像中与原始图像边缘对应的像素值；使用噪声方差区间估计的稳健局部线性近似获取描述噪声对信号强度依赖关系的表函数；计算原始图像噪声方差逐像素估计噪声图。组合双向非局部过滤通过组合NL均值过滤权重和双向过滤权重进行求平均。降低残留LF噪声使用双向过滤作为低通过滤的金字塔方案。获得最终图像后，基于过滤后图像与原始图像的结构相似性和像素相似性对最终图像和原始图像进行组合，进行过平滑校正。

Description

数字X射线图像降噪方法

技术领域

本发明涉及数字X射线图像处理领域，能够用来解决与对使用X射线等高能辐射获得的数字X射线帧进行降噪有关的问题。更具体地，本发明用来降低数字X射线帧序列（frame series）中的图像噪声。

背景技术

目前，医疗领域在数字X射线帧序列的处理方面包括各种方法：锐度增强（sharpness enhancement），解剖结构加重（anatomicalstructures emphasis），自动ROI（region of interest，感兴趣区）检测，等等。图像降噪方法被认为是提高数字X射线帧序列的质量的最重要的方法。图像降噪能够提高视觉感和低反差细节的可视性，并且能够提高其它处理算法的运算质量。因此，在数字射线照相领域，迫切需要开发一种过滤噪声的方法，使得能够保留诊断信息，并且能够以高分辨率、高帧频、良好的信号强度依赖噪声水平进行实时运算。

在开发过滤数字帧序列的算法过程中，研究人员面临两个问题。第一个问题涉及对基本上依赖于信号强度的图像噪声进行估计。第二个问题涉及在过滤质量和过滤速度之间寻求平衡。下面详细说明每个问题。

首先考虑依赖于信号的噪声的估计问题。在每个良好的数字图像过滤方法中，均使用噪声谱和振幅分布数据（amplitudedistribution data）。很大一部分降噪方法使用噪声方差作为参数。有很多作者在其出版物中提出了这样一种模型，即，将数字图像的噪声视为具有正态振幅和谱分布的白噪声的平稳随机过程。在检测中还存在降低瞬态噪声即方差值依赖于信号强度的噪声的问题，并且研究结果表明，噪声方差值越精确，降噪算法的结果越好（引证文件2、11）。

下面考虑数字X射线图像噪声的来源。探测器测量X射线辐射的衰减强度（即，穿透被检测物体的强度）。在辐照期间，每个探测器单元（detector cell）通过吸收光子平均积累

个电子。所积累的电子的数量N可以用泊松分布随机变量来模型化：

$P(N=n)=\exp (-\stackrel{\‾}{N})\frac{{\stackrel{\‾}{N}}^n}{n!},n\≥0$

所吸收的光子的随机变动称为光子噪声。在现代探测器中，光子噪声是主要的噪声源。其它噪声源包括探测器噪声、读出噪声（readout noise）、热噪声、放大器噪声（amplifier noise）、量化噪声（quantization distortion）等。上述噪声源的累积效应可以用高斯分布随机变量来模型化（引证文件2、10）。在广泛使用的模型中，对于线性电子电路，光子噪声和其它噪声源的方差线性依赖于真实的信号值（引证文件33）：

σ²(u(p))＝au(p)+b    （1）

其中，u(p)是图像中的像素p＝(i,j)的信号强度值，σ²(u(p))是依赖于信号强度的噪声方差值。

因此，数字图像噪声线性依赖于信号强度。关于信号依赖噪声的估计，有很多出版物（引证文件2、10、12、14、17、22、29）。因此，在引证文件14中，描述了一种对数字医疗帧序列逐帧进行噪声估计的非参数方法，其中，对实时运算的算法进行了加强。在引证文件10中，描述了一种数字帧噪声估计的双参数方法。从探测器获得的、未经伽马校正等非线性变换的原始数字图像的噪声可以通过在信号随机变量的基础上添加一项来模型化：

u_n(p)＝u(p)+σ²(u(p))n(p)    （2）

其中，u_n(p)是正被观测的噪声图像中的像素的信号强度值，σ²(u(p))是基于上式（1）的噪声方差对信号强度的依赖关系，n(p)∈N(0,1)是正态随机变量。有一种建立噪声方差的模型曲线的方法，该方法考虑导致曝光不足或过度曝光的传感器非线性，即，线性函数在动态范围的边界处的非线性。

上述第二个问题是在过滤运算质量和高运算速度之间提供平衡。现代X射线诊断快速处理高分辨率图像的复杂要求与高质量的过滤结果是相矛盾的。目前，有非常多的数学工具来进行数字图像的降噪。总体来说，这些降噪方法可以概括为两类过滤方法：

－变换域内的过滤

－原始图像域内的过滤

第一类过滤方法（变换域内的过滤）通过使用一定的变换来提供源图像的变换。然后，基于一定的规则例如使用适当的非线性运算符，来处理系数。该处理的目的是降低噪声的能量，以期提高图像锐度。该过滤处理结束于被重构（reconstruction）的图像的逆变换。该类方法包括大量基于傅立叶变换和小波变换的方法，例如引证文件25、28；基于复杂小波变换的方法，例如引证文件16、23；基于小曲线（curvelet）变换的方法，例如引证文件8；等等。这类过滤的运算速度和运算质量由所使用的变换和系数处理技术决定。

第二类过滤方法（原始图像域内的过滤）也包括大量方法，例如：各向异性扩散，例如引证文件31；变差算法（variationalgorithm），例如引证文件21；双向（bilateral）过滤，例如引证文件27；非局部均值（non-local means，NL-means）过滤，例如引证文件5；核回归（kernel regression），例如引证文件26；最优空间适应（optimal spatial adaptation）过滤，例如引证文件15；等等。总体上，这些技术均涉及在图像自身的像素之间或像素组（块（patch））之间进行相似性搜索，可能考虑几何数据（例如图像结构的局部方向），并且对这些像素的强度进行修正（求平均）算法，其中，修成的程度取决于像素（块）的相似程度。

引证文件：

1、Bilcu R.C.，Vehvilainen M.Fast non-local meansfor image denoising.In Digital Photography III Proc.of SPIE IS&T，Russel A.Martin，Jeffrey M.DiCarlo，andNitin Sampat，Eds.，February 2007，vol.6502.

2、Bosco A.，Battiato S.，Bruna A.R.，Rizzo R.NoiseReduction for Cfa Image Sensors Exploiting Hvs Behavior.Sensors 9(3)，1692-1713，2009.

3、Buades A.，Coll B.，Morel J.-M.A review of imagedenoising algorithms，with a new one.MultiscaleModeling and Simulation (SIAM interdisciplinaryjournal)，Vol.4(2)，pp.490-530，2005.

4、Buades A.，Coll B.，Morel J.-M.Image Denoising byNon-Local Averaging.Proc.of IEEE Int.Conf.on Acoustics，Speech，and Signal Processing，vol.2，pp.25-28，Mar.2005.

5、Buades A.，Coll B.，Morel J.-M.On image denoisingmethods.SIAM Multiscale Modeling and Simulation，4(2)：490-530，2005.

6、Buades A.，Coll B.，Morel J.-M.The staircasingeffect in neighborhood filters and its solution.IEEETransactions on Image Processing，Vol.15(6)，pp.1499-1505，2006.

7、Burt P.，Adelson E.The Laplacian pyramid as acompact image code.IEEE Trans.Commun.，vol.COM-31，pp.532-540，Apr.1983.

8、Candes E.J.，Donoho D.Curvelets，multiresolutionrepresentation，and scaling laws.In SPIE conferenceon Signal and Image Processing：Wavelet Applicationsin Signal and Image Processing VIII，San Diego，2000.

9、Dauwe A.，Goossens B.，Luong H.，Philips W.A FastNon-Local Image Denoising Algorithm.In Proc.SPIEElectronic Imaging，vol.6812，San Jose，USA，Jan 2008.

10、Foi A.，Trimeche M.，Katkovnik V.，EgiazarianK.Practical poissonian-gaussian noise modeling andfitting for single-image raw-data.Image Processing，IEEE Transactions on17，1737-1754(October 2008).

11、Foi A.Practical Denoising of Clipped orOverexposed Noisy Images.Proc.16th European SignalProcess.Conf.，EUSIPCO 2008，Lausanne，Switzerland，August 2008.

12、Forstner W.ImagePreprocessing for FeatureExtraction in Digital Intensity，Color and Range images.In Springer Lecture Notes on Earth Sciences，1998.

13、Gino M.Noise，Noise，Noise.http://www.astrophys-assist.com/educate/noise/noise.htm

14、Hensel M.，Lundt B.，Pralow T.，Grigat R.-R.Robust and Fast Estimation of Signal-Dependent Noisein Medical X-Ray Image Sequences.In Handels，H.，etal.，ed.：Bildverarbeitung fur die Medizin 2006：Algorithmen，Systeme，Anwendun-gen.Springer(2006)46-50.

15、Kervrann C.，Bourlanger J.Optimal spatialadaptation for patch-based image denoising.IEEETransactions on Image Processing，vol.15，no.10，October 2006.

16、Kingsbury N.G.The dual-tree complex wavelettransform：A new technique for shift invariance anddirectional filters.In Proc.8th IEEE DSP Workshop，Utah，Aug.9-12，1998，paper no.86.

17、Liu C.，Szeliski R.，Kang S.B.，Zitnick C.L.，Freeman W.T.Automatic estimation and removal of noisefrom a single image.IEEE Transactions on PatternAnalysis and Machine Intelligence 30(2)，299-314，2008.

18、Mahmoudi M.，Sapiro G.Fast image and videodenoising via nonlocal means of similar neighborhoods.IEEE Signal Processing Letters，12(12)：839-842，2005.

19、Peter J.，Govindan V.K.，Mathew A.T.Robustestimation approach for nonlocal-means denoising basedon structurally similar patches.Int.J.Open ProblemsCompt.Math.，vol.2，no.2，June 2009.

20、Press W.，Teukolsky S.，Vetterling W.，FlanneryB.Numerical Recipes in C：The Art of ScientificComputing，Second Edition.New Yor k：CambridgeUniversity Press，1992.

21、Rudin L.，Osher S.，Fatemi E.Nonlinear totalvariation based noise removal algorithms.Physica D，60：259-268，1992.

22、Salmeri M.，Mencattini A.，Rabottino G.，LojaconoR.Signal-dependent Noise Characterization forMammographic Images Denoising.IMEKO TC4 Symposium(IMEKOTC4′08)，Firenze，Italy，September 2008.

23、Selesnick I.W.The double-density dual-treediscrete wavelet transform.IEEE Trans.SignalProcessing，vol.52，no.5，pp.1304-1314，May 2004.

24、StarckJ.，Murtagh F.，Bijaoui A.Image Processingand Data Analysis：The Multiscale Approach.CambridgeUniversity Press，1998.

25、Strang G.，Nguyen T.Wavelets and Filter Banks.Wellesley-Cambridge Press.

26、Takeda H.，Farsiu S.，Milanfar P.Kernelregression for image processing and reconstruction.IEEE Transactions on Image Processing，vol.16，no.2，pp.349-366，February 2007.

27、Tomasi C.，Manduchi R.Bilateral Filtering forGray and Color Images.In Proc.6thInt.Conf.ComputerVision，New Delhi，India，1998，pp.839-846.

28、Vetterli M.，Kovacevic J.Wavelets and SubbandCoding.Prentice Hall，Englewood Cliffs，NJ，1995.

29、Waegli B.Investigations into the NoiseCharacteristics of Digitized Aerial Images.In：Int.Arch.For Photogr.And Remote Sensing，Vol.32-2，pages341-348，1998.

30、Wang J.，Guo Y.，Ying Y.，Liu Y.，Peng Q.Fastnon-local algorithm for image denoising.In Proc.OfIEEE International Conference on Image Processing(ICIP)，2006，pp.1429-1432.

31、Weickert J.Anisotropic Diffusion in ImageProcessing.Tuebner Stuttgart，1998.

32、ЯнеБ.Цифровая обработка изображений.М.，Техносфера，2007，583с.

33、BerndDigital Image Processing，Springer2005

发明内容

本发明提供一种包含两步的数字X射线图像降噪技术。第一步涉及信号强度依赖噪声的估计。第二步涉及降噪和一些辅助步骤，例如双向非局部过滤（bilateral nonlocal filtering）、残留LF噪声的金字塔去除（pyramidal removal）以及过平滑（over-smoothing）校正的组合。

对于数字医疗X射线图像噪声估计，与本发明最接近的是引证文件14，根据该引证文件，使用以下步骤基于原始图像来建立信号依赖噪声的估计：

-通过原始图像的LF过滤来估计真实信号，计算原始图像与其估计之间的差，从而获得噪声图像；

-在噪声图像中，通过任意方法去除与原始图像中的突变（边缘、单个“热”像素）相对应的像素值；

－将真实图像估计的强度范围分为若干个区间，对于每个区间，积累与真实图像估计的像素相对应的噪声帧像素值；

－在每个区间中，使用积累在该区间中的噪声图像像素值来计算噪声方差。

本发明使用上述引证文件中的噪声估计原理。主要区别在于：

（1）通过对噪声帧中与原始图像中的边缘相对应的像素值进行形态学提取（morphological extraction），来去除噪声帧中与原始图像中的突变相关的像素值；

（2）对噪声方差的区间估计进行稳健（robust）局部线性近似，从而生成描述噪声对信号强度的依赖关系的表函数；

（3）基于真实图像估计和计算出的表函数，计算噪声图（noisemap）作为原始图像噪声方差的逐像素估计的图。

对于数字医疗X射线图像降噪，与本发明最接近的是NL均值算法（引证文件5）。对于数字X射线图像的降噪，NL均值算法提供良好的结果。当前像素的信号强度估计通过对域求平均来进行：

$\hat{u}\left(p\right)=\underset{q\∈S(p,t)}{\mathrm{\Σ}}w(p,q)u\left(q\right)/\underset{q\∈S(p,t)}{\mathrm{\Σ}}w(p,q)---\left(3\right)$

其中，S(p,t)是以当前像素p为中心、半径为t的域（即，域的尺寸为2·t+1×2·t+1）；w(p,q)＝F(||U(p,f)-U(q,f)||_2，a，σ,k)是求平均权重；F是降函数；U(p,f)、U(q,f)是分别以像素p和q为中心、半径为f的像素块； ${\left|\right|U(p,f)-U(q,f)\left|\right|}_{2,a}=\sqrt{{\underset{h\∈S(0,f)}{\mathrm{\Σ}}G(h,a)\·(u(p+h)-u(q+h))}^2}$ 是像素块之间的高斯加权G(h,a)＝exp(-h·h^T/a²)欧几里德距离（Euclideandistance）；σ是图像噪声的均方根误差；k是控制求平均水平的参数；u(q)是原始图像像素的像素强度；

是当前像素的总估计。

在标准情况下，例如在引证文件4中，权重函数定义如下：

$F({\left|\right|U(p,f)-U(q,f)\left|\right|}_{2,a},\mathrm{\σ},k)=\exp ({-\left|\right|U(p,f)-U(q,f)\left|\right|}_{2,a}^2/k^2{\mathrm{\σ}}^2)$

假设噪声是额外的固定方差，即，噪声图像模型由公式（2）描述，其中，对于每个像素p，σ²(u(p))＝σ²＝const。引证文件中的NL均值的基本方案涉及很多方面，例如NL均值的改进、加速和实际应用（引证文件1、9、18、19、30）。

NL均值过滤非常难以计算。引证文件3描述了加速NL均值过滤运算的一个可能方式，其中，该引证文件的作者提出使用向量NL均值。

在该变种算法中，图像I(p)被分割为M个块U(p_m,f)，m＝1,M，

其中，p_m是分割块的中心像素。有重叠地进行块的分割。当前分割块的像素强度估计由如下公式进行：

$\hat{U}(p_m,f)=\underset{q\∈S(p_m,t)}{\mathrm{\Σ}}w(p_m,q)U(q,f)/\underset{q\∈S(p_m,t)}{\mathrm{\Σ}}w(p_m,q)---\left(4\right)$

其中， $w(p_m,q)=\exp ({-\left|\right|U(p_m,f)-U(q,f)\left|\right|}_2^2/k^2{\mathrm{\σ}}^2)$ 是块求平均权重； ${\left|\right|U(p_m,f)-U(q,f)\left|\right|}_2^2={\underset{h\∈S(0,f)}{\mathrm{\Σ}}(u(p_m+h)-u(q+h))}^2$ 是当前像素块和比较像素块之间的平方欧几里德距离；

是当前分割块的像素强度的最终估计。对最终块估计组合块重叠的求平均，生成作为输出项的过滤后图像。向量方案的应用能够对过滤进行加速，但处理后的图像的质量有一定的损失。

上述基本过滤方案（3）和（4）有很多缺点：

－残留噪声和特征伪像较高，它们以突变（discontinuity）形式沿着图像纹理区域的边缘出现；

－具有残留LF噪声；

－对细节过度过滤（过平滑）。

本发明的一个目的在于提供一种功能更强大的X射线图像降噪方法，具体来说，能够降低以突变形式沿着图像纹理区域的边缘出现的残留噪声和特征伪像，降低残留LF噪声，并且不对细节过度过滤（过平滑）。

为了提高数字X射线图像处理的质量，本发明提供一种数字X射线帧降噪方法，包括：获取原始图像、估计信号依赖噪声方差、组合双向非局部过滤、降低残留LF噪声、以及过平滑校正。估计信号依赖噪声方差的步骤包括：对噪声图像中与原始图像的边缘相对应的像素值进行形态学提取；通过使用噪声方差区间估计的稳健局部线性近似，来获取描述噪声对信号强度的依赖关系的表函数；计算噪声图，所述噪声图是原始图像噪声方差的逐像素估计。组合双向非局部过滤的步骤包括：通过组合NL均值过滤权重和双向过滤权重，来进行求平均。降低残留LF噪声的步骤使用金字塔方案，所述金字塔方案使用双向过滤作为低通过滤。在降低残留LF噪声的步骤后，获得最终图像，通过基于过滤后图像与原始图像的结构相似性和像素相似性对最终图像和原始图像进行组合，来进行过平滑校正。

优选地，在组合双向非局部过滤的步骤，基于使用线性低通滤波器平滑后的原始图像来计算求平均权重。

优选地，在组合双向非局部过滤的步骤，使用NL均值过滤权重和双向过滤权重分别占预定比例的组合权重函数，从而将NL均值过滤和双向过滤进行组合。

优选地，在组合双向非局部过滤的步骤，通过将NL均值过滤与双向过滤分开，来使用能够分开的双向过滤。

附图说明

以下附图示出本发明的技术方案、可能的实施方式、以及所获得的技术效果，其中：

图1示出获取X射线图像的装置。

图2示出头部的X射线图像。

图3示出表函数（附图标记5所示的实线），其描述对于图2所示的图像，正态噪声偏差对信号强度的依赖关系。

图4示出图2所示的图像的正态噪声偏差图。

图5A～5E示出本发明的方法的效果。

具体实施方式

X射线帧的获取由例如图1所示的装置进行。该装置包括X射线管1，其发射X射线束2。X射线束2进入探测器3。探测器3由闪烁屏（scintillation screen）（未图示）和光敏阵列（未图示）构成。闪烁屏与光敏阵列的有效表面光耦合。进入探测器3的X射线束2被闪烁屏转换为可见光，该可见光被探测器的传感器转换为数字形式。

在本发明中，滤波器在若干个步骤内降低噪声。以下详细说明每个所述步骤。

步骤1：对依赖于信号强度的噪声方差进行估计。在该步骤，通过对当前图像进行LF线性过滤（引证文件14），来进行原始图像估计。为了提高速度，可以使用最简单的线性过滤（例如，二项式过滤（binomial filter））。基于所获得的平滑后的图像来生成噪声图像，即，原始图像与过滤后的图像之间的差。

使用最简单的过滤所进行的图像估计是不理想的，这是因为边缘处被过度平滑。计算原始图像与平滑后的图像之间的差生成噪声图像，该噪声图像除了包括平滑区域中的噪声像素之外，还包括一定数量的与突变（例如，解剖结构边缘和非噪声像素）相对应的像素。这些像素能使所估计的方差值显著失真，因而应该从计算中去除。为此，有多种不同的方法（引证文件10、22），这些方法均对平滑后的原始图像的导数进行限值，其中，阈值由信噪局部估计（S/N localestimation）来确定。在包括大量细节的图像区域中，这种估计是不令人满意的。因此，在本发明中，当去除边缘时，使用更简单的方法来选择边缘，从而无需计算标准偏差的导数和局部估计。这种去除非噪声像素的形态学方法基本如下构成：

1、将噪声图像分为两个值，即，正、负变化的二值图像；

2、为了选择与区域相对应的边缘，使用后述的膨胀（dilatation）对所获得的图像进行形态学侵蚀（erosion）；

3、将处理后的图像合并成一个二值图像，即，原始图像的边缘图。

为了更好地保留细微结构，基于2×2小掩膜（2×2窗口）来进行侵蚀和膨胀。

当计算区间噪声方差估计时，确定估计图像的强度的最小值和最大值（强度范围的边界），选择子区间（例如，32个灰度级）。然后，对于估计图像的每个像素，找到含有该像素值的区间，并且使用适当的噪声图像像素值来计算该区间内的噪声方差估计（边缘像素已被去除）。为了计算噪声方差的区间估计，可以使用不同的公式，例如常见的无偏估计或使用中值绝对偏差（median absolute deviation）的稳健估计（引证文件14、10）。该过程生成描述噪声方差对信号强度的依赖关系的表函数。

构建边缘图时的不精确可能导致计算区间方差估计时的大的误差。因此，通过使用下述迭代非正常值去除（technique iterativeoutlier removal）技术（引证文件14），可以使这些估计更精确：对于每个区间，迭代进行去除噪声像素的绝对值超过阈值即三倍噪声标准偏差的像素和在该区间内再次计算噪声方差估计。

计算了区间噪声方差估计之后，估计噪声方差对强度的依赖关系。在参数方法中，噪声模型（1）的参数估计可以以一定的方式来构建（例如，最小二乘法、似然函数最小化、方向最优化（directionaloptimization））。还可以考虑传感器线性（引证文件10）。然而，如在引证文件14中所述，构建足以描述噪声方差对信号强度的依赖关系的参数模型因为一些因素而面临很多困难。这些因素可能包括传感器非线性以及对原始图像的非线性预处理（例如，取对数）。因此，基于区间方差估计，构建更可行、更通用的方法，即，噪声对强度的依赖关系的非参数估计。

在本发明中，使用非参数方法来构建上述依赖关系，其中，基于所获得的噪声方差的区间估计，来建立插值表函数。该表函数是基于区间方差估计的稳健局部线性近似而形成的。稳健方法的使用还可以减少非正常值（outlier）的影响（区间方差估计中大的误差），而局部近似则可以再现描述噪声对强度的依赖关系的复杂曲线。

因此，所获得的针对原始图像每个强度值的表函数与噪声方差估计相关。表的输入项可以是例如估计图像的强度。在图3中，示出了图2所示的图像的表函数（实线5）以及噪声均方差的区间估计（圆圈6）。

在实际处理中，对于参数噪声估计，可以进行变换以稳定原始图像的噪声方差（引证文件24、33）。因此，过滤信号依赖噪声的任务在于降低额外的非信号依赖噪声方差。在本发明中，当进行非参数噪声估计和构建噪声图时，使用下面的方法：基于估计图像和插值表来构建噪声图，在该噪声图中，每个像素是原始图像中的相应像素的噪声均方差估计。

噪声图提供了在实际应用中具有足够精度的逐像素噪声估计。图4示出了噪声图的一个例子，它是为图2所示的图像构建的。

该例子示出过度曝光（传感器饱和/耀斑）、在头部周围有个光斑（参见图2的附图标记4）时的噪声估计算法运算。

步骤2：组合过滤。当使用标准的过滤方案（公式（3））时，随着像素块尺寸的增加，在信噪峰值处测得的过滤质量也提高。除了质量提高之外，当使用向量NL均值（公式（4））时，随着像素块尺寸的增加，过滤速度也增加。然而，当在纹理显著的区域和强度剧烈变化处（边缘）增加像素块的尺寸时，可能发现极大的残留噪声和特征伪像（characteristic artifacts），它们的表现形式是在边缘的方向上突变。纹理区域中存在的残留噪声和特征伪像是因为降噪水平的降低而导致的，而降噪水平的降低是因为相似块的数量减少而导致的。在处理后的图像的平滑区域，可以看到大的残留噪声。在本发明中，采用组合降噪方法，该方法基于双向和NL均值的组合，旨在消除标准方法的上述不足。将这种方法称为“组合过滤算法”或“组合过滤”。组合非局部双向过滤是基于以下原理。

组合非局部双向过滤的原理。在求平均时，对于每个像素，基于两个矩阵即双向（亮度）和非局部（结构）使用其与当前像素的距离来计算权重。对于NL均值过滤，如果当前像素和比较像素在非局部均值过滤方面相似，即，它们的块在结构上相似，则结构权重占主导。相反，如果块之间相互不相似，而是当前像素和比较像素的强度相似，则增加双向权重的影响。为了降低残留噪声，使用平滑后的图像计算求平均权重，该平滑后的图像可以通过使用最简单的低通滤波器对原始图像进行过滤来获得。

标准过滤方案（3）可以修改如下：

$\hat{u}\left(p\right)=\underset{q\∈S(p,t)}{\mathrm{\Σ}}w(p,q)u\left(q\right)/\underset{q\∈S(p,t)}{\mathrm{\Σ}}w(p,q)---\left(5\right)$

其中，w(p,q)＝V(b(p,q),n(p,q))是定义组合过滤的求平均权重的函数，以满足组合过滤原理的方式来选择； $b(p,q)=G(q-p,a)\·F((\stackrel{\‾}{u}\left(p\right)-\stackrel{\‾}{u}\left(q\right)),\mathrm{\σ}\left(\stackrel{\‾}{u}\left(p\right)\right),k_b,h_b)$ 是双向过滤权重；G(s,a)＝exp(-s·s^T/a²)是高斯函数； $n(p,q)=F({\left|\right|\stackrel{\‾}{U}(p,f)-\stackrel{\‾}{U}(q,f)\left|\right|}_2,\mathrm{\σ}\left(\stackrel{\‾}{u}\left(p\right)\right),k_n,h_n)$ 是NL均值过滤权重；

分别是原始图像的像素和像素块的值；

是原始图像的强度依赖均方根误差。为了计算像素在强度和结构方面的相似性，使用连续变换阈值函数F(r,σ,k,h)＝exp(-max(r²-k·σ²,0)/h²)，其中，k、h是参数。

组合过滤的向量形式如下：

$\hat{U}(p_m,f)=(\underset{q\∈S(p_m,t)}{\mathrm{\Σ}}W(p_m,q,f)\·U(q,f))\·/\underset{q\∈S(p_m,t)}{\mathrm{\Σ}}W(p_m,q,f)---\left(6\right)$

其中，W(p_m,q,f)＝V(B(p_m,q,f),n(p_m,q))确定组合过滤的权重函数的值；B(p_m,q,f)是双向函数的权重函数的值，B(p_m,q,f)＝{b(p_m+s,q+s),s∈S(0,f)}，计算为 $b(p_m+s,q+s)=G(s,a)\·F((\stackrel{\‾}{u}(p_m+s)-\stackrel{\‾}{u}(q+s)),\mathrm{\σ}\left(\stackrel{\‾}{u}(p_m+s)\right),k_b,h_b),$ 其中G(s,a)是高斯函数； $n(p_m,q)=F({\left|\right|\stackrel{\‾}{U}(p_m,f)-\stackrel{\‾}{U}(q,f)\left|\right|}_2,\mathrm{\σ}\left(\stackrel{\‾}{u}\left(p_m\right)\right),k_n,h_n)$ 是NL均值过滤的权重，它对全部比较像素块是一样的；

是原始图像的强度依赖均方根误差；“·”是单项块乘法/除法；

是当前块的最终像素强度估计。

在上述组合过滤方案（5）、（6）中，可以由多种方式进行的权重函数V的选择是非常重要的。为了使组合过滤的速度更快，通过按一定比例使用亮度和结构权重来构建权重函数，可以生成简单但过滤质量高的组合过滤。

V(b,n)＝(1-λ)·b(p,q)+λ·n(p,q)    （7）

组合过滤的向量形式如下：

V(B,n)＝(1-λ)·B(p_m,q,f)+λ·n(p_m,q)    （8）

λ应该接近NL均值过滤权重值的最大可能值，例如，λ≈0.8。权重函数（7）或（8）使得能够将双向过滤和非局部过滤分开。当将这两个过滤分开时，组合过滤意味着具有双向过滤和非局部过滤的性能，并且其结果是这两个过滤的结果的组合。将过滤分开使得能够选择更具速度的过滤参数（例如，降低双向过滤的求平均区域），并且能够去除各块相互重叠的区域内的无用的双向过滤像素。使用向量形式的组合过滤（6）和公式（8）所示的权重函数，可以将组合过滤表示如下：

$\hat{U}(p_m,f)=$

$=(\mathrm{\λ}\·{\hat{U}}_n(p_m,f)+(1-\mathrm{\λ})\·{\hat{U}}_b(p_m,f))\·/(\mathrm{\λ}\·c\left(p_m\right)+(1-\mathrm{\λ})\·C(p_m,f))$

其中， ${\hat{U}}_n(p_m,f)=\underset{q\∈S(p_m,t)}{\mathrm{\Σ}}w(p_m,q)U(q,f)$ 是使用向量过滤进行NL均值过滤的结果（未进行权重归一化）；是双向过滤后的图像的像素块（未进行权重归一化），即， ${\hat{U}}_b(p_m,f)=\{\hat{u}\left(s\right),s\∈S(p_m,f)|\hat{u}\left(s\right)=\underset{q\∈S(s,t)}{\mathrm{\Σ}}b(s,q)u\left(q\right)\};$ c(p_m)、C(p_m,f)分别是NL均值过滤和双向过滤的权重之和。

对于除法形式的组合过滤（9），当使用可分离的双向过滤时，能够进一步提高过滤速度，从而更接近向量形式的NL均值过滤的速度。

图5A～5E示出了使用权重（8）的组合过滤方法（6）的有效性，其中，图5A示出X射线测试图案片断；图5B示出由向量形式的标准NL均值过滤（4）处理后的X射线测试图案片断；图5C示出由组合过滤处理后的X射线测试图案片断；图5D示出图5A和图5B之间的差；图5E示出图5A和图5C之间的差。组合过滤的使用能够对数字边缘处的可以看见的伪像进行降噪，而这些伪像存在于使用标准NL均值过滤时。

步骤3：L F降噪。在强度相同、噪声方差为常量的平滑图像区域中进行双向和NL均值过滤是最简单的线性相邻求平均过滤（即，移动平均（moving average））。因为使用有限个区域来求平均，因此，噪声谱中与L F相对应的成分不会被降低，因而在平滑图像区域中显著可见。

在本发明中，使用线性回归的双向过滤在图像金字塔中进行残留L F噪声的降噪（引证文件7、27）。使用如下方式构建图像金字塔：

$\left\{\begin{array}{c}{A}_0\left(p\right)=I\left(p\right),\\ A_k\left(p\right)=\mathrm{Reduce}\left(A_{k-1}\left(p\right)\right),\\ H_{k-1}\left(p\right)=A_{k-1}\left(p\right)-\mathrm{Expand}\left(A_k\left(p\right)\right),k=\stackrel{\‾}{1,n}.\end{array}\right.---\left(10\right)$

其中，I(p)是原始图像；A_k(p)、H_k(p)是图像金字塔的第k阶低通（近似）和高频（细节）图像；Reduce是过滤并2倍降低图像采样频率的运算符（比例缩小（scaling）），即，Reduce(A)＝(A*L_r)↓2，其中，L_r是线性低通滤波器，↓2是比例缩小；Expand是聚集采样（crowdingsample）和图像过滤运算符，即，Expand(A)＝(A↑2)*L_e，其中，↑2是聚集采样，L_e是低通滤波器。

对低通金字塔成分（10）估计噪声方差分布。为此，对组合过滤添加标准偏差为σ的高斯噪声。因为上述在平滑图像区域中进行过滤与相邻求平均过滤的相似性，获得依赖于金字塔阶数的低通噪声均方根误差：

$\left\{\begin{array}{c}{\mathrm{\σ}}_0=\mathrm{\σ}/(2\·t+1),\\ {\mathrm{\σ}}_k={\mathrm{\σ}}_0\·q(k,L_r),k=\stackrel{\‾}{1,n},\end{array}\right.---\left(11\right)$

其中，σ₀是由组合过滤进行处理后的图像噪声的均方根误差；σ_k是金字塔（10）的图像噪声A_k(p)的均方根误差；q(k,L_r)是均方根误差在图像A_k(p)中的衰减系数。

为了降低残留噪声，对金字塔（10）的图像A_k(p)进行双向过滤。为了降低阶梯伪像（stairstepping artifact）（引证文件6），使用如下线性回归双向过滤：

${\hat{a}}_k(i,j)=a\·i+b\·j+c---\left(12\right)$

参数(a,b,c)的选择使函数 $\underset{a,b,c}{\min }\left({\underset{(x,y)\∈S(i,j,t)}{\mathrm{\Σ}}w(x,y)(a_k(x,y)-a\·i-b\·j-c)}^2\right)$ 最小化；S(i,j,t)是金字塔（10）的图像k的当前像素的a_k(i,j)为中心、t为半径的区域；w(x,y)＝G(x-i,y-j,a)·F((a_k(i,j)-a_k(x,y)),σ_k(i,j))是求平均权重，其中，相邻像素由高斯函数G(x-i,y-j,a)＝exp(-((x-i)²+(y-j)²)/a²)描述，相邻像素强度由函数 $F((a_k(i,j)-a_k(x,y)),{\mathrm{\σ}}_k(i,j))=\exp ({-(a_k(i,j)-a_k(x,y))}^2/{\mathrm{\σ}}_k^2(i,j))$ 估计；σ_k(i,j)是噪声均方根误差，取决于金字塔近似图像的强度，由噪声图和以上关系（11）计算；

是当前像素强度的总估计。使用可分离滤波器（12）来加速过滤。

使用过滤后的金字塔近似成分来重建残留LF噪声被降噪的图像：

${\hat{A}}_{k-1}\left(p\right)={\hat{H}}_{k-1}\left(p\right)+\mathrm{Expand}\left({\hat{A}}_k\left(p\right)\right),k=\stackrel{\‾}{n,1}.---\left(13\right)$

其中，

是使用滤波器（12）所示的低通图像过滤所构建的金字塔（10）。

步骤4：过平滑校正。在本发明中，通过混合过滤后图像和其原始图像，来校正对有用的HF信息的过度平滑。混合权重基于过滤后图像与其原始图像的像素相似性和结构相似性来计算。为了计算结构相似性，计算图像块之间的欧几里德距离；为了计算像素相似性（即，强度相似性），计算相应图像的像素强度之间的距离。

令

u(p)分别表示过滤后图像和原始图像中像素p的强度值，则使用如下阈值函数来计算帧像素相似性：

$w_i\left(p\right)=\exp (-\max ({(\hat{u}\left(p\right)-u\left(p\right))}^2-k_i\·{\mathrm{\σ}}^2\left(u\left(p\right)\right),0)/h_i^2)---\left(14\right)$

其中，k_i、h_i是确定函数典型形状的参数；σ²(u(p))是像素p的噪声方差值。过滤后图像与其原始图像的结构相似性计算如下：

$w_s\left(p\right)=\exp (-\max ({\left|\right|\hat{U}\left(p\right)-U\left(p\right)\left|\right|}_2^2-k_s\·{\mathrm{\σ}}^2\left(u\left(p\right)\right),0)/h_s^2)---\left(15\right)$

其中，k_s、T_s是决定阈值函数形状的参数；

是相应图像的像素块之间的欧几里德距离。

为了计算图像最终的像素相似性和结构相似性，将所计算出的像素相似性和结构相似性按照如下方式合并：

w(p)＝w_i(p)·w_i(p)+(1-w_i(p))·w_s(p)    （16）

然后，过平滑校正包括根据如下规则将过滤后图像与原始图像合并：

${\hat{u}}_c\left(p\right)=w\left(p\right)\·\hat{u}\left(p\right)+(1-w\left(p\right))\·u\left(p\right)---\left(17\right)$

相似性测量（16）和校正（17）能够对图像中的单个突变（所谓的“突发（burst）”）进行过滤。

实施例

为了估计图像噪声方差，首先获取估计图像和噪声图像。为此，使用3×3尺寸H₁＝[1,2,1]/4,

的低通线性二项式滤波器对原始图像I(x,y)进行过滤。由此，获得平滑后的图像I_e(p)＝I*H，其中，“*”代表卷积。然后，计算噪声图像N_e(p)＝I(p)-I_e(p)。

当基于噪声图像N_e(p)从方差计算中去除边缘时，形成如下正、负变化图像：

$N_e^{-}\left(p\right)=\left\{\begin{array}{cc}1,& N_e\left(p\right)<0\\ 0,& N_e\left(p\right)\&GreaterEqual;0\end{array}\right.,$ $N_e^{+}\left(p\right)=\left\{\begin{array}{cc}1,& N_e\left(p\right)>0\\ 0,& N_e\left(p\right)\&le;0\end{array}\right.$

为了选择与边缘相对应的区域，使用2×2掩膜的后述膨胀对上述二值图像进行形态学侵蚀：

$B_e^{-}\left(p\right)={\mathrm{dilate}}_{\left[2x2\right]}\left({\mathrm{erode}}_{\left[2x2\right]}\left(N_e^{-}\left(p\right)\right)\right)$

$B_e^{+}\left(p\right)={\mathrm{dilate}}_{\left[2x2\right]}\left({\mathrm{erode}}_{\left[2x2\right]}\left(N_e^{+}\left(p\right)\right)\right)$

然后，组合这些图像得到原始图像的边缘图：

$E\left(p\right)=B_e^{-}\left(p\right)\&cap;B_e^{+}\left(p\right)$

为了计算区间噪声方差估计，获取图像I_e(p)的强度最小值I_min和最大值I_max，选择节距h_M，然后以节距h_M将强度范围分为若干个区间M_i（例如，h_M＝32）。对于每个图形像素I_e(p)，找到包含该像素值的区间，然后使用图像N_e(p)的相应像素值来计算该区间M_i的噪声方差估计σ²(i)，其中，去除边缘图的值E(p)＝1的像素。当计算区间噪声方差估计时，使用如下无偏估计：

${\mathrm{\&sigma;}}^2\left(i\right)=\underset{j=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{(N_e^i\left(j\right)-{\stackrel{\&OverBar;}{N}}_e^i)}^2/(n_i-1)$

其中，

是噪声图像N_e(p)在区间M_i中的像素的值，n_i是区间M_i中所积累的像素值的总个数，

是区间M_i中的噪声像素值的平均值。

对于每个区间M_i，使用所获得的区间噪声方差估计σ²(i)来迭代地去除噪声像素值的绝对值超过阈值即三倍噪声标准偏差的噪声像素以及然后在该区间中再次计算噪声方差估计：

$N_e^i[k+1]=\left\{N_e\left(p\right)\right|(N_e\left(p\right)\&Element;N_e^i[k\left]\right)\&cap;(N_e\left(p\right)\&le;3\mathrm{\&sigma;}\left(i\right))\}$

其中，

是第k次迭代区间M_i中的噪声像素值的集合。在后续计算中，仅使用积累了足够数量（例如，至少500个）的噪声像素的区间。此外，平均值

显著不等于零（偏差大于区间步长h_M的一半）的区间也不予考虑，这是因为在这些区间内残留的非噪声像素值占主导的可能性很大。

基于所获得的噪声方差区间估计，建立插值表函数，以估计噪声方差对强度的依赖关系。通过使用噪声方差区间估计的稳健分段线性近似，来建立该表函数。为此，对于区间M_i，选择适当的步长h_I和半径r_I（可以根据区间M_i中的像素数量n_i来选择）。在对上一步表函数所获得的值进行近似时，使用如下规则：

${\widehat{\mathrm{\&sigma;}}}^2(k\&CenterDot;h_I)=a\&CenterDot;m_k+b$

其中，k是近似数；m_k是区间M_i(k·h_I)的中心；参数a,b是基于绝对值偏差最小而计算的（引证文件20）：

$\underset{j={\mathrm{kh}}_I-r_I}{\overset{{\mathrm{kh}}_I+r_I}{\mathrm{\&Sigma;}}}|{\widehat{\mathrm{\&sigma;}}}^2\left(j\right)-a\&CenterDot;m_k-b|\&RightArrow;\underset{a,b}{\min }$

上述处理得到在整个强度范围[I_min,I_msc]进行插值的表的值即，得到所需的噪声方差对信号强度的依赖关系的插值表。

根据估计图像和插值表来生成噪声图，在该噪声图中，每个像素用来估计原始图像相应像素的噪声方差。

为了在过滤质量和过滤速度之间提供良好的平衡，使用表达式（9）进行组合双向非局部过滤，其中，λ＝0.85，求平均域的半径t＝5，像素块的半径被选择为f＝4，使一个像素重叠地选择块。为了计算像素在强度和结构上的相似性，使用连续变换阈值函数F(r,σ,k,h)＝exp(-max(r²-k·σ²,0)/h²)，其中，k、h是参数。

当强度相似性k_b＝9时，结构相似性k_n＝2.1。采用这样的方式选择h_b，使得对函数F提供预定的衰减速度，为此，对于给定的阈值kσ²例如kσ²＝10·σ²，选择F(r,σ,k,h)的函数值例如F＝0.1，从而使h_b＝9。当计算权重函数时，不考虑几何距离。以同样的方式，为NL均值过滤的权重函数选择h_n的值：对于kσ²＝3·σ²，定义F＝0.5，从而使h_n＝1。用来计算求平均权重的估计图像通过使用尺寸为H₁＝[1,2,1]/4,

的低通线性二项式滤波器对原始图像进行过滤来计算。值σ从噪声估计阶段所计算的噪声图获得。

在对残留的L F噪声进行过滤的阶段，使用阶数n＝4的金字塔图像（10），其中，作为Reduce中的低通滤波器，L_r＝[1,1]/2，

在该过滤中，LF金字塔成分中的噪声均方根误差基于q(k,L_r)＝2^-k、

而降低。当在运算符Expand进行聚集采样之后建立HF金字塔成分时，使用低通滤波器L_e＝[1,4,6,4,1]/16、

在低通过滤时，使用双向过滤，不考虑几何距离，过滤半径为5。权重函数参数与组合过滤中的双向过滤一致。当计算权重函数时，不考虑几何距离。

在过平滑校正时，为了计算结构相似性，使用半径为5的块。定义权重函数的形状的参数与组合过滤权重函数的参数相匹配。当计算像素强度相似性的权重函数时，不考虑几何距离。

Claims (4)

1.一种数字X射线帧降噪方法，包括：获取原始图像、估计信号依赖噪声方差、组合双向非局部过滤、降低残留LF噪声、以及过平滑校正，

估计信号依赖噪声方差的步骤包括：对噪声图像中与原始图像的边缘相对应的像素值进行形态学提取；通过使用噪声方差区间估计的稳健局部线性近似，来获取描述噪声对信号强度的依赖关系的表函数；计算噪声图，所述噪声图是原始图像噪声方差的逐像素估计；

组合双向非局部过滤的步骤包括：通过组合NL均值过滤权重和双向过滤权重，来进行求平均；

降低残留LF噪声的步骤使用金字塔方案，所述金字塔方案使用双向过滤作为低通过滤；

在降低残留LF噪声的步骤后，获得最终图像，通过基于过滤后图像与原始图像的结构相似性和像素相似性对最终图像和原始图像进行组合，来进行过平滑校正。

2.根据权利要求1所述的数字X射线帧降噪方法，其特征在于，在组合双向非局部过滤的步骤，基于使用线性低通滤波器平滑后的原始图像来计算求平均权重。

3.根据权利要求1所述的数字X射线帧降噪方法，其特征在于，在组合双向非局部过滤的步骤，使用NL均值过滤权重和双向过滤权重分别占预定比例的组合权重函数，从而将NL均值过滤和双向过滤进行组合。

4.根据权利要求1或2所述的数字X射线帧降噪方法，其特征在于，在组合双向非局部过滤的步骤，通过将NL均值过滤与双向过滤分开，来使用能够分开的双向过滤。

Images