基于剪切波变换的医学图像融合算法

阿都建华 王邦平 王 珂 王 艳

(成都信息工程学院软件工程学院，成都610225)

引言

医学影像技术已经是临床诊断中不可替代的技术手段，在临床诊断中得到了广泛地应用，并取得了巨大的成效。由于不同的成像设备采用的传感器原理不同，医学影像技术有多种成像模式，分别产生不同的医学图像，这些图像反映人体脏器信息也有所不同。医学图像主要包括形态图像、结构图像和功能图像。医学显微图像属于形态图像，X线图像、超声、CT、MRI 等属于结构图 像，PET、SPECT、FMRT 及EIT 等属于功能图像［1］。在临床上，将不同模态的医学图像进行适当融合，使多种模态的信息有机结合，可获得信息更为丰富的复合图像，为临床诊断和治疗提供更为可靠、准确的依据［2-8］。

近年来，由于基于多尺度几何分析的图像融合方法具有多分辨率、稀疏描述等特性，被广泛应用于图像处理领域。小波变换是最典型的多尺度分解工具，然而小波变换不能很好地表示二维图像中的线奇异性，容易导致伪吉布斯现象。剪切波是一种克服了小波变换缺点的新颖多尺度几何分析工具，由Guo 及其合作者于2007 年通过特殊形式的具有合成膨胀的仿射系统构造的一种接近最优的多维函数稀疏表示方法［9］。剪切波［9-11］具有简单的数学结构，可以通过对一个函数进行伸缩、平移、旋转而生成一组基函数，使其和多分辨分析关联起来。与图像融合中常用的其他多尺度分解工具相比，剪切波变换具有自身的优势。剪切波变换不仅具有与曲线波［12］相同的接近最优的非线性误差逼近阶，而且在频率空间中剪切波是逐层细分的，因而具有更好的表示性能。同时，剪切波在剪切过程中和轮廓波［13-14］不同，它没有方向数目和支撑基尺寸大小的限制，而且剪切波逆变换只需对正向变换中的剪切滤波器进行加和处理，不需要像轮廓波变换那样对方向滤波器逆合成，因此剪切波的实现过程比轮廓波具有更高的计算效率。

Miao 等将剪切波变换应用于图像融合，并证明了它的融合结果比其他多尺度方法捕获更多的图像细节信息［15］。郑红等将剪切波变换应用于红外与可见光图像融合，并针对红外与可见光图像的特性，采用区域显著性方法融合分解系数［16］。笔者在深入研究剪切波变换的基础上，结合医学图像的特性，提出了一种基于剪切波变换的医学图像融合算法。该算法首先通过剪切波变换，对原始图像进行多尺度、多方向的分解，得到原始图像的高频子带与低频子带，然后对低频子带的系数采用非负矩阵分解［17］(NMF)融合算法进行融合，对高频子带系数则设计了一种基于人类视觉系统特性的视觉能量对比度方法进行选择，最后通过重构融合后的高、低频子带系数获得最终的融合图像。实验证明，本研究提出的融合算法所获取的融合图像能融合不同模态医学图像的主要信息，提高了图像的临床诊断和治疗价值，取得了良好的效果。

1 基于剪切波变换的图像融合方法

1.1 剪切波变换

在论文［10］中，剪切波可定义为

令ψ ∈L2(R2)，需满足以下条件:

1)对于ξ = (ξ1，ξ2)∈(ξ2≠0)，有(ε)=(ε1，ε2)=(ε1(ε2/ε1)，其中(ε)为ψ(ε)的傅里叶变换;

2)ψ1为连续小波，∈C∞(R)，suppψ1⊂［-2，- 1/2］∪［1/2，2］;

假定

对每一个j ≥0，存在

图1 剪切波频域剖分图以及剖分子区域的几何特征。(a)剪切波的频域支撑;(b)剪切波频域剖分图Fig． 1 The tiling of the frequency by the shearlets and the geometric characteristics of the split region． (a)The tiling of the frequency by the shearlets;(b)The size of the frequency support of a shearlet

集合:

是L2(D0)∨= {f ∈L2(R2):supp ^f ⊂D0}的一个Parseval 框架。

1.2 基于Shearlet 变换的融合规则

首先对经过严格配准的原始图像进行剪切波变换，将原图分解成由剪切波系数表示的低频子带Sc=0(i，j)与一系列高频子带(i，j)，其中c = 1，2，…，N，表示图像的第c 个高频子带，k = 1，2，…，M，(i，j)为分解方向参数，表示剪切波系数所在方向的位置，从而实现图像的多尺度，多方向分解;然后将低频子带系数和高频子带系数融合，再进行逆运算，得到最后的融合图像。剪切波的具体分解步骤如下［16，18］:

1)通过非降采样塔式分解，将原始图像分解为各个尺度图像T0，T1，…，TN，其中T0表示分解后的粗尺度图像，T1～TN表示分解后的第1 细分尺度图像至第N 细分尺度图像;

2)通过快速二维傅里叶变换，将各细分尺度图像T1，…，TN变换至频域FT1，…，FTN;

3)将各细分尺度频域图像FT1，…，FTN作为剪切波滤波器组的输入，通过的多方向分解得到各高频子带频域系数，通过对的快速二维傅里叶反变换得到一系列高频子带系数(i，j);

4)基于步骤2 的结果T0，直接使S0(i，j)= T0，综合步骤4 的结果(i，j)，即得到剪切波变换的全部系数;

5)将分解后源图像对应的低频子带系数和高频子带系数按特定的融合算法进行融合处理，分别得到融合图像的低频子带系数和高频子带系数;

6)通过剪切波逆变换，得到融合图像。

图2 融合算法的流程Fig．2 The flow diagram of fusion algorithm

2 基于剪切波变换的医学图像融合算法

2.1 融合算法流程

首先对经过严格配准的原始图像进行剪切波变换，将原图分解成由剪切波系数表示的低频子带Sl0(x，y)与一系列高频方向子带Si，l(x，y)，其中i= 1，2，…，N 为图像的第i 个高频方向子带，l = 1，2，…，M 为分解方向参数，(x，y)表示剪切波系数所在方向的位置，从而实现了原始图像的多尺度、多方向分解。本研究采用NMF 算法完成低频子带系数的融合，得到融合后的低频子带系数;对于高频方向子带系数融合则采用视觉能量对比度方法完成，得到融合后的高频方向子带系数;最后通过Shearlet 逆运算，得到融合图像。本研究提出的融合算法流程如图2 所示。

2.2 低频子带系数融合规则

2.2.1 非负矩阵分解理论

非负矩阵分解(non-negative matrix factorization，NMF)是一种新的矩阵分解方法［17］，由Lee 和Seung 于1999 年在Nature 上的一篇论文中提出。在矩阵分解的过程中，该方法始终约束所有的元素，将基作为非负数存在，即要求所有分量始终为纯加性的描述，同时使分解后的所有分量也均为非负值，并且降低了矩阵的维度。非负矩阵分解(NMF)问题可描述为:已知一个非负矩阵V，求非负的n × r 矩阵W 和非负的r × m 矩阵H，使

对于 NMF 问题的求解，常用欧氏距离和Kullback-liebler 散度作为目标函数，有

因为对于一个给定矩阵V，矩阵W 和H 的最佳选择是要使得V 和WH 之间的重构误差最小。因此，NMF 的求解问题实际上是个优化问题，可以描述为

上述公式可以应用梯度下降法来求解得到结果，求解过程也是收敛的，Lee 和Seung 已经在理论上对此做出了证明［19］，并得到求解矩阵W 和矩阵H的迭代规则，即

式中，“.* ”和“. /”分别表示矩阵中各元素的相乘和相除。

2.2.2 基于非负矩阵分解的融合算法

在传感器成像的过程中，由于传感器自身的因素或外界的各种影响，常常会引入各种噪声，因此待融合的原始观测图像实际是就是客观真实世界在成像过程中引入了这些噪声而形成的。在非负矩阵分解算法中，可以假设V = WH + ε(ε 表示噪声)，此时噪声ε 会在式(12)所示的迭代算法中趋于收敛，这个过程恰好符合图像融合的过程［20-21］。因此，联系图像融合过程，如果假设观测图像为V，真实图像为W，噪声为ε，那么V 可以理解为W 与ε 之和，这样NMF 可以有效地应用于图像融合。

由NMF 算法理论可知，通过迭代运算方法能够针对原始数据矩阵V 得到一个基于部分的近似表示形式WH。其中，W 的列数(即特征基的数量r)是一个待定量，将直接决定特征子空间的维数。对于特定的数据集，隐藏在数据集内部特征空间的维数是确定的。也就是说，当选取的r 与实际数据集特征空间的维数一致时，所得到的特征空间以及特征空间的基最有意义。当r =1 时，通过迭代算法将得到唯一一个含有源数据全部特征的特征基。

由上述内容可知，NMF 与图像融合能很好地结合在一起应用。假设有k 幅来自于多传感器的大小为m × n 的观测图像f1，f2，…，fk，那么可以将观测图像逐个元素地按照行优先的方式存储到一个列向量中，这样就可以得到一个mn × k 的矩阵V。V中包含k 个列向量v1，v2，…，vk，每个列向量代表k幅观测图像中的一幅图像的信息，有

对这个观测矩阵V 进行非负矩阵分解，分解时取r =1，则可得到一个唯一的特征基W。显然，此时的W 包含了参与融合的k 幅图像的完整特征，将特征基W 还原到源图像的像素级上，即可得到比源图像效果都好的图像。

2.3 高频子带系数融合规则

2.3.1 人类视觉系统

图像对比度是指对一幅图像中灰度值反差大小的测度，反差越大意味着对比度越大，反差越小意味着对比度越小。人眼能够分辨的亮度差异所要求的最小光亮度差值ΔL 称为亮度辨别阈值，对比灵敏度是指人眼具有的分辨亮度差异的能力。亮度辨别阈值ΔL 也不是固定不变的，其大小在不同的背景亮度L 下并不相同。也就是说，即使客观的亮度是增强的，但在亮度L 的变化没有增加到L+ ΔL 之前，人眼是感受不到这些亮度变化的，只有当亮度L 变化达到或超过L + ΔL 时，人眼才能感受到这些亮度变化。因此，在研究中，采用ΔL/L 这一相对灵敏度阈值才更为合理。实验表明，在无背景光时，ΔL/L 的值近似为一个常数，其数值为0.02。

在对视觉系统的进一步研究中发现，对比度敏感门限与背景亮度的关系更接近指数规律，同时得到了一些更准确的对比度敏感函数(contrast sensitivity function)［22-25］。其中，根据文献［24］，对比度敏感门限可表示为

式中:I0为当I = 0 (背景亮度为零)时的对比度敏感门限;α 为幂函数的指数，是常数，可根据视觉生理实验获得，其值为0.6 ～0.7。

2.3.2 高频子带系数融合规则

人类视觉系统对图像局部对比度的变化十分敏感，因为对比度反映了图像内纹理、边缘等信息的变化特征，同时包含了图像的高频信息以及相对于局部图像背景的强度。鉴于此，融合方法可有选择地突出被融合图像的对比度信息，以求达到良好的视觉效果。在高频子带系数的融合算法中，笔者利用人类视觉特性来提取图像的细节信息，并根据高频子带系数的特点设计了一种视觉能量对比度方法，进行融合图像系数的选择，有

式中:j 表示第j 个高频子带，l 表示第l 个方向;Rj，l(x，y)表示像素点(x，y)的灰度与局部背景平均灰度的局部对比度;Ej，l(x，y)是以像素点(x，y)为中心的局部区域能量和;¯mj，l(x，y)为低频子带系数以(x，y)为中心的某一局部区域均值;α 为视觉常数，通常取0.6 ～0.7 间的值。

局部对比度Rj，l(x，y)表示高频子带图像中的系数与通过上节方法融合后获得的低频子带图像的局部区域灰度平均值的比值，即局部对比度。由前述人类视觉系统对比度敏感度的内容可知，其值的大小在一定程度上能反映出待融合高频子带图像的像素点融合于图像后的效果。Rj，l(x，y)的计算方式为

式中:Sj，l(x，y)为j 尺度、l 方向的 高频子带系数;为以像素点(x，y)为中心、尺寸大小为m × n 的低频子带图像的局部灰度均值，其中SL(x +i，y +j)为融合后的低频子带系数。

式中，Sj，l(x + i，y + j)为j 尺度、l 方向的高频子带系数。

局部区域能量和Ej，l(x，y)表示高频子带图像的局部区域能量和，分别计算图像f1和f2的局部区域能量值(m，n)和(m，n)，有

对于经过严格配准的源图像f1和f2，则在经过shearlet 分解后，先分别计算得到它们的视觉能量对比度，然后再选择高频子带系数，有

在高频子带系数的选择中，将选取视觉能量对比度较大者作为最后的融合图像的低频子带系数。

3 结果

为了验证上述图像融合方法的有效性，选取两组经过严格配准的医学图像进行融合实验。为了能充分比较本研究提出的融合算法，实验1 和实验2 中都采用了3 种图像融合方法作为对比方法。，包括拉普拉斯金字塔(LP)图像融合方法、基于小波变换(DWT)的图像融合方法和基于简单非下采样Contourlet 变换(NSCT)的图像融合方法。在实验中，拉普拉斯金字塔分解的层次为4 层;基于小波变换的图像融合方法和基于NSCT 的简单融合方法都采用简单的平均值方法选取低频子带系数，采用最大绝对值方法选取高频子带系数，NSCT 的多尺度分解层次为3，方向分解级数分别为4、3、2。对于所提出的算法，剪切波变换分解层数为3，2 个细分尺度分解方向数分别为6，区域大小尺寸取3 × 3 。在实验中，还采用了信息熵(information entropy，IE)、空间频率(spatial frequency，SF)、标准差(standard deviation，SD)、平均梯度(average grads，AG)和均方根交叉熵(root mean square cross entropy，RCE)作为融合图像的客观评价指标。

实验1 采用严格配准后的MRI 图像和CT 图像作为实验的源图像，如图3 所示。其中，图3(a)是CT 图像，图3(b)是MRI 图像，图像的尺寸大小均为256 像素×256 像素，且已经过严格配准。图3(c)～图3(f)分别为基于DWT 方法获取的融合图像、基于LP 方法获取的融合图像、基于NSCT 简单方法获取的融合图像和基于本研究提出方法获取的融合图像。从主观视觉角度可以明显看出，图3(c)的边缘位置明显出现了伪吉布斯现象，图3(e)则非常暗淡，MRI 图像信息明显融入较少。图3(e)相对图3(f)亮度较暗，可视性效果略差。表1 是实验1 的客观评价指标，可以看出本研究提出的方法在IE、SF、AG 三项指标中均取得了最大值，在RCE指标中取得了最小值。IE 最大，说明融合图像包含的信息量最多，融合效果最好;SF 最大，说明图像在空间域的变化最活跃，融合图像的质量最好;AG 最大，说明图像中的细节最丰富，图像清晰程度最高;RCE 最小，说明融合图像和原始图像之间的差异水平最小，即融合效果最好。SD 指标的值小于DWT方法，由于SD 表示像素点的灰度值与融合图像平均值的偏离程度，其值越大说明融合图像中细节变化越丰富，融合质量就越好，因此仅该项客观指标不是最优值。由上述分析可知，图3(f)在这4 项客观指标中取得了最好的融合结果，仅SD 一项指标的计算结果不是最优。因此，从客观评价指标同样可以看出，本研究提出方法的获取的融合图像效果最好。

图3 MRI 图像与CT 图像的融合结果比较。(a)CT 图 像;(b)MRI 图像;(c)DWT 融合图像;(d)LP 融合图像;(e)NSCT 融合图像;(f)本方法融合图像Fig． 3 Fusion results on MRI image and CT image． (a)CT image;(b)MRI image;(c)Fusion result of the DWT-based method;(d)Fusion result of the LP-based method;(e)Fusion result of the NSCT-based method;(f)Fusion result of the proposed method

表1 实验1 融合图像客观评价指标比较Tab． 1 The objective evaluation on the experiment 1

实验2 采用Brain Web［26］中提供的不同形态大脑图像作为实验图像，融合结果如图4 所示。其中，图4(a)是T1 图像，图4(b)是PD 图像，图像的尺寸大小均为216 像素×181 像素，且已经过严格配准。图4(c)～图4(f)分别为基于DWT 方法获取的融合图像、基于LP 方法获取的融合图像、基于NSCT简单方法获取的融合图像和基于本研究提出方法获取的融合图像。从主观视觉角度可以明显看出，图4(c)～图4(e)的图像效果明显较差，图像中边缘信息比较模糊，整体对比度差。图4(f)的融合质量相对较好，融合信息丰富，整体对比度高，显得更为清晰。因此，在实验2 中所提出的融合方法取得了更好的主观视觉融合效果。表2 是实验2 的客观评价指标，从中可以看出，本方法同样在IE、SD、SF

及RCE 指标中取得了最优值，仅AG 一项指标的计算结果不是最优。由实验1 中对各项客观评价指标的分析可知，图4(f)在大部分客观指标中取得了最好的融合结果。因此，从客观评价指标同样可以看出，本方法获取的融合图像效果最好。

图4 T1 图像与PD 图像的融合结果比较。(a)T1 图像;(b)PD 图像;(c)DWT 融合图像;(d)LP 融合图像;(e)NSCT 融合图像;(f)本方法融合图像Fig．4 Fusion results on T1 image and PD image．(a)T1 image;(b)PD image;(c)Fusion result of the DWT-based method;(d)Fusion result of the LPbased method;(e)Fusion result of the NSCT-based method;(f)Fusion result of the proposed method

表2 实验2 融合图像客观评价指标比较Tab． 2 the objective evaluation on the experiment 2

4 讨论和结论

由于图像传感器种类繁多，应用环境和目的各不相同，所以图像融合算法也具有一定的针对性。笔者主要针对不同模态医学图像进行融合研究，使多种模态的信息能更好地有机结合，形成一幅新的信息丰富、易于判别的融合图像，以提高临床诊断和治疗的准确性。多尺度几何分析工具在图像融合中已经广泛应用，其有效性也已经在实践中得到证明。不同的多尺度几何分析工具具有不同的特点，相应的图像融合效果也各不相同。本研究采用最新的多尺度几何分析工具——剪切波变换作为图像分解变换的数学工具，有效克服了小波变换的缺陷。此外，在分解后得到的高频和低频子带图像融合中，采用的融合规则对融合图像的质量至关重要，也是基于多尺度几何分析图像融合中的核心研究内容。

本研究以一种新颖的多尺度几何分析工具——剪切波变换为基础，提出一种新的医学图像融合算法。本算法利用剪切波变换的多分辨率和多方向性特点，将待融合的原始图像进行分解，得到低频子带图像和高频方带图像。在低频子带系数的融合中，采用非负矩阵分解方法，再利用人类视觉系统对局部对比度较为敏感的特点，提出最大视觉能量对比度方法，进行高频方向子带系数的融合，最后通过剪切波逆变换得到融合图像。实验证明，本研究提出的融合方法比传统的小波方法、拉普拉斯金字塔方法、NSCT 方法，在主观融合质量和客观评价指标方面都有较大提高。

［1］ 陈晓艳，李健楠，王化祥. 一种电阻抗图像与CT 图像融合方法研究［J］.中国生物医学工程学报，2012，30(6):892 -896.

［2］ Constantinos S，Pattichis MS，Micheli-Tzanakou E. Medical imaging fusion applications:An overview［C］ //Conference Record of the Thirty-Fifth Asilomar Conference on Signals，Systems and Computers． Pacific Grove:IEEE，2001:1263 -1267.

［3］ Das S，Kundu MK. NSCT-based multimodal medical image fusion using pulse-coupled neural network and modified spatial frequency ［J］. Medical and Biological Engineering and Computing，2012，50(10):1105 -1114.

［4］ Li Shutao， Yin Haitao， Fang Leyuan. Group-sparse representation with dictionary learning for medical image denoising and fusion［J］. IEEE Transactions on Bio-medical Engineering，2012，59(12):3450 - 3459.

［5］ Singh R，Srivastava R，Prakash O，et al． Multimodal medical image fusion in dual tree complex wavelet transform domain using maximum and average fusion rules［J］. Journal of Medical Imaging and Health Informatics，2012，2(2):168 -173.

［6］ Townsend DW，Beyer T. A combined PET/CT scanner:the path to true image fusion［J］. British journal of Radiology，2002，75(9):S24 -S30.

［7］ Yang Liu，Guo Baolong，Ni Wei. Multimodality medical image fusion based on multiscale geometric analysis of contourlet transform［J］.Neurocomputing，2008，72(1):203 -211.

［8］ 杨立才，刘延梅，刘欣，等. 基于小波包变换的医学图像融合方法［J］. 中国生物医学工程学报，2009，28(1):12 -16.

［9］ Guo Kanghui，Labate D. Optimally sparse multidimensional representation using shearlets［J］. SIAM journal on mathematical analysis，2008，39(1):298 -318.

［10］ Kutyniok G，Labate D. Resolution of the wavefront set using continuous shearlets［J］. Trans Amer Math Soc，2009，361(5):2719 -2754.

［11］ Easley G， Labate D， Lim WQ. Sparse directional image representations using the discrete shearlet transform［J］. Applied and Computational Harmonic Analysis，2008，25(1):25 -46.

［12］ Candes EJ，Donoho DL. Curvelets:A surprisingly effective nonadaptive representation for objects with edges［C］//in Curve and Surface Fitting:Saint-Malo，Nashville，USA:Vanderbilt University Press，2000:105 -120.

［13］ Do MN，Vetterli M. The contourlet transform:an efficient directional multiresolution image representation ［J］. IEEE Transactions on Image Processing，2005 ，14 (12):2091 -2106.

［14］ Cunha AL，Zhou J，Do MN. The nonsubsampled contourlet transform: theory， design， and applications ［J］. IEEE Transactions on Image Processing，2006 ，15 (10):3089 -3101.

［15］ Miao Qiguang，Shi Cheng，Xua Pengfei，et al． A novel algorithm of imagefusion using shearlets［J］. Optics Communications，2011，284(6):1540 -1547.

［16］ 郑红，郑晨，闫秀生，陈海霞. 基于剪切波变换的可见光与红外图像融合算法［J］. 仪器仪表学报，2012，33(7):1613 -1619.

［17］ Lee DD，Seung HS. Learning the parts of objects by non-negative matrix factorization［J］. Nature，1999，401(6755):788 -791.

［18］ Lim WQ. The discrete shearlet transform:A new directional transform and compactly supported shearlet frames［J］. IEEE Transactions on Image Processing，2010 ，19(5):1166 -1180.

［19］ Seung D，Lee L. Algorithms for non-negative matrix factorization［J］. Advances in Neural Information Processing Systems，2001，13(2):556 -562.

［20］ 苗启广，王宝树. 基于非负矩阵分解的多聚焦图像融合研究［J］. 光学学报，2005，25(6):755 -759.

［21］ 苗启广，王宝树. 图像融合的非负矩阵分解算法［J］. 计算机辅助设计与图形学学报，2005，17(7):2029 -2032.

［22］ Legge GE，Foley JM. Contrast masking in human vision［J］.JOSA，1980，70(12):1458 -1471.

［23］ Watson AB. DCT quantization matrices visually optimized for individual images［C］// Human Vision，Visual Processing and Digital Display IV. San Jose:SPIE，1993:202 -216.

［24］ Watson AB. Efficiency of a model human image code［J］. JOSA A，1987，4(12):2401 -2417.

［25］ Legge GE. A power law for contrast discrimination［J］. Vision Research，1981，21(4):457 -467.

［26］ http://www.bic.mni.mcgill.ca/brainweb/.