基于多向经验模式分解的医学图像融合方法

摘要

本发明涉及基于多向经验模式分解的医学图像融合方法，该方法采用多向经验模式分解对采集的不同种类的医学源图像进行多尺度多向分解，获得源图像的多尺度多向的高频分量内蕴模式分量，按照区域能量规则进行融合处理，可有效地提取各源图像的高频细节信息；对源图像的低频剩余分量采用能量贡献规则进行融合处理，最后反变换获取融合图像，融合后的图像有效提高融合图像的目视效果，避免小波、超小波、窗口经验模式分解融合算法引起的融合图像出现局部畸变或缺失的缺点，无需人为进行参数的选择，并且能够很好地提取源图像的细节信息，实现自适应地医学图像融合基于全新的多尺度分解结构，具有完全数据驱动的自适应性，具有更强的细节信息获取能力。

说明书

技术领域

本发明属于图像融合技术领域，涉及基于多向经验模式分解的医 学图像融合方法。

背景技术

图像融合是以图像为主要研究内容的信息融合技术，是将两幅或 多幅图像合成为一副图像，以获取对同一场景的更加精确、更为全面、 更为可靠的图像描述。图像融合技术通过有效地利用多幅图像间的冗 余性和互补性，使融合后的图像更适合人眼视觉系统，适合计算机的 理解、分析及后续处理的需求。

医学图像融合是图像融合中的一项重要应用。目前，医学成像设 备主要有：计算机断层扫描(Computer Tomography，CT)，磁共振成 像(Magnetic Resonance Imaging，MRI)或较新的正电子发射断层扫 描(Positron Emission Tomography，PET)。CT图像反应的是人体组 织的X线吸收系数，而骨骼组织的吸收系数最大，因此CT图像中骨 骼组织最清晰。MRI图像反应的是人体组织的质子密度，而软组织的 质子密度最高，因此MRI图像中软组织最清晰。PET是目前惟一可在 活体上显示生物分子代谢、受体及神经介质活动的新型影像技术，当 疾病早期处于分子水平变化阶段，病变区的形态结构尚未呈现异常， MRI、CT检查还不能明确诊断时，PET即可发现病灶所在。将三类 图像进行融合，可有效地提高医生的诊断效果。

目前，基于多分辨率、多尺度分解的融合算法，在图像融合中得 到广泛应用。各国研究者提出Wavelet变换、Ridgelet变换、Curvelet变 换、Contourlet变换和Bandelet变换等多种小波和超小波变换的处理 方法，就是这方面的重要研究成果。但无论是基于哪种小波，在图像 融合中都存在一个问题：融合后的图像会在局部位置出现畸变，这在 一般应用中还可接受，但在医学图像中，这就有可能造成医生的误诊， 影响病者的后续治疗。因此，工程界和数学界从未停止过探索更好的 分解算法。

1999年，美国宇航局的Norden E.Huang教授发明了经验模式分 解算法(Empirical Mode Decomposition，EMD)，能将非稳定非线性 信号按频率做自适应分解。二维经验模式分解是一维EMD分解算法 在二维平面上的推广，可用于图像数据的分析和处理，通过将原始图 像自适应的分解为有效数量的子图像，可以将图像从高频到低频的局 部窄带细节信息内蕴模式分量分解出来，剩余分量表示图像的趋势。 分解出来内蕴模式分量具有当前图像的纹理信息。但传统二维经验模 式分解有缺陷：分解得到的内蕴模式分量图像中有暗斑。因此严重影 响了传统二维经验模式分解在图像处理领域中的应用。后来出现的窗 口经验模式分解WEMD较好地解决了传统二维经验模式分解的缺陷， 又保留了传统二维经验模式分解自适应分解特性，并已在图像融合中 得到应用，但融合后的图像由于出现少量的待融合图像信息缺失还有 待进一步地提高。为进一步提高融合效果，本发明提出了多向经验模 式分解算法，其应用还处于起步阶段，少有关于图像融合的讨论，基 于多向经验模式分解算法的融合规则也少有具体实现，需深入探索。

综上，目前现有的融合技术应用于医学图像融合，还存在一些不 足：基于小波、超小波的融合图像会出现局部畸变，传统经验模式分 解方法分解得到的内蕴模式分量图像暗斑都对融合结果有很大的影 响，窗口经验模式分解的融合图中，会出现少量的细节信息缺失，并 且融合规则的好坏对融合效果也有巨大的影响。

发明内容

本发明的目的是提供一种基于多向经验模式分解的医学图像融 合方法，以达到既不出现小波、超小波融合算法融合后的图像出现畸 变，又不出现窗口经验模式分解在医学图像融合中出现细节缺失，还 无需人为进行参数的选择，并且能够很好地提取源图像的细节信息， 实现自适应地医学图像融合。

为实现上述目的，本发明的基于多向经验模式分解的医学图像融 合方法采用多向经验模式分解方法对采集到的源图像进行多尺度分 解，获得源图像的多级尺度多方向的高频内蕴模式分量，按照区域能 量规则进行融合处理，并对源图像的低频剩余分量采用能量贡献规则 进行融合处理，最后反向重构获取融合图像，该方法的步骤如下：

(1)采用多向经验模式分解对匹配好的不同种类的待融合医学 源图像进行多尺度多向分解，获得源图像的多尺度多向的高频分量内 蕴模式分量imf_ij和低频剩余分量r_i，其中i＝1，2，3，……，m，m为待融合图 像的数量，j＝1，2，3，……，n，n为分解得到的imf的级数；

(2)将相同级的待融合图像的高频分量内蕴模式分量imf_ij按照 区域能量规则进行融合处理，产生融合图像的第j级内蕴模式分量 imf_j；

(3)将待融合图像的r_i，按照整体能量贡献规则进行处理，得到 融合图像的剩余分量r；

(4)将imf_j和r反向重构并进行坐标系反变换得到融合图像。

所述的步骤(1)中对每一幅源图像进行图像高频到低频的多尺 度多方向分解，首先，分解处理的第1级内蕴模式分量是图像所含有 的最高频率分量，源数据减去第1级内蕴模式分量得到第1级剩余分 量；对第1级剩余分量再分解，得到第2级内蕴模式分量和第2级剩 余分量；以此类推，得到n级内蕴模式分量和第n级剩余分量。

所述的多向经验模式分解的处理过程包括如下步骤：

(11)以待融合图像x的中心为原点，将直角坐标变换到极坐标， 得到变换矩阵y，s＝1，s表示极径的长度；

(12)r_s0＝y_s，y_s为极坐标系中具有相同极径长度的那一行或列， j＝1，j表示分解第j级内蕴模式分量；

(13)确定r_s(j-1)的所有局部极值点，并组成极大值点集和极小值点 集；

(14)根据分解级数j；

(a)设定当前最大窗口N，初始窗口M，窗口K＝M，M数值为奇 数；

(b)以r_s(j-1)中一个数值为当前中心，如果在窗口K内的极大值点 个数和极小值点个数相等，则求取窗口K内数值均值mean，转到步 骤(d)；

(c)K＝K+2，如果K＜N，转到步骤(b)；如果K≥N，直接求取窗 口K内数值均值mean；

(d)以mean作为当前数值的局部均值，转到下一个数值，K＝M， 转到步骤(b)，直至遍历r_s(j-1)的所有数值；

(15)用得到的所有局部均值点构成r_j-1的均值线h_j-1，并计算 imf_sj＝r_s(j-1)-h_s(j-1)，r_sj＝h_s(j-1)，j＝j+1；

(16)重复步骤(13)到(15)的操作，直到规定的分解层级数n， s＝s+1；

(17)重复步骤(12)到(16)的操作，直到所有相同极径长度的 列或行都分解完毕，得到变换矩阵y的每一级内蕴模式分量imf_j和最 后一级剩余分量r。

所述的步骤(13)中，图像的局部极大值点为灰度值比周围3区 域2个相邻像素点灰度值都高的点，图像的极小值点为灰度值比周围 2个相邻像素点灰度值都低的点。

所述的步骤(14)(a)中，最大窗口N是由所有局部极大和极小 值点分别组成极大值点集S1和极小值点集S2；分别遍历极值点集S1 和S2找出紧邻当前极值点距离最近的极值点，并计算两点间的距离 P_i(i＝1，2，3，…p；p为极大值的数目)和P′_j(j＝1，2，3，…q；q为极小值的数目)，以P_i或P′_j的和均值中值较小的一个取整作为最大窗口N。当N为偶数时， 执行N＝N+1操作，初始窗口M统一定义 等于3。

所述的步骤(2)包括：

(21)计算相同级imf_ij当前像素为中心的窗口Q的能量(Q为3× 3)，形成融合规则矩阵s_ij，i＝1，2，3，……，m(m为待融合图像的数量)， j＝1，2，3，……，n(n为分解得到的imf的级数)，q表示窗口Q中的一个 元素；$s_{\mathrm{ij}}(x,y)=\underset{q\in Q}{\Sigma }\left[\mathrm{abs}\left({\mathrm{imf}}_{\mathrm{ij}}(x,y)\right)\right]$

(22)比较相同级s_ij相同位置的值，选取最大的s_ij对应的imf_ij相同 位置的像素作为融合imf_j相同位置像素的值；

(23)遍历所有像素，得到最终的融合imf_j。

所述的步骤(3)包括：

(31)计算每一个剩余分量r_i的融合权重w_i，i＝1，2，3，……，m(m为 待融合图像的数量)，

$w_i=\frac{\underset{x=1}{\overset{N}{\Sigma }}\underset{y=1}{\overset{M}{\Sigma }}{r}_i(x,y)}{\underset{i}{\overset{m}{\Sigma }}\left[\underset{x=1}{\overset{N}{\Sigma }}\underset{y=1}{\overset{M}{\Sigma }}{r}_i(x,y)\right]}$

(32)根据权重计算出r，

所述的步骤(4)反向重构得到极坐标系下的融合图像I’， 再将极坐标化为直角坐标，得到最终融合图像I。

本发明的基于多向经验模式分解的医学图像融合方法，采用多向 经验模式分解对采集到的源图像进行多尺度多方向分解，分解过程继 承了传统经验模式分解的优点：完全由数据驱动，不需要预先设置滤 波器；求取各级内蕴模式分量类似高频滤波过程，与小波，超小波相 比，可获取更为丰富的细节信息，使得融合图像的质量得到提高，与 窗口经验模式分解相比，解决了窗口经验模式分解融合算法出现少量 信息丢失的问题。获得各级尺度的高频部分，按照区域能量规则对各 级内蕴模式分量进行融合处理，可最大化的融合细节信息，减少非细 节信息的影响。对源图像的低频剩余分量采用能量贡献规则进行融合 处理，可保留各剩余分量的能量分布。本方法既不出现小波、超小波 融合算法融合后的图像出现畸变，又不出现窗口经验模式分解在医学 图像融合中出现细节缺失，还无需人为进行参数的选择，并且能够很 好地提取源图像的细节信息，实现自适应地医学图像融合基于全新的 多尺度分解结构，具有完全数据驱动的自适应性，具有更强的细节信 息获取能力；运用区域能量规则对高频的内蕴模式分量进行融合处 理，运用能量贡献规则对低频的剩余分量进行融合处理，使得融合后 图像的质量得到提高，对于医学图像融合领域的应用具有重要意义和 实用价值。

附图说明

图1是本发明方法流程示意图；

图2是传统经验模式分解和多向经验模式分解对Lena图像两级分 解结果对比图；

图3是本发明融合效果与小波和超小波的融合对比图。

图2中(a)为Lena图像，(b)为传统经验模式分解得到的第1 级内蕴模式分量，(c)为传统经验模式分解得到的第2级内蕴模式分 量，(e)为传统经验模式分解得到的剩余分量，(f)为多向经验模式 分解得到的第1级内蕴模式分量，(g)为多向经验模式分解得到的第 2级内蕴模式分量，(h)为多向经验模式分解得到的剩余分量。

图3中(a)为CT图像，(b)为MRI图像，(c)为小波融合效果， (d)为超小波curvelet融合效果，(e)为窗口经验模式分解融合效果， (f)为本发明融合效果。

具体实施方式

为了更好地理解本发明的技术方案，下面结合附图对本发明的实 施方式进行详细的说明。

如图1所示，本发明首先对源图像x₁，x₂，…x_m进行多向经验模式分 解，得到每幅源图像的n级内蕴模式分量和一个剩余分量，对各级内 蕴模式分量按照区域能量规则进行融合处理；对源图像的低频剩余分 量采用能量贡献规则进行融合处理；最后反变换获取融合图像。

本发明的具体实施如下：

1.运用MDEMD算法将匹配好的待融合源图像x₁，x₂，…x_m分别进 行相同级数的分解，得到内蕴模式函数分量imf_ij和剩余分量r_i，其中 i＝1，2，3，……，m(m为待融合图像的数量)，j＝1，2，3，……，n(n为分解得 到的imf的数量)

多向经验模式分解的处理过程包括如下步骤：

步骤1：以图像x的中心为原点，将直角坐标化变换到极坐标，得 到变换矩阵y；

步骤1：以待融合图像x的中心为原点，将直角坐标变换到极坐标， 得到变换矩阵y，s＝1，s表示极径的长度；

步骤2：r₀＝y_j(y_j为极坐标系中具有相同极径的那一行或列，j表 示极径的长度)，i＝1，j＝1；

步骤2：r_s0＝y_s，y_s为极坐标系中具有相同极径长度的那一行或列， j＝1，j表示分解第j级内蕴模式分量；

步骤3：确定r_i-1的所有局部极值点，并组成极大值点集和极小值点 集；图像的局部极大值点为灰度值比周围3区域2个相邻像素点灰度 值都高的点，图像的极小值点为灰度值比周围2个相邻像素点灰度值 都低的点。

步骤3：确定r_s(j-1)的所有局部极值点，并组成极大值点集和极小值 点集；图像的局部极大值点为灰度值比周围3区域2个相邻像素点灰 度值都高的点，图像的极小值点为灰度值比周围2个相邻像素点灰度 值都低的点。

步骤4：根据分解层数i，

步骤4：根据分解层数j，

(a)设定当前最大窗口N，初始窗口M，窗口K＝M，M数值为 奇数；

步骤4(a)中，最大窗口N的定义为：由所有局部极大和极小值 点分别组成极大值点集S1和极小值点集S2；分别遍历极值点集S1和 S2找出紧邻当前极值点距离最近的极值点，并计算两点间的距离 P_i(i＝1，2，3，…p；p为极大值的数目)和P′_j(j＝1，2，3，…q；q为极小值的数目)，以P_i或P′_j的和均值中值较小的一个取整作为最大窗口N。当N为偶数时， 执行N＝N+1操作，初始窗口M统一定义 等于3。

(b)以r_s(j-1)中一个数值为当前中心，如果在窗口K内的极大值 点个数和极小值点个数相等，则求取窗口K内数值均值mean，转到 步骤(d)；

(c)K+2，如果K＜N，转到步骤(b)；如果K≥N，直接求取窗 口K内数值均值mean；

(d)以mean作为当前数值的局部均值，转到下一个数值，K＝M， 转到步骤(b)，直至遍历r_s(j-1)的所有数值；

步骤5：用得到的所有局部均值点构成r_j-1的均值线h_j-1，并计算 imf_sj＝r_s(j-1)-h_s(j-1)，r_sj＝h_s(j-1)，j＝j+1；

步骤6：重复步骤3到步骤5的操作，直到规定的分解层级数n， s＝s+1；

步骤7：重复步骤2到步骤6的操作，直到所有相同极径长度的列 或行都分解完毕，得到变换矩阵y的每一级内蕴模式分量imf_j和最后 一级剩余分量r。

多向经验模式分解实现了图像高频到低频的多尺度自适应分解 过程。首先，分解处理的第1级内蕴模式分量是图像所含有的最高频 率分量，源数据减去第1级内蕴模式分量得到第1级剩余分量；对第 1级剩余分量再分解，得到第2级内蕴模式分量和第2级剩余分量； 以此类推，得到n级内蕴模式分量和第n级剩余分量。对Lena图像 做2级分解，得到2级内蕴模式分量和一个剩余分量，如图2所示。 可见多向经验模式分解也可很好地解决了传统经验模式分解出现暗 斑的问题。多向经验模式分解继承了传统经验模式分解的优点，它的 也基是根据信号自适应产生的，使得它具有良好的时频局部性。

2.将相同级的待融合图像的高频分量内蕴模式分量imf_ij按照区域 能量规则进行融合处理，产生融合图像的第j级内蕴模式分量imf_j；

步骤1：计算相同级imf_ij当前像素为中心的窗口Q的能量(Q为3 ×3)，形成融合规则矩阵s_ij，i＝1，2，3，……，m(m为待融合图像的数量)， j＝1，2，3，……，n(n为分解得到的imf的级数)，q表示窗口Q中的一个 元素；

$s_{\mathrm{ij}}(x,y)=\underset{q\in Q}{\Sigma }\left[\mathrm{abs}\left({\mathrm{imf}}_{\mathrm{ij}}(x,y)\right)\right]$

步骤2：比较相同级s_ij相同位置的值，选取最大的s_ij对应的imf_ij相 同位置的像素作为融合imf_j相同位置像素的值。

步骤3：遍历所有像素，得到最终的融合imf_j。

3.将待融合图像的r_i，按照整体能量贡献规则进行处理，得到融 合图像的剩余分量r。

步骤1：计算每一个剩余分量r_i的融合权重w_i，i＝1，2，3，……，m(m 为待融合图像的数量)

$w_i=\frac{\underset{x=1}{\overset{N}{\Sigma }}\underset{y=1}{\overset{M}{\Sigma }}{r}_i(x,y)}{\underset{i}{\overset{m}{\Sigma }}\left[\underset{x=1}{\overset{N}{\Sigma }}\underset{y=1}{\overset{M}{\Sigma }}{r}_i(x,y)\right]}$

步骤2：根据权重计算出r

$r=\underset{i=1}{\overset{m}{\Sigma }}(w_i*r_i)$

4.反向重构得到极坐标系下的融合图像I’，再进行MDEMD反变 换得到最终的融合图像I。

步骤1：反向重构得到极坐标融合图像I’；

$I^{,}=\underset{j=1}{\overset{n}{\Sigma }}{\mathrm{imf}}_j+r$

步骤2：将极坐标化为直角坐标，得到最终融合图像I。

为了验证本发明的有效性，使用了同一部位的CT图像和MRI 图像进行融合处理。图3中(a)、(b)分别为CT图像和MRI图像， (c)为小波融合效果，(d)为超小波curvelet融合效果，(e)为窗口 经验模式分解融合效果，(f)为本发明融合效果。对比可见，小波和 超小波融合的图像出现局部畸变，不能最优的表示各待融合图像融合 细节纹理信息，窗口经验模式分解在融合图像中出现少量的信息缺 失，如融合图像中部颜色偏白组织两侧的黑色条状带的信息就缺失 了。本发明的融合图像，细节清晰，无畸变，融合后信息保存完整， 最优地表示了各待融合图像融合细节纹理信息。为客观评价融合效 果，这里选用以下几种评价指标来进行融合结果的客观评价。

A.信息熵(Information Entropy，IE)：熵值的大小表示图像所包含 的平均信息量的多少，通过对图像信息熵的比较可以得出图像的细节 表现能力。其中p_i为灰度值相等的像素数与总的像素数的比，熵的大 小反映了图像携带的信息量的多少，熵值越大，说明携带的信息量越 大。其计算公式为：

$\mathrm{IE}=-\underset{i=0}{\overset{L-1}{\Sigma }}{p}_i{\log }_2{p}_i$

p_i＝D_i/D

B.平均梯度(Average Gradient，AG)：它反映图像F对微小细节反 差和纹理变化特征表达能力的指标，也反映了图像的清晰度，平均梯 度越大，图像越清晰。其计算公式为：

$\mathrm{AG}=\frac{1}{(M-1)(N-1)}\underset{i=1}{\overset{(M-1)(N-1)}{\Sigma }}\sqrt{[{\left(\frac{\partial F}{\partial x}\right)}^2+{\left(\frac{\partial F}{\partial y}\right)}^2]/2}$

C.互信息(Mutual Information，MI)：源图像A、B和融合图像F间 的互信息。定义为：

$\mathrm{MI}[(A,B);F]=\underset{i=0}{\overset{L-1}{\Sigma }}\underset{j=0}{\overset{L-1}{\Sigma }}\underset{k=0}{\overset{L-1}{\Sigma }}{p}_{\mathrm{ABF}}(i,j,k)\ln \frac{p_{\mathrm{ABF}}(i,j,k)}{p_{\mathrm{AB}}(i,j)p_F\left(k\right)}$

其中p_ABF(i，j，k)为图像A、B、F的归一化灰度联合直方图，p_AB(i，j)为 为图像A、B的归一化灰度联合直方图。互信息量体现了融合图像从原 始图像中提取信息的多少，互信息量愈大，则提取的信息越多.

D.相关系数(Correlation Coefficient，CC)：衡量融合图像F和理 想图像R的相关程度。定义为：

$\mathrm{CC}=\frac{\left|\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}\right[R(i,j)-{\mu }_R\left]\right[F(i,j)-{\mu }_F\left]\right|}{\sqrt{\left(\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{[R(i,j)-{\mu }_R]}^2\right)\left(\underset{i=1}{\overset{m}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{[F(i,j)-{\mu }_F]}^2\right)}}$

μ_R、μ_F分别为F和R的均值。相关系数体现了融合图像F和标准参 考图像R间相似性，相关系数愈大，则图像的接近度越好。

E.扭曲程度(Degree of Distortion，DD)：扭曲程度直接反应融合 图像的失真程度，扭曲程度愈小，表示图像的失真程度愈小。

$\mathrm{DD}=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}|F(x,y)-R(x,y)|$

由表1数据可知，虽然本文提出的融合算法的IE小于小波融合算 法和Curvelet融合算法，但这是黑色的背景区域和局部细节畸变信息 引起的。但在总体上，客观指标优于其它算法。

综上所述，本文提出的融合算法，其图像融合综合效果最好。

表1融合图像评估参数

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而 非限制。尽管实施例对本发明进行了详细说明，本领域的普通技术人 员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离 本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围 当中。