CN102136138B - 一种医学图像降噪与增强的处理方法 - Google Patents

Abstract

本发明提供一种实时医学视频图像去噪方法，应用于实时获取或已存储的医学视频图像序列，利用相邻两帧图像的高度相关性，对每一帧视频图像实施时空域的递归的滤波处理，滤波窗口的大小可以根据噪音水平自适应地来调节。本发明结合时域的运动补偿方法和空域上的降噪增强算法，针对图像当前像素点，应用运动补偿方法和形态学算子检测运动区域、噪音区域和背景区域，在噪音区域中估计图像的整体噪音水平。

Description

一种医学图像降噪与增强的处理方法

技术领域

本发明涉及对获取的医学图像进行降噪与增强的方法。

背景技术

医学超声B模式图像中通常存在大量的斑点噪声，会给超声图像质量带来明显的下降，并掩盖了某些重要组织的病变，这给医生的诊断和识别某些特定的疾病带来了较大的困难，并且会有漏诊和误诊的风险。

通常采用滤波的方法来抑制信号的高频部分，在滤去噪声的同时也会将有用的边缘信息去除，例如邻域平均法、中值滤波法都会将图像中有临床意义的边缘和细节会过平滑；由此，一些基于边缘保留的滤波算法得到迅速的发展，例如各向异性滤波(anisotropic diffusion)，BM3D，和nonlocal means等医学图像去噪技术。这给在医学图像后处理阶段提高医学图像质量提供了一个巨大的空间。

各向异性滤波器是基于物理学中的扩散方程发展起来的。Perona和Malik于1990将其引入到图像处理当中作为一种边界保留(edge-preserving)平滑滤波器，其成果发表在IEEE Transactions on Pattern Analysis and Machine Intelligence，vol.12中的名为″Scale space and edge detection using anisotropic diffusion″的文章中。他们的扩散方程如下

$I_t=\mathrm{div}(g(\▿I)\▿I).$

其中I代表图像，I_t为原图像对进化时间t的偏导数，g(·)为一个与梯度相关的0到1的标量，称作扩散系数。进化时间t也可以看成是尺度空间里的尺度，即t越大，图像越平滑。

根据显式离散化方案，图像I在进化时间t可以表达为

$I^t=I^{t-1}+\mathrm{\δ}\·\mathrm{div}(g(\▿I)\▿I),$

其中δ为进化时间步长，在显式方案中其取一个不能大于0.25的正数以保证数值稳定性。当采用隐式或半隐式离散方案时，δ可不受此限制而取较大的值。隐式和半隐式离散方案可参考相关文档。

Perona和Malik为扩散系数提出了如下两个模型。

$g(\▿I)=e^{-{\left(\right||\▿I||/K)}^2},$

$g(\▿I)=\frac{1}{1+{\left(\frac{\left|\right|\▿\left|\right|}{T}\right)}^2},$

其中T为梯度模值的阈值。根据上述两个扩散系数模型，梯度模值

与扩散系数g成反比关系，并且当

远大于K时，g趋为0，相反，当

很小时，g趋为1。

因此，在梯度模值大的图像区域，也就是边界信息强的区域，扩散停止；在梯度模值小的区域，一般为物体内部或背景，扩散增大。从而实现了在保留边界信息的同时，图像的其他区域得到平滑的效果。

应该可以看到，因为扩散系数g是一个标量，使得流量(flux)与梯度方向平行，所以Perona-Malik模型从严格的物理意义上来讲是一种非线性的各项同性(non-linear inhomogeneous isotropic)扩散模型。具体解释可见Weickert在1998由B.G.Teubner Stuttgart出版的″Anisotropic Diffusion in Image Processing″一书中的第一章。为方便起见，在本发明中仍然称Perona-Malik模型为各向异性。

基于Perona-Malik模型，后继衍生出很多不同的平滑滤波器，比较著名的有Yu和Acton在IEEE Transactions on Image Processing，vol.11，no.11上发表的″Speckle reducing anisotropic diffusion″。Yu和Acton使用了一种叫做instantaneous coefficient of variation来替代梯度值作为边缘检测指标，并使用图像中预先分割出的一块噪音区域生成扩散阈值。

在各向异性滤波的实际应用中，因为需要多次扩散以达到目标效果，速度成为阻碍其应用的一大制约因素。Acton于1998在IEEE Transact ions on Image Processing，vol.7，no.3中发表了一篇名为″Multigrid anisotropic diffusion″的文章。在其中他设计了一个基于多尺度的各项异性滤波器已达到加速和去除因低频信号丢失造成的伪像的目的。

为了实现真正意义上的各向异性滤波，Weickert于1999年在InternationalJournal of Computer Vision，vol.31，no.2中发表了一篇名为″Coherence-enhancing diffusion filtering″的文章。在文中，Weickert使用了扩散张量(diffusion tensor)，一个2x2矩阵D，代替扩散系数g。扩散方程因此可写成

$I_t=\mathrm{div}(D\·\▿I)$

使用扩散张量，使得扩散可以沿着任意方向θ进行而不仅限于梯度方向。沿任意θ方向的扩散可分解为沿梯度方向与沿切线方向分别做扩散。扩散张量D可分解为

$D=\left(\begin{array}{cc}{v}_1& v2\end{array}\right)\left(\begin{array}{cc}{g}_1& 0\\ 0& g_2\end{array}\right)\left(\begin{array}{c}{v_1}^T\\ {v_2}^T\end{array}\right).$

其中

为当前点的法线方向，

为当前点的切线方向，g₁和g₂分别为沿着法线方向和切线方向的扩散系数，定义如下。

g₁＝α

$g_2=\left\{\begin{array}{cc}\mathrm{\α}& \mathrm{if\κ}=0,\\ \mathrm{\α}+(1-\mathrm{\α})\exp \left(\frac{-T}{\mathrm{\κ}}\right)& \mathrm{else}\end{array}\right.$

其中κ为一致性检测(coherence measure)，α为一个很小的正值，Weickert建议取0.001，T为一个大于0的阈值。

一致性检测κ是一个从结构张量(structure tensor)，J，导出的非负数。结构张量J定义为

$J_{\mathrm{\ρ}}=K_{\mathrm{\ρ}}*(\▿I\⊗\▿I)$

$=\left(\begin{array}{cc}{K}_{\mathrm{\ρ}}*I_x^2& K_{\mathrm{\ρ}}*\left(I_x{I}_y\right)\\ K_{\mathrm{\ρ}}*\left(I_x{I}_y\right)& K_{\mathrm{\ρ}}*I_y^2\end{array}\right),$

I_x和I_y分别为x和y方向的梯度值，K_ρ为一个方差为ρ的高斯滤波器。一致性检测κ定义为

κ＝(μ₁-μ₂)²，

μ₁和μ₂分别为J的最大和最小的特征值。

在Weickert的方向性各向异性模型中，因为α为一很小的值，在法线方向上的扩散可以认为被一直抑制，在切线方向上的扩散强度是由一致性检测κ来确定的。当κ→0时，梯度在法线方向和切线方向为接近的值，在切线方向上的扩散也被抑制。当κ＞＞T时，此时梯度在一个方向的分量远大于另一个方向的分量，也就是最有可能是边界的地方，切线方向上的扩散系数将接近值1.

发明内容

本发明提出一种医学图像降噪与增强的处理方法。基于方向性各向异性模型与多尺度方案，针对不同的检查器官，本发明可以实现对噪音区域和结构体区域不同程度的降噪和增强效果。尽管在本发明叙述中以超声B模式图像为例，但不应将此视为对本发明的限制。本发明适用于多种医疗成像，包括但不仅限于超声，X射线，磁共振，CT，等。

本发明为解决上述技术问题所采用的技术方案为：

本发明主要是提出了一种医学图像降噪与增强的实现方法，应用于实时或已存储的医学图像或视频图像序列。在一个优选实施方式中，给定一幅医学图像，本发明首先使用拉普拉斯金字塔模型(Laplacian pyramid)将其分解为L层，其中l₀代表原图，l₁Kl_L-1由前一层的2倍降采样获得。从最高一层l_L-1开始，使用一致性检测将图像中的结构体与噪音区域区分开。对噪音区域进行N_L-1次法线和切线方向扩散；对结构体区域进行S_L-1次切线方向扩散。在实际应用中，为了达到更好的视觉效果，在对结构体区域沿切线方向进行扩散的同时，法线方向的扩散系数可设为一个相对较小的固定值，例如g₁＝0.01。当l_L-1层扩散完成后，将其升采样并重建为l_L-2层图像l_R(L-2)，并照l_L-1层所述方法进行扩散。以此类推。当l_R(1)扩散处理完成后，将其重建为l₀层图像l_R(0)，并可作为最终处理结果输出。在一些条件下，l_R(0)也可按照上述方法进行扩散滤波，并将滤波后的图像作为输出。为了达到更佳的视觉效果，可以对l_R(0)进行后处理，包括对比度增强和与原图融合。

一种医学图像降噪与增强的处理方法，其给定一幅医学图像I后包括以下步骤：

A.使用多尺度模型将输入图像I分解为L个尺度，设0为原图尺度，L-1为最小尺寸尺度；

B.从最小尺度图像开始，设当前图像为u，

a.判断图像u是否需要处理，如不需要，跳到步骤d，

b.使用分割算法将u分为噪音区域与结构体区域，

c.使用方向性各向异性滤波器，在噪音区域进行N次扩散滤波，在结构体区域进行S次平滑扩散滤波，N和S为非负整数，

d.升采样结果到下一个大尺度并使用多尺度模型进行重建，

e.将重建结果设为u并重复步骤a到e直到u达到尺度0，

f.当u达到尺度0时，判断图像u是否需要处理，如果需要，执行步骤b到c，否则直接输出u；

C.对滤波结果进行后处理或直接输出结果。

当所述的多尺度模型中的尺度数量为1时，其特征为，

A.使用分割算法将I分为噪音区域与结构体区域，

B.使用方向性各向异性滤波器，在噪音区域进行N次扩散滤波，在结构体区域进行S次平滑扩散滤波，N和S为非负整数，

C.对滤波结果进行后处理或直接输出结果。

所述分割算法可包含以下步骤：

A.对于当前图像u，进行一致性检测κ，

B.使用预设阈值或自动阈值检测算法，对κ进行噪音区域与结构体区域分割，获得的阈值作为扩散滤波阈值。

在噪音区域中，扩散沿梯度与切线方向同时进行，以达到平滑噪音区域的目的。

在结构体区域中，扩散主要沿梯度方向进行，以达到增强结构体的目的。

对滤波图像的后处理包括对比度增强。

对滤波图像的后处理包括原始输入图像融合。

在使用多尺度模型进行重建之前，对当前尺度下的残值图像也进行方向性的各项异性滤波。

所述的多尺度模型为拉普拉斯金字塔模型。

判断在某一尺度上的图像是否需要进行处理由用户预先设置。

采用上述方法，可以针对不同的检查器官设置不同的噪音区域扩散次数和结构体扩散次数，从而灵活的控制图像降噪与增强的强度。由于使用拉普拉斯金字塔的多尺度模型，从而大大降低了整体处理时间。同时，在图像重建过程中，避免了因低频信号丢失造成的伪像。

附图说明

图1是本发明一个多尺度实施例的医学图像降噪与增强方法流程示意图；

图2是本发明一个单尺度实施例的医学图像降噪与增强方法流程示意图。

具体实施方式

下面根据附图和实施例对本发明作进一步详细说明：

附图1是该发明一个优选实施例的图像降噪和增强方法流程图，该处理方法适用于大多数医疗图像的降噪和增强，例如超声B模式图像，X线透视机图像，磁共振图像，和CT图像等。本发明主要使用扩散张量以实现对噪音区域和结构体区域分别进行不同强度的扩散滤波，同时使用多尺度以提高运行速度并且降低由于低频信号丢失所造成的伪像，提高图像颗粒度。方案具体流程如下。

从步骤100开始，首先从图像获取设备或图像储存器获取一幅医学图像，I。

在步骤101中，使用拉普拉斯金字塔模型将I分解为L(L＞1)层金字塔模型，其中l₀代表原图，l₁Kl_L-1由前一层经高斯滤波后2倍降采样获得。设h₀Kh_L-2为l₀Kl_L-2所对应的残余图像(residue image)，

h_i＝l_i-K_σ*l_i，

K_σ为方差为σ的高斯滤波器。

在步骤102中，取L-1层图像为u，设i＝L-1。

在步骤103中，判断当前图像u是否需要被处理，如果不需要，跳到步骤107执行。此判断可由用户预先设定并作为参数输入。

在步骤104中，将u作为输入，计算一致性检测κ。使用大津法(Otsu)对κ自动计算阈值T，将u分为噪音区域和结构体区域。

在步骤105中，使用扩散张量D对u进行N_i次扩散。法线和切线方向上扩散系数g₁和g₂定义如下，

$g_1\left(\mathrm{\κ}\right)=g_2\left(\mathrm{\κ}\right)=1-\exp \left(\frac{-C_{m_1}}{{(\mathrm{\κ}/T)}^{m_1}}\right),$

其中m₁决定了扩散系数随着κ变化的速度。m₁值越大，扩散系数变化越快。

为一常数以使当κ＜T时，流量(flux，)上升；当κ＞T时，流量下降。当m₁＝2时，

在每次扩散完成后，需要对κ进行更新。为了节省处理时间，在实践当中，可由用户自定义κ的更新频率。

在步骤106中，使用扩散张量D对u进行S_i次扩散。法线和切线方向上扩散系数g₁和g₂定义如下，

g₁(κ)＝α

$g_2\left(\mathrm{\κ}\right)=\left\{\begin{array}{cc}\mathrm{\α}& \mathrm{if\κ}=0.\\ \mathrm{\α}+(1-\mathrm{\α})\exp \left(\frac{-T}{\mathrm{\κ}}\right)& \mathrm{else}\end{array}\right.$

在本方案中，α＝0.01。类似于步骤105，在每次扩散完成后，也需要对κ进行更新。

同样为了节省处理时间，在实践当中，可由用户自定义κ的更新频率。

在步骤107中，对处理后的u进行升采样并使用h_i-1进行i-1层图像重建得到u_r。在重建之前，可以使用步骤105和106中的扩散系数对h_i-1进行修正以达到更佳的效果。

在步骤108中，设u＝u_r，i＝i-1。

在步骤109中，判断i是否为0。如果不是，重复步骤103到步骤109。如果i为0，说明此时已达到原图尺度。在一些情况下，可对u进行一次步骤104到步骤106的扩散滤波，否则执行步骤110。

在步骤110中，可进一步对u进行后处理以达到更佳的视觉效果。在本优选方式中，使用对比度增强并且将结果与原输入图像进行融合。

在步骤111中将u作为最终结果输出。

附图2是当不使用多尺度的情况下，本发明的一个优选方式。

从步骤200开始，首先从图像获取设备或图像储存器获取一幅医学图像，I。设u＝I。

步骤201与步骤104类似，不再重复。

步骤202与步骤105类似，不再重复。

步骤203与步骤106类似，不再重复。

步骤204与步骤110类似，不再重复。

在步骤205中，输出最终处理结果。

本领域技术人员不脱离本发明的实质和精神，可以有多种变形方案实现本发明，以上所述仅为本发明较佳可行的实施例而已，并非因此局限本发明的权利范围，凡运用本发明说明书及附图内容所作的等效结构变化，均包含于本发明的权利范围之内。

Claims (7)

1.一种医学图像降噪与增强的处理方法，特征在于，其给定一幅医学图像I后包括以下步骤：

A.使用多尺度模型将输入图像I分解为L个尺度，设0为原图尺度，L-1为最小尺寸尺度；

B.从最小尺度图像开始，设当前图像为u，

a.判断图像u是否需要处理，如不需要，跳到步骤d，

b.将图像u作为输入，计算一致性检测κ，使用预设阈值或自动阈值检测算法，对κ进行噪音区域与结构体区域分割，获得的阈值作为扩散滤波阈值T，将图像u分为噪音区域和结构体区域，

c.使用方向性各向异性滤波器，在噪音区域进行N次扩散滤波，法线和切线方向上扩散系数g₁和g₂定义如下为 

其中m₁＝2，

在结构体区域进行S次平滑扩散滤波，法线和切线方向上扩散系数g₁和g₂定义为

其中α＝0.01，N和S为非负整数，

d.升采样结果到下一个大尺度并使用多尺度模型进行重建，

e.将重建结果设为u并重复步骤a到e直到u达到尺度0，

f.当u达到尺度0时，判断图像u是否需要处理，如果需要，执行步骤b到c，否则直接输出u；

C.对滤波结果进行后处理或直接输出结果。 

2.根据权利要求1所述的医学图像降噪与增强的处理方法，其特征在于当所述的多尺度模型中的尺度数量为1时，其特征为，

A.将图像I作为输入，计算一致性检测κ，使用预设阈值或自动阈值检测算法，对κ进行噪音区域与结构体区域分割，获得的阈值作为扩散滤波阈值T，将图像I分为噪音区域和结构体区域，

B.使用方向性各向异性滤波器，在噪音区域进行N次扩散滤波，法线和切线方向上扩散系数g₁和g₂定义如下为其中m₁＝2， 

在结构体区域进行S次平滑扩散滤波，法线和切线方向上扩散系数g₁和g₂定义为其中α＝0.01，N和S为非负整数，

C.对滤波结果进行后处理或直接输出结果。

3.根据权利要求1或2所述的医学图像降噪与增强的处理方法，其特征在于对滤波图像的后处理包括对比度增强。

4.根据权利要求1或2所述的医学图像降噪与增强的处理方法，其特征在于对滤波图像的后处理包括原始输入图像融合。

5.根据权利要求1所述的医学图像降噪与增强的处理方法，其特征在于在使用多尺度模型进行重建之前，对当前尺度下的殘值图像也进行方向性的各项异性滤波。

6.根据权利要求1所述的医学图像降噪与增强的处理方法，其特征在于所述的多尺度模型为拉普拉斯金字塔模型。 

7.根据权利要求1所述的医学图像降噪与增强的处理方法，其特征在于判断在某一尺度上的图像是否需要进行处理由用户预先设置。 

Images