CN1804903A - 基于尺度空间分解与重构的x射线图像均衡显示处理方法 - Google Patents

Abstract

本发明公开了一种基于尺度空间分解与重构的X射线图像均衡显示处理方法，该方法的步骤包括：首先将输入图像进行归一化和对数变换，然后选用线性或者非线性图像滤波方法进行图像的尺度空间分解，然后对分解产生的不同尺度下的图像信息进行加权，然后将加权后的信号重构，将重构的信号进行反归一化及反对数处理，得到均衡后的图像。在分解过程中，可以采用空域平滑滤波算子，如高斯滤波算子，高斯微分滤波算子进行线性尺度空间分解，或者采用中值曲率驱动方程进行非线性滤波及分解。该方法创造性地从尺度空间分解的角度出发，进行图像均衡处理，与人眼的视觉模型相符合，因而具有良好的处理效果。

Description

基于尺度空间分解与重构的X射线图像均衡显示处理方法

                             技术领域

本发明属于图像处理领域，涉及X射线图像的处理，尤其涉及X射线图像显示质量的改善，具体涉及一种基于尺度空间分解与重构X射线图像的均衡显示处理方法。

                             背景技术

X射线成像是一类广泛应用于医疗诊断和工业领域的成熟技术。从20世纪80年代起，数字化的X射线成像技术CR及DR逐渐取代了传统的胶片成像技术。虽然两种技术的空间分辨率都接近胶片，但是却具有很宽的对比度范围，其动态范围可以达到16比特。对于相近密度的组织差异，有很高的辨别程度，因此，能够捕获所有的与诊断相关的信息，然而，由于X射线成像***本身的特点，一方面，显示***的MTF性能会对图像的质量造成一定程度的恶化；另一方面，数字X射线成像***中X射线散射、电子噪声、光量子噪声等各种不利因素的影响，使得数字X射线成像显示结果呈现出组织间对比度相对过大，而细节模糊，对比度较差的特点。

对于X射线图像显示而言，普通CRT显示设备，其显示能力只有8到9比特，在有限的灰度范围内，无法完全表达CR及DR图像包含的所有信息，因此，在显示过程中，通常采用查找表的方法，将16比特的灰度信息，通过调整窗宽、窗位、变换曲线等，将有限宽度灰度范围的图像信息映射到8比特的显示范围内，供观察者观测，重复调节窗宽，窗位，用这种方法，将图像的灰度信息分次，分部分地呈现给观察者，由观察者将这些信息在大脑中理解、融合，形成对整幅图像的判断。这样势必会影响观察者对图像诊断信息的正确理解和判断，并且增大了工作复杂度，再者，人眼对灰度的分辨能力有限，也给诊断带来困难。因此，如何改善X射线图像在有限灰度范围内的显示质量显得非常重要，X射线图像均衡方法的目的在于，一方面要进行动态范围压缩，抑制不包含诊断信息的大尺度信号，同时要保持或者相对放大图像中各个尺度上与诊断密切相关的细节信息的对比度在显示设备的灰度显示范围及人眼的可视范围之内，达到整体的均衡效果，为观察者提供好的显示效果，提高图像的诊断价值，从而起到更好的辅助诊断的作用。

对于已有的相关方法，有以下三种：

1.国外专利申请(Patent No.：EP1501049；US2005018894；JP2005040590-A)所述的基于频域子带分解的方法，该方法首先在频域将图像分解，提取带通信息，对各频率子带信号乘以相应的增益系数，然后重建出增强后的图像。通过调整最高频段信号的增益系数，控制图像的锐化程度；通过调整带通信号的增益系数，控制细节对比度的扩展程度；动态范围的调整可以由抑制最低频率子带信号来实现。由这三方面的参数来相对独立地独立控制图像的锐化程度、细节对比度及动态范围。但是该方法采用的固定的三级滤波对于信号的分离过于粗略，在文献“Equalized Contrast Display Processing forDigital Radiography”(Medical Imaging 2002：proceeding of SPIE Vol.4681)中提到，针对大多数医学图像而言，6到7级的分解能够对生理解剖结构特征及细微的诊断信息进行有效地分离，从而得到更好的处理结果。

2.国外专利申请(Patent No.：US5978518)所述的基于反锐化掩模的方法，该发明采用传统的线性反锐化掩模图像增强方法，借鉴Sezan et al，“Automatic AnatomicallySelective Image Enhancement in Digital Chest Radiography”(IEEE Transactions on MedicalImaging，vol.8，No.2，Jun，1989)文献中所述根据生理解剖的特点对象素所处的不同生理区域，进行选择性增强，又融合了Snyder所提出的噪声估计方法，使得锐化的同时，尽可能地不放大噪声，并且说明该方法可以推广到多级分解框架，但是该发明并没有给出在多级分解的情况下，各级增强系数之间的关系和各级滤波器之间的关系，从而仅仅停留在对图像高频细节进行增强的基础上，而缺乏对图像各个尺度信息的整体考虑。

3.另一种是多尺度边缘提取的方法。国外专利申请(Patent No.：US5805721；EP527525)中采用的基于拉普拉斯金字塔的频域分解与重构的图像均衡方法，该框架基于DOG算子的图像边缘检测技术进行图像分解。根据Marr的边缘检测理论，信号可以用一族由不同σ下的²G算子检测到的零交叉点的平面完全恢复，而当两个高斯滤波器的尺度相差一倍频程时，DOG算子是²G算子的有效近似，因此，采用这个框架可以将图像的交流分量分解为一族不同尺度下的边缘信息，针对不同的边缘信息进行增强，达到提升高频响应的目的，通过分离并抑制直流分量，实现动态范围压缩。

                                发明内容

本发明的目的在于，提出一种基于尺度空间分解与重构的X射线图像均衡显示处理方法，该方法借鉴匹配滤波的思想，从尺度空间分解的角度出发，并针对X射线图像本身的特点构造。

本发明的方法在图像均衡过程中，利用一系列尺度由大到小的滤波器，逐级抽取图像尺度信息，在抽取的过程中，当前图像大尺度信息由该级图像滤波的结果得到，而用图像与其滤波结果之差，作为该级分解保留下的尺度相对较小的信息，作为下一级待分解对象，每一级相对小的尺度信息都能被完整地保留下来，从而将图像的尺度信息从大尺度到小尺度，一层层分解出来，然后压缩与诊断信息无关的大尺度信号对比度，达到抑制动态范围的目的，根据需要，对各级小尺度信号做适当的增强或保留，本发明对各级尺度信号的增强程度随着尺度的增大而增大，且各尺度信号处理相对独立，改善细节对比度，从而实现整个图像的均衡处理结果。并且给出了各级分解过程中，增强系数与尺度的选择方法。

为了实现上述目的，本发明采取如下技术解决方案：

一种基于尺度空间分解与重构的X射线图像均衡显示处理方法，其特征在于，该方法具体的步骤是：

a.归一化处理

将输入的X射线图像通过归一化处理及对数变换得到F₀(x，y)；

b.尺度空间分解

定义算子M_t为滤波算子，其尺度参数为t，则用尺度为t₀的滤波算子M_t0对F₀(x，y)进行滤波，则将滤波后的图像F₁(x，y)表示为：

$F_1(x,y)=M_{t_0}{F}_0(x,y)$

     D₀(x，y)＝F₀(x，y)-F₁(x，y)

重复上述步骤，用尺度为t_N-1的滤波算子M_tN-1对D_N-2(x，y)进行平滑，得到F_N(x，y)，D_N-2(x，y)与F_N(x，y)之差即D_N-1(x，y)。如下所示

$F_N(x,y)=M_{t_{N-1}}{D}_{N-2}(x,y)$

     D_N-1(x，y)＝D_N-2(x，y)-F_N(x，y)

因此有：

     F₀(x，y)＝D₀(x，y)+F₁(x，y)

     D₀(x，y)＝D₁(x，y)+F₂(x，y)

     D₁(x，y)＝D₂(x，y)+F₃(x，y)

     

     D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

所以，输入的F₀(x，y)表示为：

$F_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)+D_{N-1}(x,y);$

c.重建

每一次分解过程都可以看作将待分解图像D_i-1(x，y)分解为D_i-1(x，y)的高频部分D_i(x，y)与低频部分F_i+1(x，y)的过程，保留每次分解获得的F_i+1(x，y)，以及D_N-1(x，y)，即可恢复原始图像F₀(x，y)，恢复得到的图像表示为F₀′(x，y)，因此，加上均衡效果的重建过程可表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1}$

其中，α_i(i＝1，2，…，N+1)为均衡系数，其选择原则为0＜α_i≤α_i+1≤1；

d.反归一化及反对数变换

最后将F₀′(x，y)通过反归一化及反对数变换，及灰度范围重映射得到均衡处理后的图像。

上述图像的尺度空间分解与重构的X射线图像均衡显示处理方法，其特征在于，在线性尺度空间下，M_t为空域的平滑卷积算子，则M_tF(x，y)表示为图像与平滑滤波算子的卷积，为F(x，y)*G(x，y)(其中*代表卷积操作)，其中滤波器G(x，y)的尺度由σ控制；具体包括下列步骤：

a.将输入的X射线图像通过归一化处理及对数变换得到F₀(x，y)；

b.用尺度为σ₀的尺度函数G₀与F₀进行卷积，得到F₁，F₀与F₁之差记做D₀，如下所示；

   F₁(x，y)＝F₀(x，y)*G₀(x，y)

   D₀(x，y)＝F₀(x，y)-F₁(x，y)

重复上述步骤，用尺度为σ_i的尺度函数G_i对D_i-1进行平滑，得到F_i+1，D_i-1与F_i+1之差即D_i，如下所示：

   F_i+1(x，y)＝D_i-1(x，y)*G_i(x，y)

   D_i(x，y)＝D_i-1(x，y)-F_i+1(x，y)

因此有，

   F₀(x，y)＝D₀(x，y)+F₁(x，y)

   D₀(x，y)＝D₁(x，y)+F₂(x，y)

   D₁(x，y)＝D₂(x，y)+F₃(x，y)

           

   D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

所以，输入的F₀表示为：

$F_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)+D_{N-1}(x,y);$

c.均衡显示的重建

每一次分解过程都看作将待分解图像D_i-1分解为D_i-1的高频部分D_i与低频部分F_i+1的过程。保留每次分解获得的F_i+1，以及D_N-1，即可恢复原始图像F₀，恢复得到的图像表示为F₀′；因此，加上均衡效果的重建过程表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1}$

其中，α_i(i＝1，2，…，N+1)为均衡系数，其选择原则为0＜α_i≤α_i+1≤1；

d.最后将F₀′通过反归一化及反对数变换，及灰度范围重映射得到均衡处理后的图像。

在非线性尺度空间下，采用偏微分方程(PDE)进行尺度空间分解，定义算子M_t为MCM滤波算子，方程为：

$\frac{\∂F}{\∂t}=\left|\right|\▿F\left|\right|\mathrm{div}\left(\frac{\▿F}{\left|\right|\▿F\left|\right|}\right)$

设输入的原始图像F(x，y)是上式在t＝0时刻的解，F(0，x，y)＝F₀(x，y)为方程初值，则F₀(x，y)在任意尺度t下的平滑结果可表示为M_tF₀(x，y)；则实现步骤如下：

a.将输入的X射线图像通过归一化处理及对数变换得到F₀(x，y)；

b.用尺度为t₀的尺度函数G₀与F₀进行卷积，得到F₁，F₀与F₁之差记做D₀，如下所示：

$F_1(x,y)=M_{t_0}{F}_0(x,y)$

   D₀(x，y)＝F₀(x，y)-F₁(x，y)

重复上述步骤，用尺度为t_N-1的尺度函数M_tN-1对D_N-2进行平滑，得到F_N，D_N-2与F_N之差即D_N-1(x，y)，如下所示：

$F_N(x,y)=M_{t_{N-1}}{D}_{N-2}(x,y)$

   D_N-1(x，y)＝D_N-2(x，y)-F_N(x，y)

因此有，

   F₀(x，y)＝D₀(x，y)+F₁(x，y)

   D₀(x，y)＝D₁(x，y)+F₂(x，y)

   D₁(x，y)＝D₂(x，y)+F₃(x，y)

         

   D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

输入的原始图像F₀可以表示为

$F_0(x,y)=M_{t_0}{F}_0(x,y)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{M}_{t_i}{D}_{I-1}(x,y)+D_{N-1}(x,y);$

c.均衡显示的重建

每一次分解过程都看作将待分解图像D_i-1分解为D_i-1的高频部分D_i与低频部分F_i+1的过程。保留每次分解获得的F_i+1，以及D_N-1，即可恢复原始图像F₀，恢复得到的图像表示为F₀′；因此，加上均衡效果的重建过程表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1}$

其中，α_i(i＝1，2，…，N+1)为均衡系数，均衡系数的选择原则为0＜α_i≤α_i+1≤1；

d.最后将F₀′通过反归一化及反对数变换，及灰度范围重映射得到均衡处理后的图像；

本发明方法的优点是不仅能将传统的滤波方法应用到其中，并且能够将非线性的尺度空间分解方法纳入到其中，这种非线性的分解方法由PDE模型实现。相对于传统的空域滤波算子卷积的方法，基于PDE模型的滤波方法，在提取尺度信号的同时，能够更有效地保持边缘。将线性与非线性尺度空间分解纳入到同一图像分解与重构框架下，这种兼容性是目前别的均衡方法所没有体现出来的，并且从原理与应用两方面对这种均衡方法的有效性得到了验证，因此，具有广泛的适用范围。

通过实验证明，本发明具有良好的均衡效果，能够有效的提高X射线图像各组织间对比度差，和细节对比度，在医用和工业使用方面都能起到良好的辅助佐用。

                         附图说明

图1为图像均衡显示处理方法的分解与重构示意图；

图2为基于线性尺度空间图像均衡显示处理方法的分解与重构示意图；

图3为图像均衡显示处理采用的医用X射线胸片原始图像，其空间分辨率为3001*3004，灰度动态范围为16比特；

图4为基于高斯平滑算子的图像均衡显示处理后图像；

图5为基于高斯微分平滑算子的图像均衡显示处理后图像；

图6为基于中值曲率驱动方程的图像均衡显示处理后图像；

图7为原图与三种实验结果相同位置和大小细节区域截图，其中(a)为胸片原始图像的局部截图(b)为基于高斯平滑算子的图像均衡显示处理后图像的局部截图，(c)为基于高斯微分平滑算子的图像均衡显示处理后图像的局部截图，(d)为基于中值曲率驱动方程的图像均衡显示处理后图像的局部截图。

                       具体实施方式

本发明公开的基于尺度空间分解与重构的X射线图像均衡显示处理方法，该方法的步骤包括：首先将输入图像进行归一化和对数变换，然后选用线性或者非线性图像滤波方法进行图像的尺度空间分解，然后对分解产生的不同尺度下的图像信息进行加权，然后将加权后的信号重构，将重构的信号进行反归一化及反对数处理，得到均衡后的图像，其处理流程如图1所示。在分解过程中，可以采用空域平滑滤波算子，如高斯滤波算子，高斯微分滤波算子进行线性尺度空间分解，其处理流程如图2所示。或者采用中值曲率驱动方程进行非线性滤波及分解，处理流程参见图1。然后压缩大尺度信号对比度，达到抑制动态范围的目的，增强各级小尺度信号，实现整个图像的均衡处理结果。

图像均衡显示处理采用的医用X射线胸片原始图像如图3所示，其空间分辨率为3001*3004，灰度动态范围为16比特。

在本发明对于尺度空间分解提出了两种具体分解方法，一种是线性尺度空间分解，通过平滑滤波算子卷积实现，另一种是非线性尺度空间分解，由偏微分方程实现。由这两种分解方法分别构成两种不同的图像均衡的具体实现方法。

由线性尺度空间分解与重构构成的图像均衡显示处理方法具体步骤为：

1)将输入的X射线图像通过归一化处理及对数变换得到F₀。

2)用尺度为σ₀的尺度函数G₀与F₀进行卷积，得到F₁，F₀与F₁之差记做D₀，如下所示。

    F₁(x，y)＝F₀(x，y)*G₀(x，y)

    D₀(x，y)＝F₀(x，y)-F₁(x，y)

3)重复上述步骤，用尺度为σ_i的尺度函数G_t对D_i-1。进行平滑，得到F_i+1，D_i-1与F_i+1之差即D_i。如下所示

    F_i+1(x，y)＝D_i-1(x，y)*G_i(x，y)

    D_i(x，y)＝D_i-1(x，y)-F_i+1(x，y)

因此有，

    F₀(x，y)＝D₀(x，y)+F₁(x，y)

    D₀(x，y)＝D₁(x，y)+F₂(x，y)

    D_l(x，y)＝D₂(x，y)+F₃(x，y)

         

    D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

所以，输入的F₀可以表示为

$F_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)+D_{N-1}(x,y);$

4)均衡显示的重建方法。每一次分解过程都可以看作将待分解图像D_i-1分解为D_i-1的高频部分D_i与低频部分F_i+1，的过程。保留每次分解获得的F_i+1，以及D_N-1，即可恢复原始图像F₀，恢复得到的图像表示为F₀′。因此，加上均衡效果的重建过程可表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1}$

其中，α_i(i＝1，2，…，N+1)为均衡系数。均衡系数的选择原则为

 O＜α_i≤α_i+1≤1；

5)最后将F₀通过反归一化及反对数变换，及灰度范围重映射得到均衡处理后的图像如图4和图5。

基于线性尺度空间分解与重构的X射线图像均衡显示处理方法，可以分别采用两种平滑卷积算子。第一种是应用广泛的高斯平滑算子(即高斯滤波器)，形式如下：

$G_i(x,y,{\mathrm{\σ}}_i^2)={\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_i}e}^{-\frac{x^2+y^2}{{2\mathrm{\σ}}_i^2}}$

当采用该平滑算子时，选择分解层数为六层，即选用了五个高斯滤算子。滤波算子的尺度参数σ与均衡系数α与尺度如表1所示，对于医用X射线胸片原始图像图1均衡显示处理结果见附图4。

                      表1实施方法一参数列表

α0	α1	α2	α3	α4	α5
						0.2	0.6	0.7	0.8	0.9	1.0
σ0	σ1	σ2	σ3	σ4
						16	8	4	2	1

本发明采用的第二种算子是高斯微分平滑算子(即高斯微分滤波器)，形式如下：

$Q_i(x,y,{\mathrm{\σ}}_i^2)=(1-{\mathrm{\α}\▿}^2)G(x,y,{\mathrm{\σ}}_i^2)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_i}{[1+\frac{\mathrm{\α}}{{\mathrm{\σ}}_i^2}(2-\frac{x^2+y^2}{{\mathrm{\σ}}_i^2})]e}^{-\frac{x^2+y^2}{{2\mathrm{\σ}}_i^2}}$

与高斯平滑算子相比，高斯微分算子在平滑滤波的同时，有更好的边缘保持性能，当采用该平滑算子时，选择参数同表1，对于医用X射线胸片原始图像图1均衡显示处理结果见附图5。

由于基于高斯的空域平滑滤波是各向同性的，无法识别图像的重要信息，因此，在降低噪声的同时，使得图像的重要信息，如边缘变得模糊了。当尺度参数从小变大时，将产生边缘移位，因此，在较粗的分辨率下检测到的边缘，要得到准确的位置，必须追踪到原始的图像。从而在重建过程中，可能会引入边缘伪影，因此改进滤波技术，在滤波的同时，保护重要的边缘信息始终是关注的问题。

第二种具体实现方法具有与第一种具体实现方法类似的处理流程，所不同的是，第一种方法采用传统的平滑滤波算子卷积来获得输入图像不同尺度下的映象，而第二种方法采用偏微分方程(PDE)代替平滑滤波算子卷积进行尺度空间分解。本发明中所采用的PDE模型为中值曲率驱动方程(MCM)：

$\frac{\∂I}{\∂t}=\left|\right|\▿I\left|\right|\mathrm{div}\left(\frac{\▿I}{\left|\right|\▿I\left|\right|}\right)$

设输入的原始图像I(x，y)是上式在t＝0时刻的解I(0，x，y)＝I₀(x，y)为方程初值，定义算子M_t为MCM滤波算子，则I(x，y)在任意尺度t下的平滑结果可表示为M_tI(x，y)，以它来代替实施方法一中的卷积滤波。则输入的原始图像I₀可以表示为：

$I_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{M}_{t_i}I(x,y)+D_{N-1}(x,y)$

均衡后的输出I₀′表示为

$I_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{M}_{t_i}I(x,y)*{\mathrm{\α}}_i+D_{N-1}(x,y)*{\mathrm{\α}}_{N+1}$

其中，α_i(i＝1，2，…，N+1)为均衡系数，其选择原则为0＜α_i≤α_i+1≤1；

t_i(i＝1，2，…，N)为尺度因子，其选择原则为t₀＞t₁＞…＞t_N-1＞0。

基于非线性尺度空间分解与重构的X射线图像均衡方法下，本发明所采用的PDE模型为中值曲率驱动方程(MCM)：

$\frac{\∂I}{\∂t}=\left|\right|\▿I\left|\right|\mathrm{div}\frac{\▿I}{\left|\right|\▿I\left|\right|})$

当采用该具体实施方法时，选择参数如表2所示，对于医用X射线胸片原始图像图1均衡显示处理结果见附图6。

              表2  实施方法二参数列表

α0	α1	α2	α3	α4	α5
						0.2	0.6	0.7	0.8	0.9	1.0
t0	t1	t2	t3	t4
						256	64	16	4	1

上述三种算法的局部细节对比度分析：

Beghdadi在1989年提出引入边缘检测算子来构造对比度的定义，他提出的对比度定义为：

$C_{\mathrm{kl}}=\frac{|f_{\mathrm{kl}}-\stackrel{\‾}{E_{\mathrm{kl}}}|}{|f_{\mathrm{kl}}+\stackrel{\‾}{E_{\mathrm{kl}}|}}$

$\stackrel{\‾}{E_{\mathrm{kl}}}=\frac{\underset{(i,j)\∈W_{\mathrm{kl}}}{\mathrm{\Σ}}{\mathrm{\Δ}}_{\mathrm{ij}}\·f_{\mathrm{ij}}}{\underset{(i,j)\∈W_{\mathrm{kl}}}{\mathrm{\Σ}}{\mathrm{\Δ}}_{\mathrm{ij}}}$

    Δ_ij＝|f_ij- f|

式中，f_kl表示滑动窗口W_kl中心点的灰度值，f_ij表示边缘检测算子对应的窗口的中心点的灰度值(两个窗口不同)。Δ_ij称为像素(i，j)在滑动窗口中的边缘值。可以通过Laplacian，sobel等算子计算得到。本发明采用Laplacian算子得到Δ_ij，

表示窗口W_kl中的平均边缘灰度。

分别计算图中每点的对比度参数C_kl，设图像尺寸为M×N，则将图像的对比度参数C_M定义为：

$C_M=\underset{k=0}{\overset{M}{\mathrm{\Σ}}}\underset{l=0}{\overset{N}{\mathrm{\Σ}}}{C}_{\mathrm{kl}}/M\×N$

实验中，选择每幅图中相同位置尺寸为201×201大小的肺部纹理细节截图(图7)。滑动窗口大小选择为15×15，则实验结果如表3所示。

      表3  实验结果截图细节对比度分析

	图7(a)	图7(b)	图7(c)	图7(d)
					CM	0.0037	0.0105	0.0106	0.0105

从表3中可以看出，原始图像的细节对比度最低，某些细节未能达到人眼可识别的程度，而通过三种算法的改善，局部的细节对比度均有显著提高，实验中，细节纹理部分的平均对比度扩展为原始图像的2.8倍左右，显示出该算法有较好的局部对比度扩展能力，可以在有限的灰度显示范围内提高图像显示质量。

综上所述，本发明的X射线图像的空间尺度分解与重构均衡显示处理方法针对X射线图像自身特点，通过实际应用表明，能够兼顾图像整体对比度调整及局部对比度增强，从而有利于医用诊断及工业使用，具有实际应用价值。

Claims (6)

1.一种基于尺度空间分解与重构的X射线图像均衡显示处理方法，其特征在于，该方法具体的步骤是：

a.归一化处理

将输入的X射线图像通过归一化处理及对数变换得到F₀(x，y)；

b.尺度空间分解

重复上述步骤，用尺度为t_N-1的滤波算子M_tN-1对D_N-2(x，y)进行平滑，得到F_N(x，y)，D_N-2(x，y)与F_N(x，y)之差即D_N-1(x，y)，如下所示：

$F_N(x,y)=M_{t_{N-1}}{D}_{N-2}(x,y)$

D_N-1(x，y)＝D_N-2(x，y)-F_N(x，y)

因此有：

F₀(x，y)＝D₀(x，y)+F₁(x，y)

D₀(x，y)＝D₁(x，y)+F₂(x，y)

D₁(x，y)＝D₂(x，y)+F₃(x，y)



D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

所以，输入的F₀(x，y)表示为：

$F_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)+D_{N-1}(x,y)$

c.重建

每一次分解过程都看作将待分解图像D_i-1(x，y)分解为D_i-1(x，y)的大尺度部分D_i(x，y)与小尺度部分F_i+1(x，y)的过程，保留每次分解获得的F_i+1(x，y)，以及D_N-1(x，y)，即可恢复原始图像F₀(x，y)，恢复得到的图像表示为F₀′(x，y)，因此，加上均衡效果的重建过程表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1}$

其中，α_i(i＝1，2，…，N+1)为均衡系数；

d.反归一化及反对数变换

最后将F₀′(x，y)通过反归一化及反对数变换，及灰度范围重映射得到均衡处理后的图像。

2.如权利要求1所述的方法，其特征在于，在线性尺度空间下，M_t为空域的平滑卷积算子，则M_tF(x，y)表示为图像与空域平滑滤波算子的卷积，为F(x，y)*G(x，y)，其中*代表卷积操作，滤波器G(x，y)的尺度由参数σ控制；采用的滤波器有以下两种：

1)高斯滤波器，其表达式为：

$G_i(x,y,{\mathrm{\σ}}_i)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_i}{e}^{-\frac{x^2+y^2}{2{\mathrm{\σ}}_i^2}}$

2)高斯微分滤波器，其表达式为：

$Q_i(x,y,{\mathrm{\σ}}_i)=(1-{\mathrm{\α}\▿}^2)G(x,y,{\mathrm{\α}}_i^2)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_i}[1+\frac{\mathrm{\α}}{{\mathrm{\σ}}_i^2}(2-\frac{x^2+y^2}{{\mathrm{\σ}}_i^2})e^{-\frac{x^2+y^2}{2{\mathrm{\σ}}_i^2}}$

具体包括下列步骤：

a.将输入的X射线图像通过归一化处理及对数变换得到F₀(x，y)；

b.用尺度为σ₀的尺度函数G₀与F₀进行卷积，得到F₁，F₀与F₁之差记做D₀，如下所示：

F₁(x，y)＝F₀(x，y)*G₀(x，y)

D₀(x，y)＝F₀(x，y)-F₁(x，y)；

重复上述步骤，用尺度为σ_i的尺度函数G_i对D_i-1进行平滑，得到F_i+1，D_i-1与F_i+1之差即D_i，如下所示：

F_i+1(x，y)＝D_i-1(x，y)*G_i(x，y)

D_i(x，y)＝D_i-1(x，y)-F_i+1(x，y)；

因此有：

F₀(x，y)＝D₀(x，y)+F₁(x，y)

D₀(x，y)＝D₁(x，y)+F₂(x，y)

D₁(x，y)＝D₂(x，y)+F₃(x，y)



D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

所以，输入的F₀表示为：

$F_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)+D_{N-1}(x,y);$

c.均衡显示的重建

每一次分解过程都看作将待分解图像D_i-1分解为D_i-1的大尺度部分D_i与小尺度部分F_i+1的过程，保留每次分解获得的F_i+1，以及D_N-1，即可恢复原始图像F₀，恢复得到的图像表示为F₀′；因此，加上均衡效果的重建过程表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1}$

其中，α_i(i＝1，2，…，N+1)为均衡系数；

d.最后将F₀′通过反归一化及反对数变换，及灰度范围重映射得到均衡处理后的图像。

3.如权利要求1所述的方法，其特征在于，在非线性尺度空间下，采用偏微分方程进行尺度空间分解，定义算子M_t为MCM滤波算子，方程为：

$\frac{\∂F}{\∂t}=\left|\right|\▿F\left|\right|\mathrm{div}\left(\frac{\▿F}{\left|\right|\▿F\left|\right|}\right)$

设输入的原始图像F(x，y)是上式在t＝0时刻的解，F(0，x，y)＝F₀(x，y)为方程初值，则F₀(x，y)在任意尺度t下的平滑结果可表示为M_tF₀(x，y)；则实现步骤如下：

a.将输入的X射线图像通过归一化处理及对数变换得到F₀(x，y)；

b.用尺度为t₀的尺度函数G₀与F₀进行卷积，得到F₁，F₀与F₁之差记做D₀，如下所示：

$F_1(x,y)=M_{t_1}{F}_0(x,y)$

D₀(x，y)＝F₀(x，y)-F₁(x，y)

重复上述步骤，用尺度为t_i的尺度函数G_i对D_i-1进行平滑，得到F_i+1，D_i-1与F_i+1之差即D_i，如下所示：

F_i+1(x，y)＝M_iD_i-1(x，y)

D_i(x，y)＝D_i-1(x，y)-F_i+1(x，y)

因此有，

F₀(x，y)＝D₀(x，y)+F₁(x，y)

D₀(x，y)＝D₁(x，y)+F₂(x，y)

D₁(x，y)＝D₂(x，y)+F₃(x，y)



D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

所以，输入的F₀表示为：

$F_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)+D_{N-1}(x,y)$

c.均衡显示的重建

每一次分解过程都看作将待分解图像D_i-1分解为D_i-1的大尺度部分D_i与小尺度部分F_i+1的过程，保留每次分解获得的F_i+1，以及D_N-1，即可恢复原始图像F₀，恢复得到的图像表示为F₀′；因此，加上均衡效果的重建过程表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1};$

其中，α_i(i＝1，2，…，N+1)为均衡系数；

d.最后将F₀′通过反归一化及反对数变换，及灰度范围重映射得到均衡处理后的图像，则输入的原始图像F₀(x，y)表示为： $F_0(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{M}_{t_i}F(x,y)+D_{N-1}(x,y)$ 均衡后的输出I₀′表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{M}_{t_i}{F}_0(x,y)*{\mathrm{\α}}_i+D_{N-1}(x,y)*{\mathrm{\α}}_{N+1};$

其中，α_i(i＝1，2，…，N+1)为均衡系数。

4.如权利要求1所述的方法，其特征在于，所述的α_i(i＝1，2，…，N+1)为均衡系数，其选择原则为0＜α_i≤α_i+1≤1，α0用来控制动态范围，其余均衡系数用来控制各个尺度信号的局部对比度；尺度参数t_i(＝1，2，…，N)的选择原则为0＜t_i+1＜t_i，分解层数N不超过7。

5.如权利要求2所述的方法，其特征在于，在线性尺度空间下，所述尺度参数σ_i(i＝1，2，…，N)的选择原则为σ₀＞σ₁＞…＞σ_N-1＞0，具体选择σ₀＝2×σ₁＝4×σ₂＝…＝2^N-1×σ_N-1＞0。

6.如权利要求3所述的方法，其特征在于，在非线性尺度空间下，所述尺度参数t_i(i＝1，2，…，N)的选择原则为t₀＞t₁＞…＞t_N-1＞0，具体选择t₀＝4×t₁＝16×t₂＝…＝4^N-1×t_N-1＞0。

Images