基于引导滤波与稀疏表示的医学图像融合

王兆滨，马一鲲，崔子婧

(兰州大学信息科学与工程学院 兰州 730000)

在医学影像学中，常用的成像技术主要有计算机断层扫描(computed tomography, CT)、磁共振成像(magnetic resonance imaging, MRI)、正电子发射计算机断层扫描(positron emission computed tomography, PET)、单光子发射计算机断层扫描(single photon emission computed tomography,SPECT)等。其中，通过CT 和MRI 获得的图像可以提供解剖信息，具有较高的空间分辨率；前者可区分骨骼、血管等有密度差异的组织，后者可区分软组织但却无法提供骨骼信息。通过SPECT 和PET获得的图像能够反映人体的功能结构和代谢信息，但其空间分辨率一般较低。因此，若能将人体同一部位的不同医学图像融合，就可用一幅图像表示多种结构和功能信息[1]。医生就能根据融合图像做出更精确的诊断。基于以上原因，多模态医学图像融合已经成为图像处理领域的研究热点之一。

多模态医学图像融合是指将不同医学成像技术采集到的同一部位图像进行融合的过程[2]。获得的融合图像与单一成像技术提供的图像相比，其信息更加全面。然而，尽管不同成像技术的图像具有显著差异，这些图像间仍存有大量冗余信息。因此，在有效提取互补信息的同时尽可能剔除冗余信息是获得高质量融合图像的关键。

根据融合方法划分，图像融合算法可分为空间域和变换域算法。空间域算法根据融合规则直接处理图像每个像素，是逻辑较为简单的算法。变换域算法则是通过数学变换方法将图像转换为频域表示方式并得到相应的频域系数，随后采用合适的融合规则对系数进行融合，最后用相应逆变换得到融合图像。对于变换域算法来说，变换方法和融合规则的选择至关重要。

目前，已有多种融合算法应用于多模态医学图像融合领域，如稀疏表示、多尺度变换、神经网络等。然而，这些算法在图像质量和细节保存方面仍有较多缺陷。因此，本文结合引导滤波器与自适应稀疏表示，提出了一种新的多模态医学图像融合算法，用于提升融合图像的清晰度。

1 自适应稀疏表示模型

由于不同图像块的结构差异较大，传统稀疏表示算法通常需要学习一个高度冗余的字典来满足图像重建的要求。然而，这种高冗余字典不仅会使生成的图像产生视觉伪影，而且还会花费较长的运算时间。与其相比，自适应稀疏表示模型不仅构建了一个冗余字典，而且还构建了一组紧凑的子字典并在融合过程中自适应地选择这些子字典进行运算。该模型由文献[3]首次提出并应用于图像融合和去噪。此外，自适应联合稀疏模型(adaptive joint sparse model, JSM)[4]也被提出，其同样是基于自适应稀疏表示的概念。随后，文献[5]根据图像的相似结构将图像分解为相似模型、平滑模型和细节模型。其中，对细节模型使用自适应稀疏表示进行融合，提高了融合效率。文献[6]将自适应稀疏表示与修正后的空间频率结合用于多模态医学图像融合。

随后，自适应稀疏表示模型广泛应用于多尺度变换融合算法中。文献[7]为更好地提取图像边缘和方向信息，利用多个滤波器对输入图像进行分解，采用自适应稀疏表示模型对子带进行融合。文献[8] 使用非下采样剪切波变换(non-subsampled shear wave transform, NSST)分解图像，利用自适应稀疏表示模型融合高频子带，得到了较好的融合结果。文献[9]也使用了NSST 对输入图像进行分解，但自适应稀疏表示模型被用于低频子带的融合；而高频子带的空间频率被用于神经元的反馈输入，并且激励一个脉冲耦合神经网络(pulse coupled neural network, PCNN)进行融合，也得到了较好的结果。文献[10] 则利用自适应的卷积稀疏表示模型融合遥感图像。文献[11] 使用拉普拉斯金字塔分解输入图像，使用自适应稀疏表示融合了每一层子带，通过实验证明了该方法应用于医学图像融合的有效性。

在自适应字典的构建中，使用的训练集是由从多幅高质量图像中随机采样出的大量图像块构成的。同时，为保证训练集中的图像块都具有丰富的边缘结构信息，训练集需要删除强度方差小于给定阈值的图像块[12-14]。因此假设最终的训练集为P={p1,p2,···,pM}，即有M个满足以上条件的图像块。然后，对集合P中的图像块根据其梯度主导方向进行分类，构建出一组子词典训练集。子字典训练集的构建采用方向分配方法[15]，通过计算每个图像块的梯度方向直方图来决定主导方向。如对于图像块pm(i,j)∈P，其水平梯度Gx(i,j)和垂直梯度Gy(i,j)通过Sobel 算子计算得出：

随后，根据pm中所有像素的梯度幅值和方向来绘制pm的方向直方图。假设将360°平均划分为K份，用于构建K个子字典，如图1a 所示为K=4。其次，根据pm中每个像素的梯度方向将所有像素划分为K个方向。同时，将属于同一方向的所有像素的梯度幅值相加作为该方向在直方图中的统计量，由此得到pm的方向直方图。最后，选择直方图中统计量最大的方向作为pm的主导方向。由此，所有图像块根据其主导方向划分为K个图像块子集{Pk|k=1,2,···,K}∈P。每个子集的图像块都有相同的主导方向因而都具有相似的结构。根据这些子集训练出的子字典将会更紧凑。

图1 自适应稀疏表示模型子词典图解

此外，模型最终需要构建出K+1 个字典D={D0,D1,···,DK}，其中字典D0如图1b 所示，该字典是由集合P中所有图像块训练得到的，其余子字典{Dk|k=1,2,···,K}则是通过其相应子集学习得出的，如图1c～图1f 所示。

在实际融合过程中，将待融合的输入图像切分成与pm大小一致的图像块，并自适应地选择合适子字典进行融合。子字典的选择过程如图2 所示。设km表示为切分后的某一输入图像块应被分配到第km个子集，km按以下公式得到：

图2 自适应词典的选择过程

式中，θmax=max{θ1,θ2,···,θK}， θK表示该输入图像块的方向直方图中第K类的统计量，θmax为K个统计量的峰值；k*=argmax{θk|k=1,2,···,K}是θmax的指标，若km=0，则表示输入的是不规则图像块，将采用字典D0进行融合。

本文应用K-SVD(K-means singular value decomposition)算法对字典进行训练，训练集是通过对30 组多模态医学图像进行采样得到的10 万个大小为8×8 的图像块。通过实验确定了参数K的最佳值为6，最终本文训练出7 个字典，每张字典的大小均为64×256。

2 引导滤波模型

引导滤波器可较好地保持边缘信息，同时具有较低的复杂度。由于其良好的保边平滑性，大量研究将其应用于图像处理领域。文献[16] 为在融合中挖掘出更多有效信息，提出基于目标提取和引导滤波增强的图像融合算法。文献[17] 引入感知因子对引导滤波器进行改进，对原融合权重进行修正，明显改善了融合图像的边缘特性。此外，为更好保存有效信息以及提高空间连续性，文献[18]通过结合复剪切波变换(complex shearlet transform,CST)与引导滤波器的优点提出了新的图像融合方案。文献[19]使用改进的PCNN 模型与引导滤波结合，更好地匹配了人眼视觉感知。文献[20]将引导滤波器与小波变换结合，克服了引导滤波器对噪声图像处理的缺陷，并提出了新颖的权值映射融合规则。文献[21] 用引导滤波器和显著性映射分别计算输入图像分解后的权重图，得到较好的融合结果。

此外，引导滤波器也可以看作是一种图像的多尺度变换法，并且可以结合图像的像素值和空间信息来构造权重图。文献[22] 基于这种理念结合随机游走模型用于融合多聚焦图像，其客观评价结果表明该方法可以得出很好的融合效果。因此，本文算法利用引导滤波器对输入图像进行分解，从而得到细节层和基础层。下面将介绍本文使用的引导滤波器模型。

引导滤波器可看作是一种具有局部窗口的线性滤波器，定义为：

式中，G为 引导图像，用于指导输入图像I进行滤波，得到输出结果O；Wij(G)为I在(i,j)处的像素点权重值，其根据G在(i,j)处的局部窗口wk中的所有像素计算得出。

Wij(G)的计算过程如下：

式中， μk和 σk分别为局部窗口wk中的像素平均值和方差；Gi与Gj为两个相邻的像素点的像素值； ε为惩罚值。通过上式可得出，当像素点Gi与Gj在边界时，(Gi-μk)(Gj-μk)为异号，否则为同号。异号时，其权重较低，远小于同号时的权重，因此图像平坦区域的像素将会被添加一个较大的权重值，从而削弱了滤波器平滑模糊的效果，起到保持边缘细节的作用。此外，惩罚值ε同样显著影响该滤波器的滤波效果。当ε较小时，其滤波效果与上述相同，起到保持边缘细节的作用；而当ε较大时，引导滤波器就可以看作是一个均值滤波器，其平滑效果会更加显著，但却起不到保边平滑的效果。

此外，也可以从线性滤波器的角度来理解引导滤波器：

该式表明引导图像的梯度变化与结果图的梯度变化呈线性关系。因此，一般来说，引导图像的作用就是在提示输出图像不同输入区域的位置和梯度信息，以便更好地区分和保存输入图像的平滑和边缘区域。

3 基于引导滤波与自适应稀疏表示的多模态医学图像融合

本文结合自适应稀疏表示模型和引导滤波器，提出一种新的多模态医学图像融合算法，该算法的流程图如图3 所示。

图3 本文算法流程

1) 图像分解。使用高斯滤波器将输入图像A和B分解为细节层图像和基础层图像，其过程如式(12)～式(14)所示：

式中，G表 示高斯函数，使用大小为3×3的高斯核。

通过高斯滤波器得到基础层B1和B2后，将输入图像与其基础层相减得到细节层D1和D2：

2) 结合拉普拉斯滤波器和显著性区域构建基础层的权值图。首先使用大小为3×3的拉普拉斯滤波器L 作用于输入图像A和B，获得其高通图像H1和H2：

对得到的高斯图像H1和H2的绝对值进行显著性区域计算，构建显著性图像S1和S2，并得到其相对应的初始权值图P1和P2：

式中，Y是大小为(2τg+1)(2σg+1)的高斯低通滤波器，本实验中 τg和 σg的值设置为5。

若仅使用P1和P2对图像进行融合会导致融合结果的拼接感过强，丢失大量边缘信息，从而使得融合图像会不符合人类的视觉系统。因此，将P1和P2作为输入图像A和B在引导滤波器中的引导图就可以计算出基础层的最终权值图：

4) 融合细节层。利用自适应稀疏表示模型对细节层D1和D2进行融合，从而构建出融合细节层FD。由于医学图像的大小均为256×256，因此不需要图像配准过程。通过计算细节层图像块的方向直方图，然后根据式(5)自适应地选择不同字典对细节层D1和D2的图像块进行稀疏表示：

式中，Dk1和Dk2为自适应挑选的字典； α1和 α2分别为细节层D1和D2的稀疏系数。随后，对稀疏系数进行融合。在本文算法中选择依据稀疏系数的方差进行融合。因此，设稀疏系数 α1和 α2的均值分别为μ1和 μ2，则 α1和 α2的方差为：

根据稀疏系数的方差大小和稀疏矩阵的稀疏度d1=‖α1‖0、d2=‖α2‖0来指定稀疏系数的融合规则：

最后，根据融合后的稀疏系数进行细节层重构，得到融合细节层：

5) 计算融合结果图。将融合后的基础层图像与融合细节层图像相加，得到融合图像F：

4 实验结果与分析

为验证提出方法的有效性，选择3 组多模态脑部病变医学图像进行对比实验。这些图片由网站(http://www.med.harvard.edu/aanlib/home.html)获取。实验使用的3 组CT 和MRI 图像为：肉瘤患者的第17 组脑切片、脑膜瘤患者的第17 组脑切片和急性中风言语停止患者的第18 组脑切片。所有的CT 和MRI 图像的大小均为256×256。

本文算法与拉普拉斯金字塔(laplacian-pyramid,LP)、离散小波变换(discrete wavelet transform,DWT)[27]、轮廓波变换(contourlet transform, CVT)[28,29]、非下采样轮廓波变换(nonsubsampled contourlet,NSCT)[30]、稀疏表示(sparse representation, SR)和自适应稀疏表示(adaptive sparse representation, ASR)等方法进行了比较。实验结果如图4～图6 所示。

图4 肉瘤患者的融合结果

图6 急性中风患者的融合结果

4.1 融合结果的指标评价

所有算法的融合结果通过以下8 个指标进行评价：标准方差(standard deviation, SD)、平均梯度(average gradient, AG)、熵(entropy, E)[23]、空间频率(spatial frequency, SF)[24]、互信息(mutual information,MI)[25]、融合质量QAB/F、融合损失LAB/F和融合伪影NAB/F。其中QAB/F、LAB/F和NAB/F这3 个指标是由文献[26]提出，具有互补信息，且三者之和为1。

SD、 AG、E、 SF仅由融合图像F计算得出。因此，设F的图像大小为M×N，Fi,j表示F在(i,j)处的像素值，则标准差 SD反应了F的像素离散程度；SD越大表明像素值分布越分散，越小则表明分布更集中：

式中，Pi为像素值为i的像素在F中出现的概率。

空间频率 SF反映了灰度值的变化率，可用于评估图像的细节及纹理的多少。 SF值越大说明F的纹理边缘信息更丰富。 SF由横向频率、纵向频率、主对角线频率和副对角线频率组成：

式中， RF是F中所有处于同一行的相邻两像素差值的平方平均数；同理， CF、MDF、SDF分别是F中所有处于同一列、同一主对角线、同一副对角线的相邻两像素差值的平方平均数。

此外， MI、QAB/F、LAB/F、NAB/F由输入图像A和B以及融合图像F联合计算得出。其中，互信息MI表示F从A和B中获取的信息总量。 MI越大，表明融合图像从输入图像中获取的信息量越多，融合效果越好：

式中， |W|为滑动窗口 ω滑动的次数；s(A|ω)和s(B|ω)为A和B在滑动窗口 ω内的显著性特征；SSIM(A,F|ω)和SSIM(B,F|ω)为F在 ω处与A和B的结构相似度。

融合损失LAB/F是对融合过程中丢失信息的度量，即表示输出中的有价值信息在融合图像中丢失的程度。因此，其值越低，丢失的信息就会越少：

融合伪影NAB/F用于评价引入融合图像的无用信息的数量。这些信息在任何输入图像中都没有相对应的特征。因此，融合伪影本质上就是错误的信息，其值越小融合效果越好：

不同算法结果的评价指标数值由表1～表3 给出。所有指标中的最优数值加粗表示。从表1 可以看出，本文方法在E、 MI、QAB/F和NAB/F指标上的结果是最优，同时在 AG和 SF指标上取得了次优。在表2 中，本文算法在 SD、E、 SF、QAB/F和LAB/F指标上取得最优结果；同时，在 AG和 MI指标上取得次优。最后，在表3 中，本文算法在E、S F、QAB/F和LAB/F指标上取得最优结果；同时，在 AG和 MI指标上取得次优。

表1 第一组实验融合结果的客观评价指标对比

表2 第二组实验融合结果的客观评价指标对比

表3 第三组实验融合结果的客观评价指标对比

从整体来看，本文算法在E、QAB/F上均取得了最优结果，且在 MI和LAB/F指标上也可获得最优或次优结果。这证明了本文算法相比于其他算法减少了输入图像的信息损失，保存了最多的输入信息，同时具有优良的纹理保留和细节信息的能力，因而在 AG和 SF指标上多次取得了最优或次优结果。最后，对于 SD和NAB/F指标，尽管本文算法分别存在一组最优结果，但其整体表现不佳，这意味着该算法的融合结果存在一定的伪影干扰，且融合图像的像素分布较为集中。

此外，其他6 种算法在某些指标上的表现也较为优异，甚至可以取得超越本文算法的结果，如表1 和表3 所示，LP 算法在 SD指标上取得了最优结果。表2 中的结果也表明LP 算法取得了次优结果。而SR 算法和ASR 算法在 AG、 MI、 SF这3 个指标上表现较为优异，在3 组实验中多次取得最优和次优结果。最后，DWT 算法在NAB/F指标上的表现较为出色，这表明其融合结果有较少的伪影。

4.2 融合结果的视觉评价

为了进一步证明本文所提算法的有效性，对图4～图6 中所有算法融合结果的部分细节和边界区域进行了放大处理，以便在视觉评价方面对比本文算法与其他算法。所有算法的放大结果如图7～图9 所示，其中融合图像中的方框为选中放大的图像区域，而放大后的图像在原图像的右侧展示并由对应的箭头连接。

图7 第一组实验融合结果的部分细节和边界区域放大对比

图9 第三组实验融合结果的部分细节和边界区域放大对比

从图7a～图9a 可以看出，LP 算法的融合结果存在亮度过高的情况，导致该算法会丢失MRI 图像的部分边缘细节。而这种情况同样出现在CVT算法中，如图9c 所示。

对于DWT 算法，从表2 可以看出，其在NAB/F指标上取得了最优结果。然而，从图8b 可以看出，其相对应融合图像的纹理信息和图像对比度与其他算法相比明显较差。这意味着该算法保存的信息较少，从表2 和表3 中也可看出，DWT 算法在与信息保存相关的QAB/F指标上的数值结果是最差的。

图8 第二组实验融合结果的部分细节和边界区域放大对比

对于SR 算法和ASR 算法，尽管在指标评价分析中，两者的融合结果在与评价细节和纹理信息保存相关的 AG、 MI、 SF指标上都能取得优异的数值，如表1～表3 所示；但在视觉分析中，从图7e～图9e 以及图7f～图9f 中可以看出，这两种算法融合结果的放大区域与其他算法相比更加模糊，这是由于图像重建的图像块效应导致的。而对于本文算法，尽管也会受到影响，但与SR 算法和ASR 算法相比有明显的提升。

最后，从整体来看，NSCT 算法和本文所提算法在纹理和边界信息都取得了较优异的视觉效果。但是通过比较图5b 以及图8a～图8f 可以看出，6 种经典算法都存在不同程度的边界损失情况，而通过图8g 可以看出，本文算法与其他6 种算法相比保存了最多的边界信息。因而，与NSCT 算法相比，本文算法更具鲁棒性。

图5 脑膜瘤患者的融合结果

4.3 算法的时间复杂度评价

在算法的时间复杂度方面，与传统的变换域融合算法(LP、DWT、CVT、NSCT) 相比，基于稀疏表示的算法(SR、ASR、本文算法) 更为复杂。对于ASR 算法来说，每个待融合的图像块都要根据其梯度特征自适应地选择最适合的子字典，因而该算法会比SR 算法更加耗时。而在本文提出的算法中，ASR 模型仅用于细节层的融合。因此，为评价SR、ASR 和本文算法的时间复杂度，在Inter Core i7-10875H@2.30 GHz 八核处理器、16.0 GB内存及NVIDIA GeForce RTX 2060 6 GB 显卡的硬件环境下进行实验。每个算法的程序以前文所述的3 组实验图像为输入重复运行10 次，计算每次的运行时间，并取均值作为最终结果来评价该算法的时间复杂度。

3 种算法的平均运行时间如表4 所示。其中，ASR 算法的运行时间在3 组实验中均为最长，而本文算法的运行时间在所有实验中均是最短。除此之外，通过对比SR 和ASR 算法的3 组实验可以看出，不同的输入图像对其运行时间影响较大，尤其是在第三组实验中，SR 和ASR 算法的运行时间均高于前两组实验。而对于本文算法来说，其运行时间的变化较为稳定，受输入图像影响较小。

表4 SR、ASR 和本文算法的平均运行时间s

综上所述，与SR 和ASR 算法相比，本文算法的时间复杂度更低，且不受输入图像影响；然而，与LP、DWT、CVT、NSCT 算法相比，该算法的时间复杂度仍有较大改进空间。

5 结 束 语

本文提出了一种基于自适应稀疏表示和引导滤波的多模态图像融合算法。该算法将待融合图像分解为细节层和基础层；利用基于引导滤波器的加权融合方法融合基础层；同时，利用自适应稀疏表示模型融合细节层。实验结果表明，本文算法在纹理和边界信息保存上均取得了优异结果；同时，时间复杂度也优于基于稀疏表示的方法。然而，本文算法在融合图像的对比度、伪影以及块效应上存在一定缺陷，需要进一步改进；算法的时间复杂度也需要进一步降低。