1 引言

在工业产品的疵病检测、内部结构分析等质量检测工作中,数字X射线成像检测方法一直以来都占据重要的地位。但是,受工件结构和材料等方面的约束,在X射线透照方向上,大的有效厚度差异造成单一能量参数透照图像中过曝光和欠曝光现象共存,工件的结构信息无法得到全面的反映^[1-2]。为实现大厚度比目标的检测,研究人员提出分区透照、厚度补偿、变能量透照等方法^[3]。其中,变能量X射线检测方法是通过递增X管电压来获取不同厚度范围的透照子图,并通过子图的加权融合得到最终图像,具有一定的先进性^[4-6]。虽然该方法能够在单幅图像中完整保留目标结构信息,但是加权融合过程导致图像灰度范围远远超出显示设备的能力,融合图像无法完整显示。

针对上述高动态数字图像和低动态显示设备间的矛盾,需要对高动态图像进行色阶映射以满足观察者的主观体验。目前,色阶映射中最简单的方式是线性映射,该算法计算量小,在图像动态范围与显示设备动态范围相近时能得到令人满意的效果^[7-8]。研究者分析了箝位伽马映射、对数映射和指数函数映射算子,与线性映射算子相比,上述算子能够得到较好的效果,但是当图像动态范围与显示设备动态范围相差较大时效果一般。研究者提出了基于直方图调整技术的色阶映射算子,该方法与图像增强领域的直方图均衡化方法类似,认为高动态范围图像的直方图呈稀疏分布,但是融合的高动态范围X射线数字图像的直方图分布不符合上述假设 ^[9]。

因此,本文针对变能量X射线子图融合的高动态范围图像的色阶映射问题,开展了多分辨梯度信息增强算法的研究。该算法充分利用X射线重点表征结构变化信息的特点,将高动态范围图像进行多分辨分解,以获取不同变化程度的结构信息,并对其进行增强。通过融合多分辨梯度增强图像,利用恢复算法获取特定色阶范围的灰度图像,实现图像的可靠显示。

2 变能量X射线高动态范围成像的基本原理

X射线在穿透物体的过程中遵循比尔衰减定律,由于X射线辐射能谱为连续能谱,则在X射线与物质相互作用过程中,X射线的衰减特性可以表示为

I=∑jI0(Ej)exp[-μm(Ej)ρd],(1)

式中:I₀为入射射线的初始强度;I为透射射线的强度;μ_m为特定能量X射线下物质的质量衰减系数,该系数是光子能量E_j的函数;j为能量抽样编号;ρ为物质密度;d为射线在穿透物质过程中经过的路径长度。

由射线与物质相互作用的机理可知,通过增加射线源管电压数值可以有效增加X射线的能量,从而增加射线对目标的透照能力,即有效透照厚度随管电压数值的增加而增加。在工件有效透照厚度范围确定的情况下,通过管电压的递增调节,可以获取不同厚度工件的有效子图 GUi(x,y)(i=1,…,n),其中U_i为管电压值,(x,y)为图像像素点的位置,i为管电压的取值编号,n为选定的管电压的数量。通过融合不同能量下的有效子图,可获得具有完整信息的工件X射线数字图像G(x,y)。

G(x,y)=Fuse[GU1(x,y),GU2(x,y),…,GUn-1(x,y),GUn(x,y)]=K1GU1(x,y)+K2GU2(x,y)+…+KnGUn(x,y),(2)

式中:Fuse(·)为图像融合算子;K₁,…,K_n为子图加权系数^[10-12]。

由(2)式可知,在变能量X射线成像过程中,最终的融合图像是各电压参数下对应透照子图的有效区域的加权拼接结果。由于(2)式中权系数通常大于1,因此融合图像的数字动态范围得到极大的扩展,最大灰度和最小灰度的比值通常超过10000∶1。上述结果类似于具有超高动态范围的探测器探测到的结果。对于具有超大灰度分布范围的数字图像,常规显示器无法有效显示目标整体的结构信息。因此,必须对其进行色阶映射,以匹配显示设备的动态范围。

3 基于多分辨梯度增强的色阶映射算法

3.1 基于梯度增强的色阶映射算法

由第2节内容可知,变能量X射线成像获取的融合图像具有超大的数字动态范围,无法与显示器的动态范围匹配。从实际应用的角度出发,在X射线检测领域中,人们更关注目标结构的信息。目标结构信息在图像中表现为灰度的变化,因此利用上述X射线成像检测的特点,本文提出基于图像梯度信息的增强方案,即在对具有高动态范围的融合图像进行亮度范围压缩的过程中,尽可能地保留灰度变化信息。

按照上述思想,对变能量X射线成像获取的高动态范围图像进行处理。假设原始图像表示为G(x,y),ÑG(x,y)为图像对应的梯度域信号。定义梯度域加权函数为φ(x,y),梯度域增强过程为

G'(x,y)=∇G(x,y)·φ(x,y),(3)

式中:G'(x,y)为梯度域信号加权后的图像。为避免图像的失真,在对梯度域信号进行加权处理时,仅对其数值进行处理,梯度方向信息则完全保留。

加权函数φ(x,y)的设计原则是使梯度图像中的较大梯度得到一定程度的衰减,较小梯度得到一定程度的增强,从而保证增强后的梯度信息能够尽可能均匀分布并得到完整的保留。最终将增强的梯度图像恢复为一定色阶范围内的灰度图像。但是,对于数字图像来说,其梯度必须为一保守场,即图像G(x,y)的梯度必须满足

∂2[∇G(x,y)]∂x∂y=∂2[∇G(x,y)]∂y∂x。(4)

但是,经过增强的梯度图像不再满足(4)式。因此,最终恢复的灰度图像I(x,y)需要采用优化算法进行求解。估计的恢复图像必须满足

I^(x,y)=argminI‖∇I(x,y)-G'(x,y)‖2,(5)

式中: I^(x,y)为最终恢复图像I(x,y)的估计值。即恢复图像的梯度与增强处理后的梯度图像逼近。在优化过程中也可加入正则项,以降低噪声的影响。

3.2 多分辨梯度的融合优化

对于X射线透照图像来说,目标结构的变化程度表现为不同程度的灰度变化,通过单一分辨率图像无法有效提取并增强所有的结构信息。因此,本文利用图像的多分辨特性,通过高斯金字塔下采样获取不同分辨率层级的图像,进而通过不同分辨率层级梯度信息,构建对应层级的加权函数,并对不同层级梯度图进行增强。

图像金字塔是一种利用多分辨率来表征图像的结构。在该结构中,一系列图像均来源于同一张原始图,且图像分辨率由大到小呈金字塔形状排列。

金字塔的底部是高分辨率,顶部是低分辨率。当向金字塔上层移动时,图像的尺寸和分辨率均呈降低的趋势。高斯金字塔(Gaussian pyramid)是通过高斯平滑和亚采样处理后获得的一系列下采样图像,也就是说第K层高斯金字塔通过平滑、亚采样处理后就可以获得第K+1层高斯图像,高斯金字塔包含了一系列低通滤波器,其截止频率从上往下呈2倍递增。拉普拉斯金字塔(Laplacian pyramid)可以认为是残差金字塔,可记录高斯金字塔每一级下采样后再上采样的图像信息与下采样前的图像信息的差异。拉普拉斯金字塔与高斯金字塔搭配进行上采样,能够完整地恢复出每一层级下采样前的图像。金字塔的图像如图1所示。

图 1. 金字塔的图像

Fig. 1. Image of pyramid

下载图片 查看所有图片

首先,通过高斯金字塔下采样操作构造原始图像G(x,y)的金字塔,其包含全分辨率的原始图像G(x,y)和多层级较低分辨率的图像G_i(x,y)。其中,每个层级分辨率图像的梯度可以采用中心差分表示:

∇Gi=Gi(x+1,y)-Gi(x-1,y)2i+1,Gi(x,y+1)-Gi(x,y-1)2i+1。(6)

定义该分辨率层级上对应的加权函数为φ_i(x,y):

φi(x,y)=‖∇Gi(x,y)‖G-β',(7)

式中: G-为平均梯度模值;β'为调制因子,决定梯度加权的程度。

多分辨梯度融合的思想是针对各分辨率层级的图像进行加权函数的构造,即得到φ_n(x,y),φ_n-1(x,y),…,φ₀(x,y),其中φ_n(x,y)对应最小分辨率层级,φ₀(x,y)对应最大分辨率层级。针对较低分辨率层级的加权函数,使用高斯金字塔和拉普拉斯金字塔进行上采样操作,并将结果与相邻较高分辨率层级的加权函数相乘,以此类推可以获取多分辨融合的加权函数。即

φ'i(x,y)=Lup[φi(x,y)]φ'i-1(x,y)=φ'i(x,y)·φi-1(x,y),(8)

式中:L_up(·)为上采样操作符。

将多分辨融合的加权函数作用到原始分辨率梯度图上,可获得最终增强的梯度图像,进而通过(5)式对灰度图像进行恢复。

对恢复后的灰度图像采用结构相似性(SSIM)判据进行量化评价。SSIM是一种全参考的图像质量评价指标 ,它分别从亮度l(x,y)、对比度c(x,y)、结构s(x,y)三个方面评价样本X和样本Y的相似性。

l(x,y)=2μXμY+C1μX2+μY2+C1,(9)c(x,y)=2σXσY+C2σX2+σY2+C2,(10)s(x,y)=σXY+C3σXσY+C3,(11)

式中:μ_X、μ_Y分别为图像X和Y的均值;σ_X、σ_Y分别为图像X和Y的方差;σ_XY为图像X和Y的协方差;C₁、C₂、C₃为常数,为了避免分母为0的情况,通常取C₁= (K1L)2,C₂= (K2L)2,C₃=C₂/2,一般K₁=0.01,K₂=0.03,L=255。则SSIM指数函数为

M(x,y)=[l(x,y)]α[c(x,y)]β[s(x,y)]γ,(12)

式中:α,β,γ为调整模型中三个成分所占的权重。

SSIM的取值范围为[0,1],值越大,表示图像失真越小。为了得到简化形式,将α,β,γ设为1,可以得到

M(x,y)=(2μXμY+C1)(2σXσY+C2)(μX2+μY2+C1)(σX2+σY2+C2)。(13)

4 实验及数据分析

为验证本文算法的可行性,设计实验对具有大厚度比的目标进行变能量X射线扫描成像。

实验工件的实物如图2所示,其由尼龙、钢、塑料三种材料组成,三种材料的密度和原子序数有很大差异,故透照方向上的有效厚度存在差异。变能量X射线成像系统中射线源为上海GE公司型号为ISOVOLT-450的X射线机,实验过程中固定管电流为2 mA,管电压为40~450 kV可调。探测器为上海Varian公司的paxscan-2520平板成像器件,物理动态范围为60 dB。

图 2. 实验工件实物图

Fig. 2. Physical drawing of experimental workpiece

下载图片 查看所有图片

在变能量X射线成像实验过程中,通过调节管电压来实现不同厚度段的有效透照,提取子图中的有效区域,并将其融合以获取高数字动态范围融合图像。图3所示为不同管电压下的透照图像。

图 3. 不同管电压下的透照图像。 (a) 50 kV;(b) 70 kV;(c) 90 kV;(d) 100 kV;(e) 120 kV;(f) 140 kV

Fig. 3. Trans-illumination images under different tube voltages. (a) 50 kV; (b) 70 kV; (c) 90 kV;(d) 100 kV; (e) 120 kV; (f) 140 kV

下载图片 查看所有图片

可以看出,由于工件的有效透照厚度差异偏大,特定的管电压只能够保证特定厚度段内的区域得到有效透照,较薄区域出现过曝光现象,而厚度较大区域则无法反映内部结构信息。因此,通过提取并融合不同管电压透照图中的有效区域,可以获取反映工件内部结构信息的高动态范围融合图像,图像的灰度范围为0~16442。为实现融合图像的可靠显示,利用线性和非线性色阶映射算法对高动态范围图像进行增强处理,得到的结果如图4所示。

图4(a)是对图像进行全局线性灰度压缩的映射结果,映射后灰度范围为0~255,图像中暗区的内部结构信息无法得到高对比度显示。图4(b)是对图像进行对数变换映射的结果,映射后灰度范围为0~255,工件的暗区结构信息得到良好的显示,但是原图像中亮区的部分信息丢失严重。图4(c)是对图像进行伽马值为0.3的映射结果,图像中亮区部分的结构信息得到保留,但是暗区的细节信息对比度较差。图4(d)是对图像中不同材质区域进行分区,区域内灰度分布比较集中,因此分区线性拉伸处理结果可以很好地表征区域内的结构信息,但是也存在明显的边界。

图 4. 不同色阶映射法处理变能量融合图像后的结果。(a)灰度线性映射的显示结果;(b)灰度对数映射的显示结果;(c)灰度伽玛校正的显示结果;(d)分区灰度映射的显示结果

Fig. 4. Different color gradation map processing results of variable energy fusion image. (a) Display result of grayscale linear mapping; (b) display result of grayscale logarithmic mapping; (c) display result of grayscale gamma correction; (d) display result of partition grayscale mapping

下载图片 查看所有图片

按照本文方法,对融合的图像进行多分辨梯度增强处理。在实验过程中,对全分辨图像进行下采样操作,下采样截止频率从上往下呈2倍递增,根据所得到的下采样截止频率,依次获取4个分辨率层

级图像。针对具体的分辨率层级图像,开展梯度函数的构造并对其进行增强,融合多分辨率增强梯度信息,获取最终增强的水平和垂直全分辨率梯度图,如图5所示。

图 5. 高动态范围图像的增强显示处理结果。(a)原始水平方向梯度图;(b)增强的水平方向梯度图;(c)原始垂直方向梯度图;(d)增强的垂直方向梯度图

Fig. 5. Enhancement display processing results of high dynamic range image. (a) Original horizontal gradient diagram; (b) enhanced horizontal gradient diagram; (c) original vertical gradient diagram; (d) enhanced vertical gradient diagram

下载图片 查看所有图片

对比原始全分辨率梯度图和增强后的梯度图,可以看出,增强后的梯度图含有更多的结构信息,不同材质区域的结构信息均得到了增强。利用增强后的全分辨率梯度图,进行灰度图像的恢复,得到的结果如图6所示。通过有效控制灰度范围,最终恢复的图像可以与常规显示设备匹配,且不同有效厚度区域的目标结构信息均得到了增强,能够反映工件的内部结构信息。

图 6. 所提方法获取的增强图像

Fig. 6. Enhanced image obtained by proposed method

下载图片 查看所有图片

为定量分析不同处理方法的效果,本文利用局部透照方法获取目标物体不同区域的透照图像,实验中准直器窗口位置以图4(d)中的区域划分为参考进行设置。虽然局部透照方法无法获取目标的完整信息,但是能够保证特定区域内结构信息的完整,因此局部透照图像可以作为全局透照图像的对比参考图像。局部透照的分区示意图如图7所示。

图 7. 目标的局部透照分区示意图

Fig. 7. Schematic of local trans-illumination partition of target

下载图片 查看所有图片

利用结构相似性判据SSIM对不同方法获取的结果进行比较,结果如表1所示。

从定量评价结果可以看出,在线性映射、对数映射、幂指数映射处理结果中,局部区域的SSIM值偏低,无法有效表征该区域的结构信息。分段线性处理方法表现良好,但是整体显示结果存在明显的边界。本文方法处理结果略次于分段线性处理方法, 但是各区域评价值比较均匀,可以完整表征目标工件的整体结构信息。

5 结论

针对变能量X射线成像获取的高动态范围数字图像的增强显示问题开展了研究,采用图像梯度非线性增强方法,在增强梯度的基础上通过优化迭代算法对灰度级图像进行恢复,使其能够与常规显示设备的动态范围匹配,以得到增强显示的效果。同时结合图像的多分辨特性,利用重采样技术构造图像金字塔,分别实现不同分辨率层级的梯度增强,并将不同分辨率层级梯度信息进行融合,从而有效保留并增强不同变化程度的目标结构信息。