改进的直方图均衡化与双边滤波的X射线图像对比度增强算法

赵爱玲，张鹏程，杨一鸣，陈 燕，刘 祎，桂志国

(1.中北大学 信息与通信工程学院，山西 太原 030051；2.生物医学成像与影像大数据山西重点实验室(中北大学)，山西 太原 030051)

0 引 言

X射线成像广泛应用于医学诊断、工业及航天等领域[1].但X射线检测厚薄不均工业构件时，易出现图像对比度低，像素分布不均匀，且在获取、传输过程中易受到噪声污染，噪声会覆盖图像的细节，导致图像清晰度低，一些细节很难被观察和分析的问题，因此，需要对X射线图像增强，以便于后续的应用[2-3].

图像增强方法有灰度变换、直方图均衡化(Histogram Equalization，HE)和Retinex算法等.灰度变换虽然可以提高对比度但容易导致图像的一些细节缺失；Retinex方法使增强后的图像曝光严重[4]；HE算法的本质是图像的灰度级经过某种变换，使结果的灰度级分布较广且呈现均匀分布，提高了图像的对比度，但是传统直方图均衡(Traditional Histogram Equalization，THE)方法针对不同的图像会造成增强效果不明显或者过度增强，因此有必要对其进行优化.

子直方图均衡算法解决了传统直方图均衡化的亮度保持特性问题，其代表算法是BBHE算法[5]和DSIHE算法[6].递归子直方图均衡化技术的代表算法是RMSHE算法[7]和RSIHE算法[8]，该技术受参数制约且图像的增强效果不会随着递归参数的增加而越来越好.直方图修正技术解决了传统直方图均衡化的过度增强问题，修正的直方图均衡化技术有直方图剪切技术和直方图频率加权技术两种[9].直方图剪切技术的代表算法为BHEP[10]算法.直方图频率加权技术的典型算法为 WTHE[11].AGCWD[12]算法是对直方图频率加权技术的改进，通过Gamma函数来调整累计分布函数.

近年来，一些研究人员常常将直方图均衡化和灰度变换、滤波等进行结合以提高图像的增强效果.万智萍[13]提出一种对直方图截取和平滑再对其均衡化的红外图像增强算法.董丽丽[14]提出将Laplace算法处理后与直方图均衡化后的结果图按比例融合的图像增强算法.何文[15]提出一种用改进的高通滤波器和直方图均衡化对X射线图像增强的算法.上述提到的几种方法对于X射线图像来说，图像的增强效果不好而且有的方法参数太多，不容易调整到视觉效果最好的状态.

为了解决上述不足，本文提出改进的直方均衡化与双边滤波的X射线图像增强算法.首先采用改进的直方图均衡化算法增强图像的对比度，然后采用双边滤波对增强后的图像去噪，最后采用改进的基于摄影的全局色调映射方法对图像的动态范围压缩，使图像清晰地显示在显示设备上.

1 本文X射线图像增强算法

基于改进的直方图均衡化与双边滤波的X射线图像增强算法过程如图 1 所示.

图 1 本文X射线图像增强算法流程图

由图 1 可知，所提增强算法主要有3个阶段：1)基于改进的直方图均衡化算法增强图像对比度；2)基于双边滤波的噪声抑制；3)改进的基于摄影的全局色调映射.

1.1 基于改进的直方图均衡化算法

THE算法易使图像的灰度级缺失从而使图像清晰度降低，因此，本文对THE算法进行改进，提出一种图像亮度保持的子直方图均衡化算法，其工作步骤如下：

1)采集X射线图像，并利用图像低照度区域的阈值LT和图像高照度区域的阈值LH，将直方图分为3个子直方图.

2)分别计算3个子直方图的概率密度函数，并对概率密度函数进行修正.

3)通过调整后的概率密度函数与概率密度函数的总和的比值得到子直方图的修正累计分布函数.

4)用伽马矫正和修正累计分布函数去修正输出函数曲线.

首先，利用key[16-17]值计算出图像低照度区域的阈值LT和图像高照度区域的阈值LH，从而将直方图分为3个子直方图L1，L2，L3.图像亮度小于LT时，为低照度子直方图；大于LH时为高照度子直方图；介于LT和LH之间的亮度区域为中照度子直方图.其具体的实现过程为

(1)

(2)

式中：Lav为对数亮度平均值；L(x)为图像在点x处的亮度值；w为图像中的像素总数；δ为一个较小的值，是为了避免像素值为0时带来的运算异常；Lmax为图像的最大亮度值；Lmin为最小亮度值.将图像归一化后根据式(3)、式(4)得到图像低、高照度区域的阈值.

LT=Lmaxn-(0.9+0.1key)(Lmaxn-Lminn),

(3)

LH=Lminn-(0.6+0.4key)(Lmaxn-Lminn),

(4)

式中：Lmaxn,Lminn分别为图像归一化后的最大值和最小值.

L=L1∪L2∪L3,

(5)

L1={L(i,j)≤LT,∀L(i,j)∈L},

(6)

L2={LT

(7)

L3={L(i,j)>LH,∀L(i,j)∈L},

(8)

式中：L1，L2，L3分别为得到的3个子直方图.每个子直方图的概率密度函数为

(9)

下面对概率密度函数调整，可以避免边缘效应的产生.

(10)

式中：pdfmaxi和pdfmini分别代表3个子直方图中概率密度函数的最大值和最小值；αi是调节参数.因此，修正的累计分布函数为

(11)

(12)

式中：∑pdfwi是子直方图的修正的概率密度函数的总和.

一般来说，X射线图像像素的分布特点为：图像的前景(细节)像素主要分布在灰度级较小的区间内，而背景像素主要分布在灰度级较高的区间.

根据X射线图像像素分布特点，为得到一个理想的图像增强效果，本文运用伽马矫正方法针对3个子直方图采取与之相适应的γ系数使子直方图的像素级得到不同程度的拉伸.在低照度子直方图中，像素分布紧密且数量较多，因此需将图像像素值的动态范围进行扩展；在中照度子直方图中，像素分布并不是非常集中和密集且数量较少，因此，应把处于中间段的像素值基本保持不变的像素映射到直方图均衡化图像中；对于高照度子直方图，由于并不包含图像信息，因此可对图像像素进行压缩.

经过上述分析，伽马校正参数的选择是非常重要的，但通过参数γ调整图像的对比度时容易丢掉一些细节信息，图像的增强效果差.因为轻微地改变参数γ就会对结果产生很大的变化，而用累计分布函数cdf作为自适应参数时，其修正过程中图像强度缓慢增长，所以采用修正累计分布函数和Gamma函数去修正输出函数曲线.其表达式为

Ti(L)=Lmaxi(L/Lmaxi)γi=Lmaxi(L/Lmaxi)1-cdfwi(L),

(13)

式中：L为灰度级，i=1，2，3；Lmaxi分别为3个子直方图灰度级的最大值；Ti(L)为经过伽马矫正子直方图均衡化算法之后的图像的像素值.

1.2 基于双边滤波的噪声抑制

X射线图像在获取、传输以及格式转换过程中，都会受到噪声的影响[18]，且经过改进的直方图均衡化增强算法后，噪声也被相应放大且噪声会覆盖图像的部分细节，因此，本文采用双边滤波对图像进行降噪以更好地抑制噪声.

双边滤波是一种非线性滤波器且可以抑制X射线图像中的噪声[19].双边滤波器既考虑了空间距离又考虑了像素值间的相似性，其表达式为

(14)

(15)

(16)

(17)

式中：Ts，Tp为经过改进的直方图均衡化后的图像T(T=T1∪T2∪T3)在点s，p处的像素值；Js为像素点s经过滤波后的像素值.Ω为像素点s的邻域范围；μ(s)为点s处的归一化系数；σr为控制空间邻近度因子的衰减程度；σd为控制亮度相似度因子的衰减程度.本算法中，σr取值为5，σd值为3，滤波器的窗口大小为N=3.

1.3 改进的基于摄影的全局色调映射

X射线图像的动态范围很大且每个像素值可能需要多个字节进行保存，而一般的图像显示设备只能显示8 bit灰度图像或者24 bit彩色图像，故若直接将X射线图像显示在显示设备上时需对图像像素值压缩或者剪切，但这样会造成图像信息以及细节的丢失[20].基于摄影的全局色调映射[21]可以看作是灰度拉伸的过程，该方法只是从整体对灰度范围压缩，而忽略了灰度级之间的关系，导致图像的局部对比度降低，因此，本文对色调映射算法进行改进，改进后的算法如式(18)所示.该算法综合考虑了图像的全局对比度和局部对比度，且在拉伸过程中减小高亮度区域的像素值的分布范围，同时扩大低亮度区域的像素值的分布范围，进而达到了增强图像细节的效果.

(18)

式中：Lla(i,j)为图像的像素值与其均值的比值；Vl(i,j)为局部对比算子，在局部邻域里对图像的对比度进行拉伸，使图像可以看到更多的细节；Vg(i,j)是全局对比算子，全局对比算子可以增强较暗区域的灰度级，增强了图像的对比度.

(19)

(20)

式中：J(x)是经过滤波后图像的像素值；Lavg是图像像素值的平均值；δ是一个极小的值，用来避免像素值为0时带来的运算异常.

(21)

(22)

式中：m为滤波的窗口大小，本文设置窗口大小为3*3；ω为滤波模板.

(23)

全局对比算子Vg(i,j)反映了图像像素值与像素均值之间的差异.β是全局对比算子的增强系数，此处β取值为0.6，图像的全局对比度随着β值的增大而减小，当β值增大时，在灰度级较低的区域可以看到更多的细节信息，反之当β值减小时，图像的很多细节信息不能很清晰地显示出来.

2 结果分析

实验环境：Inter(R)Core(TM)i7-7700CPU@3.6GHz，内存为16.0 GB.操作系统为Windows10，采用Matlab2016a仿真实现.为了说明本文算法的有效性及其通用性，本文通过两组不同的数据进行算法的验证.

第一组实验参数为：X射线能量为95 kV、150 μA，采用平板探测器接收数据，探测器型号为迈普3521，探测器大小为2 352*2 944，探元大小为0.2 mm，故图像大小为2 352* 2944，像素为2个字节.

第二组实验参数为：X射线能量为110 kV、150 μA，采用平板探测器进行接收数据，探测器型号为PaxScan2520D，探测器大小为1 024*1 024，探元大小为0.2 mm，故图像大小为1 024*1 024，像素为2个字节.

2.1 主观分析

为验证本文算法的有效性，分别利用本文提出的算法、董丽丽[14]提出的基于边缘信息融合的细节增强算法(Histogram Equalization with Edge Fusion, HEEF)以及何文[15]提出的基于高频强调滤波的医学X光图像增强算法对X射线图像增强，其结果如图 2 和图 3 所示.

图 3 第二组测试图各算法的结果图及对应的直方图

由图 2 和图 3 的结果图及对应的直方图可知：本文提出的方法有良好的视觉效果且可以看到更多的细节，图像边缘没有模糊且没有伪影产生，图像增强效果比较自然.由表 1 的局部区域图可知：对于区域1和2，从HEEF算法、何文提出的算法以及本文的算法所得结果中都可以清晰地看到图像的细节.对于区域3，HEEF算法在图像的边缘有光晕及伪影产生，何文所提的算法导致图像的灰度级缺失,而本文算法将图像细节清晰地显示出来，并且没有伪影生成.由表 2 的局部区域图可知：本文所提的算法可以将细节清晰地显示出来，如区域1中显示的字母“X O”和区域2中的数字“16”对比度最高，何文算法次之.在区域3中“DC”两个字母在这三种方法中都较为模糊，但本文方法对比度较高.

表 1 第一组测试图各算法结果图的局部区域

表 2 第二组测试图各算法结果图的局部区域

表 3 各算法的评价指标

故综合来看，本文所提的算法在视觉效果上相对于其他两种方法更好.

2.2 客观分析

为了验证本文算法的客观性，这里采用信息熵、对比度和方差三种无参照图像的评价指标对图像质量进行评价.

由表 3 可知：本文所提算法的对比度和方差相对来说比较大，表明本文提出的算法可以有效提高图像的对比度.本文所提算法的信息熵最大，说明该算法可以看到更多的细节信息.因此，与其他几种算法相比，本文提出的算法可以提高图像的对比度和清晰度.

综上所述，本文算法对于X射线图像增强，在主观视觉上，与其他算法相比，可以看到更多的细节信息且提高了图像的对比度，在客观评价中，其各项评价指标总体来说优于其他算法.

3 结 论

高动态范围X射线图像增强是比较新的研究课题，本文结合直方图均衡化、双边滤波以及色调映射，综合考虑了图像的全局和局部信息，组合成一个能应用于低对比度、低照度的工业X射线增强算法.用方差、信息熵和对比度三个评价指标对本文提出的算法的有效性进行了验证.实验表明，与其他方法相比，本文提出的算法可以更有效地提高图像的对比度和清晰度，使处理后的图像更加符合人眼的视觉效果.