CN109816617B - 基于导向滤波和图论显著性的多模态医学图像融合方法 - Google Patents

Abstract

本发明请求保护一种基于导向滤波和图论显著性的多模态医学图像融合方法，主要用于融合多模态医学图像，包括步骤：S1得到待融合图像的平滑子带图像和细节子带图像；S2利用图论显著性方法构造融合规则重组平滑子带系数和细节子带系数；S3重构融合后的平滑子带系数和细节子带系数得到融合后的图像。本发明基于导向滤波算法的分解方法时间复杂度较低，保证融合图像且在分解过程中能有效的减少噪声的干扰，增加方法的鲁棒性；同时利用图论的显著性方法来提取平滑子带图像的亮度、颜色和方向特征，能保留灰度图像的亮度信息和伪彩色图像的颜色信息，有效解决其他融合方法融合图像亮度信息丢失和色彩失真问题，保证融合图像的质量。

Description

基于导向滤波和图论显著性的多模态医学图像融合方法

技术领域

本发明属于图像处理技术领域，用于融合多个模态的医学图像，主要是灰色图像(MRI)和伪彩色图像(PET)的融合，即MRI-PET融合和灰度图像(MRI)和灰度图像(CT)的融合即MRI-CT的融合，特别涉及利用导向滤波算法和图论显著性算法。

背景技术

目前，图像融合是图像处理研究领域的一个研究热点，主要的研究方法可概括为两类：基于空域的图像融合方法和基于频域的图像融合方法。

基于空域的图像融合方法是直接对原图像构造融合规则来得到融合图像，常用的基于空域的图像融合方法包括简单的加权平均，基于PCA变换的融合方法和基于IHS变换的融合融合方法。该类方法的优点在于操作简单，时间复杂度低，但会丢失大量的边缘轮廓等重要的细节信息。基于频域的图像融合方法主要包括三大类：基于塔型变换的图像融合方法、基于小波变换的图像融合方法和基于几何变换的图像融合方法。与基于空域的融合方法相比，基于频域的图像融合方法能在多个尺度上提取显著性特征，使得融合后的图像能尽可能多的保留原图像的信息。基于塔型变换的图像融合方法主要将原始图像分解为不同尺度的高频和低频图像，优点在于能提取更多的图像信息，但高频图像不具备方向性；基于小波变换的图像融合方法能在水平、垂直、对角三个方向上进行分解，但不能很好的表示边角信息；基于轮廓波变换的图像融合方法能在多个方向上分解原图像，但不具备平移不变性，会产生伪影。为了解决多尺度分解方法产生伪影和时间复杂度高的问题，李树涛等人提出了基于边缘保持滤波器分解的图像融合方法，具有时间复杂度低，能很好的保留原图像的边缘轮廓信息的优点。

目前图像融合方法众多，虽都各具优点，但也存在很多问题，总结起来主要主要包括(1)时间复杂度高；(2)亮度信息丢失；(3)色彩失真；(4)边缘信息丢失；(5)抗干扰能力差。导向滤波算法具有平滑去噪的能力，图论显著性方法能提取图像的亮度、颜色、方向信息，使得融合后的图像能保留原图像的亮度信息和颜色信息，并具有一定的抗噪能力，因此将两种算法结合起来用于图像融合，能取得很好的融合结果。

发明内容

本发明旨在解决以上现有技术亮度信息丢失，色彩失真和抗噪能力差的问题。提出了一种有效平滑噪声，并将图像分解为多个尺度，能在多个尺度上提取特征的基于导向滤波和图论显著性的多模态医学图像融合方法，主要用于融合灰度图像(MRI)和伪彩色图像(PET)，即一种基于导向滤波和图论显著性的多模态医学图像融合方法。本发明的技术方案如下：

一种基于导向滤波和图论显著性的多模态医学图像融合方法，其包括以下步骤：

S1、首先利用导向滤波算法对原始输入图像A和B进行L次不同的平滑操作，平滑半径为r＝(2¹,2²,...2^L)，得到一系列平滑子带图像C_i和一系列细节子带图像D_i，D_i＝I-C_i,其中i表示第i个尺度。

S2、其次，利用图论显著性构造融合规则，主要包含三个步骤：活跃度测量步骤、系数分组步骤及系数重组步骤；

S3、最后利用求和方法对子带图像系数进行重构，得到融合后的图像F，F＝F^C+F^D，F^C表示融合后的平滑子带图像F^D表示融合后的细节子带图像。

进一步地，所述步骤S1的导向滤波算法假设引导图像I和输出图像Q之间是一个线性模型:

式(1)中，ω_k表示以i为中心，大小为(2r+1)×(2r+1)的窗口，(a_k,b_k)是ω_k中的两个系数。我们把输入图像P作为引导图像I，要求输出图像Q尽可能多的保留原图像的信息，即输出图像Q和输入图像P的误差尽可能的小，误差表达如下:

式(2)中，∈是一个正则化系数，防止求得的a_k太大。通过最小二乘法求得：

式(3)和式(4)中，μ_k和

分别表示ω_k中I的均值和方差，|ω|表示ω_k中的像素总和，/>

表示P在ω_k中的均值。我们利用导向滤波算法进行图像分解时，利用滑动窗口对原始输入图像A和B进行L次不同的平滑操作，平滑半径为r＝(2¹,2²,2³)，得到一系列平滑子带图像C_i。再通过D_i＝I-C_i得到一系列细节子带图像D_i，其中i表示第i个尺度，i＝1,2,3。由于在平滑滤波时(a_k,b_k)会被多个窗口ω_k包含，所以我们最终求得输出图像即平滑子带图像为：

式(5)中，

|ω|表示ω_k中的像素总和，ω_i表示以i为中心，大小为(2r+1)×(2r+1)的窗口。

进一步地，所述步骤S2的活跃度测量步骤：由于A和B的平滑子带图像系数包含原图像的大量结构、亮度和颜色信息，利用图论显著性方法来提取这些子带图像的亮度、颜色和方向显著性特征，且每幅子带图像对应生成一幅显著图

其中/>

和/>

表示A和B的平滑子带图像的显著图，GBVS()表示图论显著性方法获取显著图，/>

和/>

表示A和B的平滑子带图像系数。对于A和B的细节子带图像系数，我们利用信息熵函数/>

来衡量每幅子带图像所携带的信息大小，其中p_I(i)表示输入图像的像素值，生成显著图/>

EN表示信息熵函数，/>

和/>

表示A和B的细节子带系数的显著图，/>

和/>

表示A和B的细节子带图像系数。

进一步地，所述步骤S2的活跃度测量步骤中利用图论显著性方法(GBVS)获取显著图的具体步骤为：

(1)特征提取和特征图生成：利用高斯采样方法在图像的r、g、b三个通道上构建图像的颜色、亮度和方向的金字塔，然后用高斯金字塔和中心-周围差值算法计算图像的亮度、颜色和方向特征，并将中心像素的尺度设为c∈{2,3,4}，周围像素的尺度设为s＝c+δ,δ∈{3,4}，以便于能在多个尺度上提取特征，

表示点到点的相减。得到亮度特征为：

式(6)中，I＝(r+g+b)/3,r、g、b表示图像的红、绿、蓝三个通道系数，I(c)和I(s)表示在c和s尺度上求解亮度特征I，I(c,s)表示在c,s尺度上的亮度特征。

通过模拟皮层中的“颜色双竞争”系统表示：在接受野中心，神经元被一种颜色激活而禁止另外一种颜色，而在周围区域则是相反的。在主视觉皮层中，红/绿，绿/红，蓝/黄，黄/蓝颜色对都存在这种空间和色彩竞争关系。因此，产生RG(c,s)和BY(c,s)特征图来反映双竞争关系。RG(c,s)和BY(c,s)表达式如下：

式(7)和式(8)中，R＝r-(g+b)/2，G＝g-(r+b)/2，B＝b-(r+g)/2，Y＝(r+g)/2-|r-g|/2-b，r、g、b表示图像的红、绿、蓝三个通道系数，R(c)、G(c)、B(c)、Y(c)表示在尺度c上求解R、G、B、Y值，R(s)，G(s)，B(s)，Y(s)表示在尺度s上求解R、G、B、Y值。

采用方向Gabor金字塔得到局部方向信息，得到方向特征图:

O(c,θ)表示在尺度c和角度为θ时的方向特征，O(s,θ)表示在尺度s和角度为θ时的方向特征。

(2)特征显著性计算和特征显著图生成：研究学者认为一个点与其相邻的点的差异越大，则其显著性越大，所以基于图论显著性方法(GBVS)在得到了多个亮度、颜色和方向特征图后，利用图论的思想和马尔科夫链算法对这些特征图进行显著性计算和特征图组合生成显著图。假设对于特征图M中的点p和点q，两点的不相似性为：

式(9)中,M(p)和M(q)分别表示p和q的特征值。

假设特征图M中的所有像素点构成一个全连接图G_A，图中任意一个顶点与其他n-1个顶点相连，G_A中点p到点q的有向边的权值为：

式(10)中

(p_x,p_y),(q_x,q_y)表示点p和点q的空间坐标，δ表示随机参数。

在G_A中定义马尔科夫链，并将边上的权重归一化到[0,1],使得对于任意一个顶点，作为始点的所有边的权重之和为1。图论显著性方法将G_A中的节点类比为马尔科夫链的状态，有向边的权重类比为转移概率，马尔科夫链模型为：

式(11)中a_i和a_j表示状态，p_ji表示状态转移概率，n表示n个状态，最后能得到稳定状态的马尔科夫链表示特征显著性状态，得到特征显著图S(a₁,a₂,…,a_n)。

进一步地，所述步骤S2中的系数分组步骤中，由于显著图的值越大代表所携带的重要信息越多，通过分别比较同一幅原图像的平滑子带和细节子带图的显著图的值的大小来生成决策图D，且决策图的值由0或1组成，构建方式如下：

式(12)中

和/>

分别代表图像A的平滑子带图像和细节子带图像的决策图，/>

和/>

表示A图像的平滑子带图像和细节子带图像在第i个尺度的显著图，n表示n个尺度，本文中n＝3，式(13)中/>

和/>

分别代表图像B的平滑子带图像和细节子带图像的决策图，

和/>

表示B图像的平滑子带图像和细节子带图像在第i个尺度的显著图，n表示n个尺度，本文中n＝3。

进一步地，所述步骤S2中的系数重组步骤主要利用加权平均的方法对子带系数进行重组，加权平均方法的思想是把融合系数看做是输入图像与权重图的线性乘积之和，把系数分组步骤中得到的决策图D作为权重图，最后分别对原图像A和B的多个尺度的平滑子带系数和细节子带系数进行重组得到A和B融合的平滑图像系数和细节图像系数的方法如下：

式(14)和式(15)中，n表示n个尺度，本文中n＝3，

和/>

分别代表图像A的融合后的平滑子带图像系数和细节子带图像系数，/>

和/>

表示A图像在第i个尺度上的平滑子带系数和细节子带系数，/>

和/>

表示B图像在第i个尺度上的平滑子带系数和细节子带系数,/>

表示在图像A的平滑子带图像和细节子带图像的决策图，/>

和/>

表示在图像B的平滑子带图像和细节子带图像的决策图，/>

的/>

和分别代表图像B的融合后的平滑子带图像系数和细节子带图像系数；

最后利用平均加权的方法将A和B的融合后的平滑子带系数和细节子带系数进行融合得到最终融合图像F的平滑子带图像和细节子带图像系数，表达式如下：

其中F^C和F^D分别表示融合后的平滑子带图像系数和细节子带图像系数。

本发明的优点及有益效果如下：

由于基于多尺度分解的融合方法能在多个尺度上提取特征进行融合，且人眼对图像的显著性特征更为敏感，本发明提出一种基于导向滤波算法和图论显著性算法的多模态图像融合方法，利用导向滤波算法将原图像分解为一系列平滑子带图像和细节子带图像，有利于我们在多个尺度上提取特征，能保证融合图像尽可能多的保留原图像的重要信息，并且与其他多尺度分解方法相比，基于导向滤波算法的分解方法时间复杂度较低，保证融合图像且在分解过程中能有效的减少噪声的干扰，增加方法的鲁棒性；同时利用图论的显著性方法来提取平滑子带图像的亮度、颜色和方向特征，能保留灰度图像的亮度信息和伪彩色图像的颜色信息，有效解决其他融合方法融合图像亮度信息丢失和色彩失真问题，保证融合图像的质量。

附图说明

图1是本发明提供优选实施例方法的流程图

图2为输入图像A(MRI)示例

图3为输入图像B(CT)示例

图4为融合后的图像F示例

图5为本专利多模态医学图像的融合实验结果和其他6种方法实验结果

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、详细地描述。所描述的实施例仅仅是本发明的一部分实施例。

本发明解决上述技术问题的技术方案是：

本发明的方法流程如图1所示。

本发明涉及了输入原图像示例如图2、图3，融合后图像的示例如图4，本专利图像的融合实验结果和其他6种方法实验结果图示例为图5。

具体步骤：

S1、首先利用导向滤波算法对原始输入图像A和B进行L次不同的平滑操作，平滑半径为r＝(2¹,2²,...2^l)，得到一系列平滑子带图像C_i和一系列细节子带图像D_i(D_i＝I-C_i)，本文令L＝3；

S2、其次，利用图论显著性构造融合规则，主要包含三个步骤：a)活跃度测量、b)系数分组、c)系数重组。

活跃度测量步骤：由于A和B的平滑子带图像包含原图像的大量结构、亮度和颜色信息，利用图论显著性方法来提取这些子带图像的亮度、颜色和方向显著性特征，且每幅子带图像对应生成一幅显著图

其中/>

和/>

表示A和B的平滑子带图像的显著图，/>

和/>

表示A和B的平滑子带图像系数。对于A和B的细节子带图像，利用信息熵函数/>

来衡量每幅子带图像所携带的信息大小，其中p_I(i)表示输入图像的像素值生成显著图S，

系数分组：由于显著图的值越大代表所携带的重要信息越多，通过分别比较同一幅原图像的平滑子带和细节子带图的显著图的值的大小来生成决策图D，且决策图的值由0或1组成，构建方式如下：

式(12)中,

和/>

分别代表图像A的平滑子带图像和细节子带图像的决策图，

和/>

表示A图像的平滑子带图像和细节子带图像在第i个尺度的显著图，n表示n个尺度，本文中n＝3，式(13)中,/>

和/>

分别代表图像B的平滑子带图像和细节子带图像的决策图，/>

和/>

表示B图像的平滑子带图像和细节子带图像在第i个尺度的显著图，n表示n个尺度，本文中n＝3。

系数重组：利用加权平均的方法对子带系数进行重组，加权平均方法的思想是把融合系数看做是输入图像与权重图的线性乘积之和，本方法把系数重组步骤中得到的决策图D作为权重图，最后分别对原图像A和B的多个尺度的平滑子带系数和细节子带系数进行重组得到A和B融合的平滑图像系数和细节图像系数的方法如下：

式(14)中的

和/>

分别代表灰度图像A的融合后的平滑子带图像系数和细节子带图像系数，式(15)中/>

的/>

和分别代表伪彩色图像B的融合后的平滑子带图像系数和细节子带图像系数。

最后利用平均加权的方法将A和B的融合后的平滑子带系数和细节子带系数进行融合得到最终融合图像F的平滑子带图像和细节子带图像系数，表达式如下：

S3、最后利用求和方法对子带图像系数进行重构，得到融合后的图像F(F＝F^C+F^D)。

为了评估本算法的性能，本发明采用哈佛大学公开数据集中的两种病例作为实验图像集，在本实验中主要采用6种基于频域的多尺度分解方法进行对比，分别是基于小波变换的图像融合方法(DWT),基于梯度金字塔变换的图像融合(GRP),基于导向滤波的图像融合方法(GFF),基于非下采样轮廓波变换的图像融合方法(NSCT),基于稀疏表示的图像融合方法(SR)和基于并行显著性的图像融合方法。

本实验选择了哈佛大学公开数据集的四类病例的40张图像和一种加噪图像的10张图像来进行灰度图像(MRI)和伪彩色图像(PET)即MRI-PET的融合和灰度图像(MRI)和灰度图像(CT)即MRI-CT的融合，对图像的融合质量进行定量评价，主要采用SSIM、MI、IWSSIM、QAC、VSI、GSM作为客观评价指标。实验结果对比表如表1，表2和表3所示。

表1 MRI-PET图像融合客观评价指标

表2 MRI-CT图像融合客观评价指标

/>

表3加噪MRI-PET图像融合客观评价指标

以上这些实施例应理解为仅用于说明本发明而不用于限制本发明的保护范围。在阅读了本发明的记载的内容之后，技术人员可以对本发明作各种改动或修改，这些等效变化和修饰同样落入本发明权利要求所限定的范围。

Claims (2)

1.一种基于导向滤波和图论显著性的多模态医学图像融合方法，其特征在于，包括以下步骤：

S1、首先利用导向滤波算法对原始输入图像A和B进行L次不同的平滑操作，平滑半径为r＝(2¹,2²,...2^L)，得到一系列平滑子带图像C_i和一系列细节子带图像D_i，D_i＝I-C_i,其中i表示第i个尺度；I表示亮度特征；

S2、其次，利用图论显著性构造融合规则，主要包含三个步骤：活跃度测量步骤、系数分组步骤及系数重组步骤；

S3、最后利用求和方法对子带图像系数进行重构，得到融合后的图像F，F＝F^C+F^D，F^C表示融合后的平滑子带图像，F^D表示融合后的细节子带图像；

所述步骤S2的活跃度测量步骤具体包括：

由于A和B的平滑子带图像系数包含原图像的大量结构、亮度和颜色信息，利用图论显著性方法来提取这些子带图像的亮度、颜色和方向显著性特征，且每幅子带图像对应生成一幅显著图

其中/>

和/>

表示A和B的平滑子带图像的显著图，GBVS()表示图论显著性方法获取显著图，/>

和/>

表示A和B的平滑子带图像系数，对于A和B的细节子带图像系数，利用信息熵函数

来衡量每幅子带图像所携带的信息大小，其中p_I(i)表示输入图像的像素值，生成显著图/>

EN表示信息熵函数，/>

和/>

表示A和B的细节子带系数的显著图，/>

和/>

表示A和B的细节子带图像系数；

所述步骤S2的活跃度测量步骤中利用图论显著性方法(GBVS)获取显著图的具体步骤为：

(1)特征提取和特征图生成：利用高斯采样方法在图像的r、g、b三个通道上构建图像的颜色、亮度和方向的金字塔，然后用高斯金字塔和中心-周围差值算法计算图像的亮度、颜色和方向特征，并将中心像素的尺度设为c∈{2,3,4}，周围像素的尺度设为s＝c+δ，δ∈{3，4}，以便于能在多个尺度上提取特征，

表示点到点的相减，得到亮度特征为：

式中，I＝(r+g+b)/3，r、g、b表示图像的红、绿、蓝三个通道系数，I(c)和I(s)表示在c和s尺度上求解亮度特征I，I(c，s)表示在c，s尺度上的亮度特征；

通过模拟皮层中的“颜色双竞争”系统表示：在接受野中心，神经元被一种颜色激活而禁止另外一种颜色，而在周围区域则是相反的，在主视觉皮层中，红/绿，绿/红，蓝/黄，黄/蓝颜色对都存在这种空间和色彩竞争关系，因此，产生RG(c，s)和BY(c，s)特征图来反映双竞争关系，RG(c，s)和BY(c，s)表达式如下：

/>

式中R＝r-(g+b)/2，G＝g-(r+b)/2，B＝b-(r+g)/2，Y＝(r+g)/2-|r-g|/2-b，r、g、b表示图像的红、绿、蓝三个通道系数，R(c)、G(c)、B(c)、Y(c)表示在尺度c上求解R、G、B、Y值，R(s)，G(s)，B(s)，Y(s)表示在尺度s上求解R、G、B、Y值；

采用方向Gabor金字塔得到局部方向信息，得到方向特征图：

O(c，θ)表示在尺度c和角度为θ时的方向特征，O(s，θ)表示在尺度s和角度为θ时的方向特征；

(2)特征显著性计算和特征显著图生成：研究学者认为一个点与其相邻的点的差异越大，则其显著性越大，所以基于图论显著性方法GBVS在得到了多个亮度、颜色和方向特征图后，利用图论的思想和马尔科夫链算法对这些特征图进行显著性计算和特征图组合生成显著图，假设对于特征图M中的点p和点q，两点的不相似性为：

式中M(p)和M(q)分别表示点p和q的特征值；

假设特征图M中的所有像素点构成一个全连接图G_A，图中任意一个顶点与其他n-1个顶点相连，G_A中点p到点q的有向边的权值为：

式中

表示点p和点q的空间坐标，σ表示随机参数；

在G_A中定义马尔科夫链，并将边上的权重归一化到[0,1],使得对于任意一个顶点，作为始点的所有边的权重之和为1；图论显著性方法将G_A中的节点类比为马尔科夫链的状态，有向边的权重类比为转移概率，马尔科夫链模型为：

式中a_i和a_j表示状态，p_ji表示状态转移概率，n表示n个状态，最后能得到稳定状态的马尔科夫链表示特征显著性状态，得到特征显著图S(a₁,a₂,...,a_n)；

所述步骤S2中的系数分组步骤中，由于显著图的值越大代表所携带的重要信息越多，通过分别比较同一幅原图像的平滑子带和细节子带图的显著图的值的大小来生成决策图D，且决策图的值由0或1组成，构建方式如下：

式中

和/>

分别代表图像A的平滑子带图像和细节子带图像的决策图，/>

和/>

表示A图像的平滑子带图像和细节子带图像在第i个尺度的显著图，n表示n个尺度，n＝3，式中

和/>

分别代表图像B的平滑子带图像和细节子带图像的决策图，/>

和/>

表示B图像的平滑子带图像和细节子带图像在第i个尺度的显著图，n表示n个尺度；/>

所述步骤S2中的系数重组步骤主要利用加权平均的方法对子带系数进行重组，加权平均方法的思想是把融合系数看做是输入图像与权重图的线性乘积之和，把系数分组步骤中得到的决策图D作为权重图，最后分别对原图像A和B的多个尺度的平滑子带系数和细节子带系数进行重组得到A和B融合的平滑图像系数和细节图像系数的方法如下：

式中，n表示n个尺度，

和/>

分别代表图像A的融合后的平滑子带图像系数和细节子带图像系数，/>

和/>

表示A图像在第i个尺度上的平滑子带系数和细节子带系数，/>

和

表示B图像在第i个尺度上的平滑子带系数和细节子带系数,/>

表示在图像A的平滑子带图像和细节子带图像的决策图，/>

和/>

表示在图像B的平滑子带图像和细节子带图像的决策图，/>

的/>

和分别代表图像B的融合后的平滑子带图像系数和细节子带图像系数；

最后利用平均加权的方法将A和B的融合后的平滑子带系数和细节子带系数进行融合得到最终融合图像F的平滑子带图像和细节子带图像系数，表达式如下：

其中F^C和F^D分别表示融合后的平滑子带图像系数和细节子带图像系数。

2.根据权利要求1所述的基于导向滤波和图论显著性的多模态医学图像融合方法，其特征在于，所述步骤S1的导向滤波算法假设引导图像I和输出图像Q之间是一个线性模型:

式中ω_k表示以i为中心，大小为(2r+1)×(2r+1)的窗口，(a_k,b_k)是ω_k中的两个系数，把输入图像P作为引导图像I，要求输出图像Q尽可能多的保留原图像的信息，即输出图像Q和输入图像P的误差尽可能的小，误差表达如下:

式中∈是一个正则化系数，防止求得的a_k太大，通过最小二乘法求得：

式中，μ_k和

分别表示ω_k中I的均值和方差，|ω|表示ω_k中的像素总和，ω_k表示以i为中心，/>

表示P在ω_k中的均值；/>

利用导向滤波算法进行图像分解时，利用滑动窗口对原始输入图像A和B进行L次不同的平滑操作，平滑半径为r＝(2¹,2²,2³)，得到一系列平滑子带图像C_i，再通过D_i＝I-C_i得到一系列细节子带图像D_i，其中i表示第i个尺度，i＝1,2,3，由于在平滑滤波时(a_k,b_k)会被多个窗口ω_k包含，最终求得输出图像即平滑子带图像为：

式中，

|ω|表示ω_k中的像素总和ω_k表示以i为中心，大小为(2r+1)×(2r+1)的窗口。/>

Images