CN101930599B - 一种医学图像增强方法及系统 - Google Patents

Abstract

本发明属于医学图像处理技术领域，具体为一种边缘保留的医学图像增强方法及系统。本发明应用基于梯度适应的滤波器将图像平滑，来估计图像的每个像素点的梯度能量和噪音能量，即计算一个滤波误差能量函数。然后，应用动态规划的办法，将所得到的误差能量函数再量化到四个不同的级别。将所得到量化后的误差能量和图像的纹理特征节后，构造一组量化的上下文。最后依据不同的量化上下文，应用回归分析的方法，构造不同的参数的滤波器，从而实现对医学图像的增强处理。

Description

一种医学图像增强方法及系统

技术领域

本发明属于医学图像处理技术领域，具体涉及一种对医学图像进行去噪处理的方法及系统。

背景技术

图像去噪技术面临最主要的挑战就是在图像噪音被消除的同时，一些重要的图像细节或边缘信息也被弱化。这使得边缘保留滤波器的设计成为医学图像处理一个研究热点。传统的边缘保留滤波器应用数学形态学算子，可以在去图像加性或乘性噪音的同时，又使边缘信息得到增强，如中值滤波器等。然而，它们依然无法克服一些弱边界被弱化的问题。为了解决这一问题，本发明应用信息论中的上下文量化技术和自适应回归分析方法，设计一个非局部的边缘保留滤波器，在保持重要的边缘信息同时，图像的噪音得到最大限度的去除。

医学图像由于临床诊断的要求通常是图像的空间和灰度分辨率较高，然而其图像质量经常由于成像过程所产生的噪音污染而下降，例如超声图像、低剂量的CT图像及X-ray图像。因此，在将图像发送给临床医生做诊断前，通常医学图像需要降噪或增强等预处理。通常的增强办法是在一个固定大小的局部窗口内，应用一个空间滤波器平滑输入的噪音图像。但是这样的处理会严重破坏图像本身所具有的纹理特征和边缘细节，而且图像的对比度会下降。

为了解决这一实际问题，近年来很多边缘保留的滤波器得到广泛地研究并应用到医学图像增强的软件开发上。例如，中值滤波器、双边滤波器，各项异性滤波器等被广泛地应用到医疗影像设备中，但是他们大部分只对高斯噪音或乘性噪音比较有效，对解决具有混合型噪音的医学图像增强问题效果并不理想，例如MRI图像的增强问题。

本发明主要涉及设计具有一个边缘保留的、自适应的滤波器，在去除图像噪音的同时，图像的边缘细节和纹理同时得到增强。不像传统的噪音估计方法，应用信息论中的上下文量化技术来估计噪音不依赖于噪音信号和图像信号本身，这样的估计方法有很强的鲁棒性，适用于所有的噪音模型未知的去噪问题。

综上所述，现有图像增强技术都无法适用于大部分医学图像去噪问题，不同的医学图像去噪问题需要不同的方法来解决，这给医学图像增强预处理的软件开发带来很大的难度和开发成本。

发明内容

本发明的目的在于提供一种在去除图像噪音的同时，图像的边缘细节和纹理同时得到保留的医学图像增强处理方法及系统。

本发明考虑到信息论中的上下文量化技术来解决医学图像噪音模型的估计问题。首先应用基于梯度的滤波器将图像平滑，再根据得到的平滑图像，计算滤波误差能量函数，也就是估计每个像素点上的高频信息，包括像素点周围的边缘信息和噪音。然后，应用动态规划的办法，将所得到的误差能量函数量化到四个不同的级别。随后，将所得到的量化后的误差能量和图像的纹理特征后，构造一组量化的上下文。最后依据不同的量化上下文，应用回归分析的方法针对不同的上下文模型，构造不同的参数的滤波器，实现一个自适应的滤波器。

对于给定的一幅图像I，本发明提出的医学图像增强处理方法，包括以下步骤：

A.设计一种基于梯度的图像滤波器或预测器(Gradient Adjusted Predictor)，将该滤波器作用到图像I，得到滤波图像

B.基于得到的滤波图像对每个像素点计算误差估计能量函数Δ；

C.根据得到的误差估计能量函数Δ，构造一个4级量化器Q(Δ)∈{0，1，2，3}；

D.根据量化后的Q(Δ)、原始图像I和滤波图像

构造一个上下文量化器Q(C)；

E.针对每个量化后的上下文C_Q，应用回归分析(Regression Analysis)方法求解公式(1)中的滤波系数波器系数b_k和α，在一个菱形的窗口构造一个滤波器f(x|C_Q)；

$\underset{k=1}{\overset{12}{\mathrm{\Σ}}}{b}_k{x}_k+\mathrm{\α}=y---\left(1\right)$

F.应用滤波器f(x|C_Q)对图像I进行去噪滤波。

本发明中，所述A步骤中设计滤波器，包括以下步骤：

A1.对当前图像I计算其水平方向梯度值d_h和垂直方向梯度值d_v：d_h和d_v分别为：

d_h＝(|I(i，j-2)-I(i，j-1)|+|I(i-1，j-1)-I(i-1，j)|

+|I(i-1，j)-I(i-1，j+1)|+|I(i，j+1)-I(i，j+2)|                    (2)

+|I(i+1，j)-I(i+1，j+1)|+|I(i+1，j-1)-I(i+1，j)|)/2

d_v＝(|I(i-1，j-1)-I(i，j-1)|+|I(i-2，j)-I(i-1，j)|

+|I(i-2，j+1)-I(i-1，j+1)|+|I(i-1，j+1)-I(i，j+1)|

+|I(i+1，j)-I(i+2，j)|+|I(i+1，j-1)-I(i+2，j-1)|)/2

其中这里，I(i-2，j-2)、I(i-2，j-1)…I(i+2，j+1)、I(i+2，j+2)是图像I在以(i，j)为中心的5×5邻域中(如图2所示)的像素点灰度值。

A2.根据得到的水平梯度和垂直梯度，按公式(3)设计一个基于梯度的滤波器，并作用于图像I，由此得到滤波图像

如果d_v(i，j)-d_h(i，j)＞C₁(经过点(i，j)存在水平方向强边缘：也就是d_v(i，j)-d_h(i，j)＞C₁，且有C₁＞C₂＞C₃，垂直方向的梯度值远大于水平方向的梯度值，则存在一个水平方向的强边缘)

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)}{2};$

否则如果d_v(i，j)-d_h(i，j)＜-C₁(经过点(i，j)存在垂直方向强边缘)

$\hat{I}(i,j)=\frac{I(i-1,j)+I(i+1,j)}{2};$

否则

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)+I(i-1,j)+I(i+1,j)}{4}$

$+\frac{I(i-1,j+1)-I(i-1,j-1)}{8}+\frac{I(i+1,j-1)-I(i+1,j+1)}{8};$

如果d_v(i，j)-d_h(i，j)＞C₂(经过点(i，j)存在水平方向边缘，也就是d_v(i，j)-d_h(i，j)＞C₂，且有C₁＞C₂＞C₃，垂直方向的梯度值大于水平方向的梯度值，则存在一个水平方向的边缘)

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i,j-1)+I(i,j+1)}{4};$

否则，如果d_v(i，j)-d_h(i，j)＞C₃(经过点(i，j)存在水平方向弱边缘，也就是d_v(i，j)-d_h(i，j)＞C₂，且有C₁＞C₂＞C₃，垂直方向的梯度值略微大于水平方向的梯度值，则存在一个较弱水平方向的边缘)

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i,j-1)+I(i,j+1)}{8};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₂(经过点(i，j)存在垂直方向边缘)

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i-1,j)+I(i+1,j)}{4};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₃(经过点(i，j)存在垂直方向弱边缘)

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i-1,j)+I(i+1,j)}{8};$

判断结束

判断结束；

本发明中，所述B步骤包括以下步骤：

B1.基于图像I和滤波图像

计算滤波误差图像

B2.根据得到的水平梯度和垂直梯度：d_h和d_v以及滤波误差图像g，计算误差能量值：

Δ＝d_h+d_v+(|g(i-1，j)|+|g(i+1，j)|)/2                          (4)

+(|g(i，j-1)|+|g(i，j+1)|)/2

本发明中，所述C步骤包括以下步骤：

C1.基于所得到的误差能量Δ值和滤波误差图像g，应用动态规划算法使得以下条件熵的值最小：

$-\underset{d=0}{\overset{3}{\mathrm{\&Sigma;}}}\underset{q_d\&le;\mathrm{\&Delta;}<q_{d+1}}{\mathrm{\&Sigma;}}p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})\log p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})---\left(5\right)$

p(g|q_d≤Δ＜q_d+1)是当前像素点在图像g相对于误差能量Δ得到的条件概率，q_d(d＝0，...，4)是个整数，且有0＝q₁＜q₂＜q₃＜q₄＝∞，构成了误差能量Δ的量化区间：

$\mathrm{\&Delta;}\&Element;{\&cup;}_{d=0}^3[q_d,q_{d+1})---\left(6\right)$

C2.将误差能量Δ，也就是当前像素点的误差能量函数，量化到正整数域中的四个区间里，即Q(Δ)为：

$Q\left(\mathrm{\&Delta;}\right)=\left\{\begin{array}{cc}0,& q_0\&le;\mathrm{\&Delta;}<q_1\\ 1,& q_1\&le;\mathrm{\&Delta;}<q_2\\ 2,& q_2\&le;\mathrm{\&Delta;}<q_3\\ 3,& q_3\&le;\mathrm{\&Delta;}<q_4\end{array}\right.---\left(7\right)$

Q(Δ)∈{0，1，2，3}，这里Q(Δ)也可称作标量量化器(Scalar Quantizer)，{0，1，2，3}是个包含4个整数的集合。

本发明中，所述D步骤包括以下步骤：

D1.如果当前像素点的坐标值为x＝(i，j)，取一个3×3的图像邻域，选取邻域中的像素点来抽取图像纹理特征：

x₀＝(i-1，j-1)

x₁＝(i-1，j)

x₂＝(i-1，j+1)

x₃＝(i，j-1)                        (8)

x₄＝(i，j+1)

x₅＝(i+1，j-1)

x₆＝(i+1，j)

x₇＝(i+1，j+1)

D2.在以上3×3的图像邻域内，抽取纹理特征S＝s₇s₆s₅s₄s₃s₂s₁s₀，

其中：

$s_k=\left\{\begin{array}{cc}0& I\left(x_k\right)\&GreaterEqual;\hat{I}(i,j)\\ 1& I\left(x_k\right)<\hat{I}(i,j)\end{array}\right.0\&le;k\&le;7---\left(9\right)$

D3.纹理特征S为大小8位，量化的误差能量Q(Δ)为0～3的数字，大小为2位的数字，将Q(Δ)和S做二进制位数合并，生成一个大小10位的数字，可以把它看作一个量化的上下文C_Q，也就是得到一组1024个数字，即0～1023，把这1024个数字标记成一个量化的上下文集合C＝{C_Q|C_Q＝0，...1023}

本发明中，所述E步骤包括以下步骤：

E1.在图像I中，针对每个像素点y，在其周围，取一个如下的菱形邻域R(y)，即：

		x1
						x2	x3	x4
x5	x6	y	x7	x8
						x9	x10	x11
		x12

(10)

E2.根据D步骤所得到1024个量化的上下文集合C，得到像素点集满足Y＝{yC_Q(y)∈C，C_Q(y)＝C_Q}，且有R(Y)＝{R(y)|C_Q(y)∈C，C_Q(y)＝C_Q}，设对所有R(Y)，像素点灰度值y和其菱形邻域内中的其它12个像素点灰度值x_k都满足如下关系：

$\underset{k=1}{\overset{12}{\mathrm{\&Sigma;}}}{b}_k{x}_k+\mathrm{\&alpha;}=y$

则基于得到的R(Y)，通过回归分析我们可将系数b_k和α估计出来，估计到的系数b_k和α也就是f(x|C_Q)滤波器的滤波系数。

基于上述医学图像增强处理方法，本发明的医学图像增强处理系统包括：一种基于梯度的图像滤波器或预测器(Gradient Adjusted Predictor)，一个4级量化器Q(Δ)∈{0，1，2，3}，一个上下文量化器Q(C)，

一个滤波器f(x|C_Q)。其中，各个部分的功能、作用见上面相应部分的描述。

附图说明

图1是本发明实施例二维医学图像增强的流程示意图。

图2是本发明实施例基于梯度适应的滤波器的窗口大小。

图3是本发明实施例抽取当前像素点纹理特征的邻域。

图4是本发明实施例空域滤波器所采用的菱形邻域。

图5是本发明实施例对胸片X光图像一块区域的增强处理对比。其中，(a)含有噪音的胸部X线图像的区域，(b)增强处理后的胸部X线图像的区域。

图6是本发明实施例对MRI脑部图像的增强处理对比。其中，(a)含有噪音的脑部MRI图像，(b)增强处理后的脑部MRI图像。

具体实施方式

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

图1是本发明涉及的二维医学图像增强的流程示意图，它描述了从输入初始噪音图像62，然后计算梯度图像，并基于得到的水平梯度和垂直梯度设计自适应的滤波器64，应用滤波器得到一个初始平滑图像

步骤66结合原始图像I和初始平滑图像

估计每个像素点的高频信息量或滤波误差能量函数Δ；步骤68根据最小条件熵准则将误差能量函数量化成1-4级；步骤70基于初始噪音图像I和初始平滑图像

在每个像素点3×3邻域上抽取一个大小8位的图像纹理特征。步骤72则将纹理特征和量化的误差能量函数结合成一组量化的上下文特征，然后步骤74根据这组量化的上下文特征，应用回归分析来估计在每组上下特征下的滤波器系数，最后步骤76应用得到的1024组滤波器和相应步骤72得到的上下文特征将图像增强。具体步骤如下：

1.计算噪音图像I的水平梯度和垂直梯度，基于水平和垂直梯度设计一个基于梯度适应的(Gradient Adjusted Predictor)滤波器，将该滤波器作用到图像I，得到滤波图像

2.基于得到的滤波图像

对每个像素点计算误差估计能量函数Δ；

3.根据得到的误差估计能量函数Δ，利用求解最小条件熵，设计一个误差能量函数的4级量化器Q(Δ)∈{0，1，2，3}；

4.根据量化后的Q(Δ)、原始图像I和滤波图像构造一个上下文量化器Q(C)；

5.针对每个量化后的的上下文C_Q，应用回归分析，在一个菱形的窗口构造一个滤波器f(x|C_Q)；

6.应用滤波器f(x|C_Q)对图像I进行去噪滤波。

上述处理程序包含三大部分，图像噪音的估计，图像的上下文模型的选择，根据量化的上下文设计滤波器。图像噪音的估计通常是设计一个自适应滤波器的必要过程，不同的噪音模型，估计的方法也不同。本发明区别于其他技术最关键的一点，估计出的噪音量并不直接应用于滤波器的设计，而是用来进一步地描述和量化图像像素点周围的特征的上下文模型，根据信息论中通用信源编码理论可以推断出：这样的方法好处在于，在不知道图像中的噪音模型情况下，就可以找到一个和噪音本身非常近似的统计模型。步骤1包含了水平梯度和垂直梯度的计算以及基于梯度适应的滤波器的设计：

(1)对当前图像I计算其水平方向梯度值和垂直方向梯度值：d_h和d_v分别为：

d_h＝(|I_(i，j-2)-I_(i，j-1)|+|I_(i-1，j-1)-I_(i-1，j)|+|I_(i-1，j)-I_(i-1，j+1)|

+|I_(i，j+1)-I_(i，j+2)|+|I_(i+1，j)-I_(i+1，j+1)|+|I_(i+1，j-1)-I_(i+1，j)|)/2

d_v＝(|I_(i-1，j-1)-I_(i，j-1)|+|I_(i-2，j)-I_(i-1，j)|+|I_(i-2，j+1)-I_(i-1，j+1)|

+|I_(i-1，j+1)-I_(i，j+1)|+|I_(i+1，j)-I_(i+2，j)|+|I_(i+1，j-1)-I_(i+2，j-1)|)/2

(2)根据得到的水平梯度和垂直梯度，设计一个基于梯度的滤波器并作用于图像I，由此得到滤波图像

如果d_v(i，j)-d_h(i，j)＞C₁(经过点(i，j)存在水平方向强边缘)

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)}{2};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₁(经过点(i，j)存在垂直方向强边缘)

$\hat{I}(i,j)=\frac{I(i-1,j)+I(i+1,j)}{2};$

否则

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)+I(i-1,j)+I(i+1,j)}{4}$

$+\frac{I(i-1,j+1)-I(i-1,j-1)}{8}+\frac{I(i+1,j-1)-I(i+1,j+1)}{8};$

如果d_v(i，j)-d_h(i，j)＞C₂(经过点(i，j)存在水平方向边缘)

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i,j-1)+I(i,j+1)}{4};$

否则，如果d_v(i，j)-d_h(i，j)＞C₃(经过点(i，j)存在水平方向弱边缘)

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i,j-1)+I(i,j+1)}{8};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₂(经过点(i，j)存在垂直方向边缘)

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i-1,j)+I(i+1,j)}{4};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₃(经过点(i，j)存在垂直方向弱边缘)

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i-1,j)+I(i+1,j)}{8};$

判断结束

判断结束；

梯度适应的滤波器是根据不同边缘强度和方向来设计一个不同平滑强度的滤波器，其实它也是一种边缘保留的滤波器并且属于各项异性滤波器，C₁，C₂和C₃为该滤波器的预设参数，分别来检测是否存在水平或垂直方向的强边缘、边缘和弱边缘。例如，可取C₁＝80，C₂＝32和C₃＝8。图2说明了该边缘保留滤波器所用的图像滤波窗口大小或像素点邻域。

通常来说，图像的噪音可以应用一个自适应滤波器来估计，换句话说就是

g也可以看成是一定程度上的噪音图像，但同时也包含边缘信息和纹理特征，如果将当前像素点在噪音图像g四邻域点噪音值的绝对值和梯度图像相加，我们可以近似地看成一个能量函数Δ：

Δ＝d_h+d_v+(|g_i-1，j|+|g_i+1，j|)/2+(|g_i，j-1|+|g_i，j+1|)/2

得到能量函数Δ以后，我们需要对噪音图像g的统计分布进一步估计，这样才可以近似描述噪音模型，根据贝叶斯法则和通用信源编码理论，对图像信号g最为近似的估计就是用以一个近似的模型求解最小化的条件熵值，也就是最小求解Kullback-Leiber距离：

$-\underset{d=0}{\overset{3}{\mathrm{\&Sigma;}}}\underset{q_d\&le;\mathrm{\&Delta;}<q_{d+1}}{\mathrm{\&Sigma;}}p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})\log p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})$

也就是求解一个量化器Q(Δ)∈{0，1，2，3}，由于条件熵是一个严格的凸函数，故存在唯一的极小值，可以通过动态规划的办法求解。

梯度和量化的噪音模型实际上并不能完全描述图像的上下文模型，需要加入图像的纹理特征，也就是处理滤波后的图像和周围像素点在原始图像之间的关系，步骤70在以上3×3的邻域内，抽取纹理特征S＝s₇s₆s₅s₄s₃s₂s₁s₀：

$s_k=\left\{\begin{array}{cc}0& I\left(x_k\right)\&GreaterEqual;\hat{I}(i,j)\\ 1& I\left(x_k\right)<\hat{I}(i,j)\end{array}\right.0\&le;k\&le;7$

图3说明了抽取该纹理特征所在的3×3的图像邻域且有：

x₀＝(i-1，j-1)

x₁＝(i-1，j)

x₂＝(i-1，j+1)

x₃＝(i，j-1)

x₄＝(i，j+1)

x₅＝(i+1，j-1)

x₆＝(i+1，j)

x₇＝(i+1，j+1)

步骤72将上述得到的纹理信息大小为一个大小为8位(0～255)的正整数，故纹理特征共有256个，这些特征再和步骤68得到的量化误差能量函数融合起来，共生成1024个关联的上下文模型或参数。这样得到的上下文模型能最大程度地描述图像信号模型和噪音信号模型，不同的模型可以应用不同的滤波系数来处理。

步骤74是根据72得到上下文关联模型，应用统计的方法将图像中所有像素点都归类成不同的统计模型或上下文模型，也就是将每个像素点归类成1024个模型中的一类，因此我们会得到1024组数据，每一组数据都包含若干个当前像素点y的灰度值以及它的一个菱形邻域中周围像素点的值，图4说明了本发明所采用的菱形邻域。对每一组菱形邻域的数据，可以有如下假定：

$\underset{k=1}{\overset{12}{\mathrm{\&Sigma;}}}{b}_k{x}_k+\mathrm{\&alpha;}=y$

b_k是估计的滤波器系数和α是估计出的近似噪音，x_k是当前像素点y菱形邻域内的像素点灰度值。应用基于最小二乘法的回归分析方法可以估计出滤波器系数b_k和模型噪音α的值。

在基于所有的1024组模型的滤波器参数估计完成以后，结合步骤72所得到的量化的上下文模型，步骤76应用1024组滤波器对图像进行滤波处理，完成本医学图像增强方法。

图5-6是本发明的实施例对数字X线胸片的局部区域和脑部MRI图像的增强处理结果，从处理前后的对比结果来看，无论是处理后的数字X线胸片的局部区域还是脑部MRI图像，噪音水平明显下降，但图像有效区域中纹理、细节边缘和对比度都得到增强，说明本发明实现了一种有效的边缘保留的医学图像增强技术。

Claims (2)

1.一种医学图像增强方法，其特征在于，给定一幅图像I，增强的具体步骤如下：

步骤A.设计一种基于梯度的图像滤波器，将该滤波器作用到图像I，得到滤波图像

步骤B.基于得到的滤波图像

对每个像素点计算误差估计能量值Δ；

步骤C.根据得到的误差估计能量值Δ，构造一个4级量化器Q(Δ)∈{0，1，2，3}；

步骤D.根据量化后的Q(Δ)、原始图像I和滤波图像构造一个上下文量化器Q(C)；

步骤E.针对每个量化后的上下文C_Q，在一个菱形的窗口构造一个滤波器f(x|C_Q)；

步骤F.应用滤波器f(x|C_Q)对图像I进行去噪滤波；

其中，所述步骤A中，设计滤波器步骤为：

A1.对当前图像I计算其水平方向梯度值d_h＝d_h(i，j)和垂直方向梯度值d_v＝d_v(i，j)：

d_h＝(|I(i，j-2)-I(i，j-1)|+|I(i-1，j-1)-I(i-1，j)|

+|I(i-1，j)-I(i-1，j+1)|+|I(i，j+1)-I(i，j+2)|

+|I(i+1，j)-I(i+1，j+1)|+|I(i+1，j-1)-I(i+1，j)|)/2

                                                        (2)

d_v＝(|I(i-1，j-1)-I(i，j-1)|+|I(i-2，j)-I(i-1，j)|

+|I(i-2，j+1)-I(i-1，j+1)|+|I(i-1，j+1)-I(i，j+1)|

+|I(i+1，j)-I(i+2，j)|+|I(i+1，j-1)-I(i+2，j-1)|)/2

这里，I(i-2，j-2)、I(i-2，j-1)…I(i+2，j+1)、I(i+2，j+2)是图像I在以(i，j)为中心的5×5邻域中的像素点灰度值；

A2.根据得到的水平梯度和垂直梯度，按公式(3)设计一个基于梯度的图像滤波器：

如果d_v(i，j)-d_h(i，j)＞C₁，即经过点(i，j)存在水平方向强边缘；

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)}{2};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₁，即经过点(i，j)存在垂直方向强边缘；

$\hat{I}(i,j)=\frac{I(i-1,j)+I(i+1,j)}{2};$

否则，

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)+I(i-1,j)+I(i+1,j)}{4}$

$+\frac{I(i-1,j+1)-I(i-1,j-1)}{8}+\frac{I(i+1,j-1)-I(i+1,j+1)}{8};$

如果d_v(i，j)-d_h(i，j)＞C₂，即经过点(i，j)存在水平方向边缘；

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i,j-1)+I(i,j+1)}{4};---\left(3\right)$

否则，如果d_v(i，j)-d_h(i，j)＞C₃，即经过点(i，j)存在水平方向弱边缘；

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i,j-1)+I(i,j+1)}{8};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₂，经过点(i，j)存在垂直方向边缘；

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i-1,j)+I(i+1,j)}{4};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₃，即经过点(i，j)存在垂直方向弱边缘；

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i-1,j)+I(i+1,j)}{8};$

其中，C₁、C₂和C₃为该基于梯度的图像滤波器的预设参数，分别用来检测是否存在水平或垂直方向的强边缘、边缘和弱边缘；

所述步骤B中对每个像素点计算误差估计能量值Δ的步骤为：

B1.基于图像I和滤波图像

计算滤波误差图像

B2.根据得到的水平梯度d_h和垂直梯度d_v以及滤波误差图像g，按式(4)计算误差估计能量值Δ：

Δ＝d_h+d_v+(|g(i-1，j)|+|g(i+1，j)|)/2

+(|g(i，j-1)|+|g(i，j+1)|)/2                (4)

所述步骤C中，求4级量化器Q(Δ)∈{0，1，2，3}的步骤为：

C1.基于所得到的误差估计能量值Δ和滤波误差图像g，应用动态规划算法使得以下条件熵的值最小：

$-\underset{d=0}{\overset{3}{\mathrm{\&Sigma;}}}\underset{q_d\&le;\mathrm{\&Delta;}<q_{d+1}}{\mathrm{\&Sigma;}}p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})\log p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})---\left(5\right)$

p(g|q_d≤Δ＜q_d+1)是当前像素点在滤波误差图像g相对于误差估计能量值Δ得到的条件概率，q_d是个整数，d＝0，…，4，且有0＝q₀＜q₁＜q₂＜q₃＜q₄＜∞，构成误差估计能量值Δ的量化区间：

$\mathrm{\&Delta;}\&Element;{\&cup;}_{d=0}^3[q_d,q_{d+1})---\left(6\right)$

C2.将误差估计能量值Δ，量化到正整数域中的四个区间里，即Q(Δ)为：

$Q\left(\mathrm{\&Delta;}\right)=\left\{\begin{array}{cc}0,& q_0\&le;\mathrm{\&Delta;}{q}_1\\ 1,& q_1\&le;\mathrm{\&Delta;}<q_2\\ 2,& q_2\&le;\mathrm{\&Delta;}<q_3\\ 3,& q_3\&le;\mathrm{\&Delta;}<q_4\end{array}\right.---\left(7\right)$

Q(Δ)∈{0，1，2，3}，Q(Δ)也称作标量量化器，{0，1，2，3}是4个整数的集合；

所述步骤D中，构造上下文量化器Q(C)步骤为：

D1.设当前像素点的坐标值为x＝(i，j)，取式(8)所示的3×3的图像邻域：

x₀＝(i-1，j-1)

x₁＝(i-1，j)

x₂＝(i-1，j+1)

x₃＝(i，j-1)

x₄＝(i，j+1)

x₅＝(i+1，j-1)                                   (8)

x₆＝(i+1，j)

x₇＝(i+1，j+1)

D2.在以上3×3的图像邻域内，抽取纹理特征S＝s₇s₆s₅s₄s₃s₂s₁s₀，

其中：

$s_k=\left\{\begin{array}{cc}0& I\left(x_k\right)\&GreaterEqual;\hat{I}(i,j)\\ 1& I\left(x_k\right)<\hat{I}(i,j)\end{array}\right.0\&le;k\&le;7---\left(9\right)$

D3.纹理特征S为大小8位，4级量化器Q(Δ)为取值0～3的、大小为2位的数字，将Q(Δ)和S做二进制位数合并，生成一个大小为10位的数字，把它看作一个量化的上下文C_Q，得到一组1024个数字，即0～1023，把这1024个数字标记成一个量化的上下文集合C＝{C_Q|C_Q＝0，…1023}；

所述步骤E中，构造滤波器f(x|C_Q)的步骤为：

E1.在图像I中，针对每个像素点y，在其周围，选择一个如下的菱形邻域R(y)：

E2.根据步骤D所得到的量化的上下文集合C，得到像素点集满足Y＝{y|C_Q(y)∈C，C_Q(y)＝C_Q}，且有R(Y)＝{R(y)|C_Q(y)∈C，C_Q(y)＝C_Q}，设对所有R(Y)，像素点灰度值y和其菱形邻域内中的其它12个像素点灰度值x_k都满足如下关系：

$\underset{k=1}{\overset{12}{\mathrm{\&Sigma;}}}{b}_k{x}_k+\mathrm{\&alpha;}=y---\left(1\right)$

则通过回归分析方法，估计出系数b_k和α，该系数b_k和α即为滤波器f(x|C_Q)的滤波系数。

2.一种医学图像增强系统，其特征在于包括：一种基于梯度的图像滤波器、一个4级量化器Q(Δ)∈{0，1，2，3}、一个上下文量化器Q(C)和一个滤波器f(x|C_Q)；其中：

所述基于梯度的图像滤波器的构造如下：

A1.对当前图像I计算其水平方向梯度值d_h＝d_h(i，j)和垂直方向梯度值d_v＝d_v(i，j)：

d_h＝(|I(i，j-2)-I(i，j-1)|+|I(i-1，j-1)-I(i-1，j)|

+|I(i-1，j)-I(i-1，j+1)|+|I(i，j+1)-I(i，j+2)|

+|I(i+1，j)-I(i+1，j+1)|+|I(i+1，j-1)-I(i+1，j)|)/2

                                                         (2)

d_v＝(|I(i-1，j-1)-I(i，j-1)|+|I(i-2，j)-I(i-1，j)|

+|I(i-2，j+1)-I(i-1，j+1)|+|I(i-1，j+1)-I(i，j+1)|

+|I(i+1，j)-I(i+2，j)|+|I(i+1，j-1)-I(i+2，j-1)|)/2

这里，I(i-2，j-2)、I(i-2，j-1)…I(i+2，j+1)、I(i+2，j+2)是图像I在以(i，j)为中心的5×5邻域中的像素点灰度值；

A2.根据得到的水平梯度和垂直梯度，按公式(3)设计得到基于梯度的图像滤波器：

如果d_v(i，j)-d_h(i，j)＞C₁，即经过点(i，j)存在水平方向强边缘；

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)}{2};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₁，经过点(i，j)存在垂直方向强边缘；

$\hat{I}(i,j)=\frac{I(i-1,j)+I(i+1,j)}{2};$

否则

$\hat{I}(i,j)=\frac{I(i,j-1)+I(i,j+1)+I(i-1,j)+I(i+1,j)}{4}$

$+\frac{I(i-1,j+1)-I(i-1,j-1)}{8}+\frac{I(i+1,j-1)-I(i+1,j+1)}{8};$

如果d_v(i，j)-d_h(i，j)＞C₂，即经过点(i，j)存在水平方向边缘；

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i,j-1)+I(i,j+1)}{4};---\left(3\right)$

否则，如果d_v(i，j)-d_h(i，j)＞C₃即经过点(i，j)存在水平方向弱边缘；

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i,j-1)+I(i,j+1)}{8};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₂，经过点(i，j)存在垂直方向边缘；

$\hat{I}(i,j)=\frac{\hat{I}(i,j)}{2}+\frac{I(i-1,j)+I(i+1,j)}{4};$

否则，如果d_v(i，j)-d_h(i，j)＜-C₃，即经过点(i，j)存在垂直方向弱边缘；

$\hat{I}(i,j)=\frac{3\hat{I}(i,j)}{4}+\frac{I(i-1,j)+I(i+1,j)}{8};$

其中，C₁、C₂和C₃为该基于梯度的图像滤波器的预设参数，分别用来检测是否存在水平或垂直方向的强边缘、边缘和弱边缘；

所述4级量化器Q(Δ)∈{0，1，2，3}的构造如下：

B1.基于图像I和滤波图像

计算滤波误差图像

B2.根据得到的水平梯度d_h和垂直梯度d_v以及滤波误差图像g，按式(4)计算误差估计能量值Δ：

Δ＝d_h+d_v+(|g(i-1，j)|+|g(i+1，j)|)/2

+(|g(i，j-1)|+|g(i，j+1)|)/2                       (4)

C1.基于所得到的误差估计能量值Δ和滤波误差图像g，应用动态规划算法使得以下条件熵的值最小：

$-\underset{d=0}{\overset{3}{\mathrm{\&Sigma;}}}\underset{q_d\&le;\mathrm{\&Delta;}{<q}_{d+1}}{\mathrm{\&Sigma;}}p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})\log p\left(g\right|q_d\&le;\mathrm{\&Delta;}<q_{d+1})---\left(5\right)$

p(g|q_d≤Δ＜q_d+1)是当前像素点在滤波误差图像g相对于误差估计能量值Δ得到的条件概率，q_d是个整数，d＝0，…，4，且有0＝q₀＜q₁＜q₂＜q₃＜q₄＜∞，构成误差估计能量值Δ的量化区间：

$\mathrm{\&Delta;}\&Element;{\&cup;}_{d=0}^3[q_d,q_{d+1})---\left(6\right)$

C2.将误差估计能量值Δ，也就是当前像素点的误差能量函数，量化到正整数域中的四个区间里，即Q(Δ)为：

$Q\left(\mathrm{\&Delta;}\right)=\left\{\begin{array}{cc}0,& q_0\&le;\mathrm{\&Delta;}{q}_1\\ 1,& q_1\&le;\mathrm{\&Delta;}<q_2\\ 2,& q_2\&le;\mathrm{\&Delta;}<q_3\\ 3,& q_3\&le;\mathrm{\&Delta;}<q_4\end{array}\right.---\left(7\right)$

Q(Δ)∈{0，1，2，3}，Q(Δ)也称作标量量化器，{0，1，2，3}是4个整数的集合；所述上下文量化器Q(C)的构造如下：

D1.设当前像素点的坐标值为x＝(i，j)，取式(8)所示的3×3的图像邻域：

x₀＝(i-1，j-1)

x₁＝(i-1，j)

x₂＝(i-1，j+1)

x₃＝(i，j-1)                          (8)

x₄＝(i，j+1)

x₅＝(i+1，j-1)

x₆＝(i+1，j)

x₇＝(i+1，j+1)

D2.在以上3×3的图像邻域内，抽取纹理特征S＝s₇s₆s₅s₄s₃s₂s₁s₀，

其中：

$s_k=\left\{\begin{array}{cc}0& I\left(x_k\right)\&GreaterEqual;\hat{I}(i,j)\\ 1& I\left(x_k\right)<\hat{I}(i,j)\end{array}\right.0\&le;k\&le;7---\left(9\right)$

D3.纹理特征S为大小8位，4级量化器Q(Δ)为取值0～3的、大小为2位的数字，将Q(Δ)和S做二进制位数合并，生成一个大小为10位的数字，把它看作一个量化的上下文C_Q，得到一组1024个数字，即0～1023，把这1024个数字标记成一个量化的上下文集合C＝{C_Q|C_Q＝0，…1023}；

所述滤波器f(x|C_Q)的构造如下：

E1.在图像I中，针对每个像素点y，在其周围，取式(10)所示的菱形邻域R(y)：

E2.根据步骤D3所得到的量化的上下文集合C，得到像素点集满足Y＝{y|C_Q(y)∈C，C_Q(y)＝C_Q}，且有R(Y)＝{R(y)|C_Q(y)∈C，C_Q(y)＝C_Q}，设对所有R(Y)，像素点灰度值y和其菱形邻域内中的其它12个像素点灰度值x_k都满足如下关系：

$\underset{k=1}{\overset{12}{\mathrm{\&Sigma;}}}{b}_k{x}_k+\mathrm{\&alpha;}=y---\left(1\right)$

则通过回归分析方法，估计出系数b_k和α，该系数b_k和α即为滤波器f(x|C_Q)的滤波系数。

Images