CN103049895A - 基于平移不变剪切波变换的多模态医学图像融合方法 - Google Patents

Abstract

本发明公开了一种基于平移不变剪切波变换的多模态医学图像融合方法，包括以下步骤：1）准备待融合的两幅多模态医学图像，利用平移不变剪切波变换分别将两幅图像都分解为低频近似图像和高频细节图像，再将高频细节图像分解为不同的方向子带；2）分别融合两幅图像的低频近似图像和高频细节图像分解的各个方向子带；3）对步骤2）融合后的低频近似图像和高频细节图像分解的各个方向子带，利用逆平移剪切波变换得到融合图像。本发明采用了平移不变剪切波变换的图像融合方法，能方便、经济、高效率地实现多模态医学图像数据的融合，并能充分显示并捕捉到不同模态图像的内部隐含的细节部位的形态信息和功能信息，从而满足医学运用的精确性要求。

Description

基于平移不变剪切波变换的多模态医学图像融合方法

技术领域

本发明涉及一种图像融合方法，尤其是一种基于平移不变剪切波变换的多模态医学图像融合方法，属于医学图像处理及应用技术领域。

背景技术

多模态医学图像融合在临床诊断中扮演着重要的角色，它广泛应用于图像引导手术、图像引导放射治疗、非侵入性诊断以及治疗的计划的制定中。例如，在非侵入性诊断及治疗计划中，医生常常需要利用不同模态的医学图像来看清病人体内的病变结构和细节信息来准确地确定肿瘤的相对位置以及观察病变肿瘤对放射性射线的病理反应。因此，多模态医学图像融合技术是现代医学可视化技术的关键环节，广泛应用于现代医学临床诊断中。

目前，多模态医学图像融合方法大致可分为三种策略：替换方法、算术方法、多尺度分解法。替换方法，如彩色空间变换法，可能导致图像光谱的畸变；算术方法，如贝叶斯估计法，易导致图像对比度降低而丢失重要的细节信息。通常，图像中的轮廓、边缘等特征出现在不同大小的尺度上，人类的视觉系统研究表明物理图像经过视网膜处理变成神经图像是在不同的频率通道处理的，图像特征能在不同的尺度和频率上引起人类视觉系统的敏感。因此相比于上述两种方法，基于多尺度分解的多模态医学图像融合方法更符合人类视觉系统的工作原理。

基于多尺度分解的多模态医学图像融合方法的核心问题之一是多尺度分解工具的选择。在现有的发明中，如中国专利号为20081023245.3公开的基于多尺度几何分析的SAR图像融合方法、中国专利号为200910230339公开的医学超声的基波和谐波图像融合方法、中国专利号为201010300277公开的基于图像相似度的图像融合方法等，它们或选取小波（wavelet）或曲线波（curvelet）或轮廓波（contourlet）作为多尺度变换分解工具，但是它们都有其自身的缺点。例如，二维离散的小波变换只能将图像分解为水平、垂直、对角三个高频子带，这使其只能捕获有限的方向信息，这表明基于小波变换的融合方法由于小波不能很好的表示图像中的线、边缘、轮廓等高维信息而会导致伪影现象。

剪切波是为了克服传统小波的缺点而提出的一种新颖的多尺度几何分析工具，相比在上文提到的在图像融合中常用的其他多尺度分解工具，如曲线波（Curvelet）,轮廓波变换(Contourlet)等，剪切波变换具有其独特的优势：从近似理论的观点来看，剪切波形成一个具有的小波形式的多尺度多方向的紧支集结构，它是图像信号(如边缘)的真正的二维稀疏表示，目前只有曲线波具有类似的性质。然而，曲线波的实现过程过于复杂，并且不在图像的多分辨分析框架内实现，不利于快速算法的实现。相比轮廓波，虽然剪切波具有相似的实现过程，但是剪切波在剪切过程中没有方向数目和支撑基尺寸大小的限制。此外，剪切波逆变换只需对正向变换中的剪切滤波器进行加和处理，而不需要像轮廓波对方向滤波器逆合成，这说明剪切波的实现过程具有更高的计算效率。目前，剪切波变换已经被应用于图像融合中，如文献《Anovel algorithm ofimage fusion using shearlets》（Optics Communications,2010,284(6):1540–1547.）。然而，在上述文献中，剪切波变换的离散化过程采用了下采样策略实现，由于其不具备平移不变的性质而在图像中的奇异处易产生伪吉布斯现象。平移不变性对于图像奇异处的特征信息的提取至关重要，它严重影响融合后图像的质量。

基于多尺度分解的多模态医学图像融合方法的另一个核心问题是低频系数和高频系数的融合规则的选取。目前，在融合中最常用的融合规则是对低频系数取平均值，高频系数采用对应系数绝对值最大策略或者是基于区域能量最大策略。多模态医学图像融合的目的是尽可能的获取形态图像（如CT、MRI）的解剖细节信息和功能图像（如PET、SPECT）的新陈代谢的功能信息。此外，医学图像中的特征信息比并不只是图像像素的简单叠加，图像中的特征变化具有实际的诊断意义。对待融合图像进行多尺度分解后，低频系数表示图像的近似逼近，如果在医学图像的融合中仍旧对低频系数取均值就会增加原来不存在的形态信息和丢失一些已经存在的功能信息。因此在本文中采用基于区域系数绝对值和权重的低频系数融合策略。高频系数代表图像的细节信息，此外，基于人工神经网络的融合策略越来越受到学者们的重视，如中国专利号为20061002815.3公开的发明专利：基于Contourlet变换和改进型脉冲耦合神经网络的图像融合方法，它利用了脉冲耦合神经网络作为融合规则。但是这些应用中神经网络包含许多参数，并且需要专家依靠经验去设置，这限制了它们的使用范围。然而，传统的融合策略，包括基于对应系数绝对值最大策略或者基于区域能量最大策略的融合规则，它们的共同特点是所有的运算都是在每一个高频子带单独进行的。实际研究表明，不同尺度的高频子带间存在很强的相关性。因此，传统的融合规则没有充分考虑这种相关性而造成了不同模态图像的高频子带相关性的丢失。

发明内容

本发明的目的是为了解决上述现有技术的缺陷，提供一种可在普通硬件条件下实现的基于平移不变剪切波变换的多模态医学图像融合方法。

本发明的目的可以通过采取如下技术方案达到：

基于平移不变剪切波变换的多模态医学图像融合方法，其特征在于包括以下步骤：

1）准备待融合的两幅多模态医学图像，利用平移不变剪切波变换分别将两幅图像都分解为低频近似图像和高频细节图像，再将高频细节图像分解为不同的方向子带；

2）分别融合两幅图像的低频近似图像和高频细节图像分解的各个方向子带：

2.1）对低频近似图像，采用基于区域系数绝对值和权重的融合策略进行融合；

2.2）对高频细节图像各个方向子带，采用基于子带标准差和概率密度函数权重的融合规则进行融合：

a）构建高频子带系数的隐马尔科夫树HMT模型，利用HMT模型训练各子带系数，得到每个子带的标准差和概率密度函数；

b）采用基于子带标准差和概率密度函数权重的融合规则确定融合图像的各个高频子带系数。

3）对步骤2）融合后的低频近似图像和高频细节图像分解的各个方向子带，利用逆平移剪切波变换得到融合图像。

作为一种优选方案，所述待融合的两幅多模态医学图像若为彩色图像，在步骤1）对图像分解之前，先进行IHS变换，获得图像的Intensity分量，再对Intensity分量利用平移不变剪切波变换进行分解；在步骤3）中利用逆平移剪切波变换先得到Intensity分量，再进行逆IHS变换得到融合的彩色图像。

作为一种优选方案，所述步骤1）所述利用平移不变剪切波变换分解图像，具体如下：

1.1）利用非下采样金字塔策略将其中一幅多模态医学图像f^j分解成为低频近似图像f^j+1和高频细节图像g^j+1，其中j表示图像分解的尺度，j≥1；

1.2）针对高频细节图像g^j+1，构建Meyer窗口进行多尺度剖分：

a）在伪极化网格生成剪切滤波器窗口W；

b）将W从伪极化网格系统映射回到笛卡尔坐标系统，生成新的剪切滤波器W_new；

c）计算细节图像的傅里叶变换，产生矩阵Fg^j+1；

d）将W_new作用到Fg^j+1，获得方向子带；

1.3）对各个方向子带作用逆傅里叶变换得到平移不变剪切波系数；

1.4）对另一幅多模态医学图像采用上述步骤进行分解。

作为一种优选方案，步骤2.1）所述基于区域系数绝对值和权重的融合规则如下：

a）计算低频系数在3×3邻域的绝对值的和：

$S_f(i,j)=\underset{i,j}{\mathrm{\Σ}}\left|C_f(i,j)\right|,f=A,B---\left(1\right)$

其中，C_f(i,j)表示图像f在(i,j)处的低频系数；

b）计算系数的权重：

${\mathrm{\ω}}_A=\frac{S_A(i,j)}{S_A(i,j)+S_B(i,j)},$ ${\mathrm{\ω}}_B=\frac{S_B(i,j)}{S_A(i,j)+S_B(i,j)};---\left(2\right)$

c）计算融合图像在(i,j)处的低频系数：

C_F(i,j)=ω_A×C_A+ω_B×C_B    （3）

从而完成图像的低频系数的融合。

作为一种优选方案，所述步骤2.2）中的步骤a）具体如下：

首先，定义各尺度系数间的关系：对每一个系数X，定义与其相邻的8个子带系数称为X的邻域系数NX；定义与X在上一粗糙尺度相同空间位置的系数称为X的父系数PX；定义与X在下一加细尺度相同空间位置的系数为X的子系数CX；每个父系数有4子系数，不同尺度的子带形成四叉树结构；

然后，利用两状态的高斯混合模型对各个高频子带系数建模：令a^kl _mn表示父系数n处于状态l的时候子系数m处于状态k的概率，其中状态l,k=1，2；不同尺度系数间通过状态转移矩阵A_mn进行关联：

$A_{\mathrm{mn}}=\left\{\begin{array}{cc}{a^{11}}_{\mathrm{mn}}& {a^{12}}_{\mathrm{mn}}\\ {a^{21}}_{\mathrm{mn}}& {a^{22}}_{\mathrm{mn}}\end{array}\right\}---\left(4\right)$

再次，每个子带系数c_i通过其概率密度函数唯一确定：

$f_{c_i}\left(c_i\right)=\underset{s=1}{\overset{s=2}{\mathrm{\Σ}}}{p}_i\left(s\right)\×f\left(c_i\right|s_i=s)---\left(5\right)$

其中，

是系数c_i分别处于s=2的大状态、s=1的小状态时的概率密度函数，p_i(s)是c_i处于状态s时的概率，μ_i,s表示均值，σ_i,s表示标准差；

最后，利用最大期望算法训练模型，求解参数。

作为一种优选方案，所述步骤2.2）中的步骤b）具体如下：

首先，归一化高频子带系数C_μ ^l,k(i,j)：

${C_{\mathrm{\μ}}}^{l,k}(i,j)=\frac{{C_{\mathrm{\μ}}}^{l,k}(i,j)\×f_{\mathrm{\μ}}\left({C_{\mathrm{\μ}}}^{l,k}(i,j)\right)}{\left|f_{\mathrm{\μ}}\left({C_{\mathrm{\μ}}}^{l,k}(i,j)\right)\right|}---\left(6\right)$

然后，利用求得的概率密度函数定义子带系数光滑权重因子ω：

$\mathrm{\ω}=\left\{\begin{array}{cc}\frac{f_A\left({C_A}^{l,k}(i,j)\right)}{f_B\left({C_B}^{l,k}(i,j)\right)}& f_A\left({C_A}^{l,k}(i,j)\right)\≤f_B\left({C_B}^{l,k}(i,j)\right)\\ \frac{f_B\left({C_B}^{l,k}(i,j)\right)}{f_A\left({C_A}^{l,k}(i,j)\right)}& f_A\left({C_A}^{l,k}(i,j)\right)>f_B\left({C_B}^{l,k}(i,j)\right)\end{array}\right.---\left(7\right)$

最后，融合图像在子带位置(i,j)处的系数通过如下方式计算：

${C_F}^{l,k}(i,j)=\left\{\begin{array}{cc}\frac{{C_A}^{l,k}(i,j)\×{{\mathrm{\σ}}_A}^{l,k}(i,j)+{C_B}^{l,k}(i,j)\×{{\mathrm{\σ}}_B}^{l,k}(i,j)\×\mathrm{\ω}}{{{\mathrm{\σ}}_A}^{l,k}(i,j)+{{\mathrm{\σ}}_B}^{l,k}(i,j)\×\mathrm{\ω}}& \mathrm{\ω}\≤1\\ \frac{{C_A}^{l,k}(i,j)\×{{\mathrm{\σ}}_A}^{l,k}(i,j)\×\mathrm{\ω}+{C_B}^{l,k}(i,j)\×{{\mathrm{\σ}}_B}^{l,k}(i,j)}{{{\mathrm{\σ}}_A}^{l,k}(i,j)\×\mathrm{\ω}+{{\mathrm{\σ}}_B}^{l,k}(i,j)}& \mathrm{\ω}>1\end{array}\right.---\left(8\right)$

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

本发明相对于现有技术具有如下的有益效果：

1、本发明采用了平移不变剪切波变换的图像融合方法，相比传统的小波变换只能将待融合图像分解为低频近似图像和高频的水平、垂直和对角三个方向子带，可以将图像分解为更多方向子带，这为在融合图像中捕获并包含原始的待图像中更多的方向子带信息提供了可能。

2、本发明的图像融合方法采用的平移不变剪切波相比于其它常用的图像多尺度变换工具，如曲线波等，具有实现过程简单、计算效率高的优点，相比于普通的剪切波变换，能有效克服在图像奇异处由于缺乏平移不变性而易产生的伪吉布斯现象得缺点，因而，在融合图像中能捕获更多原始图像的特征信息，得到质量更高的融合图像。

3、本发明的图像融合方法对低频近似图像，采用基于区域绝对值和权重的融合规则，能避免图像对比度的降低；对高频子带，充分考虑了各个子带间的相关性，提出利用隐马尔科夫树模型来获取这种相关性，进而提出基于高频子带系数的概率密度函数和标准差权重的高频子带系数融合规则；可以有效克服传统的基于对应系数绝对值最大或者区域能量最大的融合规则子带相关系丢失的缺点。此外，相比传统的基于神经网络的融合规则需要人工设置许多参数的缺点，本发明的图像融合方法采用的高频子带融合规则仅需用计算机自动对各高频子带进行建模训练，成功避免了融合过程中过多的人工经验的干预，因此，融合结果更加真实和可靠。

4、本发明的图像融合方法能方便、经济、高效率地实现多模态医学图像数据的融合，并能充分显示并捕捉到不同模态图像的内部隐含的细节部位的形态信息和功能信息，从而满足医学运用的精确性要求。

附图说明

图1为本发明图像融合方法流程示意图。

图2a为传统的采用小波变换分解结果示意图

图2b为本发明采用平移不变剪切波分解结果示意图。

图3为本发明平移不变剪切波不同尺度系数间的关系图。

图4为本发明隐马尔科夫树模型中一个四叉树结构示意图。

图5a为实施例1待融合的MR图像。

图5b为实施例1待融合的SPECT图像。

图5c为传统的采用小波变换对图5a和图5b的融合结果示意图。

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

图6a为对图5a的MR图像的一个高频子带利用高斯混合模型拟合曲线图。

图6b为对图5b的SPECT图像的一个高频子带利用高斯混合模型拟合曲线图。

图7a为实施例2待融合的MR图像。

图7b为实施例2待融合的PET图像。

图7c为传统的采用小波变换对图7a和图7b的融合结果示意图。

图7d为本发明对图7a和图7b的融合结果示意图。

图8a为对图7a的MR图像的一个高频子带利用高斯混合模型拟合曲线示意图。

图8b为对图7b的PET图像的一个高频子带利用高斯混合模型拟合曲线示意图。

具体实施方式

本发明的实验平台为Matlab7.1.1，计算机主机配置：Intel（R）处理器，CPU主频3.0GHz,内存2.0GB。本实验平台用本发明方法对二维多模态医学图像数据进行处理，最终得到融合图像。

实施例1：MRI―SPECT融合

本实施例的图像融合方法首先考虑了如何从原始图像中获取更多的方向信息：高频子带的分解以及保持平移不变性；其次考虑了如何将提取的低频子带和高频子带信息尽可能多的转移到融合图像中：低频子带融合避免降低图像的对比度；高频子带融合避免各个子带间的相关信息的丢失和过多人工经验的干预。

如图1所示，包括如下步骤：

步骤一：对图5a和图5b所示的待融合的SPECT图像、MRI图像进行IHS变换，获得图像的Intensity分量，然后对其进行平移不变剪切波分解。执行过程包括两步：多尺度剖分和方向局部化。

1）多尺度剖分：

平移不变剪切波变换的多尺度剖分是通过非下采样的拉普拉斯金字塔滤波器组实现的，采用“9-7”滤波器组，每一级需要对上级中所采用的滤波器按矩阵 $D=2I=\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right]$ 进行上采样，图像经过k级非下采样塔式滤波，可以得到k+1个与源图像大小相同的子带图像，这样就避免了普通剪切波实现中下采样操作而丢失平移不变特性的缺点，在图像奇异处能够提取更多的特征信息，避免伪吉布斯现象的产生。

利用非下采样金字塔策略将SPECT图像、MRI图像的其中一幅图像f^j分解成为低频近似图像f^j+1和高频细节图像g^j+1，其中j表示图像分解的尺度，j≥1；

2）方向局部化：

平移不变剪切波变换的方向局部化通过剪切波滤波器实现，对多尺度剖分得到的第j+1个高频细节图像g^j+1，执行如下过程：

2.1）针对高频细节图像g^j+1，构建Meyer窗口进行多尺度剖分；

a）在伪极化网格生成剪切滤波器窗口W；

b）将W从伪极化网格系统映射回到笛卡尔坐标系统，生成新的剪切滤波器W_new；

c）计算细节图像的傅里叶变换，产生矩阵Fg^j+1；

d）将W_new作用到Fg^j+1，获得方向子带；

2.2）对各个方向子带作用逆傅里叶变换得到平移不变剪切波系数。

同理，对SPECT图像、MRI图像中的另一幅图像进行多尺度剖分和方向局部化。

采用平移不变剪切波变换分解的结果如图2b所示。而图2a是传统的采用小波变换分解结果。

步骤二：低频近似图像融合和高频图像各个方向子带融合：

1）对低频近似图像融合，采用基于区域系数绝对值和权重的融合规则：

1.1）计算低频系数在3×3邻域的绝对值的和：

f=A,B，其中，C_f(i,j)表示图像f在(i,j)处的低频系数；

1.2）计算系数的权重：

${\mathrm{\ω}}_A=\frac{S_A(i,j)}{S_A(i,j)+S_B(i,j)},$ ${\mathrm{\ω}}_B=\frac{S_B(i,j)}{S_A(i,j)+S_B(i,j)}$

1.3）计算融合图像在(i,j)处的低频系数：

C_F(i,j)=ω_A×C_A+ω_B×C_B

2）对高频细节图像各个方向子带，采用基于子带标准差和概率密度函数权重的融合规则进行融合：

2.1）对待融合图像的高频子带系数，构建其隐马尔科夫树（HMT）模型，并用最大期望（Expectation-Maximization，简称EM）算法求取其各个子带的标准差和概率密度函数，具体如下：

首先，定义各尺度系数间的关系：对每一个系数X，定义与其相邻的8个子带系数称为X的邻域系数NX；定义与X在上一粗糙尺度相同空间位置的系数称为X的父系数PX；定义与X在下一加细尺度相同空间位置的系数为X的子系数CX；每个父系数有4子系数，不同尺度的子带形成四叉树结构，如图3和图4所示；

然后，利用两状态的高斯混合模型，对MR图像的一个高频子带利用高斯混合模型拟合曲线如图6a所示，对SPECT图像的一个高频子带利用高斯混合模型拟合曲线如图6b所示，对各个高频子带系数建模：令a^kl _mn表示父系数n处于状态l的时候子系数m处于状态k的概率，其中状态l,k=1，2；不同尺度系数间通过状态转移矩阵A_mn进行关联：

$A_{\mathrm{mn}}=\left\{\begin{array}{cc}{a^{11}}_{\mathrm{mn}}& {a^{12}}_{\mathrm{mn}}\\ {a^{21}}_{\mathrm{mn}}& {a^{22}}_{\mathrm{mn}}\end{array}\right\}$

再次，每个子带系数c_i通过其概率密度函数

唯一确定：

$f_{c_i}\left(c_i\right)=\underset{s=1}{\overset{s=2}{\mathrm{\Σ}}}{p}_i\left(s\right)\×f\left(c_i\right|s_i=s)$

其中，

是系数c_i分别处于s=2的大状态、s=1的小状态时的概率密度函数，p_i(s)是c_i处于状态s时的概率，μ_i,s表示均值，σ_i,s表示标准差；

最后，利用最大期望算法训练模型，求解参数。

2.3）归一化高频子带系数C_μ ^l,k(i,j):

${C_{\mathrm{\μ}}}^{l,k}(i,j)=\frac{{C_{\mathrm{\μ}}}^{l,k}(i,j)\×f_{\mathrm{\μ}}\left({C_{\mathrm{\μ}}}^{l,k}(i,j)\right)}{\left|f_{\mathrm{\μ}}\left({C_{\mathrm{\μ}}}^{l,k}(i,j)\right)\right|}$

2.4）利用求得的概率密度函数定义子带系数光滑权重因子ω:

$\mathrm{\ω}=\left\{\begin{array}{cc}\frac{f_A\left({C_A}^{l,k}(i,j)\right)}{f_B\left({C_B}^{l,k}(i,j)\right)}& f_A\left({C_A}^{l,k}(i,j)\right)\≤f_B\left({C_B}^{l,k}(i,j)\right)\\ \frac{f_B\left({C_B}^{l,k}(i,j)\right)}{f_A\left({C_A}^{l,k}(i,j)\right)}& f_A\left({C_A}^{l,k}(i,j)\right)>f_B\left({C_B}^{l,k}(i,j)\right)\end{array}\right.$

2.5）融合图像在子带位置(i,j)处的系数通过如下方式计算:

${C_F}^{l,k}(i,j)=\left\{\begin{array}{cc}\frac{{C_A}^{l,k}(i,j)\×{{\mathrm{\σ}}_A}^{l,k}(i,j)+{C_B}^{l,k}(i,j)\×{{\mathrm{\σ}}_B}^{l,k}(i,j)\×\mathrm{\ω}}{{{\mathrm{\σ}}_A}^{l,k}(i,j)+{{\mathrm{\σ}}_B}^{l,k}(i,j)\×\mathrm{\ω}}& \mathrm{\ω}\≤1\\ \frac{{C_A}^{l,k}(i,j)\×{{\mathrm{\σ}}_A}^{l,k}(i,j)\×\mathrm{\ω}+{C_B}^{l,k}(i,j)\×{{\mathrm{\σ}}_B}^{l,k}(i,j)}{{{\mathrm{\σ}}_A}^{l,k}(i,j)\×\mathrm{\ω}+{{\mathrm{\σ}}_B}^{l,k}(i,j)}& \mathrm{\ω}>1\end{array}\right.$

步骤三：利用逆平移不变剪切波变换得到融合的intensity分量，利用逆IHS变得到融合的彩色图像，融合结果如图5d所示。而图5c为传统的利用小波变换融合MR图像和SPECT图像的融合结果。

实施例2：MRI―PET融合

用本发明提供的方法能在普通硬件条件下实现头部MRI图像与SPECT图像的高质量融合，MRI图像、PET图像原始分辨率为256×256，如图7a和图7b所示。

步骤一：对待融合的MRI图像、PET图像进行IHS变换，获得图像的Intensity分量，然后对其进行平移不变剪切波分解。执行过程包括两步：多尺度剖分和方向局部化。整个过程同实施例1。

步骤二：低频近似图像融合和高频图像各个方向子带融合：

1）对低频近似图像融合，采用基于区域系数绝对值和权重的融合规则；

2）对高频细节图像各个方向子带，采用基于子带标准差和概率密度函数权重的融合规则进行融合，本过程中，对MR图像的一个高频子带利用高斯混合模型拟合曲线如图8a所示，对PET图像的一个高频子带利用高斯混合模型拟合曲线如图8b所示。其余同实施例1。

步骤三：利用逆平移不变剪切波变换得到融合的Intensity分量，利用逆IHS变换得到融合的彩色图像，融合结果如图7d所示。而图7c为传统的利用小波变换融合MR图像和PET图像的融合结果。

以上所述，仅为本发明优选的实施例，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明所公开的范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都属于本发明的保护范围。

Claims (6)

1.基于平移不变剪切波变换的多模态医学图像融合方法，其特征在于包括以下步骤：

1）准备待融合的两幅多模态医学图像，利用平移不变剪切波变换分别将两幅图像都分解为低频近似图像和高频细节图像，再将高频细节图像分解为不同的方向子带；

2）分别融合两幅图像的低频近似图像和高频细节图像分解的各个方向子带：

2.1）对低频近似图像，采用基于区域系数绝对值和权重的融合策略进行融合；

2.2）对高频细节图像各个方向子带，采用基于子带标准差和概率密度函数权重的融合规则进行融合：

a）构建高频子带系数的隐马尔科夫树HMT模型，利用HMT模型训练各子带系数，得到每个子带的标准差和概率密度函数；

b）采用基于子带标准差和概率密度函数权重的融合规则确定融合图像的各个高频子带系数；

3）对步骤2）融合后的低频近似图像和高频细节图像分解的各个方向子带，利用逆平移剪切波变换得到融合图像。

2.根据权利要求1所述的基于平移不变剪切波变换的多模态医学图像融合方法，其特征在于：所述待融合的两幅多模态医学图像若为彩色图像，在步骤1）对图像分解之前，先进行IHS变换，获得图像的Intensity分量，再对Intensity分量利用平移不变剪切波变换进行分解；在步骤3）中利用逆平移剪切波变换先得到Intensity分量，再进行逆IHS变换得到融合的彩色图像。

3.根据权利要求1或2所述的基于平移不变剪切波变换的多模态医学图像融合方法，其特征在于：所述步骤1）所述利用平移不变剪切波变换分解图像，具体如下：

1.1）利用非下采样金字塔策略将其中一幅多模态医学图像f^j分解成为低频近似图像f^j+1和高频细节图像g^j+1，其中j表示图像分解的尺度，j≥1；

1.2）针对高频细节图像g^j+1，构建Meyer窗口进行多尺度剖分：

a）在伪极化网格生成剪切滤波器窗口W；

b）将W从伪极化网格系统映射回到笛卡尔坐标系统，生成新的剪切滤波器W_new；

c）计算细节图像的傅里叶变换，产生矩阵Fg^j+1；

d）将W_new作用到Fg^j+1，获得方向子带；

1.3）对各个方向子带作用逆傅里叶变换得到平移不变剪切波系数；

1.4）对另一幅多模态医学图像采用上述步骤进行分解。

4.根据权利要求1或2所述的基于平移不变剪切波变换的多模态医学图像融合方法，其特征在于：步骤2.1）所述基于区域系数绝对值和权重的融合规则如下：

a）计算低频系数在3×3邻域的绝对值的和：

$S_f(i,j)=\underset{i,j}{\mathrm{\Σ}}\left|C_f(i,j)\right|,f=A,B---\left(1\right)$

其中，C_f(i,j)表示图像f在(i,j)处的低频系数；

b）计算系数的权重：

${\mathrm{\ω}}_A=\frac{S_A(i,j)}{S_A(i,j)+S_B(i,j)},$ ${\mathrm{\ω}}_B=\frac{S_B(i,j)}{S_A(i,j)+S_B(i,j)}---\left(2\right)$

c）计算融合图像在(i,j)处的低频系数：

C_F(i,j)=ω_A×C_A+ω_B×C_B    （3）

从而完成图像的低频系数的融合。

5.根据权利要求1或2所述的基于平移不变剪切波变换的多模态医学图像融合方法，其特征在于：所述步骤2.2）中的步骤a）具体如下：

首先，定义各尺度系数间的关系：对每一个系数X，定义与其相邻的8个子带系数称为X的邻域系数NX；定义与X在上一粗糙尺度相同空间位置的系数称为X的父系数PX；定义与X在下一加细尺度相同空间位置的系数为X的子系数CX；每个父系数有4子系数，不同尺度的子带形成四叉树结构；

然后，利用两状态的高斯混合模型对各个高频子带系数建模：令a^kl _mn表示父系数n处于状态l的时候子系数m处于状态k的概率，其中状态l,k=1，2；不同尺度系数间通过状态转移矩阵A_mn进行关联：

$A_{\mathrm{mn}}=\left\{\begin{array}{cc}{a^{11}}_{\mathrm{mn}}& {a^{12}}_{\mathrm{mn}}\\ {a^{21}}_{\mathrm{mn}}& {a^{22}}_{\mathrm{mn}}\end{array}\right\}---\left(4\right)$

再次，每个子带系数c_i通过其概率密度函数

唯一确定：

$f_{c_i}\left(c_i\right)=\underset{s=1}{\overset{s=2}{\mathrm{\Σ}}}{p}_i\left(s\right)\×f\left(c_i\right|s_i=s)---\left(5\right)$

其中，

是系数c_i分别处于s=2的大状态、s=1的小状态时的概率密度函数，p_i(s)是c_i处于状态s时的概率，μ_i,s表示均值，σ_i,s表示标准差；

最后，利用最大期望算法训练模型，求解参数。

6.根据权利要求1或2所述的基于平移不变剪切波变换的多模态医学图像融合方法，其特征在于：所述步骤2.2）中的步骤b）具体如下：

首先，归一化高频子带系数C_μ ^l,k(i,j)：

${C_{\mathrm{\μ}}}^{l,k}(i,j)=\frac{{C_{\mathrm{\μ}}}^{l,k}(i,j)\×f_{\mathrm{\μ}}\left({C_{\mathrm{\μ}}}^{l,k}(i,j)\right)}{\left|f_{\mathrm{\μ}}\left({C_{\mathrm{\μ}}}^{l,k}(i,j)\right)\right|}---\left(6\right)$

然后，利用求得的概率密度函数定义子带系数光滑权重因子ω：

$\mathrm{\ω}=\left\{\begin{array}{cc}\frac{f_A\left({C_A}^{l,k}(i,j)\right)}{f_B\left({C_B}^{l,k}(i,j)\right)}& f_A\left({C_A}^{l,k}(i,j)\right)\≤f_B\left({C_B}^{l,k}(i,j)\right)\\ \frac{f_B\left({C_B}^{l,k}(i,j)\right)}{f_A\left({C_A}^{l,k}(i,j)\right)}& f_A\left({C_A}^{l,k}(i,j)\right)>f_B\left({C_B}^{l,k}(i,j)\right)\end{array}\right.---\left(7\right)$

最后，融合图像在子带位置(i,j)处的系数通过如下方式计算：

${C_F}^{l,k}(i,j)=\left\{\begin{array}{cc}\frac{{C_A}^{l,k}(i,j)\×{{\mathrm{\σ}}_A}^{l,k}(i,j)+{C_B}^{l,k}(i,j)\×{{\mathrm{\σ}}_B}^{l,k}(i,j)\×\mathrm{\ω}}{{{\mathrm{\σ}}_A}^{l,k}(i,j)+{{\mathrm{\σ}}_B}^{l,k}(i,j)\×\mathrm{\ω}}& \mathrm{\ω}\≤1\\ \frac{{C_A}^{l,k}(i,j)\×{{\mathrm{\σ}}_A}^{l,k}(i,j)\×\mathrm{\ω}+{C_B}^{l,k}(i,j)\×{{\mathrm{\σ}}_B}^{l,k}(i,j)}{{{\mathrm{\σ}}_A}^{l,k}(i,j)\×\mathrm{\ω}+{{\mathrm{\σ}}_B}^{l,k}(i,j)}& \mathrm{\ω}>1\end{array}\right.---\left(8\right)$

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

Images