CN106504221A - 基于四元数小波变换上下文结构的医学图像融合新方法 - Google Patents

Abstract

本发明公开了一种基于四元数小波变换QWT(Quaternion Wavelet Transform)上下文结构的医学图像融合方法，主要解决医学图像融合时亮度、对比度等的变化以及原始图像信息丢失的问题。其实现步骤是：1)分别对待融合图像进行QWT处理得到待融合图像对应的QWT系数；2)对待融合图像的QWT系数的低频部分根据绝对值取大的融合规则进行融合，得到低频的融合系数；3)对待融合图像的QWT系数的高频部分求取高频子带系数对应的清晰度系数；4)基于上下文统计模型与清晰度系数设计融合规则，获得融合后的高频系数；5)根据得到的低频和高频的融合系数，利用QWT逆变换得到融合图像。本发明能充分整合医学图像的结构和功能信息，有效保护图像细节，改善视觉效果，相比传统的融合方法极大的提高了融合图像的质量。

Description

基于四元数小波变换上下文结构的医学图像融合新方法

技术领域

本发明涉及一种基于四元数小波变换上下文结构的医学图像融合新方法，是医学图像处理技术领域的一项融合方法，在临床医学诊断和治疗中有着广泛的应用。

背景技术

作为当今国内外研究的热点，医学图像融合技术有着很强的实用性。医学成像有很多，主要包括超声、X射线、磁共振MRI、计算机断层成像CT、正电子发射断层成像PET、单光子发射断层成像SPECT、红外线等。根据各医学图像体现的内容，可将这些医学图像大致分为两类：功能图像SPECT、PET等和解剖结构图CT、MRI、B超等。这两类图像各有其优缺点：解剖图像主要提供脏器的解剖信息，虽然有分辨率高的优点但是脏器的功能情况却无法反映出来；功能图像分辨率较差，但它能提供解剖图像所不能提供的脏器功能代谢信息。各种医学图像从各种不同的角度反映了人体的信息，单单从某一种特定的图像中是无法得到最全面的人体信息的。只凭借一种图像医生是无法从中精确地判断出患者的病症的，要想精确诊断，中间就必须有综合多种模态的医学图像的过程。

20世纪80年代研究人员逐渐开始关注到医学图像融合技术。医学中往往有多个成像设备，这些成像设备对同一个病灶区域成像之后得到多个医学图像，综合这些医学图像的信息获取互补信息，能为诊断和治疗提供更准确的依据。医学图像融合技术的目的就是结合各种医学图像的优势，减少医学图像的不确定性，让医生对病灶有更直观的观察，进而利于其做出准确的诊断。

上世纪80年代，医学图像融合技术的思想被提出，当时仅采用加权平均等简单算法，效果很不理想。到了90年代，拉普拉斯及高斯等金字塔分析方法、多分辨率分析、小波变换等基于频域的分解工具的提出及应用使得医学图像融合开始成为图像融合领域的前沿课题，并取得了很大进展。但是利用这些方法得到的融合图像多少都会有亮度、对比度等的变化以及原始图像信息的丢失，而四元数小波变换具有丰富的变换信息、近似平移不变性以及时频局部化特性，在处理图像时有着其他工具不具备的优越性，故本发明将四元数小波变换运用到医学图像融合中。

四元数小波变换QWT是一种新兴的有效的图像解析工具，它的独特之处使其在能很好 的用于图像融合领域。在图像融合领域，虽然有很多基于频域的图像融合算法，但大多遵循同一个处理模式。即待融合图像经过多尺度变换得到变换系数的低频子带与高频子带，利用低频融合规则对各个待融合图像的低频子带进行处理得到融合的低频子带，同时利用高频融合规则对各个待融合图像的高频子带进行处理得到融合的高频子带，然后通过多尺度变换对应的逆变换将子带系数从频率域转换回空间域，从而得到融合图像。为了便于描述，此处以两幅待融合图像的融合为例。同时，由于多尺度变换工具的多样性，低频和高频融合规则的多样性，尽管图像融合算法的融合框架类似，但图像融合的算法却有很多种，研究人员往往都是根据自己研究的侧重点选择融合算法。

根据上述分析，本发明提出了一种基于四元数小波变换上下文结构的医学图像融合新方法。该方法通过QWT变换将图像由空间域转换到多尺度域，得到图像的QWT系数。对于QWT低频系数融合，采用对应位置系数绝对值取大的融合规则；对于QWT高频系数融合，进行局部开窗计算清晰度值，然后采用基于上下文的融合规则。

发明内容

本发明的目的是针对上述现有技术的不足，提出一种基于四元数小波变换上下文结构的医学图像融合方法，以解决已有的多尺度医学图像融合方法所得融合图像有亮度、对比度等的变化以及原始图像信息的丢失的问题，并充分整合不同尺度医学图像的结构信息和功能信息，有效保护图像细节，增强图像对比度和边缘轮廓，改善其视觉效果，提高融合图像的质量。

本发明解决其技术问题所采用的技术方案如下：

一种基于四元数小波变换上下文结构的医学图像融合方法，其特征在于包括以下步骤：

1)通过QWT多尺度工具将图像A、B由空间域转换到多尺度域，得到图像A、B的QWT低频部分L和高频部分H；

2)对A、B的低频子带系数采用绝对值取大的融合规则进行融合，得到融合后的低频系数；

3)计算A、B的高频子带系数对应的清晰度系数；

4)基于上下文统计模型与清晰度系数设计融合规则，获得融合后的高频系数；

5)根据得到的低频和高频的融合系数，利用QWT逆变换得到融合图像F。

作为一种优选方案，所述的基于四元数小波变换上下文结构的医学图像融合方法，步骤2)低频部分根据绝对值取大的融合规则进行融合，具体描述见公式：

由于QWT系数的低频系数包括四个子带，故上式中n的取值范围为1,2,3,4。当待融合图像A的子带系数的绝对值不小于待融合图像B的子带系数则融合系数 取A的对应系数；反之，融合系数取B的对应系数。

作为一种优选方案，步骤3)高频子带系数或B,求取其对应的清晰度系数 首先，基于当前高频子带系数设置局部窗体w，然后计算局部窗体的清晰度值CW。CW即为子带系数的清晰度系数其中，j为分解尺度，n为第k,k＝1,2,3个分解方向的第n,n＝1,2,3,4,个子带。将CW执行QWT分解，由变换系数根据四元数代数计算得到各子带(LL、LH、HL、HH子带)对应的幅值与相位信息，其中低频LL子带水平方向、垂直方向、对角方向的相位分别记为PhaseLL_hor、PhaseLL_ver、PhaseLL_dia，高频LH子带水平方向的相位记为PhaseLH_hor，高频HL子带垂直方向的相位记为PhaseHL_ver，高频HH子带对角方向的相位记为PhaseHH_dia；

按下述公式计算CW：

其中：GradL_hor(i′,j′)、GradL_ver(i′,j′)和GradL_dia(i′,j′)分别代表低频幅值系数LAM(i′,j′)的水平、垂直以及对角方向的梯度值。

作为一种优选方案，步骤4)在获得清晰度系数后，采用基于上下文的融合规则实现融合规则。

四元数变换系数的分布不仅受其所在的尺度及方向的影响，还受到该系数在相邻尺度、相邻位置的系数的影响。一般情况下，变换系数的分布有其持续性、相似性与传递性。同一 尺度下，变换系数大的地方相邻的系数一般会比较大，变换系数小的地方相邻的系数就会比较小,变换系数的这一关系我们称之为邻域关系。同样，变换系数大的地方其上一尺度与下一尺度对应位置的系数也会比较大，反之亦然，变换系数的这一关系我们称之为父子关系。四元数小波变换系数的邻域关系和父子关系统称为系数的上下文关系。

QWT系数上下文设计如下式所示：

其中，代表四元数小波变换高频子带系数的活动测度，代表第n个子带的系数均值，此处，本发明采用子带系数对应的清晰度系数作为活动测度。N是当前系数的邻域，ni,nj是邻域的系数下标。是基于上下文的活动测度值，[·]^↓2表示以2为降采样因子的降采样过程，erf()为Sigmoid函数。Sigmoid函数一般用在神经网络作为其激活函数，用以引入非线性关系。erf()定义为：

最后，高频子带的融合系数由下式得到：

从而，完成图像各个高频方向子带系数的融合。

本发明与现有技术相比，其显著优点：(1)利用四元数小波变换的分解优势，能够更为精确的解释图像。(2)本发明方法所提出的基于四元数小波域的上下文融合规则，能够更加有效的鉴别子带系数的特性得到更准确的融合结果，克服了现有技术对四元数小波分解系数相关性考虑不够全面的问题。(3)本发明方法与现有技术相比融合效果更好。

附图说明

图1为本发明基于四元数小波变换上下文结构的图像融合方法的流程图。

图2为本发明计算高频系数清晰度的示意图。

图3为本发明四元数小波变换的系数的结构图。

图4为本发明高频系数处理过程的示意图。

图5a为本发明一个实施例的待融合MRI图像。

图5b为本发明一个实施例的待融合CT图像。

图5c为基于平均加权AVG的方法对图5a和图5b的融合结果。

图5d为基于小波变换DWT的方法对图5a和图5b的融合结果。

图5e为基于梯度金子塔变换GRA的方法对图5a和图5b的融合结果。

图5f为基于对比度金子塔变换CON的方法对图5a和图5b的融合结果。

图5g为基于四元数小波变换的平均加权融合方法QWT-Avg-Max对图5a和图5b的融合结果。

图5h为本发明对图5a和图5b的融合结果。

图6a为本发明一个实施例的待融合T1-MR图像。

图6b为本发明一个实施例的待融合MRA图像。

图6c为基于平均加权AVG的方法对图6a和图6b的融合结果。

图6d为基于小波变换DWT的方法对图6a和图6b的融合结果。

图6e为基于梯度金子塔变换GRA的方法对图6a和图6b的融合结果。

图6f为基于对比度金子塔变换CON的方法对图6a和图6b的融合结果。

图6g为基于四元数小波变换平均加权QWT-Avg-Max的方法对图6a和图6b的融合结果。

图6h为本发明对图6a和图6b的融合结果。

具体实施方式

下面对本发明的两个实施例MRI和CT医学图像以及T1-MR和MRA医学图像结合附图作详细说明。本实施例在以本发明技术方案为前提下进行，如图1所示，详细的实施方式和具体的操作步骤如下：

步骤1，准备待融合的两幅M×N(本实施例中M＝256,N＝256)源图像A和B，利用QWT分别将两幅图像分解为低频部分L与高频部分H。QWT的分解层数设置为两层j＝2，方向个数为3。

步骤2，对A、B的QWT系数的低频系数Z＝A,B，根据绝对值取大的融合规则进行融合，得到低频的融合系数具体描述见公式：

由于QWT系数的低频系数包括四个子带，故上式中n的取值范围为1,2,3,4。当待融合图像A的子带系数的绝对值不小于待融合图像B的子带系数则融合系数 取A的对应系数；反之，融合系数取B的对应系数。

步骤3，高频子带系数z＝A或B,求取其对应的清晰度系数首先，基于当前高频子带系数设置局部窗体w,w＝4×4，然后计算局部窗体的清 晰度值CW。CW即为子带系数的清晰度系数其中，j为分解尺度，n为第k,k＝1,2,3个分解方向的第n,n＝1,2,3,4,个子带。将CW执行QWT分解，分解尺度为1，按下述公式计算CW：

其中：PhaseLH_hor代表四元数小波变换后LH子带的水平方向的相位信息：PhaseHL_ver代表四元数小波变换后HL子带的垂直方向的相位信息，PhaseHH_dia代表四元数小波变换后HH子带的对角方向的相位信息。GradL_hor(i′,j′)、GradL_ver(i′,j′)和GradL_dia(i′,j′)分别代表低频幅值系数LAM(i′,j′)的水平、垂直以及对角方向的梯度值。

步骤4，在获得清晰度系数后，采用基于上下文的融合规则实现融合规则。

四元数变换系数的分布不仅受其所在的尺度及方向的影响，还受到该系数在相邻尺度、相邻位置的系数的影响。一般情况下，变换系数的分布有其持续性、相似性与传递性。同一尺度下，变换系数大的地方相邻的系数一般会比较大，变换系数小的地方相邻的系数就会比较小,变换系数的这一关系我们称之为邻域关系。同样，变换系数大的地方其上一尺度与下一尺度对应位置的系数也会比较大，反之亦然，变换系数的这一关系我们称之为父子关系。四元数小波变换系数的邻域关系和父子关系统称为系数的上下文关系。

QWT系数上下文设计如下式所示：

其中，代表四元数小波变换高频子带系数的活动测度，此处，本发明采用子带系数对应的清晰度系数作为活动测度。是基于上下文的活动测度值，[·]^↓2表示以2为降采样因子的降采样过程，erf()为Sigmoid函数。Sigmoid函数一般用在神经网络作为其激活函数，用以引入非线性关系。erf()定义为：

最后，高频子带的融合系数由下式得到：

从而，完成图像各个高频方向子带系数的融合。

步骤5，根据得到的高频和低频的融合系数，利用QWT逆变换得到融合图像。

本发明的效果可以通过以下实验结果作进一步说明：

1、实验条件与方法

硬件平台：处理器Intel Core i3-4150，CPU主频2.4GHz，内存4GB；

软件平台：计算机采用64位的Windows 7专业版，Matlab R2012b；实验中采用两组已配准的源图像，即MRI和CT医学图像以及T1-MR和MRA医学图像，图像大小均为256×256，bmp格式。

实验时采用现有的五种融合方法作为对比方法，其中：

方法1为基于平均加权AVG的融合方法；

方法2为基于小波变换DWT的融合方法；

方法3为基于梯度金子塔变换GRA的融合方法；

方法4为基于对比度金子塔变换CON的融合方法；

方法5为基于四元数小波变换平均加权QWT-Avg-Max的融合方法。

2、仿真内容

仿真一：遵循本发明的技术方案，对MRI和CT医学图像(见图5(a)和图5(b))进行融合，图5(c)-图5(h)为本发明方法和对比方法的仿真实验结果。从主观视觉效果来看，方法1的融合图像丢失了部分信息，而且对比度较低；方法2的融合图像中MRI图像的信息整体亮度降低，产生了一定的失真；方法3的融合图像细节丢失且也存在亮度下降的问题；方法4的融合图像中CT图像信息基本融合且对比度大大升高，但是MRI图像信息丢失严重，亮度也大大降低；方法5的融合图像存在许多伪影；本发明的融合图像尽可能的保留了两幅原始图像的有用信息，而且在亮度、对比度等评价上依然保持着原始图像的标准，没有产生失真或者扭曲等现象。

仿真二：遵循本发明的技术方案，对T1-MR和MRA医学图像(见图6(a)和图6(b))进行融合，图6(c)-图6(h)为本发明方法和对比方法的仿真实验结果。从主观视觉效果来看，方法1的融合图像对比度大大降低，细节模糊；方法2的融合图像对比度降低的同时某些细节部分的信息丢失了；方法3的融合图像与方法2的融合图像结果相似，但是在关键信息上亮度 的降低更为严重；方法4的融合图像产生了很多噪声点，导致融合图像不均匀，无法准确反映两幅原始图像的信息；方法5的融合图像反映了两幅原始图像的基本信息，在亮度和对比度等度量上也表现不错，但是存在细小特征的丢失，而且图像整体感觉比较模糊。本发明的融合图像较好的保留了原图像的细节和有用信息，整体效果理想。

将本发明的融合结果与对比方法的融合结果进行客观指标评价。

表1给出了MRI和CT医学图像融合结果的客观评价指标，

表2给出了T1-MR和MRA医学图像融合结果的客观评价指标，

其中，最优数据以加粗的形式表示。

表1.MRI和CT医学图像融合结果的客观评价指标

融合方法	STD	EN	AG	SF	MI	Qabf
							方法1	35.5610	5.9196	3.9869	10.4552	5.2315	0.4281
方法2	47.2734	5.9969	6.8740	18.7529	2.0803	0.6215
							方法3	41.8812	6.2197	5.4719	14.3150	2.9434	0.6878
方法4	44.5781	4.6789	4.8582	21.6927	2.0441	0.2877
							方法5	35.7910	5.9881	4.8078	12.2589	3.2392	0.3802
本发明	62.9844	6.7932	7.2336	18.1123	5.1287	0.7478

表2.T1-MR和MRA医学图像融合结果的客观评价指标

融合方法	STD	EN	AG	SF	MI	Qabf
							方法1	47.2944	5.1260	5.4013	13.0686	3.6555	0.3784
方法2	55.7298	5.3645	9.0421	22.2881	2.8047	0.4577
							方法3	48.7541	5.7278	7.2559	17.5448	3.2371	0.4995
方法4	55.3919	4.7830	11.4372	29.5195	3.3181	0.2738
							方法5	72.7948	5.6723	8.2632	19.5600	3.6193	0.4493
本发明	73.1765	5.7502	8.1832	19.9434	3.7897	0.5198

上述评价指标中，标准差STD越大，证明图像内部反差越大，能获取到的信息更多；熵EN越大则代表融合图像融入待融合图像的信息越多；平均梯度AG和空间频率SF越大，则融合图像的细节信息越完整；而互信息MI与边缘转换率Qabf的值越大，就说明原始图像与融合图像之间的像素相似度以及结构相似度越高，融合图像的效果越好。

从表1可以看出，本发明的客观指标除了SF和MI外其他指标都优于对比方法。从表2可以看出，本发明的客观指标除了AG和SF外其他指标都优于对比方法。

从各仿真实验的融合结果可以看出，无论是主观感知还是客观比较本发明的融合图像都具有较好的图像质量，全局清晰，含有丰富的边缘和细节信息。仿真实验检验了本发明的有效性。

Claims (5)

1.一种基于四元数小波变换上下文结构的医学图像融合方法，其特征是：首先对待融合图像通过QWT将图像由空间域转换到频率域得到图像的QWT系数，对于QWT低频系数，采用对应位置系数绝对值取大的融合规则；对于QWT高频系数，进行局部开窗计算清晰度值，然后采用基于上下文的融合规则；最后根据得到的低频和高频的融合系数，利用QWT逆变换得到融合图像。

2.根据权利要求1所述的基于四元数小波变换上下文结构的医学图像融合方法，其特征在于包括以下具体步骤：

1)通过QWT将图像A、B由空间域转换到频率域得到图像A、B对应的QWT系数；

2)对A、B的QWT系数的低频部分根据绝对值取大的融合规则进行融合，得到低频的融合系数；

3)计算A、B的高频子带系数对应的清晰度系数；

4)基于上下文统计模型与清晰度系数设计融合规则，获得融合后的高频系数；

5)根据得到的低频和高频的融合系数，利用QWT逆变换得到融合图像。

3.根据权利要求2所述的基于四元数小波变换上下文结构的医学图像融合方法，其特征在于：步骤2)低频部分根据绝对值取大的融合规则进行融合，具体描述见公式：

$C_{F,L}^n(x,y)=\left\{\begin{array}{c}{C}_{A,L}^n(x,y),if\left|C_{A,L}^n(x,y)\right|\&GreaterEqual;\left|C_{B,L}^n(x,y)\right|\\ C_{B,L}^n(x,y),if\left|C_{A,L}^n(x,y)\right|<\left|C_{B,L}^n(x,y)\right|\end{array}\right.,$

由于QWT系数的低频系数包括四个子带，故上式中n的取值范围为1,2,3,4。当待融合图像A的子带系数的绝对值不小于待融合图像B的子带系数则融合系数取A的对应系数；反之，融合系数取B的对应系数。

4.根据权利要求2所述的基于四元数小波变换上下文结构的医学图像融合方法，其特征在于：步骤3)计算高频子带系数对应的清晰度系数。高频子带系数z＝A或B,求取其对应的清晰度系数首先，基于当前高频子带系数设置局部窗体w，然后计算局部窗体的清晰度值CW。CW即为子带系数的清晰度系数其中，j为分解尺度，n为第k,k＝1,2,3个分解方向的第n,n＝1,2,3,4,个子带。将CW执行QWT分解，由变换系数根据四元数代数计算得到各子带(LL、LH、HL、HH子带)对应的幅值与相位信息，其中低频LL子带水平方向、垂直方向、对角方向的相位分别记为PhaseLL_hor、PhaseLL_ver、PhaseLL_dia，高频LH子带水平方向的相位记为PhaseLH_hor，高频HL子带垂直方向的相位记为PhaseHL_ver，高频HH子带对角方向的相位记为PhaseHH_dia；

按下述公式计算CW：

CW＝∑_{i′,j′∈w}|GradL_hor(i′,j′)×PhaseLL_hor(i′,j′)×PhaseLH_hor(i′,j′)|+∑_{i′,j′∈w}|GradL_ver(i′,j′)×PhaseLL_ver(i′,j′)×PhaseHL_ver(i′,j′)|+∑_{i′,j′∈w}|GradL_dia(i′,j′)×PhaseLL_dia(i′,j′)×PhaseHH_dia(i′,j′)|；

$\left\{\begin{array}{c}{\mathrm{GradL}}_{hor}(i^{\&prime;},j^{\&prime;})=LAM(i^{\&prime;},j^{\&prime;}+1)-LAM(i^{\&prime;},j^{\&prime;})\\ {\mathrm{GradL}}_{ver}(i^{\&prime;},j^{\&prime;})=LAM(i^{\&prime;}+1,j^{\&prime;})-LAM(i^{\&prime;},j^{\&prime;})\\ {\mathrm{GradL}}_{dia}(i^{\&prime;},j^{\&prime;})=LAM(i^{\&prime;}+1,j^{\&prime;}+1)-LAM(i^{\&prime;},j^{\&prime;})\end{array}\right.;$

其中：GradL_hor(i′,j′)、GradL_ver(i′,j′)和GradL_dia(i′,j′)分别代表低频幅值系数LAM(i′,j′)的水平、垂直以及对角方向的梯度值。

5.根据权利要求2所述的基于四元数小波变换上下文结构的医学图像融合方法，其特征在于：步骤4)对于QWT高频系数获得清晰度系数之后，采用基于上下文的融合规则实现融合规则。上下文的构建过程如下：

四元数变换系数的分布不仅受其所在的尺度及方向的影响，还受到该系数在相邻尺度、相邻位置的系数的影响。一般情况下，变换系数的分布有其持续性、相似性与传递性。同一尺度下，变换系数大的地方相邻的系数一般会比较大，变换系数小的地方相邻的系数就会比较小,变换系数的这一关系我们称之为邻域关系。同样，变换系数大的地方其上一尺度与下一尺度对应位置的系数也会比较大，反之亦然，变换系数的这一关系我们称之为父子关系。四元数小波变换系数的邻域关系和父子关系统称为系数的上下文关；。

QWT系数上下文设计如下式所示：

${\mathrm{SM}}_{k,n}^j(x,y)=\left\{\begin{array}{c}\left(erf\right({\max }_{ni,nj\&Element;N}\left(M_{k,n}^j\right(x_{ni},y_{nj})-{\stackrel{\&OverBar;}{M}}_{k,n}^j)),,j=1;\\ \max \{\left(erf\right({\max }_{ni,nj\&Element;N}\left(M_{k,n}^j\left(x_{ni},y_{nj}\right)-{\stackrel{\&OverBar;}{M}}_{k,n}^j\right)),{\&lsqb;{\mathrm{SM}}_{k,n}^{j-1}(x,y)\&rsqb;}^{\&DownArrow;2}\},j\&Element;\&lsqb;2,J\&rsqb;\end{array}\right.;$

其中，代表四元数小波变换高频子带系数的活动测度，代表第n个子带的系数均值，此处，本发明采用子带系数对应的清晰度系数作为活动测度。N是当前系数的邻域，ni,nj是邻域的系数下标。是基于上下文的活动测度值，[·]^↓2表示以2为降采样因子的降采样过程，erf()为Sigmoid函数。Sigmoid函数一般用在神经网络作为其激活函数，用以引入非线性关系。erf()定义为：

$erf\left(\tilde{x}\right)=\frac{1}{1+e^{-\left(\tilde{x}\right)}};$

最后，高频子带的融合系数由下式得到：

$C_{F,H}^{j,k,n}(x,y)=\left\{\begin{array}{c}{C}_{A,H}^{j,k,n}(x,y),if{\mathrm{SM}}_{A,H}^{j,k,n}(x,y)\&GreaterEqual;{\mathrm{SM}}_{B,H}^{j,k,n}(x,y)\\ C_{B,H}^{j,k,n}(x,y),if{\mathrm{SM}}_{A,H}^{j,k,n}(x,y)<{\mathrm{SM}}_{B,H}^{j,k,n}(x,y)\end{array}\right.;$

从而，完成图像各个高频方向子带系数的融合。