基于多尺度边缘保持分解与PCNN 的医学图像融合

郭淑娟，高 媛，秦品乐，王丽芳

（中北大学大数据学院山西省生物医学成像与影像大数据重点实验室，太原 030051）

0 概述

图像融合是指将不同模态的相关图像合并为一幅新图像。在此过程中，要求生成的融合图像具有良好的清晰度，并使不同模态的图像特征尽可能融合到一幅图像中。图像融合在多模态医学图像中的应用是一个热门的研究方向。随着现代医学成像技术的发展，结合不同的医学成像方式能够获得器官和组织的全面信息，例如计算机断层扫描（Computed Tomography，CT）图像仅显示高密度组织，磁共振（Magnetic Resonance，MR）图像仅显示高水分组织，正电子发射断层扫描（Positron Emission Tomography，PET）图像仅反映不同人体之间组织强度的差异［1］，而将这些医学图像的互补特征融合到一幅图像中，能够提高医生的诊断效率和准确性。

图像融合技术可分为空间域方法和变换域方法两类，其中，变换域方法可以更好地表征图像特征，弥补空间域方法在细节提取方面的不足。在变换域方法中，多尺度分析方法通过将源图像映射至不同的尺度空间提取潜在的重要信息，符合人眼视觉系统的生理机制，被认为是一种主流的图像融合方法［2］。现有的多尺度分析方法主要有基于离散小波变换［3］、基于剪切波变换［4］、基于轮廓波变换［5］、基于非下采样轮廓波变换（Non-Subsampled Contourlet Transform，NSCT）［6］和基于非下采样剪切波变换（Non-Subsampled Shearlet Transform，NSST）［7］的算法，其中，基于NSCT 和NSST 的算法具有平移不变性，能够很好地抑制伪吉布斯现象。

脉冲耦合神经网络（Pulse Coupled Neural Network，PCNN）是一种通过研究哺乳动物视觉皮层而获得的具有仿生机制的单层神经网络模型［8］，其不需要学习或者训练，且具有全局耦合性和脉冲发射特性，能够提取图像的全局特征并且增强图像细节，有利于图像的实时处理。为提高融合图像的视觉感知效果并充分利用图像的空间信息，一些多尺度分析与PCNN 相结合的图像融合方法相继被提出［9］，如NSCT-SF-PCNN［10］和NSST-PAPCNN［11］。然而，传统的多尺度分析方法大多采用线性滤波器，由于无法保留边缘从而导致分解阶段的强边缘处出现模糊，使得边缘处产生光晕［12］。

为解决该问题，研究人员将保边滤波器引入图像融合领域。保边滤波器具有平移不变性、空间一致性和边缘保护性能，而且计算效率较高。在此方面：文献［13］提出一种多尺度双边滤波器对源图像进行分解；文献［14］提出一种基于引导滤波的图像融合方案，其对源图像进行两尺度分解，通过引导滤波优化权重实现了空间一致性；文献［15］提出以多尺度极值滤波对源图像进行多尺度分解，采用对比度融合规则使融合图像更符合人眼视觉感知特性；文献［16］提出一种加权最小二乘滤波与引导滤波相结合的新方案，与只基于引导滤波的方法相比，该方案可以获得更好的效果。尽管上述方法在一定程度上保护了边缘信息，但是基础层主要包含粗尺度结构信息，直接对其融合可能在之后的融合过程中丢失或模糊保边滤波保留的显著的边缘特征［17］。

为保留图像完整清晰的边缘信息，提高融合图像的视觉感知效果，本文提出一种新的医学图像融合方法。将多尺度边缘保持分解方法引入图像分解过程，对已配准的源图像进行分解得到低频层和高频层。在此基础上，利用改进的空间频率以及区域能量激励PCNN 对两者进行融合，并通过逆变换获得最终的融合图像。

1 相关理论

1.1 加权最小二乘滤波与高斯滤波

文献［18］提出一种基于加权最小二乘（Weighted Least Square，WLS）优化框架的非线性滤波器。给定一个输入图像g，WLS 滤波器的目的是使滤波后的平滑图像u尽可能与输入图像g近似，则滤波后的图像u可以表示为以下目标函数：

其中，下标p表示像素的空间位置，（up-gp）2用以保证滤波后的平滑图像u更接近输入图像g，是一个正则化项，通过求u的偏导数来实现平滑，λ是正则化参数，用以调整两项之间的平衡，ax，p和ay，p代表平滑因子。ax，p和ay，p的计算公式如下：

其中，l是输入图像g的对数亮度通道，指数α确定对g的梯度的敏感度，而ε是一个非常小的常数，通常为0.000 1，可防止在g恒定的区域被零除。

高斯滤波是一种线性滤波，由于没有考虑相邻像素的影响，其对图像滤波时去除高频信息的同时也模糊了边缘，因此在图像的平滑和去噪中被广泛应用。高斯滤波公式可表示为：

1.2 非下采样方向滤波器组

文献［19］提出一种方向滤波器组（Directional Filter Bank，DFB）。DFB 将图像划分为具有方向性的子带，捕捉图像几何特征，但是其中的下采样操作使其缺乏平移不变性。非下采样方向滤波器组（Non-subsampled Directional Filter Bank，NDFB）［20］与DBF 类似，可以对图像进行多方向分解，提取图像方向细节特征，但是其进行上采样操作，因而能够获得位移不变的方向扩展。NDFB 是由扇形滤波器和棋盘滤波器构成的四方向滤波器组，基本结构如图1 所示。

图1 NDFB 基本结构Fig.1 Basic structure of NDFB

1.3 脉冲耦合神经网络

图2 PCNN 的简化模型Fig.2 Simplified model of PCNN

脉冲耦合神经网络（PCNN）是一个M×N的二维网络，其中每一个神经元对应于图像的一个特定的像素。每个神经元由三部分组成，分别是接收域、调制域和脉冲产生器。一种PCNN 的简化模型和数学表达式如图2 和式（5）所示。其中：下标L表示图像解等级；Fij(n)和Lij(n)分别是PCNN 的馈送输入和链接输入，处于PCNN 模型的接收域，Fij(n)通常接收的是图像像素的归一化灰度值，Lij(n)则通过突触权重与8 个相邻神经元的先前放电状态相关；aL和VL是Lij(n)的时间衰减常数和振幅增益；Uij(n)代表内部活动，处于PCNN 模型的调制域，通过Fij(n)和Lij(n)进行非线性调制得到，其中参数β是连接强度。

脉冲发生域控制脉冲的发生：将Uij与动态阈值θij进行比较，若Uij＞θij，Yij=1 即神经元点火，产生脉冲输出；反之Yij=0，aθ和Vθ是θij的时间衰减常数和振幅增益。将上述过程多次迭代直到满足设定的迭代条件，把此时由神经元点火形成的点火映射图作为输出结果。

2 本文医学图像融合方法

对已配准的源图像进行图像融合的主要步骤包括图像分解、图像融合和图像重构。本文在图像分解过程中引入多尺度边缘保持分解方法，对已配准的源图像进行分解得到低频层和高频层，而在图像融合过程中针对得到的低频层和高频层不同的特点，分别采用不同的融合规则进行融合，在图像重构中获得最终的融合的图像。融合过程如图3所示。

图3 医学图像融合过程Fig.3 Process of medical image fusion

2.1 多尺度边缘保持分解

WLS 滤波器可以有效地在模糊和锐化之间进行最佳折衷，更好地保留边缘，解决双边滤波在边缘处容易产生梯度反转的问题，弥补双边滤波无法很好地在多尺度上提取到细节信息的不足，从而更准确地提取多尺度信息。由于高斯滤波器只考虑像素的空间分布而没有考虑邻近像素的影响，因此具有对图像进行滤波时能够剔除细节信息、模糊边缘的特点。结合两种滤波的特性，本文对源图像进行多尺度分解。分解过程由以下步骤组成：

1）对输入图像进行WLS 滤波分解得到细节层和基础层。使用WLS 滤波对输入图像进行保边平滑处理得到基础层，细节层则由输入图像与基础图像做差得到。基础层按照式（1）计算。设Cinput是输入图像，则细节层表示为：

其中，CD_layer和CB_layer是经过WLS 滤波分解得到的细节层和基础层。由于WLS 滤波的保边特性，因此基础层含有强边缘信息。

2）为避免光谱损失，使用高斯滤波器再次分解基础层提取边缘信息，得到低频层和边缘层。低频层图像依据式（4）计算，则边缘层计算公式为：

其中，CL_layer和CE_layer是基础层经过高斯滤波分解得到的低频层和边缘层。基于WLS 滤波和高斯滤波的分解单元结构如图4 所示。

图4 基于WLS 滤波和高斯滤波的分解单元结构Fig.4 Decomposition unit structure based on WLS filtering and Gaussian filtering

3）将边缘层CE_layer和细节层CD_layer进行叠加，构建图像的高频层CH_layer。同时为提取高频层系数中存在的大量几何特征，对高频层使用NSDFB 进行方向分析，捕获高频层系数不同方向的特征。边缘保持分解框架如图5 所示。

图5 边缘保持分解框架Fig.5 Frame of edge-preserving decomposition

将得到的低频层CL_layer作为下一级分解的输入图像，经过多级分解实现多尺度分解。通过上述分解模式最终将源图像分解为一个低频层和一系列不同方向的高频层方向子带。将源图像A 和源图像B通过多尺度分解得到的低频层表示为和，不同方向的高频层方向子带表示为和，分别表示源图像A 和源图像B 在第i级别第d方向上的高频层。

2.2 图像融合规则

2.2.1 低频层融合规则

低频层图像是源图像去除高频信息之后的近似图像，包含图像大量的能量信息，其中子带相邻系数之间存在区域相关性［21］。区域能量将邻域内每个中心像素的能量表征为其自身和邻域内像素值的平方和，考虑了图像局部特征，减少了融合图像中的不连续性，边缘信息在该融合模式下能够很好地得以保留。因此，本文将区域能量作为PCNN的外部激励融合低频子带。首先计算低频层的区域能量作为PCNN 的外部激励，然后使用表征图像清晰度的平均梯度作为PCNN 的连接强度，最后通过PCNN 输出的点火图进行融合。低频层融合过程如下：

其中，W是一个3×3 的滑动窗口。为体现低频层子带的平滑特性，选取W为：

其中，η是一个大于0 的常数，用于调节链接强度的值，本文设置为0.2(i，j)是计算所得的平均梯度。的计算公式如下：

3）将上述计算的区域能量作为PCNN 模型的外部激励，平均梯度构建的链接系数作为PCNN 模型的链接强度，通过PCNN 模型得到低频子带相应的点火图

2.2.2 高频层融合规则

高频层反映了图像的纹理细节和边缘信息，对最终的融合结果有着至关重要的影响。改进的空间频率能够反映图像的边缘和细节信息，同时表示图像灰度值的变化，因此，本文采用改进的空间频率激PCNN 对高频层进行融合。首先计算每一级不同方向的高频层方向子带改进的空间频率作为PCNN 的外部刺激，然后使用表征图像清晰度的平均梯度作为PCNN 的连接强度，最后通过PCNN 输出的点火图进行融合。高频层融合过程如下：

3）将上述计算的改进的空间频率作为PCNN 模型的外部激励，将平均梯度构建的链接系数作为PCNN 模型的链接强度，通过PCNN 模型得到高频方向子带相应的点火图

通过NDFB 逆变换，可以得到第i层的高频层融合结果

2.3 图像重构

将得到的高频层融合结果和低频层融合结果重构为最终的融合结果：

3 实验结果与分析

3.1 实验参数设置

本文实验数据选取自哈佛大学医学院的公开数据，配准的医学图像是大小为256 像素×256 像素、颜色深度为8 bit 的灰度图。

在图像分解阶段，基于WLS 分解的参数主要包括平滑参数λ和边缘保留参数α。λ是一个正则化参数，用于保持数据项和平滑度项之间的平衡，增加λ可产生更平滑的图像。α通过非线性缩放梯度来确定对梯度的敏感程度，增加α可保留更清晰的边缘。本文参考文献［16］，设置λi=[0.075，0.6，4.8]，α=0.6，其中，i=1，2，…，L。基于高斯滤波分解的主要参数是用于控制图像平滑效果的高斯分布的标准差σ，第1 级设置为σ1=1.1，其他级设置为σi+1=kσi，其中，k是相邻级之间的分解比例因子，设置为k=2。图像分解等级设置为L=1，2，3，方向分析的分解水平分别设置为第1 层方向数为16 方向，第2 层和第3 层分解分别为8 方向和4 方向，采用的方向滤波器为“vk”，将第L级的得到的低频图像作为下一级分解的输入图像以实现多尺度分解。

在图像融合阶段，本文参考文献［9］，将用于融合图像的PCNN 的参数设置为：aL=10，aθ=0.815，VL=1.0，Vθ=10，n=200

3.2 实验结果分析

3.2.1 不同融合方法的比较

本文采用多发性脑梗塞MR-T1 和MR-T2 图像、脑弓形虫CT 和MR 图像、正常脑部脑CT 和MR 图像进行3 组仿真实验，在64 位Windows10 操作系统、Matlab2014b 平台上进行仿真。选取以下5 种融合方法作为对比方法：1）基于NSCT 的图像融合方法（NSCT）；2）基于多尺度变换和稀疏表示的图像融合方法（MST-SR）［22］；3）基于NSCT域下空间频率激励PCNN的图像融合方法（NSCT-SF-PCNN）［10］；4）基于NSST 与自适应PCNN 的图像融合方法（NSST-PAPCNN）［11］；5）基于引导滤波的图像融合方法（GFF）［14］。NSCT采用高频子带区域能量取大、低频子带加权平均的融合规则，MST-SR 的多尺度分解方法选取为LP，分解等级为4。6 种方法的实验结果如图6～图8 所示。

图6 不同方法的多发性脑梗塞MR-T1/MR-T2图像融合结果Fig.6 Fusion results for MR-T1/MR-T2 image of multiple cerebral infarction by different methods

图7 不同方法的脑弓形虫CT/MR 图像融合结果Fig.7 Fusion results for CT/MR image of toxoplasma gondii by different methods

图8 不同方法的正常脑部CT/MR 图像融合结果Fig.8 Fusion results for CT/MR image of normal brain by different methods

由实验结果可以看出：NSCT 和NSCT-SF-PCNN的融合结果丢失了骨质信息，同时NSCT 的融合图像整体亮度偏低；NSCT-SF-PCNN 在第3 组实验中的融合结果MR 信息严重缺失并且出现边缘模糊的现象；NSST-PAPCNN 的融合结果虽然较好地保存了骨质信息，但是在边缘处产生了边缘模糊和光晕现象；MST-SR 能够较好地保存骨质信息和边缘信息，但是在第2 组实验中MST-SR 融合图像MR信息有所缺失，存在对比度下降和细节丢失的问题；GFF 的方法能够有效保护边缘信息，但是第1 组和第2 组实验结果存在明显的细节丢失现象。

为更客观地评价图像的融合结果，选取边缘评价因子QAB/F、平均梯度、互信息、标准差和空间频率5 种指标进行客观评价。3 组实验的客观评价指标如表1～表3 所示，其中加粗数据表示最优值。可以看出：在5 种评价指标中，各组实验中本文方法的QAB/F均达到最大值，验证了其在边缘保持方面的可行性和有效性；此外，除第1 组实验的空间频率和第2 组、第3 组实验的标准差之外，本文方法其余各实验指标均达到最大值，表明其能较好地保留原图像信息，体现图像的纹理细节。

表1 不同方法的多发性脑梗塞MR-T1/MR-T2 图像融合评价指标Table 1 Evaluation indexes of MR-T1/MR-T2 image fusion of multiple cerebral infarction by different methods

表2 不同方法的脑弓形虫CT/MRI 图像融合评价指标Table 2 Evaluation indexes of CT/MRI image fusion of toxoplasma gondii by different methods

表3 不同方法的正常脑部CT/MRI 图像融合评价指标Table 3 Evaluation indexes of CT/MRI image fusion of normal brain by different methods

通过对以上主观视觉和客观指标的分析可知，本文方法能够很好地保留边缘信息，所得图像具有较高的对比度和丰富的细节信息，较对比方法能够得到更好的视觉效果。

3.2.2 不同滤波器的比较

在本文方法中使用不同的保边滤波器进行实验，选取引导滤波器（Guided Filter，GF）和双边滤波器（Bilateral Filter，BF）与上文使用的WLS 滤波器进行对比。双边滤波器参数设置为δr=2，4，8，δs=0.1，0.2，0.4，引导滤波器参数设置为r=2，4，8，ε=0.12，0.22，0.42。表4～表6 列出了使用不同滤波器的客观定量指标，其中，加粗数据表示最优值。可以看出，使用WLS 滤波器效果较使用其他滤波器更好。

表4 本文方法不同滤波器下的多发性脑梗塞MR-T1/MR-T2图像融合评价指标Table 4 Evaluation indexes of MR-T1/MR-T2 image fusion of multiple cerebral infarction by the proposed method using different filters

表5 本文方法不同滤波器下的脑弓形虫CT/MR 图像融合评价指标Table 5 Evaluation indexes of CT/MR image fusion oftoxoplasma gondii by the proposed method using different filters

表6 本文方法不同滤波器下的正常脑部CT/MR 图像融合评价指标Table 6 Evaluation indexes of CT/MR image fusion of normal brain by the proposed method using different filters

4 结束语

传统基于多尺度分析的图像融合方法因无法保留边缘特征而导致边缘处产生光晕。针对该问题，本文提出一种多尺度边缘保持分解与PCNN 结合的医学图像融合方法。在分解过程保护边缘的同时进行多尺度多方向分解，并采用PCNN 作为融合规则，从而改善融合图像的视觉感知效果。实验结果表明，该方法能够突出医学图像的边缘轮廓并增强图像细节，有利于将更多的显著特征从源图像分离并转移到融合图像中。本文方法使用了大量参数，而参数设置将影响融合结果的质量，因此，下一步将从参数自适应角度着手对其进行优化。