基于简单射线散射模型的X射线图像增强算法

张文璐，刘 祎，张鹏程，桂志国

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

0 引 言

射线检测技术通常应用于无损检测(Nondestructive Tesing,NDT)领域，对工业零件、材料、设备以及各种样本进行安全性检测[1-3].然而在传统的射线成像技术中，由于物体本身的结构与散射的X射线叠加，所得图像的可见度通常会受到限制.为了克服这一问题所广泛使用的方法有：抗散射网格技术，空气间隙技术和使用测量、数学物理建模或两者结合的散射校正方法[4-9].然而，这些方法需要特殊的设备和额外的工作来研究X射线的散射特性.另外，在一种情况下有用的方法可能在另一种情况下不起作用，比如，抗散射网格被认为是改善图像对比度必不可少的工具，但由于散射X射线的高透射率，它们在高能X射线检测系统中是没有用的.

X射线图像常用的增强方法有灰度变换法、直方图均衡化法(Histogram Equalization，HE)和Retinex算法等.灰度变换法虽然可以提高对比度但容易导致图像一些细节的缺失；Jobson等[10]提出了单尺度、多尺度等多种形式的Retinex算法，该算法得到的像素值是一定邻域内像素的加权平均值，但增强效果随尺度变化较大，会使增强后的图像整体对比度偏低，曝光严重；传统直方图均衡化(Traditional Histogram Equalization，THE)算法的本质是图像的灰度级经过某种变换，使灰度级呈现均匀分布且分布较广，提高了图像的对比度，但是针对不同的图像会出现增强效果不明显或者过度增强现象；Sherrier等[11]提出了局部自适应的直方图均衡化增强方法，避免了当图像直方图出现多波分布时的过度增强现象，但直方图均衡化过程中会合并一些灰度级较少的灰度值，使得图像的一些细节信息被丢掉.

K.S.Kim等[12]将简单射线散射模型应用于X射线图像的增强，此散射模型与大气光散射模型基本一致，大气光散射模型被广泛应用于单幅图像去雾.He等[13]基于大气光散射模型提出的暗通道先验(DCP)被认为是最经典的单幅图像去雾的方法之一，通过观察大量户外无雾图像得出结论：在绝大部分非天空的区域里，大部分像素至少会有一个颜色通道拥有非常低的亮度值且趋近于零.以此为基础，Peng等[14]，Nair等[15]和Huang等[16]改进了DCP，以获得更优异的去雾效果.此外，Li和Zheng[17]通过使用全局导向滤波(G-GIF)来复原有雾图像.Berman等[18-19]提出了非局部先验(NLP)，描述了从场景亮度到观测图像亮度值的变化，以估计透射率图和大气光值.Meng等[20]引入了新的边界约束和上下文正则化来计算图像的透射率.

本文针对以上问题提出一种基于简单X射线散射模型的图像增强算法来提高单一射线图像的可见度.首先，根据简单散射模型复原X射线图像，在复原过程中加入了新的中值算子，可以有效抑制射线图像中的脉冲噪声成分；同时，对细化之前的透射率图像应用了导向滤波，使灰度值跳跃的边缘部分的复原效果得到明显改善，并使光晕伪影现象最小化；然后，使用一种灰度优化的方法对复原图像进行后处理，在提高对比度的同时保留图像中的小细节.最后，通过实验来验证该算法对于低对比度、低亮度射线图像的可行性.

1 简单射线散射模型

1.1 简单射线散射模型

图1 为射线图像增强所用的简单射线散射模型的原理图.图中，I(x,y)是观察物的射线图像，μstru和μsurr分别为物体结构和周围介质的线性衰减系数，l为射线路径.

图1 简单射线散射模型原理图Fig.1 Schematic diagram of simple radiographic scattering model

在传统射线成像技术中，与Koschmieder模型[21]类似，射线图像I(x,y)可以用式(1)简单地表示为

I(x,y)=t·J(x,y)+(1-t)·A，

(1)

(2)

(3)

式中：t(x,y)表示物体周围介质的直接透射率；J(x,y)表示物体原本的结构也就是未被散射影响的图像；A为射线散射强度.

如式(1)所示，如果要复原可见度差的射线图像，就需要根据简单射线散射模型估计出射线散射强度A和物体周围介质的直接透射率t(x,y)，复原后的射线图像J(x,y)可表示为

(4)

1.2 射线散射强度A

利用图像去雾算法中He等[13]提出的暗通道先验(DCP)来消减逆问题的难度.如果为灰度图像，那么DCP表示图像每个局部块中的最小值，换句话说就是对图像进行最小值滤波.

射线散射强度A可以由DCP来估算.首先，从输入图像I(x,y)的DCP图像中选取亮度值前0.1%的像素，根据这些像素的位置在输入图像的I(x,y)中找到相对应的像素，再选取其中亮度值最高的像素点作为散射强度A.

1.3 透射率t(x,y)的估算

式(4)中的未知数比已知数更重要，所以，单一图像的复原存在很强的欠约束问题.Meng等[20]提出基于边界约束和上下文正则化的去雾算法，对式(4)中的透射率t(x,y)应用了更多约束，通过获得物体周围介质透射率的固有边界，并将其与基于上下文正则化和L1范数的加权函数相结合，就转化为了一个优化问题来求解未知的直接透射率.

在几何上，一个被散射影响的像素会被“推”向散射强度A(见图2)，得到像素I(x,y).所以，可以通过逆推这个过程来获得没有被散射影响的像素J(x).考虑到图像的亮度有上下边界，故

C0≤J(x)≤C1,∀x∈Ω,

(5)

式中：C0和C1是两个与输入图像亮度相关的常量(本文中C0=100，C1=300).对任意一个像素x都可以计算出相应的边界约束点Jb(x1)和Jb(x2)，构成如图2 所示的亮度方块.

图2 亮度方块Fig.2 Brightness cube

上述给出的J(x)的必要条件相当于又对t(x)施加了边界约束，t(x)的边界约束为

0≤tb(x)≤t(x)≤1,

(6)

式中：tb(x)是t(x)的下边界.

(7)

(8)

式中：ωx是以x为中心的局部块.

为解决复原图像出现的光晕伪影现象，Meng[20]引入了一个基于上下文正则化和L1范数的加权函数W(x,y).

W(x,y)(t(y)-t(x))≈0,

(9)

式中：x和y是两个相邻的像素.

通过计算像素的像素值差来构造加权函数，基于两个像素颜色向量平方差所构造的加权函数为

Wj(i)=e-∑c∈{r,g,b}|(Dj⊗Ic)i|2/2σ2,

(10)

式中：σ是一个规定的参数；Dj是一组高阶微分滤波器，为了保护图像的边缘和角落，它由8个Kirsch算子和1个拉普拉斯算子组成.

因为L1范数比L2范数对异常值普遍具有更好的鲁棒性，所以使用L1范数的积分代替L2范数的积分来进行正则化.同时，再引入一组高阶微分滤波器，对整个图像区域整合上下文约束会得到对t(x)的上下文正则化为

∑j∈ω‖Wj∘(Dj⊗t)‖1,

(11)

式中：ω是一个索引集；∘表示逐元素乘法运算符；⊗表示卷积算子；Wj是权重矩阵.

(12)

(13)

(14)

2 基于简单射线散射模型的改进算法

由式(13)求出细化后的透射率图像之后，就可以根据式(4)得到Meng[20]的算法所复原出的图像.但此时复原出的图像对比度偏低，噪声明显，在像素值跳跃的边缘部分复原效果不自然.

2.1 改进的边界约束图像

(15)

2.2 导向滤波

(16)

一般滤波过程可描述为

Wlp(Ig)=

(17)

导向滤波通过代价函数E(ak,bk)来寻找使输出和输入差值最小的系数(ak,bk).

(18)

式中：(ak,bk)为线性系数，在一个窗口ωk内假设为常数，ε为正则化参数，以防止ak过大.窗口大小通常由其半径r定义，r是中心像素到外部像素的像素距离.由于使用正方形窗口，因此总窗口大小为(2r+1)×(2r+1).式(18)的解如式(19)和式(20)所示.

(19)

(20)

(21)

图3 Meng的算法改进前后的比较Fig.3 Comparison of Meng’s algorithm before and after improvement

2.3 灰度优化

为解决复原图像整体对比度偏低的问题，使用一种灰度优化的方法对图像进行后处理，这种算法可以在提高图像对比度的同时保留图像细节.完整算法为

(22)

式中：J(x,y)为输入图像；JG(x,y)为灰度优化后的输出图像；max为输入图像的最大像素值；min为输入图像的最小像素值.

图4 本文算法简化流程图Fig.4 A simplified flowchart of the proposed algorithm

3 实验结果及分析

3.1 本文算法图像增强结果

图5 给出了X射线图像各阶段增强过程的图像，包括初始图像、边界约束图、复原图像和灰度优化后的最终增强结果图像.如图5(a)和(d)所示，增强后的图像整体对比度有显著提高，零件的大部分结构更加清晰可见，圈内所示的零件结构在增强后的图像中更容易识别，视觉上也更加立体.

图5 本文所提算法的图像增强过程Fig.5 The image enhancement process of the proposed algorithm

3.2 各算法对比的定性分析

本文分别采用CLAHE算法[24]、Meng的算法和本文算法对四种构件的X射线图像进行图像增强，结果如图6 所示.由图6 可知：1)对于各原始图像，由于本次实验的构件薄厚不均匀，X射线能量过小则无法穿透构件厚的部分，看不到构件的结构；X射线能量过大则会导致构件薄的部分过度曝光，丢失构件细节；适中的X射线能量可以兼顾构件的薄、厚部分，但会导致所成图像可见度低、对比度低.2)经过CLAHE算法增强后图像整体亮度、对比度有所提升，但光晕伪影现象严重，模糊了图像的边缘和细节部分.3)经过Meng的算法处理后，图像整体亮度明显提高，构件整体结构清晰可见，局部细节也被凸显出来，但图像1灰度值跳跃的边缘部分效果不自然，图像2背景被改变，且图像对比度偏低，亮度跳跃的细节部分会出现轻微过曝光现象，图像3有严重的噪声，图像4整体亮度、对比度都偏低.4)经过本文所提算法增强，图像整体的亮度、对比度都有明显提升，构件的整体结构清晰可见且视觉上更加立体，构件的局部细节也得到了充分的展示且层次感清晰，图像的过曝光现象比Meng的算法增强的图像明显减少.

图6 本文算法与其他图像增强算法结果的对比Fig.6 The results of this algorithm are compared with other image enhancement algorithms

3.3 各算法对比的定量分析

1)空间频率

空间频率[25]是指观察者的每一度视角内，图像的亮暗做正弦调制的周数，单位为周/(°).空间频率不同，图像中的噪声对视觉的影响也不同；空间频率越大，图像噪声对人眼视觉的影响就越小，图像细节凸显得越好，这与人眼主观感受是一致的，换言之，对于人眼视觉而言，空间频率越大，图像质量越好.图像的空间频率可由式(25)求得.

(23)

(24)

(25)

式中：m,n分别为图像的宽和高；fr和fc分别为图像的行空间频率和列空间频率；I表示要计算的图像；(i,j)为图像像素位置索引；fs为图像整体的空间频率.表1 为图6 所示的X射线图像增强后的空间频率.

表1 不同增强算法处理后的X射线图像的空间频率Tab.1 Spatial frequency of X-ray images processed by different enhancement algorithms

2)平均梯度

平均梯度[26]可以敏感地反映出图像对微小细节反差的凸显能力，在图像中，某一个方向的灰度级变化率越大，则它的梯度值也就越大.所以，平均梯度反映了图像的清晰度和细节变化，平均梯度越大说明图像细节越清晰，图像质量越好.平均梯度的表达式为

(26)

表2 不同增强算法处理后的X射线图像的平均梯度Tab.2 Average gradient of X-ray image processed by different enhancement algorithms

将表1 和表2 的数据与图7 进行对比分析.由表1的数据可以看出，3种增强方法的空间频率指数最大的均是本文算法，与人眼视觉效果基本一致；由表2的数据可以看出，1，2和4号图经过3种增强算法处理后，平均梯度值最大的均为本文算法，虽然3号图的处理结果平均梯度值最大，但与本文算法的平均梯度值相差很小，与人眼视觉效果基本一致，证明了本文算法的优越性.

4 结 论

本文针对单一X射线图像的增强提出了一种基于边界约束和上下文正则化的图像去雾改进算法，并将其应用于X射线图像.该算法加入了新的中值算子，同时对细化之前的透射率图像应用了导向滤波，最后用一种灰度优化的方法对复原图像进行后处理.用空间频率和平均梯度两个评价指标对本文提出算法的有效性进行了验证.结果表明，与其他方法相比，本文提出的算法可以更有效地提高图像的对比度和清晰度，图像的细节和结构更细致，并且使光晕伪影现象最小化，噪声明显减少，视觉效果更好.