融合NSST和稀疏表示的PET和MRI图像融合

康家银，陆 武，张文娟

1(江苏海洋大学 电子工程学院，江苏 连云港 222005) 2(连云港市第一人民医院 核医学科，江苏 连云港 222061) 3(江苏海洋大学 计算机工程学院，江苏 连云港 222005)

1 引 言

正电子发射断层(Positron Emission Tomography，PET)是一种核医学成像技术，也是当今最先进的无创伤性诊断手段之一[1].磁共振成像(Magnetic Resonance Imaging，MRI)是一种能够呈现人体解剖信息的结构型成像[2].PET图像具有能够反映人体不同组织的新陈代谢情况的颜色信息，但其空间分辨率低，结构信息不明显；MRI图像的空间分辨率高、软组织结构信息强，但缺乏反应人体新陈代谢功能的颜色信息[3,4].通过影像融合技术可以将PET图像和MRI图像各自的优势有效结合，从而使得融合后的图像在保留PET图像颜色信息的同时，又具有MRI图像明显的空间结构信息，以便融合后的图像将更有助于临床中对疾病的诊断和治疗[5].

近年来，针对PET和MRI图像的融合，研究者们已经提出了不同的策略和算法.其中一类是基于变换和替换策略的IHS(Intensity-Hue-Saturation)方法[5,6].这类方法首先将PET变换到IHS空间中；然后用匹配后的MRI(将MRI与PET在IHS空间中的I分量进行直方图匹配)去替换PET在IHS空间中的I分量；最后经反变换得到融合后的PET图像.这类方法通常会得到一幅空间信息强(分辨率高)的PET图像，但会使得在信息替代的过程中损失PET图像的颜色信息.另一类用于PET和MRI图像融合的方法是基于多分辨率分析(Multi-Resolution Analysis,MRA)的方法，如小波变换等[7,8].这类方法首先将PET和MRI图像进多尺度分解(分解成系数和变换基)；然后将分解后的系数根据某种融合规则进行融合；最后将融合后的系数和变换基进行反变换，从而得到最终的融合图像.这类方法可以有效的保留PET图像的色彩信息，但在提升PET图像的空间结构信息方面具有不可避免的局限性.第三类用于PET和MRI图像融合的方法是基于稀疏表示(Sparse Representation，SR)的方法[9,10].这类方法首先求出PET和MRI图像各自的稀疏表示系数；然后将二者的稀疏系数采用某种融合规则进行融合；再利用融合后的稀疏系数和用于重构的过完备字典重构出最终的融合图像.稀疏表示在图像融合上取得了引人注目的效果，但目前用于PET和MRI图像融合的稀疏表示模型中,用于图像重构的字典的训练/构造是针对整幅图像进行的，即从整幅图像中提取图像块来构造全局字典，从而导致计算复杂度的增加.此外，由于构造全局字典时没有考虑图像中块与块之间的结构相似性，致使求解得到的稀疏系数不够准确，从而难以准确的重建出目标图像[11].

非下采样剪切波变换(Non-Subsampled Shearlet Transform，NSST)[12]，一种先进的多分辨率、多尺度分析方法.除了具有其它一些多分辨率分析方法具有的局部性、多尺度和多方向性特点外，NSST还具有平移不变性、计算复杂性低等特点.因此，NSST被成功的应用于图像融合中[13-15].

根据以上分析，本文结合NSST和稀疏表示各自的优点，提出一种基于NSST和稀疏表示的PET和MRI图像融合方法和框架.首先，利用NSST将PET和MRI图像分解成高频和低频成分.然后，使用梯度值较大法将对应的高频部分进行融合；基于SR将对应的低频部分进行融合，即利用SR获取能表征低频成分的稀疏系数，再对获取的稀疏系数进行融合，并获得融合的低频成分.最后，利用NSST逆变换得到PET与MRI的融合图像.该方法不仅能够有效的保持PET图像中的颜色信息，而且还能够很好的集成MRI图像中丰富的结构信息.

2 相关工作

2.1非下采样剪切波变换

不同于之前提出的多分辨率分析(MRA)方法(如轮廓波变换[16]、曲线波变换[17]等)，非下采样剪切波变换(NSST)是一种新颖的多分辨率、多尺度分析方法.NSST除具有其它MRA方法的特点，如局部性、多尺度和各向异性等外，还具有计算复杂度更低、实现效率更高、分解方向数不受限制等优点[13].

图1 NSST分解示意图Fig.1 Diagram of NSST decomposition

2.2 图像的稀疏表示理论

近年来，稀疏表示(SR)成为研究的热点并已成功的应用到图像处理和计算机视觉的各个方面，如超分辨率重建[18]，图像重构[19]，多源/多传感器/多聚焦图像的融合[20-22]等.图像的稀疏表示模型是基于一种假设，即一幅图像可以被事先定义好的过完备字典中的“原子”稀疏性的表示[23].稀疏表示使得图像的信息主要集中在少数“原子”上，而这些“原子”反映了图像的主要特征和内部结构.图像的稀疏表示模型可表示为：

(1)

其中x∈n是大小为的图像块(以列向量形式表示)；α∈N是稀疏表示系数；‖·‖0是l0范数，表示向量α中非零元素的个数；D∈n×N(n≪N)是过完备字典；ε是预先设定的逼近误差容限，且ε≥0.

式(1)是一个非凸函数，所以上述稀疏表示最优化问题，即(1)的求解是一个NP-难问题[24].因此需要采用特殊的方法对(1)进行求解.其中一类常用的方法是用l1范数代替l0范数，从而将(1)中的非凸优化问题转化为l1范数的凸优化问题，于是便可得到如下的基追踪(Basis Pursuit，BP)优化问题：

(2)

3 本文算法

3.1 算法框架

基于前面引入的NSST和SR基本原理，本文提出的PET和MRI图像融合基本框架如图2所示(假设源PET和MRI图像已经进行了配准)，其实现步骤为：

2)将PET和MRI高频子带采用梯度值最大规则进行融合，具体的融合规则将在随后的小节中详述.

3)对于二者的低频子带利用SR进行融合.类似的，融合规则的具体实现方法将在后面的小节中给出.

4)最后，将整合后的高、低频子带进行NSST反变换，从而得到最终的融合图像.

图2 PET和MRI图像融合框架Fig.2 Image fusion framework for PET and MRI images

3.2 基于梯度值法的高频子带融合

图像经过NSST分解后，边缘、轮廓等许多源图像的细节信息包含在高频子带中.高频子带融合的目的就是从源图像中尽可能多的保留这些有用的细节信息.传统的基于“绝对值取大”的高频子带融合规则仅仅考虑了子带图像中“像素”自身，而没有考虑该像素与其周围像素间的空间关系.因此简单的采用传统的基于“绝对值取大”的高频子带融合规则容易导致源图像中的细节信息的丢失和引入额外的人为 “噪声”.为此，本文采用如下的基于“梯度值取大”的高频子带融合规则：

(3)

其中IPH(i,j)，IMH(i,j)和IFH(i,j)分别为PET，MRI和融合的高频子带图像中(i,j)处的值；grad(·)为求PET和MRI高频子带中(i,j)处梯度值的操作，本文采用Sobel算子求梯度值，其形式如图3所示.

图3 Soble梯度算子Fig.3 Sobel gradient operator

3.3 基于SR的低频成分融合

图像经过NSST分解后，源图像的信息主要集中在低频子带图像中.因此，低频子带的融合非常关键.为此，本文基于稀疏表示进行PET和MRI低频子带图像的融合.

3.3.1 自适应局部字典的构造

在图像的稀疏表示模型中，字典(无论是用于稀疏系数求解的编码字典还是用于图像表示的重构字典)的构造至关重要.经典的稀疏表示中，字典的构造或学习都是基于整幅图像的，即构造或学习得到的是针对整幅图像的全局字典，从而难以体现图像中的局部相似性.此外，全局字典的引入也势必会增加计算复杂度.因此，本文采用如下的策略构造能够反应图像中块与块之间结构相似性的局部自适应字典.

图4 局部字典构造示意图Fig.4 Diagram of local dictionary construction

3.3.2 低频子带系数的融合

本文提出的基于稀疏表示的PET和MRI图像低频子带融合框架如图5所示，分为四个阶段(见图5中的①，②，③，④)，具体的实现过程如下：

第一阶段：联合字典的构建

第二阶段：稀疏编码

第三阶段：稀疏系数的融合

图5 基于SR的低频子带融合示意图Fig.5 Diagram of low frequency component fusion based on SR

第四阶段：稀疏重构

重复上述1)到7)步，直到将低频子带(IPL和IML)中所有的“像素点”(共N对)处理完毕.

4 实验结果与分析

为了验证本文所提出的用于PET和MRI图像融合的算法框架的性能，我们用一组公开的数据集进行了相关的实验.

4.1 实验数据介绍

本文用于实验的PET和MRI图像来自于哈佛大学医学院的“全脑”(The Whole Brain)数据集1.本文选用了30组PET和MRI图像进行了融合实验，图像大小均为256×256.

4.2 算法参数设置

本文实验中算法的参数设置如下：

NSST分解级数：L=3；

搜索窗口滑动步长：6；

“最”相关的图像块个数：s=60.

4.3实验结果分析

本文选取了30组数据进行实验，并将本文方法得到的结果与其它方法的结果进行了比较.其中用于比较的方法分别来自于文献[5](记为IHS-Retina)、文献[6](记为IHS-Gabor)和文献[10](记为CT-SR).图6和图7为两组数据(分别以Coronal和Transaxial两个视角呈现)的融合结果.

图6 PET和MRI图像的融合结果(Coronal视角)

图7 PET和MRI图像的融合结果(Transaxial视角)

由图6和图7可知，就保留MRI图像的结构信息而言，算法NSST(融合结果分别见图6(e)和图7(e))和本文算法(融合结果分别见图6(h)和图7(h))明显好于其它几种方法，见图6和图7的各个子图中用矩形框 "标注"的区域；而在其它方法中，IHS-Retina方法的效果最不理想，即模糊化程度最严重(见图6(c)中的椭圆形区域和图7(c)中的圆形区域).此外，相较于图6(e)和图7(e)，图6(h)和图7(h)的结构信息更接近原始的MRI图像，见图6和图7的各个子图中用椭圆形框 "标注"的区域.

1http://www.med.harvard.edu/aanlib/

另外，由图6和图7可知，在保留PET图像颜色信息方面，本文方法略逊于IHS-Retina方法，与其它方法相当.此外，将各融合结果与原始MRI图像(图6(b)、图7(b))的脑壳(Skull)部分进行对照，可见唯有本文方法得到的融合结果与原始MRI图像最接近(灰度一致)，而其它方法得到的结果中均有不同程度的伪彩色.

综合以上两方面的因素可知，就视觉效果而言，相对于其它几种方法，本文方法取得了更好的融合效果，即在较好的保持PET图像颜色信息的同时，很好的融合了MRI图像的结构信息.

为了进一步评价融合效果，本文采用一些客观评价指标对融合结果进行客观、定量的评价和分析.本文用于评价融合效果的指标包括信息熵(Information Entropy，IE)[25]、平均梯度(Average Gradient，AG)[25]、光谱差异性(Spectral Discrepancy、SD)[5]和Qabf[26].各个评价指标的具体形式如下：

信息熵IE主要用于衡量融合后图像中所含信息量的大小.IE越大，表明图像中包含的信息越丰富.本文中，IE越大，表明融合后，PET图像包含的信息(色彩信息和结构信息)越丰富.IE定义如下：

(4)

其中P(i)表示灰度值为i(i=0,1,…,L-1)的像素出现的概率.

平均梯度AG反应的是融合图像的清晰度(或空间分辨率)的高低.AG越大，表明图像层次越多，越清晰.本文中，AG越大，表明融合后，PET图像层次越丰富，空间分辨率越高.AG定义如下：

k=R,G,B

(5)

光谱差异性SD衡量的是融合图像的光谱(颜色)特性；具体的，SD反应了融合图像与源PET图像的颜色差异性.SD越小，表明融合图像的光谱特性越好.本文中，SD越小，表明融合后，PET图像的色彩信息和原始PET图像越接近.SD定义如下：

(6)

其中Fk(x,y)和Ok(x,y)分别为融合后的图像与输入的PET图像在(x,y)处的值；R，G，B为PET图像的三个颜色分量.

Qabf主要用于评价融合过程中源图像的边缘信息保持程度.Qabf值越大，表明融合图像与源图像间的边缘相似度越高，或源图像的边缘信息损失越小.本文中，Qabf值越大，主要表明融合后的PET图像中包含了更多的MRI图像的边缘信息(结构信息).Qabf定义如下：

(7)

利用不同方法进行PET和MRI图像融合后得到的量化指标如表1所示(其中表1中的结果是30组数据的平均值).

表1 不同方法融合PET和MRI图像的结果
Table 1 Fusion results of PET and MRI images using different methods

融合方法IEAGSDQabfIHS-Retina5.02210.1175.9070.565IHS-Gabor5.59710.9136.4310.635NSST5.66311.0856.2740.647SR5.69711.5066.1370.612CT-SR5.74110.4027.5510.563本文方法5.83612.0516.0290.681

注：加粗者为最优结果.

由表1可知，就信息熵(IE)、平均梯度(AG)和Qabf而言，本文方法的性能均优于其它方法.就光谱差异性(SD)而言，本文方法的性能略逊色于IHS-Retina，但优于其它四种方法.综上可知，就定量指标而言，相对于其它几种方法，本文提出的方法取得了更好的融合效果.

4.4 算法参数L对融合性能的影响

本文提出的算法框架中涉及到一个重要的参数，即NSST的分解层级L.本小节结合4.3节中给出的四个客观量化指标(IE、AG、SD和Qabf)分析不同的L对融合性能的影响.图8中，(a)、(b)、(c)和(d)分别为不同分解层级L对应的IE、AG、SD和Qabf.

图8 分解层级L对融合性能的影响

由图8(a)和图8(d)可知，随着分解层级的加大，融合性能会随之提升，但当L>3时，性能提升便不明显.同样的，由图8(b)和图8(c)可知，随着分解层级的加大，融合性能会随之提升，但是提升不明显.此外，分解层级越大，运算耗时便越长.因此，综合考虑上述各个因素，最终选择了L=3的分解层级.

5 结 论

本文提出了一种用于PET和MRI图像融合的方法，即利用非下采样剪切波变换(NSST)和稀疏表示(SR)的PET和MRI图像融合方法.该方法充分的结合了NSST和SR各自的优点.利用公开的脑部数据集进行实验的结果表明，本文提出的方法取得了很好的融合效果，即在有效保留PET图像“颜色”信息的同时，很好的保持了MRI图像的空间结构信息.此外，与其它已有的用于PET和MRI图像融合的方法相比，本文的方法无论在视觉效果还是在客观的定量评价中，都取得了最好的图像融合效果.