0 引言

在医学领域，医生需要对同时具有高空间和高光谱信息的单幅图像进行研究和分析，以便对疾病进行准确诊断和治疗^[1]。这种类型的信息仅从单模态图像中无法获取，例如：计算机断层成像(Computed Tomography, CT)能够捕捉人体的骨结构，具有较高的分辨率，而磁共振成像(Magnetic Resonance, MR)能够捕捉人体器官的软组织如肌肉、软骨、脂肪等细节信息。因此，将CT和MR图像的互补信息相融合以获取更全面丰富的图像信息，可为临床诊断和辅助治疗提供有效帮助^[2]。

目前应用于脑部医学图像融合领域比较经典的方法是基于多尺度变换的方法：离散小波变换(Discrete Wavelet Transform, DWT)^[3]、平稳小波变换(Stationary Wavelet Transform, SWT)^[4]、双树复小波变换(Dual-Tree Complex Wavelet Transform, DTCWT)^[5]、拉普拉斯金字塔(Laplacian Pyramid, LP)^[6]、非下采样Contourlet变换(Non-Subsampled Contourlet Transform, NSCT)^[7]。基于多尺度变换的方法能够很好地提取图像的显著特征，并且计算效率较高，但是多尺度分解水平不易确定，分解层数太少则不能提取足够多的细节信息，分解层数太多则融合过程对图像误配准更加敏感，分解层数需要在细节提取能力和对误配准鲁棒性之间做出折中。另外，传统的融合策略也会导致融合图像对比度丢失。随着压缩感知^[8]的兴起，基于稀疏表示的方法^[9-10]被广泛用于图像融合领域，并取得了极佳的融合效果。Yang等^[11]使用冗余的离散余弦变换(Discrete Cosine Transform, DCT)字典稀疏表示源图像，并使用“选择最大”的规则融合稀疏系数。DCT字典是一种隐式字典，易于快速实现，但是其表示能力有限。Aharon等^[12]提出K奇异值分解(K-means-based Singular Value Decomposition, K-SVD)算法用于从训练图像中学习字典。与DCT字典相比，学习的字典是一种自适应于源图像的显式字典，具有较强的表示能力。学习的字典中将仅从自然图像中采样训练得到的字典称为单字典，单字典可以表示任意一幅与训练样本类别相同的自然图像，但对于结构复杂的脑部医学图像，使用单字典既要表示CT图像又要表示MR图像，难以得到精确的稀疏表示系数。Ophir等^[13]提出小波域上的多尺度字典学习方法，即在小波域上对所有子带分别使用K-SVD算法训练得到所有子带对应的子字典。多尺度字典有效地将解析字典和学习的字典的优势相结合，能够捕捉图像在不同尺度和不同方向上包含的不同特征。但是所有子带的子字典也是单字典，使用子字典对所有子带进行稀疏表示仍然难以得到精确的稀疏表示系数，并且分离的字典学习时间效率较低。Yu等^[14]提出基于联合稀疏表示的图像融合方法兼具去噪功能。这种方法是对待融合的已配准源图像本身进行字典学习，根据联合稀疏模型的第一种模型(Joint Sparsity Model-1, JSM-1) 提取待融合图像的公共特征和各自特征，再分别组合并重构得到融合图像。这种方法由于是对待融合源图像本身训练字典所以适用于脑部医学图像，可以得到精确的稀疏表示系数；但是对于每对待融合的源图像都需要训练字典，时间效率低，缺乏灵活性。

针对上述问题，本文提出一种适用于脑部CT/MR图像融合的方法。不同于单字典的样本采集方式，该方法首先将高质量的脑部CT/MR图像对作为训练集，使用改进的K-SVD算法联合训练分别得到耦合的CT字典和MR字典，再结合空间域的方法融合CT和MR字典得到融合字典。融合字典很好地保留了CT和MR训练图像的特征，使用这样一个融合字典就可以同时准确地表示CT和MR源图像并得到精确的稀疏表示系数。其次，在耦合的特征空间下使用改进的K-SVD算法联合训练字典较其他几种经典字典学习方法提高了字典学习的时间效率。最后，在融合阶段将源图像列向量的均值去掉，这样只保留了图像纹理信息而不是图像的强度值，再使用不同的融合规则组合重构各部分得到融合图像。实验结果表明本文提出的方法在主观视觉和客观参数上性能较优。

1 稀疏表示和字典学习

稀疏表示理论的基本假设是：任何给定的自然信号都可以表示为字典中原子的稀疏线性组合^[15]。对给定的信号x∈Rⁿ，信号x的稀疏表示可以描述为x=Dα，其中矩阵D=(d₁, d₂, …, d_m)∈R^n×m是一个字典，每个列向量d_i∈Rⁿ(i=1, 2, …, m)是字典D的一个原子，向量α∈R^m是含有一些非零元素的稀疏系数。稀疏系数α的计算式如下：

$\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{\alpha }}} {\left\| {\mathit{\boldsymbol{\alpha }}} \right\|_0}$

(1)

s.t. ‖x-Dα‖₂²≤ε

式中：向量α的稀疏水平是由l₀范数衡量的；‖α‖₀表示向量α中非零元素的个数；ε表示允许偏差的精度。问题(1) 的求解过程称为“稀疏编码”。

在稀疏表示中，字典D的选择是非常关键的。构造字典的两种主要方法是：解析方法和基于学习的方法。与解析方法相比，基于学习的方法可以提高字典的性能和灵活性，应用范围更广。

具体地，令{x_i}(i=1, 2, …, Num)是一个包含Num个样本的训练集，把这些样本排列成列组成一个矩阵X=(x₁, x₂, …, x_Num)∈R^n×Num。字典学习就是从X中学习一个字典D∈R^n×m，则X=DA，A=(α₁, α₂, …, α_Num)∈R^m×Num是稀疏系数矩阵。最有代表性的字典学习算法是最优方向法(Method of Optimal Directions, MOD)^[16]和K-SVD算法^[12]。这两种学习算法都是通过最小化稀疏逼近误差‖X-DA‖_F²学习字典，对应的优化问题如下：

$\mathop {{\rm{min}}}\limits_{{\mathit{\boldsymbol{D}}}, {\mathit{\boldsymbol{A}}}} \left\| {{\mathit{\boldsymbol{X}}} - {\mathit{\boldsymbol{DA}}}} \right\|_\rm{F}^2\;\;$

(2)

${\rm{s}}{\rm{.t}}{\rm{.}}\;\;\forall i\;{\left\| {{{\mathit{\boldsymbol{\alpha }}}_i}} \right\|_0} \le \tau $

式中：符号‖·‖_F表示F范数；τ是稀疏矩阵A的稀疏度。MOD和K-SVD都是迭代算法，通过交替地执行稀疏编码和字典更新这两个步骤得到问题(2) 的解。

2 改进的耦合字典学习和融合字典

精确的稀疏表示是图像融合成功的关键。经典的模型是从自然图像当中训练一个超完备字典来描述源图像。与自然图像相比，医学图像具有灰度不均匀、对比度低等特点，使用这样的超完备字典并不能精确地表示CT和MR图像。基于此，本文提出使用改进的K-SVD算法联合训练分别得到耦合的CT字典和MR字典，耦合的字典学习能够更高效完整地提取CT和MR图像的特征，然后再结合空间域的方法融合CT字典和MR字典得到融合字典。由于融合字典同时包含了CT和MR图像的特征，所以使用融合字典表示这两种源图像能够得到更精确的稀疏表示系数，从而提高脑部医学图像融合的质量。本文将提出的基于改进耦合字典学习的脑部CT/MR图像融合方法简记为ICDL (Improved Coupled Dictionary Learning)。图 1为ICDL方法的流程。

2.1 改进的耦合字典学习

本文使用的训练集是8个已经配准的CT和MR图像对，这些图像对来自哈佛大学医学院^[17]，如图 2所示。从训练集中采样得到向量对{X^C, X^R}，定义X^C$ \buildrel \Delta \over = $(x₁^C, x₂^C, ...x_n^C)∈R^d×n为n个采样的CT图像向量组成的矩阵，X^R$ \buildrel \Delta \over = $(x₁^R, x₂^R, ...x_n^R)∈R^d×N作为对应的n个采样的MR图像向量组成的矩阵。根据文献[18]和^[19]，本文的目标是在耦合的特征空间下训练两个字典D^C, D^R∈R^d×N，使得CT和MR训练图像在相应的字典作用下具有相同的稀疏表示系数A∈R^N×n。图 3为ICDL方法的字典学习过程。

传统的耦合字典学习问题可表示为：

$\mathop {{\rm{min}}}\limits_{{{\mathit{\boldsymbol{D}}}^{\rm{C}}}, {{\mathit{\boldsymbol{D}}}^{\rm{R}}}, {\mathit{\boldsymbol{A}}}} \left\| {{{\mathit{\boldsymbol{X}}}^{\rm{C}}} - {{\mathit{\boldsymbol{D}}}^{\rm{C}}}{\mathit{\boldsymbol{A}}}} \right\|_2^2 + \left\| {{{\mathit{\boldsymbol{X}}}^{\rm{R}}} - {{\mathit{\boldsymbol{D}}}^{\rm{R}}}{\mathit{\boldsymbol{A}}}} \right\|_2^2\;$

(3)

s.t. ‖A‖₀≤τ, ‖d_i^C‖₂²≤1, ‖d_i^R‖₂²≤1, ∀1≤i≤N其中：A是X^C和X^R的联合稀疏编码；d_i^C和d_i^R分别是D^C和D^R字典原子的第i列；τ是稀疏矩阵A的稀疏度。

本文的耦合字典训练使用改进的K-SVD算法^[20]，该算法是在字典学习代价函数式(3) 中加入支撑完整的先验信息^[21]，交替地更新D^C、D^R和A，对应的训练优化问题如下：

$\mathop {{\rm{min}}}\limits_{{{\mathit{\boldsymbol{D}}}^{\rm{C}}}, {{\mathit{\boldsymbol{D}}}^{\rm{R}}}, {\mathit{\boldsymbol{A}}}} \left\| {{{\mathit{\boldsymbol{X}}}^{\rm{C}}} - {{\mathit{\boldsymbol{D}}}^{\rm{C}}}{\mathit{\boldsymbol{A}}}} \right\|_2^2 + \left\| {{{\mathit{\boldsymbol{X}}}^{\rm{R}}} - {{\mathit{\boldsymbol{D}}}^{\rm{R}}}{\mathit{\boldsymbol{A}}}} \right\|_2^2$

(4)

s.t. ‖A‖₀≤τ, A⊙M=0

式中：⊙代表点乘；掩膜矩阵M由元素0和1组成，定义为M={|A|=0}，等价于如果A (i, j)=0则M (i, j)=1，否则为0。因此A⊙M=0能使A中所有零项保持完备。引入辅助变量：

${\mathit{\boldsymbol{\bar X}}} = \left( \begin{array}{l} {{\mathit{\boldsymbol{X}}}^{\rm{C}}}\\ {{\mathit{\boldsymbol{X}}}^{\rm{R}}} \end{array} \right)$

(5)

${\mathit{\boldsymbol{\bar D}}} = \left( \begin{array}{l} {{\mathit{\boldsymbol{D}}}^{\rm{C}}}\\ {{\mathit{\boldsymbol{D}}}^{\rm{R}}} \end{array} \right)$

(6)

则问题式(4) 可以等价地转化为：

$\begin{array}{l} \quad \quad \mathop {{\rm{min}}}\limits_{{\mathit{\boldsymbol{\bar D}}}, {\mathit{\boldsymbol{A}}}} \left\| {{\mathit{\boldsymbol{\bar X}}} - {\mathit{\boldsymbol{\bar DA}}}} \right\|_2^2\\ {\rm{s}}{\rm{.t}}.\quad \;{\left\| {\mathit{\boldsymbol{A}}} \right\|_0} \le \tau , \;\;{\mathit{\boldsymbol{A}}} \odot {\mathit{\boldsymbol{M}}} = 0 \end{array}$

(7)

式(7) 的求解过程分为稀疏编码和字典更新两个步骤。

首先，随机矩阵初始化字典D₀^C和D₀^R。通过求解式(8) 来实现对系数矩阵A的更新：

$\begin{array}{l} \quad \quad {\mathit{\boldsymbol{A}}} = \mathop {{\rm{arg \ min}}}\limits_{\mathit{\boldsymbol{A}}} \left\| {{\mathit{\boldsymbol{\bar X}}} - {\mathit{\boldsymbol{\bar DA}}}} \right\|_2^2\\ {\rm{s}}{\rm{.t}}.\;\;{\mathit{\boldsymbol{A}}} \odot {\mathit{\boldsymbol{M}}} = 0 \end{array}$

(8)

分别对系数矩阵A中每一列的非零元素进行处理，而保持零元素完备，则式(8) 可以转换为式(9)：

${{\mathit{\boldsymbol{\alpha }}}_i} = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_{{{\mathit{\boldsymbol{\alpha }}}_i}} \left\| {{{{\mathit{\boldsymbol{\bar x}}}}_i} - {{{\mathit{\boldsymbol{\bar D}}}}_i}{{\mathit{\boldsymbol{\alpha }}}_i}} \right\|_2^2$

(9)

式中:D_i是D对应A的非零支集的子矩阵; α_i是A第i列的非零部分。问题(9) 由系数重用正交匹配追踪(Coefficient Reuse Orthogonal Matching Pursuit, CoefROMP)算法^[19]求解，由此可得到更新的稀疏系数矩阵A。

其次，在字典更新阶段，式(7) 的优化问题可以转化为式(10)：

$\begin{array}{l} \quad \quad \left\{ {{\mathit{\boldsymbol{\bar D}}}, {\mathit{\boldsymbol{A}}}} \right\} = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_{{\mathit{\boldsymbol{\bar D}}}, {\mathit{\boldsymbol{A}}}} \left\| {{\mathit{\boldsymbol{\bar X}}} - {\mathit{\boldsymbol{\bar DA}}}} \right\|_2^2\\ {\rm{s}}{\rm{.t}}.\;\;{\mathit{\boldsymbol{A}}} \odot {\mathit{\boldsymbol{M}}} = 0 \end{array}$

(10)

则式(10) 的补偿项可以写为：

$\begin{array}{l} \left\| {{\mathit{\boldsymbol{\bar X}}} - {\mathit{\boldsymbol{\bar DA}}}} \right\|_2^2 = \left\| {{\mathit{\boldsymbol{\bar X}}} - \sum\limits_{j = 1}^N {{{{\mathit{\boldsymbol{\bar d}}}}_j}} {\mathit{\boldsymbol{\alpha }}}_j^T} \right\|_2^2 = \\ \;\;\;\;\;\;\;\;\;\;\left\| {\left( {{\mathit{\boldsymbol{\bar X}}} - \sum\limits_{j \ne k} {{{{\mathit{\boldsymbol{\bar d}}}}_j}{\mathit{\boldsymbol{\alpha }}}_j^T} } \right) \odot \left( {{1_d} \cdot {\mathit{\boldsymbol{m}}}_j^T} \right) - {{{\mathit{\boldsymbol{\bar d}}}}_k}{\mathit{\boldsymbol{\alpha }}}_k^T} \right\|_2^2 = \\ \;\;\;\;\;\;\;\;\;\;\left\| {{{\mathit{\boldsymbol{E}}}_k} - {{{\mathit{\boldsymbol{\bar d}}}}_k}{\mathit{\boldsymbol{\alpha }}}_k^T} \right\|_2^2 \end{array}$

(11)

式中:d_k表示字典D中待更新的第k列; α_k^T表示稀疏系数矩阵A的第k行; m_j^T表示掩膜矩阵M的第j行用来保证α_k^T中的零元素在正确的位置。掩膜矩阵1_d·m_j^T是将行向量m_j^T复制d次得到尺寸为d×n的秩为1的矩阵，利用掩膜矩阵1_d·m_j^T可以有效地去除${\mathit{\boldsymbol{\bar X}}} - \sum\limits_{j \ne k} {{{{\mathit{\boldsymbol{\bar d}}}}_j}} {\mathit{\boldsymbol{a}}}_j^T$中那些未用到第k个原子所对应样本的列。对误差矩阵E_k进行奇异值分解(SVD)得到E_k=UΔV^T，使用矩阵U的第一列更新字典D中的原子d_k，同时将稀疏系数矩阵A中的α_k^T更新为矩阵V的第一列与Δ(1, 1) 的乘积。

最后，循环执行稀疏编码和字典更新这两个阶段，直至达到预设的迭代次数为止，输出一对耦合的D^C和D^R字典。由于字典更新阶段同时更新字典和稀疏表示系数的非零元素，使得字典的表示误差更小且字典的收敛速度更快。在稀疏编码阶段，考虑到每次迭代时都忽略前一次迭代的表示，CoefROMP算法提出利用上次迭代的稀疏表示残差信息进行系数更新，从而更快地得到所要求问题的解^[21]。

2.2 融合字典

为构造一个同时包含CT和MR两种图像特征的字典，使得这个字典能够同时精确地表示CT和MR源图像。考虑到2.1节中CT和MR训练图像在一对耦合字典的作用下具有相同的稀疏表示系数，因此将CT和MR字典中的原子作为训练图像的特征，再结合空间域的方法融合CT和MR字典得到融合字典。

具体地，令D^C, D^R∈R^d×N，L^C(n)和L^R(n)(n=1, 2, …, N)分别代表CT字典和MR字典的第n个原子的特征指标。由于脑部CT和MR图像是对应于人体同一部位由不同成像设备获取的图像，因此两者之间一定存在着公共特征和各自特征。本文提出：将特征指标相差较大的原子看作各自特征，使用“选择最大”规则融合；特征指标相差较小的原子看作公共特征，使用“平均”的规则融合。融合字典的计算式表示如下：

$\begin{array}{l} {\mathit{\boldsymbol{D}}_{\rm{F}}}\left( n \right) = \\ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{D}}}^{\rm{C}}}\left( n \right),} & {\begin{array}{*{20}{c}} {} & {} \end{array}\begin{array}{*{20}{c}} {} & {} \end{array}\begin{array}{*{20}{c}} {} & {} \end{array}{L^{\rm{C}}}\left( n \right) > {L^{\rm{R}}}\left( n \right)且} \end{array}\\ \quad \quad \left| {{L^{\rm{C}}}\left( n \right) - {L^{\rm{R}}}\left( n \right)} \right|/\left| {{L^{\rm{C}}}\left( n \right) + {L^{\rm{R}}}\left( n \right)} \right| \ge \lambda \\ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{D}}}^R}\left( n \right)} & {\begin{array}{*{20}{c}} {} & {} \end{array}\begin{array}{*{20}{c}} {} & {} \end{array}\begin{array}{*{20}{c}} {} & {} \end{array}{L^{\rm{C}}}\left( n \right) \le {L^{\rm{R}}}\left( n \right)且} \end{array}\\ \quad \quad \left| {{L^C}\left( n \right) - {L^{\rm{R}}}\left( n \right)} \right|/\left| {{L^{\rm{C}}}\left( n \right) + {L^{\rm{R}}}\left( n \right)} \right| \ge \lambda \\ \begin{array}{*{20}{c}} {\left( {{{\mathit{\boldsymbol{D}}}^{\rm{C}}}\left( n \right) + {{\mathit{\boldsymbol{D}}}^{\rm{R}}}\left( n \right)} \right)/2,\;\;} & {其他} \end{array} \end{array} \right. \end{array}$

(12)

式中：设λ=0.25，D_F∈R^d×N表示融合字典。根据医学图像的物理特性，使用信息熵^[22]作为特征指标。这种方法将稀疏域和空间域的方法结合起来，考虑了医学图像的物理特性计算字典原子的特征指标，与稀疏域的方法相比，具有更加明确的物理意义。

3 脑部CT/MR图像融合过程

不失一般性，设已经配准的脑部CT/MR源图像I_C, I_R∈R_MN。图 4为脑部CT/MR图像融合流程，融合过程如下：

1) 预处理阶段。使用步长为1的滑动窗把源图像I_C、I_R分为$\sqrt m \times \sqrt m $大小的图像块，然后将这些提取的图像块编纂成m维列向量。对于每幅源图像I_C、I_R，都有$\left( {M - \sqrt m + 1} \right)\left( {N - \sqrt m + 1} \right)$个图像块。将两幅源图像中的第j个图像块拉成的列向量记为和x_C^j和x_R^j，并减去各自的均值：

${\mathit{\boldsymbol{\hat x}}}_C^j = {\mathit{\boldsymbol{x}}}_C^j - m_C^j \cdot 1$

(13)

${\mathit{\boldsymbol{\hat x}}}_R^j = {\mathit{\boldsymbol{x}}}_R^j - m_R^j \cdot 1$

(14)

其中：m_C^j和m_R^j分别是x_C^j和x_R^j中所有元素的均值；1表示一个全1的m维列向量。

2) 融合阶段。使用CoefROMP算法^[17]求解$\left\{ {{\mathit{\boldsymbol{\hat x}}}_{\rm{C}}^j, {\mathit{\boldsymbol{\hat x}}}_{\rm{R}}^j} \right\}$的稀疏表示系数，计算式如下：

$\begin{array}{l} \quad {\mathit{\boldsymbol{\alpha }}}_{\rm{C}}^j = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_{\mathit{\boldsymbol{\alpha }}} {\left\| {\mathit{\boldsymbol{\alpha }}} \right\|_0}\;\;\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;{\left\| {{\mathit{\boldsymbol{\hat x}}}_{\rm{C}}^j - {{\mathit{\boldsymbol{D}}}_{\rm{F}}}{\mathit{\boldsymbol{\alpha }}}} \right\|_2} ＜ \varepsilon \end{array}$

(15)

$\begin{array}{l} \quad {\mathit{\boldsymbol{\alpha }}}_{\rm{R}}^j = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_{\mathit{\boldsymbol{\alpha }}} {\left\| {\mathit{\boldsymbol{\alpha }}} \right\|_0}\;\;\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;{\left\| {{\mathit{\boldsymbol{\hat x}}}_{\rm{R}}^j - {{\mathit{\boldsymbol{D}}}_{\rm{F}}}{\mathit{\boldsymbol{\alpha }}}} \right\|_2} ＜ \varepsilon \end{array}$

(16)

其中：D_F是由2.2节提出的融合字典。

将稀疏系数的l₂范数作为源图像的活跃度测量，则稀疏系数α_C^j和α_R^j通过以下的融合规则融合：

${\mathit{\boldsymbol{\alpha }}}_{\rm{F}}^j = \left\{ \begin{array}{l} {\mathit{\boldsymbol{\alpha }}}_{\rm{C}}^j & {\left\| {{\mathit{\boldsymbol{\alpha }}}_{\rm{C}}^j} \right\|_2} \ge {\left\| {{\mathit{\boldsymbol{\alpha }}}_{\rm{R}}^j} \right\|_2}, \\ {\mathit{\boldsymbol{\alpha }}}_{\rm{R}}^j, & 其它 \end{array} \right.$

(17)

均值m_C^j和m_R^j使用“加权平均”规则融合：

$m_{\rm{F}}^j = w \cdot m_{\rm{C}}^j + \left( {1 - w} \right) \cdot m_{\rm{R}}^j$

(18)

其中，w=1/(1+exp(-β(m_C^j-m_R^j))), β>0。则x_C^j和x_R^j的融合结果为：

$\mathit{\boldsymbol{x}}_{\rm{F}}^j = {{\mathit{\boldsymbol{D}}}_{\rm{F}}}{\mathit{\boldsymbol{\alpha }}}_{\rm{F}}^j + m_{\rm{F}}^j \cdot 1$

(19)

3) 重建阶段。对所有的图像块都执行上述两个步骤以得到所有图像块的融合结果。对于每个块向量x_F^j，通过反滑窗的过程重塑成$\sqrt m \times \sqrt m $的图像块并放回到对应的像素位置，再对重复像素取平均得到最终的融合图像I_F。

4 实验结果与分析

本文的实验环境为：32位Windows7操作系统、Matlab R2013a、Intel i3-2350 2.3 GHz处理器、4 GB运行内存。为验证本文方法的有效性，选取三组已经配准的脑部CT/ MR图像进行融合，分别为正常脑部CT/MR (图 5(a)、(b))、脑萎缩CT/MR (图 6(a)、(b))、脑肿瘤CT/MR (图 7(a)、(b))，图片大小均为256×256。选取的对比算法有：离散小波变换(DWT)^[3]、平稳小波变换(SWT)^[4]、非下采样Contourlet变换(NSCT)^[7]、传统基于选择最大的稀疏表示方法(Sparse Representation “choose-Max”-based, SRM)^[11]、基于K-SVD的稀疏表示方法(Sparse Representation based on K-SVD, SRK)^[12]、基于多尺度字典学习的方法(Multi-scale Dictionary Learning, MDL)^[13]。

基于多尺度变换的方法中，对于DWT和SWT方法，分解水平都设为3, 小波基分别设为“db6”和“bior1.1”。NSCT方法使用“9—7”金字塔滤波器和“c—d”方向滤波器，分解水平设为{2², 2², 2³, 2⁴}。基于稀疏表示的方法中滑动步长为1，图像块大小均为8×8，字典大小均为64×256，误差ε=0.01，稀疏度τ=6。ICDL方法使用改进的K-SVD算法，执行6个字典更新周期(Dictionary Update Cycles, DUC)和30次迭代。

由图 5~7可以看出，DWT方法的融合图像边缘纹理模糊，图像信息失真且存在块效应；与DWT方法相比，SWT和NSCT方法的融合质量相对较好，图像的亮度、对比度、清晰度有了很大的提升，但仍存在边缘亮度失真，软组织和病灶区域存在伪影的问题；SRM和SRK方法较基于多尺度变换的方法，图像的骨组织和软组织更加清晰，伪影也有所减少，能很好地识别病灶区域；MDL方法与SRM和SRK方法相比，能保留更多的细节信息，图像质量取得进一步的改善，但仍有部分伪影存在；本文提出的ICDL方法在图像的亮度、对比度、清晰度和细节的保持度上都优于其他方法，融合图像没有伪影，骨组织、软组织和病灶区域显示清晰，有助于医生诊断。

为进一步验证本文方法的有效性，使用互信息(Mutual Information, MI)^[23]、基于梯度Q^AB/F[24]、基于相位一致Q_p^[25]和基于通用图像质量指标Q_w^[26]来对融合图像进行客观评价。MI反映融合图像从源图像提取的信息量的多少，提取的信息量越多则融合图像的质量越好；Q^AB/F反映融合图像对源图像边缘特性的保留情况，取值范围在[0, 1]，越接近1表示融合图像的边缘越清晰；Q_p反映融合图像多种类型的显著特征，如边缘和角点等；Q_w反映融合图像在系数相关性、光照和对比度方面与源图像的关联，符合人类视觉系统的特点。

表 1~3为三组实验条件下不同融合方法的客观评价指标。

从表 1~3可以看出，DWT方法在MI、Q^AB/F、Q_p和Q_w四项指标上均不理想，主要是因为离散小波变换具有平移变化性、走样和缺乏方向性等缺陷；SWT和NSCT方法在上述四项指标上均明显优于DWT方法是由于它们具有平移不变性；SRM、SRK和MDL方法在MI指标上优于三层分解水平的多尺度变换的方法，而在其他三项指标上存在一些低于多尺度变换方法的情况，这可能是由于字典的表示能力不足以提取足够的细节信息和l₁范数取极大的融合规则导致的图像的灰度不连续效应；经计算，本文提出的ICDL方法的MI、Q^AB/F、Q_p和Q_w四项指标在上述三组实验条件下的均值分别为4.113 3、0.713 1、0.463 6和0.762 5。与其他几种方法相比：本文提出的ICDL方法在MI指标上的值最大，表明该方法的融合图像从源图像提取的信息最丰富；同时在Q^AB/F、Q_p和Q_w三项指标上对应的值也最大，表明ICDL方法对图像误配准处理得较好，图像融合质量较高，其主客观评价一致。

为验证本文提出的字典学习方法的时间效率，分别计算不同的图像融合方法的字典学习阶段在相同条件下进行10次实验所消耗时间的平均值。SRK、MDL、CDL和ICDL方法的字典训练时间均值分别为4.61 min、42.59 min、81.27 min和5.96 min。SRK方法的字典学习中训练数据是从50幅自然图像中提取5万个8×8图像块；MDL方法的字典学习中训练数据是在小波域中分别对每个子带提取5万个8×8图像块；传统的耦合字典学习(Coupled Dictionary Learning, CDL)方法和本文提出的ICDL方法的字典学习中训练数据均是从图 2的训练集中提取5万对8×8的图像块。下面对四种字典学习方法的时间复杂度进行分析。假设训练图像的尺寸为N×N，提取N²个大小为$\sqrt n \times \sqrt n $的训练图像块。SRK方法的时间复杂度为O (N²)；MDL方法的时间复杂度为O ((3S+1) N²)，其中S为小波变换的分解水平；CDL方法的时间复杂度为O (N⁴)，这是因为两个字典交替地更新影响了字典学习的效率；ICDL方法的时间复杂度为O (N²)，这是由于字典更新阶段同时更新字典和稀疏表示系数的非零元素而使零元素保持完备，降低了矩阵的维数，同时使用联合学习的策略，进一步提高了字典学习的时间效率。

5 结语

针对目前基于字典学习的多模态脑部医学图像融合中，使用单字典难以得到精确的稀疏表示进而影响脑部医学图像融合的质量，以及字典训练时间过长的问题，提出了一种改进耦合字典学习的脑部CT/MR图像融合方法。该方法分别对正常脑部、脑萎缩和脑肿瘤三组脑部医学图像进行了多次实验，结果表明本文提出的ICDL方法与基于多尺度变换的方法、传统基于选择最大的稀疏表示方法、基于K-SVD的稀疏表示方法以及多尺度字典学习的方法相比，不仅提高了脑部医学图像融合的质量，而且有效降低了字典训练的时间，能为临床医疗诊断提供有效帮助。但仍存在以下需要进一步研究的问题：

1) 本文提出的融合方法在空间域中使用滑窗技术易使源图像的细节被模糊。因此，进一步研究适用于稀疏表示方法的新的图像融合规则是下一步的工作重点。

2) 本文是在已配准的脑部医学图像的基础上进行的融合研究，对未精确配准的脑部医学图像进行有针对性的图像融合的方法也需要进一步的研究。

3) 本文主要实现了脑部CT和MR医学图像的高质量融合，但如何对脑部解剖图像和功能图像如MR与正电子发射断层成像(Positron Emission Tomography, PET)、MR与单光子发射计算机断层成像(Single Photon Emission Computed Tomography, SPECT)等进行高质量融合，也将是下一步继续研究的内容。