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可以表示为以下目标函数：

$ u=\underset{u}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left({\left({u}_{p}-{g}_{p}\right)}^{2}+\lambda \left({a}_{x, p}\left(g\right){\left(\frac{\partial u}{\partial x}\right)}_{p}^{2}+ \\ \quad\quad {a}_{y, p}\left(g\right){\left(\frac{\partial u}{\partial y}\right)}_{p}^{2}\right)\right) $

(1)

其中，下标$ p $表示像素的空间位置，（u_p-g_p）²用以保证滤波后的平滑图像u更接近输入图像g，$ \lambda \left({a}_{x, p}\left(g\right){\left(\frac{\partial u}{\partial x}\right)}_{p}^{2}+{a}_{y, p}\left(g\right){\left(\frac{\partial u}{\partial y}\right)}_{p}^{2}\right) $是一个正则化项，通过求u的偏导数来实现平滑，$ \lambda $是正则化参数，用以调整两项之间的平衡，$ {a}_{x, p} $和$ {a}_{y, p} $代表平滑因子。$ {a}_{x, p} $和$ {a}_{y, p} $的计算公式如下：

$ {a}_{x, p}\left(g\right)={\left({\left|\frac{\partial l}{\partial x}\left(p\right)\right|}^{\alpha }+ε \right)}^{-1} $

(2)

$ {a}_{y, p}\left(g\right)={\left({\left|\frac{\partial l}{\partial y}\left(p\right)\right|}^{\alpha }+ε \right)}^{-1} $

(3)

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

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

$ G{\left(I\right)}_{p}=\frac{1}{{K}_{p}}\sum\limits _{q\in \Omega }{g}_{\sigma }\left(‖p-q‖\right){I}_{q} $

(4)

其中，$ {K}_{p}=\sum\limits _{q\in \Omega }{\mathrm{g}}_{\sigma }\left(‖p-q‖\right) $为归一化系数，$ {g}_{\sigma } $是高斯函数，表示为$ {g}_{\sigma }=\frac{1}{\sqrt{2\mathrm{\pi }}\sigma }\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{{x}^{2}}{2{\sigma }^{2}}\right\} $。

1.2 非下采样方向滤波器组

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

1.3 脉冲耦合神经网络

脉冲耦合神经网络（PCNN）是一个M×N的二维网络，其中每一个神经元对应于图像的一个特定的像素。每个神经元由三部分组成，分别是接收域、调制域和脉冲产生器。一种PCNN的简化模型和数学表达式如图 2和式（5）所示。

$ \left\{\begin{array}{l}{F}_{ij}\left(n\right)={S}_{ij}\\ {L}_{ij}\left(n\right)={\mathrm{e}}^{-{\alpha }_{L}}{L}_{ij}(n-1)+{V}_{L}\sum\limits _{k, l}{W}_{ijkl}{Y}_{kl}(n-1)\\ {U}_{ij}\left(n\right)={F}_{ij}\left(n\right)(1+\beta {L}_{ij}(n\left)\right)\\ {\theta }_{ij}\left(n\right)={\mathrm{e}}^{-{\alpha }_{\theta }}{\theta }_{ij}(n-1)+{V}_{\theta }{Y}_{ij}(n-1)\\ {Y}_{ij}\left(n\right)=\left\{\begin{array}{l}1, {U}_{ij}\left(n\right)>{\theta }_{ij}\left(n\right)\\ 0, \mathrm{其}\mathrm{他}\end{array}\right.\\ {T}_{ij}\left(n\right)={T}_{ij}(n-1)+{Y}_{ij}\left(n\right)\end{array}\right. $

(5)

其中：下标L表示图像解等级；$ {F}_{ij}\left(n\right) $和$ {L}_{ij}\left(n\right) $分别是PCNN的馈送输入和链接输入，处于PCNN模型的接收域，$ {F}_{ij}\left(n\right) $通常接收的是图像像素的归一化灰度值，$ {L}_{ij}\left(n\right) $则通过突触权重与8个相邻神经元的先前放电状态相关；$ {a}_{L} $和$ {V}_{L} $是$ {L}_{ij}\left(n\right) $的时间衰减常数和振幅增益；$ {U}_{ij}\left(n\right) $代表内部活动，处于PCNN模型的调制域，通过$ {F}_{ij}\left(n\right) $和$ {L}_{ij}\left(n\right) $进行非线性调制得到，其中参数$ \beta $是连接强度。

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

2 本文医学图像融合方法

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

2.1 多尺度边缘保持分解

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

1）对输入图像进行WLS滤波分解得到细节层和基础层。使用WLS滤波对输入图像进行保边平滑处理得到基础层，细节层则由输入图像与基础图像做差得到。基础层按照式（1）计算。设$ {C}_{\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}} $是输入图像，则细节层表示为：

$ {C}_{\mathrm{D}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}={C}_{\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}}-{C}_{\mathrm{B}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $

(6)

其中，$ {C}_{\mathrm{D}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $和$ {C}_{\mathrm{B}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $是经过WLS滤波分解得到的细节层和基础层。由于WLS滤波的保边特性，因此基础层含有强边缘信息。

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

$ {C}_{\mathrm{E}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}={C}_{\mathrm{B}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}-{C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $

(7)

其中，$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $和$ {C}_{\mathrm{E}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $是基础层经过高斯滤波分解得到的低频层和边缘层。基于WLS滤波和高斯滤波的分解单元结构如图 4所示。

3）将边缘层$ {C}_{\mathrm{E}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $和细节层$ {C}_{\mathrm{D}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $进行叠加，构建图像的高频层$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $。同时为提取高频层系数中存在的大量几何特征，对高频层使用NSDFB进行方向分析，捕获高频层系数不同方向的特征。边缘保持分解框架如图 5所示。

将得到的低频层$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}} $作为下一级分解的输入图像，经过多级分解实现多尺度分解。通过上述分解模式最终将源图像分解为一个低频层和一系列不同方向的高频层方向子带。将源图像A和源图像B通过多尺度分解得到的低频层表示为$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}}\left(i, j\right) $和$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}}\left(i, j\right) $，不同方向的高频层方向子带表示为$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}, i, d}(i, j) $和$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}, i, d}(i, j) $，分别表示源图像A和源图像B在第i级别第d方向上的高频层。

2.2 图像融合规则 2.2.1 低频层融合规则

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

1）根据式（8）计算$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}}\left(i, j\right) $和$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}}\left(i, j\right) $的区域能量，将其作为PCNN的外部激励。

$ \mathrm{L}\mathrm{E}\left(x, y\right)=\sum\limits _{a}\sum\limits _{b}{\left[C\left(i+a, j+b\right)\right]}^{2}W\left(a, b\right) $

(8)

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

$ W=\frac{1}{9}\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right] $

2）计算$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}}\left(i, j\right) $和$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}}\left(i, j\right) $的平均梯度，构建PCNN链接强度系数，如式（9）所示：

$ {\beta }_{ij}=\frac{1}{1+\mathrm{e}\mathrm{x}\mathrm{p}(-\eta \overline{g}(i, j\left)\right)} $

(9)

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

$ \overline{g}(i, j)=\frac{1}{9}{\sum\limits _{a=-1}^{1}\sum\limits _{b=-1}^{1}\left\{\frac{\left[{g}_{1}(i+a, j+b)+{g}_{2}(i+a, j+b)\right]}{2}\right\}}^{\frac{1}{2}} $

(10)

$ {g}_{1}(i, j)={\left[C(i, j)-C(i+1, j)\right]}^{2} $

(11)

$ {g}_{2}(i, j)={\left[C(i, j)-C(i, j+1)\right]}^{2} $

(12)

3）将上述计算的区域能量作为PCNN模型的外部激励，平均梯度构建的链接系数作为PCNN模型的链接强度，通过PCNN模型得到低频子带相应的点火图$ {T}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}}(i, j) $和$ {T}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}}(i, j) $。

4）根据计算得到的点火图可以得到低频层的融合结果$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{F}} $：

$ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{F}}(i, j)=\left\{\begin{array}{l}{C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}}(i, j), {T}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}}(i, j)\ge {T}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}}(i, j)\\ {C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}}(i, j), {T}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}}(i, j)<{T}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}}(i, j)\end{array}\right. $

(13)

2.2.2 高频层融合规则

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

1）根据式（14）~式（15）计算$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}, i, d}(i, j) $和$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}, i, d}(i, j) $改进的空间频率，将计算所得的改进的空间频率D_f作为PCNN的外部激励。

$ \mathrm{N}\mathrm{M}\mathrm{S}\mathrm{F}=\sqrt{\frac{1}{M\left(N-1\right)}\left[\sum\limits _{i=1}^{M}\sum\limits _{j=2}^{N}{\left(C\left(i, j\right)-C\left(i, j-1\right)\right)}^{2}\right]+\frac{1}{N\left(M-1\right)}\left[\sum\limits _{i=2}^{M}\sum\limits _{j=1}^{N}{\left(C\left(i, j\right)-C\left(i-1, j\right)\right)}^{2}\right]+{D}_{\mathrm{f}}} $

(14)

$ {D}_{\mathrm{f}}={\left[\sqrt{\frac{1}{\left(M-1\right)\left(N-1\right)}\left[\sum\limits _{i=2}^{M}\sum\limits _{j=2}^{N}{\left(C\left(i, j\right)-C\left(i-1, j-1\right)\right)}^{2}\right]}+\sqrt{\frac{1}{\left(M-1\right)\left(N-1\right)}\left[\sum\limits _{i=2}^{M}\sum\limits _{j=2}^{N}{\left(C\left(i-1, j\right)-C\left(i, j-1\right)\right)}^{2}\right]}\right]}^{2} $

(15)

2）根据式（10）计算$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}, i, d}(i, j) $和$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}, i, d}(i, j) $的PCNN链接强度。

3）将上述计算的改进的空间频率作为PCNN模型的外部激励，将平均梯度构建的链接系数作为PCNN模型的链接强度，通过PCNN模型得到高频方向子带相应的点火图$ {T}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}, i, d}(i, j) $和$ {T}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}, i, d}(i, j) $。

4）根据计算得到的点火图可以得到高频方向子带的融合结果$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{F}, i, d} $：

$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{F}, i, d}(i, j)=\left\{\begin{array}{l}{C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}, i, d}(i, j), {T}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}, i, d}(i, j)\ge {T}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}, i, d}(i, j)\\ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}, i, d}(i, j), {T}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{A}, i, d}(i, j)<{T}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{B}, i, d}(i, j)\end{array}\right. $

(16)

通过NDFB逆变换，可以得到第i层的高频层融合结果$ {C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{F}, i} $。

2.3 图像重构

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

$ {F}_{\mathrm{s}}={C}_{\mathrm{L}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{F}}+\sum\limits _{i=1}^{L}{C}_{\mathrm{H}\_\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}^{\mathrm{F}, i} $

(17)

3 实验结果与分析 3.1 实验参数设置

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

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

在图像融合阶段，本文参考文献[9]，将用于融合图像的PCNN的参数设置为：$ {a}_{L}=10 $，$ {a}_{\theta }=0.815 $，$ {V}_{L}=1.0 $，$ {V}_{\theta }=10 $，$ n=200 $，$ {{\mathit{\boldsymbol{W}}}}_{ijkl}=\left[\begin{array}{ccc}0.707\mathrm{ }1& 1& 0.707\mathrm{ }1\\ 1& 0& 1\\ 0.707\mathrm{ }1& 1& 0.707\mathrm{ }0\end{array}\right] $。

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所示。

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

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

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

3.2.2 不同滤波器的比较

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

4 结束语

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