CN106023126A - 一种基于区域对比的医学图像融合方法 - Google Patents

Abstract

本发明公开了一种基于区域对比的医学图像融合方法，属于医学图像处理领域。该方法包括以下步骤：a、获取已配准的多模医学图像，根据多模医学图像的大小制定分块策略，按照从上到下、从左到右的原则对图像进行分块操作；b、计算每一像素在不同分块下的对比度和分块图像的直方图矩阵，并以直方图矩阵作为系数，生成每一像素的区域对比度矩阵；c、应用典型相关分析方法计算每一像素区域对比度矩阵的相关系数，求解像素区域对比度矩阵最相关时的系数矩阵；d、以求解得到的系数矩阵作为权值，生成融合图像。与现有技术相比，本发明的基于区域对比的医学图像融合方法具有融合效果好、执行效率高等特点，具有很好的推广应用价值。

Description

一种基于区域对比的医学图像融合方法

技术领域

本发明涉及医学图像处理领域，具体地说是一种基于区域对比的医学图像融合方法。

背景技术

随着成像技术的发展，医学图像已成为现代医学诊断和治疗的重要手段。目前，临床中常用的医学图像主要分为解剖图像和功能影像两种。解剖图像(如CT、MR等)能清晰的提供器官的解剖结构、形态特征，并且具有较高的分辨率；功能影像(如PET、fMRI等)可以准确提供器官的功能和新陈代谢信息，但是其分辨率较低。单独一种模态的图像只能提供患者某些方面的特征，不能反映其全部病理信息。因此，在临床实践中，往往需要有机整合多种影像中的互补信息，融合多种模态的医学图像，为精确诊断和精准治疗提供依据。

现有技术中，多模医学图像融合方法按照处理策略可以分为像素级、特征级、决策级三个层次，像素级融合方法主要有空域融合、变换域融合、智能融合等。空域方法通过对像素值本身的加权平均或形态学处理进行融合，其计算简单，融合效果较差。变换域方法首先对医学图像进行金字塔、小波或多尺度几何变换，提取变换域子带系数，针对子带系数制定不同的融合规则，最后通过反变换生成融合图像，其融合效果较好，但计算复杂度高，实时性差。智能融合方法通过神经网络、统计学习、模糊集等算法迭代计算融合系数，生成融合图像，其融合效果好坏往往依赖于学习算法的设计，且同样存在计算复杂度高，实时性差等问题。总的来讲，现有的像素级融合方法都是基于多模医学图像的结构、形态、颜色、时频信息等特征制定融合策略，缺乏对像素空间位置关系及对比度等因素的分析；同时，现有方法一般采取加权平均或绝对值取大作为融合规则，易使融合结果陷入全局最优，导致局部细节丢失或模糊，难以满足实际应用需求。

发明内容

本发明的技术任务是针对上述现有技术的不足，提供一种应用像素空间位置关系及对比度制定融合规则的基于区域对比的医学图像融合方法。

本发明的技术任务是按以下方式实现的：一种基于区域对比的医学图像融合方法，其特点是包括以下步骤：

a、获取已配准的多模医学图像，根据多模医学图像的大小制定分块策略，按照从上到下、从左到右的原则对图像进行分块操作；

b、计算每一像素在不同分块下的对比度和分块图像的直方图矩阵，并以直方图矩阵作为系数，生成每一像素的区域对比度矩阵；

c、应用典型相关分析方法计算每一像素区域对比度矩阵的相关系数，求解像素区域对比度矩阵最相关时的系数矩阵；

d、以求解得到的系数矩阵作为权值，生成融合图像。

作为优选，步骤a中优选采用以下方式对图像进行分块操作：图像分块操作分t次进行，每次分块尺寸分别为多模医学图像大小的1，1/4，1/8…1/2^t倍。t为大于等于2的正整数，定义配准后多模医学图像A和B大小均为m×n(单位：像素)，则t的最大取值为即

步骤b中采用以下方法计算每一像素区域对比度矩阵：

定义多模医学图像A和B的分块分别表示为BA_i和BB_i，i表示第i次分块，i的取值为大于等于2且小于等于t的正整数，分块内像素点(x,y)的对比度C(x,y)_i为：

$CA{(x,y)}_i=\left|\frac{A(x,y)-{\mathrm{MA}}_i}{max\left({\mathrm{BA}}_i\right)-min\left({\mathrm{BA}}_i\right)}\right|$

$CB{(x,y)}_i=\left|\frac{B(x,y)-{\mathrm{MB}}_i}{max\left({\mathrm{BB}}_i\right)-\min \left({\mathrm{BB}}_i\right)}\right|$

分块BA_i和BB_i的直方图矩阵HA_i、HB_i，则每一像素区域对比度矩阵C(x,y)为：

CA(x,y)＝[HA(x,y)₁×CA(x,y)₁…HA(x,y)_i×CA(x,y)_i…HA(x,y)_t×CA(x,y)_t]′

CB(x,y)＝[HB(x,y)₁×CB(x,y)₁…HB(x,y)_i×CB(x,y)_i…HB(x,y)_t×CB(x,y)_t]′。

步骤c利用典型相关分析求解区域对比度矩阵的系数矩阵的方法包括：

首先，生成待融合图像像素点的对比度矩阵向量Z：

Z＝[CA(x,y)CB(x,y)]

计算其协方差矩阵Σ：

$\mathrm{\Σ}=\left[\begin{array}{cc}{\mathrm{\Σ}}_{11}& {\mathrm{\Σ}}_{12}\\ {\mathrm{\Σ}}_{21}& {\mathrm{\Σ}}_{22}\end{array}\right]$

其中，∑₁₁＝Cov(CA(x,y),CA(x,y))

∑₁₂＝∑₂₁＝Cov(CA(x,y),CB(x,y))

∑₂₂＝Cov(CB(x,y),CB(x,y))

其次，令U＝A₁CA(x,y)，V＝B₁CB(x,y)作为两组变量的线性组合，则有，

Cov(U,V)＝A₁∑₁₂B₁'

其中，A₁、B₁为相关系数矩阵，即需要求解的系数矩阵,B₁'为B₁的转置矩阵。

步骤d的融合结果为：

C(x,y)＝w₁A(x,y)+w₂B(x,y)

其中，w₁、w₂为A₁、B₁前r维系数，r为CA(x,y)和CB(x,y)的秩。

与现有技术相比，本发明的基于区域对比的医学图像融合方法具有以下突出的有益效果：

(一)结合了区域相似性、像素空间位置关系等因素作为融合策略，利用相关系数制定融合规则，使每一像素点融合结果与原多模图像的相关性最高；

(二)采用特有的分块方式，兼顾像素的全局及局部对比度，避免了融合结果的局部最优。

附图说明

附图1是本发明基于区域对比的医学图像融合方法的流程图；

附图2是本发明提供的一例颅脑CT图像实例；

附图3是本发明提供的一例颅脑MR图像实例；

附图4是基于PCA方法对图2和图3融合结果实例；

附图5是基于DCT变换方法对图2和图3融合结果实例；

附图6是基于NSCT方法对图2和图3融合结果实例；

附图7是基于LPSR方法对图2和图3融合结果实例；

附图8是基于GFF方法对图2和图3融合结果实例；

附图9是基于本发明方法对图2和图3融合结果实例。

具体实施方式

参照说明书附图以具体实施例对本发明的基于区域对比的医学图像融合方法作以下详细地说明。

【实施例】

如附图1所示，本发明的基于区域对比的医学图像融合方法的具体实施步骤为：

(1)获取已配准多模态医学图像图2和图3，图2为颅脑CT图像，图3为MR图像，大小均为512×512，分别将其表示为A和B。按照从左到右、从上到下的原则，对图像进行t次不同大小的分块操作，即t∈[2,10]。t值可以根据实际需要自由选取，为了便于计算，本实施例中t＝5。

设置分块尺寸分别为原始图像大小的1，1/4，1/8…1/2^t倍。将A和B的每一分块表示为BA_i和BB_i，i表示第i次分块操作，即大小为原始图像的1/2ⁱ倍时的分块表示。

(2)计算分块BA_i和BB_i的直方图矩阵HA_i、HB_i，直方图矩阵反映了每一像素在不同的分快内的概率分布。之后，根据不同分块内的像素均值MA_i、MB_i，计算分块内像素点(x,y)的对比度，

$CA{(x,y)}_i=\left|\frac{A(x,y)-{\mathrm{MA}}_i}{max\left({\mathrm{BA}}_i\right)-min\left({\mathrm{BA}}_i\right)}\right|$

$CB{(x,y)}_i=\left|\frac{B(x,y)-{\mathrm{MB}}_i}{max\left({\mathrm{BB}}_i\right)-\min \left({\mathrm{BB}}_i\right)}\right|$

按照分块大小进行排列，以直方图矩阵HA_i、HB_i作为系数，生成每一像素的区域对比度矩阵。

CA(x,y)＝[HA(x,y)₁×CA(x,y)₁…HA(x,y)_i×CA(x,y)_i…HA(x,y)_t×CA(x,y)_t]′

CB(x,y)＝[HB(x,y)₁×CB(x,y)₁…HB(x,y)_i×CB(x,y)_i…HB(x,y)_t×CB(x,y)_t]′

(3)利用典型相关分析方法计算每一像素点对比度矩阵的相关系数。

首先，生成待融合图像像素点的对比度矩阵向量

Z＝[CA(x,y)CB(x,y)]

计算其协方差矩阵

$\mathrm{\Σ}=\left[\begin{array}{cc}{\mathrm{\Σ}}_{11}& {\mathrm{\Σ}}_{12}\\ {\mathrm{\Σ}}_{21}& {\mathrm{\Σ}}_{22}\end{array}\right]$

其中，∑₁₁＝Cov(CA(x,y),CA(x,y))

∑₁₂＝∑₂₁＝Cov(CA(x,y),CB(x,y))

∑₂₂＝Cov(CB(x,y),CB(x,y))

其次，令U＝A₁CA(x,y)，V＝B₁CB(x,y)作为两组变量的线性组合，A₁、B₁为相关系数矩阵。则有，

Cov(U,V)＝A₁∑₁₂B₁'

最后，根据Cov(U,V)，求解得到A₁、B₁。这里求解的方法有很多种，但结果总是唯一，本实施例应用Lagrange乘数法进行求解。

(4)以相关系数矩阵A₁、B₁作为权值，生成融合图像。

C(x,y)＝w₁A(x,y)+w₂B(x,y)

其中，w₁、w₂为A₁、B₁前r维系数，r为CA(x,y)和CB(x,y)的秩，本实施例中r＝1。

【实验及结果】

下面结合实施例对本发明融合效果进行说明。图2和图3为已完成配准的颅脑CT和MR影像，图像大小为512×512。图4、图5、图6、图7、图8分别为利用主成分分析(PCA)、离散余弦变换(DCT)、非下采样轮廓波变换(NSCT)、稀疏拉普拉斯变换(LPSR)、边缘保持滤波(GF)等方法融合后生成的实例图像，图9为利用本发明方法生成的融合图像。从实施例中可以看出，图4、图6、图8融合结果各丢失了部分颅内细节信息，图5、图7颅内软组织细节模糊、对比度低，图9具有较高的对比度，颅内软组织细节清晰。

表1比较了各组融合方法的性能，选取弧度/角度相似性(Q_AB|F)、联合熵(Entropy)及运行时间(Time)作为参数，实验平台为普通电脑，酷睿i5四核处理器，8G内存，Win7 64位操作系统，MATLAB R2014A运行环境。

	PCA	DCT	NSCT	LPSR	GF	本发明
							QAB|F	0.5245	0.5531	0.6104	0.6387	0.6423	0.6478
Entropy	2.5742	2.9824	3.1695	3.8616	3.6239	4.7803
							Time(s)	0.3416	0.1375	1.3010	0.0510	0.1713	0.2041

表1 各组融合方法性能参数

从表1中可以看出，与上述多种算法相比，本发明方法联合熵、弧度/角度相似性明显提高，具有较好的融合效果；同时，计算时间较小，可满足实时性需求，具有较高的执行效率。

综上所述，上述实施方式仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，进行的无创新性修改和替换，都应涵盖在本发明的保护范围之内。

Claims (5)

1.一种基于区域对比的医学图像融合方法，其特征在于包括以下步骤：

a、获取已配准的多模医学图像，根据多模医学图像的大小制定分块策略，按照从上到下、从左到右的原则对图像进行分块操作；

b、计算每一像素在不同分块下的对比度和分块图像的直方图矩阵，并以直方图矩阵作为系数，生成每一像素的区域对比度矩阵；

c、应用典型相关分析方法计算每一像素区域对比度矩阵的相关系数，求解像素区域对比度矩阵最相关时的系数矩阵；

d、以求解得到的系数矩阵作为权值，生成融合图像。

2.根据权利要求1所述的基于区域对比的医学图像融合方法，其特征在于步骤a中图像分块操作分t次进行，每次分块尺寸分别为多模医学图像大小的1，1/4，1/8…1/2^t倍，t为大于等于2的正整数，定义配准后多模医学图像A和B大小均为m×n，则t的最大取值为

3.根据权利要求1或2所述的基于区域对比的医学图像融合方法，其特征在于步骤b中采用以下方法计算每一像素区域对比度矩阵：

定义多模医学图像A和B的分块分别表示为BA_i和BB_i，i表示第i次分块，i的取值为大于等于2且小于等于t的正整数，分块内像素点(x,y)的对比度C(x,y)_i为：

$CA{(x,y)}_i=\left|\frac{A(x,y)-{\mathrm{MA}}_i}{max\left({\mathrm{BA}}_i\right)-min\left({\mathrm{BA}}_i\right)}\right|$

$CB{(x,y)}_i=\left|\frac{B(x,y)-{\mathrm{MB}}_i}{max\left({\mathrm{BB}}_i\right)-\min \left({\mathrm{BB}}_i\right)}\right|$

分块BA_i和BB_i的直方图矩阵HA_i、HB_i，则每一像素区域对比度矩阵C(x,y)为：

CA(x,y)＝[HA(x,y)₁×CA(x,y)₁ … HA(x,y)_i×CA(x,y)_i … HA(x,y)_t×CA(x,y)_t]′

CB(x,y)＝[HB(x,y)₁×CB(x,y)₁ … HB(x,y)_i×CB(x,y)_i … HB(x,y)_t×CB(x,y)_t]′。

4.根据权利要求3所述的基于区域对比的医学图像融合方法，其特征在于步骤c利用典型相关分析求解区域对比度矩阵的系数矩阵的方法为：

首先，生成待融合图像像素点的对比度矩阵向量Z：

Z＝[CA(x,y) CB(x,y)]

计算其协方差矩阵Σ：

$\mathrm{\Σ}=\left[\begin{array}{cc}{\mathrm{\Σ}}_{11}& {\mathrm{\Σ}}_{12}\\ {\mathrm{\Σ}}_{21}& {\mathrm{\Σ}}_{22}\end{array}\right]$

其中，Σ₁₁＝Cov(CA(x,y),CA(x,y))

Σ₁₂＝Σ₂₁＝Cov(CA(x,y),CB(x,y))

Σ₂₂＝Cov(CB(x,y),CB(x,y))

其次，令U＝A₁CA(x,y)，V＝B₁CB(x,y)作为两组变量的线性组合，则有，

Cov(U,V)＝A₁Σ₁₂B₁'

其中，A₁、B₁为需要求解的系数矩阵,B₁'为B₁的转置矩阵。

5.根据权利要求4所述的基于区域对比的医学图像融合方法，其特征在于步骤d的融合结果为：

C(x,y)＝w₁A(x,y)+w₂B(x,y)

其中，w₁、w₂为A₁、B₁前r维系数，r为CA(x,y)和CB(x,y)的秩。