Study on fusion method of X-ray images of ancient bronze mirrors based on lifting wavelet

文物在出土后，表面会有较厚的锈层或者其它附着物，局部甚至全部掩盖文物信息，以至于通过直接观察很难完整地显现信息，而在未知的情况下对文物进行保护或修复处理，容易对其造成伤害。X射线成像技术作为5种常用的无损检测技术之一，可以在文物检测与文物信息重建中广泛应用，在我国文物研究中，随着文物保护工作的日益完善与更加科学化的保护研究，X射线检测技术越来越受到重视和普及，使得在文物保护方面也得到了技术支持与推进^[1]。X光片上表现的信息较少地受到表面腐蚀物的干扰，它所呈现的文物面貌比肉眼直接观察更为精准。X光射入文物后，会与文物发生相互作用，被文物材质吸收或被其散射。不同材质和厚度的文物对X光的反应不同，从而使X光射线的强度在射穿文物后发生不同程度的衰减^[2]。材料中所含元素的相对原子数越大，材料厚度越厚，对X射线吸收越强。因铜镜的边缘区和纹饰区的厚度不同，最佳X光投射强度也不同。因此，分别使用边缘区和纹饰区的最佳投射强度，会得出两幅X光片，两幅X光分别集中了边缘区和纹饰区的详细信息，图像融合技术可良好地融合铜镜的详细信息。

图像融合是通过计算机技术和图像处理等技术，将多源通道收集的相同目标图像，提取每个通道的数据，并且集成有利信息用来形成高质量图像的处理过程^[3-4]。小波变换因其良好的时频局部性和多分辨性，已成为非平稳信号分析的常用工具，被广泛应用于包括图像融合在内的图像处理之中^[5]。近几十年里，图像融合研究取得了长足的发展，出现了许多的融合方法和融合策略，而基于多分辨率分析的图像融合方法考虑到人眼视觉系统的敏感特性，将图像分解到不同尺度空间，有针对性地突出图像的重要特征和细节信息^[6]。实际上，大多都是可见光与多传感器的多源图像融合，如CHEN等人使用小波变换的方法融合了红外与X光的图像融合^[7]。ZHANG等人使用稀疏表示对遥感图像进行融合^[8]。X光多源图像融合仅在医学影像领域有所应用与研究，CHEN等人使用光线跟踪配准图像算法融合多源X光图像的信息融合^[9]。YANG等人使用提升小波变换融合医学上的计算机断层扫描(computed tomography, CT)与核磁共振成像(magnatic resonance imaging, MRI)图像等^[10]。但X光多源图像融合在文物信息重建方面尤为缺乏。

图像融合可以根据信息流的不同分为不同的融合等级，分别为像素级图像融合、特征级图像融合和符号决策级图像融合^[11-13]。

图像在像素级的融合是最基本的图像融合，它将不同的物理参量进行融合。根据融合策略直接在源图上进行处理，组合源图像中的像素得到融合图像^[14]。但由于它对像素进行操作，计算机必须处理大量数据，并且处理需要很长时间，因此合并后的图像不能及时显示出来，无法实现实时处理。

特征级图像融合首先根据图像的纹理、对比度、形状等进行特征提取，然后融合源图像的特征域，最后将这些特征融合到同一图像中。特征级融合对图像信息进行了压缩，实时性有所提高。因图像特征被提取为融合信息，所以丢失了许多细节特征。

符号决策级图像融合是对图像信息更高的抽象融合，融合方法是从源图像中抽取特征和分类，处理得到各个决策的权值与置信度，并将这些转化成全局最佳决策，每个适应特征的选择取决于实际情况中的不同因素。符号决策级图像融合的计算量是最小的，且对前两种方法具有很强的依赖性。

近20多年，已经提出了许多用于构造小波及其滤波器组的方法，用于工程应用领域和信号处理，如正交小波、半正交小波和双正交小波等。

这些方法大多数在频域中构建，并且过程较复杂。1995年，贝尔实验室的SWELDENS博士提出了一种不依赖于Fourier变换的第2代小波构造方法, 利用时域提升模式构造小波。与塔式算法相比，第2代小波变换方法是一种更有效的小波变换实现方法。它的优势有以下4点：(1)不需要进行Fourier变换，变换过程直接在时域内完成，同时继承了第1代小波多分辨率的特性，具有很强的适应性与灵活的结构设计能力；(2)灵活性和通用性强，小波函数的改善是使用简单的小波函数改进小波基的特征来实现的；(3)它不再是一种特定小波函数通过尺度和平移得到，适合于不等间隔采样的构造问题；(4)存储量相对较小，执行效率高、且算法简单、运算速度快，可以分析任意长度的信号。

1998年，DAUBECHIES和SWELDENS证明任何具有有限冲击响应滤波器的离散小波变换都可以通过一系列简单的多步提升步骤来解决。即可以使用第1代小波变换算法实现的所有离散小波变换都可以使用第2代小波方法来实现。这个结论建立了离散小波变换与提升小波之间的联系。目前，第2代小波具有多种构造方法，包括基于提升方法的第2代小波构造、冗余第2代小波构造、自适应第2代小波构造、自适应冗余第2代小波构造等等。在众多领域得到广泛应用。

提升小波分解主要包含3个步骤：分解、对偶提升和主要提升^[16-17]。

分解：在符合完全重构要求的前提下，将原始输入信号x分解成偶数序列和奇数序列，即$x_{{\rm{e}}}(n)=x(2 n), x_{0}(n)=x(2 n+1) $。

对偶提升：基于初始信号的相关性，用偶数序列预测奇数序列，预测误差$d(n)=x_{\mathrm{o}}(n)-P\left[x_{\mathrm{e}}(n)\right] $，其中, P[]为预测算子，d(n)反映原始信号x中的高频信息。

主要提升：使用d(n)更新偶数序列x_e(n)，即c(n)=x_e(n)+U[d(n)]，其中, U[]为更新算子，c(n)为原始输入信号x的低频信息。

在X光图像融合算法中，首先要确定融合规则。古铜镜X光的源图像经过提升小波分解后可得到X光图像的低频分量(A_l, B_l)和高频分量(A_h, B_h)，需采用不同的规则分别对低频分量和高频分量进行融合(F_h, F_l)。目前，低频分量的融合方法主要有单像素融合，一般采用加权平均的方法进行融合。高频分量采用区域融合，包含有区域方差，模极大值等方法进行融合^[18-21]。单像素点融合方法忽略了单点像素与周围像素的相关性，会使融合后的图像过渡不自然，但其融合速度快。加权平均的融合规则会降低图像前景与背景之间的差异，使得融合后图像的对比度降低。传统方法中高频分量大多使用区域特征来确定融合规则，增加了融合后图像的对比度浮动。融合框架如图 1所示。

Figure 1.  Image fusion frame based on lifting wavelet

X光图像富有大量的能量信息，低频子带继承了X光图像的大部分能量。加权平均准则和取系数极大值等融合规则在融合过程中会使得融合图像的对比度减弱，使得图像的信息丢失。区域能量与区域方差相结合可以很好地避免这种情况，因为X光图像包含了大部分能量信息，低频分量包含了X光图像的能量信息，使用区域能量作为系数可以有效地保留X光图像的信息。区域方差反映了图像灰度的离散程度。极大地保留了图像目标特征。因此本文中提出使用区域方差与区域能量相结合，综合两者可更好地表达源图像的信息。

首先，计算低频子带以点(i, j)为中心的区域方差，区域方差定义如下：

式中，V_{l, M(N)}(i, j)为图像M, N低频子带, 以(i, j)为中心，区域大小为S×T的方差，本文中的S×T的取值为3×3；C_{l, M(N)}为图像第i行、第j列的低频系数；$ \bar{C}_{1, M(N)}$为该部分的系数平均值。

其次，图像以第i行、第j列为中心的低频子带区域能量定义如下：

融合图像的低频子带系数以区域能量与区域方差来确定，当图像M的区域能量和区域方差都大于图像N时，选取图像M的低频子带系数作为融合图像系数；当图像M的区域能量和区域方差都小于图像N时，选取图像N的低频子带系数作为融合系数；当图像M的区域能量和区域方差都不大于或都不小于图像N时，本文中以区域能量的系数作为加权融合，图像M、N的区域能量占比值如下：

最后，低频子带融合后的系数可表示为：

图像边缘信息与图像纹理信息是高频子带融合的关键，图像的边缘及纹理含有大量的信息，X光图像边缘细节包含了纹饰区的大部分信息。X光图像的边缘细节是X光图像的重要特征。人类视觉特征主要以边缘及轮廓认识图像, 包含了图像的细节特性。目前大多数高频融合规则都是取图像特征求取理论最优解，而古铜镜X光根据镜缘区和纹饰区的材质与厚度已取得相应区域的最优解。因此使用空间频率在源图像上选取最优解，避免了融合时取特征权值后，降低或增加图像对比度和信息丢失的可能性。空间频率反映了图像的总体活动与清晰度。空间频率包含了行频率与列频率, 空间频率可定义为：

式中，f_r表示空间频率的行频率，f_c表示空间频率的列频率，对于X×Y大小的图像，f_r与f_c的定义为：

式中，X, Y表示总行数和总列数; f(i, j)表示位置(i, j)处的灰度值。

定义空间频率比$ \theta = {f_{{\rm{s}}, M, k}}\left( {{2^l}, i, j} \right)/{f_{{\rm{s}}, N, k}}\left( {{2^l}, i, j} \right)$。${f_{{\rm{s}}, M, k}}\left( {{2^l}, i, j} \right) $与$ {f_{{\rm{s}}, N, k}}\left( {{2^l}, i, j} \right)$分别表示在第k个融合子带、在分辨率为2^l下, 点(i, j)的空间频率。按以下准则初步确定融合后的子带数据：

式中，$ {D_k}\left( {{2^l}, i, j} \right)$表示在第k个融合子带、在分辨率为2^l下, 点(i, j)的融合值，${{D_{M, k}}\left( {{2^l}, i, j} \right)} $与$ {{D_{N, k}}\left( {{2^l}, i, j} \right)}$分别表示源图像M、源图像N在第k个融合子带，点(i, j)位置的融合值。

邻域选取S×T(本文中取3×3)大小的验证窗口，中心像素点的系数使用邻域8个像素点的系数来源窗口来验证并修正。修正的方式如下：目标像素点融合后的图像来源于图像M，且邻域8个像素超过1/2来源于N时，将中心像素替换为N在该位置的系数; 中心像素来源于图像N，且周围像素超过1/2来源于M时，将该中心像素位置的系数修改为M。

本实验中采用3组6幅古铜镜X光图像作为实验数据，采用的图像数据为陕西省文物保护研究院采集的古铜镜X光图像。实验中采用MATLAB仿真，平台为Intel i5-2400 3.1GHz四核处理器、4G内存。采集到的3组6幅古铜镜图像如图 2所示。图 2中，每组古铜镜X光图像都是采用边缘区和纹饰区的最佳投射强度采集的X光图像。

Figure 2.  X-ray images of ancient bronze mirrors

为了能更加客观地评价融合图像的性能，实验对比了除本文中算法以外其它3种不同的算法，分别为加权平均+模极大值、加权平均+邻域方差、传统小波变换以及本文中的方法。实验结果如图 3、图 4和图 5所示。

Figure 3.  Fusion effect of ancient bronze mirrorⅠ

Figure 4.  Fusion effect of ancient bronze mirror Ⅱ

Figure 5.  Fusion effect of ancient bronze mirror Ⅲ

融合得到的效果图一般采用图像信息熵、平均梯度与标准差来对融合图像进行客观评价比较。信息熵U反映了图像中平均信息量的多少，图像的信息熵越大，说明图像的保留的信息就越丰富，其定义为^[7]：

式中，P_i表示某一灰度值在图像中出现的概率，L为图像总灰度级数。

平均梯度V表现了图像的边界与影线两侧的灰度情况变化差异，也是图像的清晰度，图像的平均梯度越大，表明图像越清晰。平均梯度V的定义为：

式中, g(m, n)表示图像在第m行、第n列的灰度值。

标准差W可以反映图像的灰度均值的状况，标准差越大，表现出的视觉效果越好，标准差的定义为：

式中，G(m, n)表示点(m, n)的灰度值，${\bar G} $表示点(m, n)的平均灰度。

实验结果如表 1、表 2和表 3所示。

从表中可以看出，本文中的算法在图像的信息熵、平均梯度和标准差上都略高于其它3种算法。分析表中数据可知，在3组铜镜实验中, 相对于其它3种算法，本文中的算法信息熵平均提升了5.76%，平均梯度平均提升了28.70%，标准差平均提升了7.70%，由此可知，本文中算法的融合图像在信息量清晰度视觉效果方面都优于其它算法。综合分析和对照数据可知：本文中提出的算法有效地保留了源图像的信息，对于边缘的传递效果更优秀。

文物研究和保护首要工作是完整科学地认识和分析文物信息，X光射线成像作为一种特殊的文物投影，它所呈现的文物面貌比肉眼观察更为精准，使用提升小波将源图像分别分解为低频子图和高频子图，并分别对其低频子图采用区域方差与能量加权、高频子图采用空间频率来确定融合规则。实验结果表明，该算法具有良好的融合效果，为文物保护工作者提供了更准确可靠的信息, 在文物信息重建领域有着广阔的应用前景。