CN103985105A - 基于统计建模的Contourlet域多模态医学图像融合方法 - Google Patents

Abstract

本发明公开了一种基于统计建模的Contourlet域多模态医学图像融合方法，主要解决医学图像融合时空间分辨率和光谱信息难以均衡的问题。其实现步骤是：1)待融合图像进行IHS变换，得到亮度、色调和饱和度；2)分别对亮度分量执行Contourlet变换，并采用EM算法估计高频子带的上下文隐马尔科夫模型CHMM参数；3)低频子带采用区域绝对值和取大的融合规则，高频子带基于CHMM和改进的脉冲耦合神经网络M-PCNN设计融合规则；4)融合后的高、低频系数执行Contourlet逆变换重构新的亮度分量；5)利用IHS逆变换获得融合图像。本发明能充分整合医学图像的结构和功能信息，有效保护图像细节，改善视觉效果，相比传统的融合方法极大地提高了融合图像的质量。

Description

基于统计建模的Contourlet域多模态医学图像融合方法

技术领域

本发明涉及一种基于统计建模的Contourlet域多模态医学图像融合方法，是医学图像处理技术领域的一项融合方法，在临床医学诊断和治疗中有广泛地应用。

背景技术

作为图像融合领域的一个研究分支和研究重点，随着医学影像技术的快速发展，多模态医学图像融合已成为当前国内外的研究热点。不同的医学影像可以提供人体相关脏器和组织的不同信息，如具有较高分辨率的CT和MRI图像提供脏器的解剖形态结构信息，而具有较低分辨率的SPECT和PET图像则提供脏器的功能代谢和血液流动信息。为了弥补单一模态图像信息的不足，为医务工作者提供患者各个器官更为清晰、可靠、准确、全面的描述，有必要通过图像融合方法将不同模态、信息互补的医学影像进行适当的集成，以便了解病变组织或器官的综合信息，从而做出准确的诊断或制订合适的治疗方案。多模态医学图像融合方法是通过处理冗余信息、整合互补信息来提高图像的分辨率，提升视觉效果，增强可识别性，增加图像的准确性，目前已被实际应用于病理分析(Pathologic Analysis)、临床诊断和治疗(Clinical Diagnoses and Therapy)中，且有广泛的研究前景和发展空间。

多模态医学图像融合方法也可大致分为空间域和变换域两类。空间域融合方法是对两幅图像在空间坐标下直接进行运算和叠加，简单直观，但往往不能满足医学应用精确性的要求。变换域融合方法有彩色空间变换法、基于IHS(Intensity-Hue-Saturation--亮度、色度、饱和度)的融合方法、基于主成份分析PCA(Principal Component Analysis)的融合方法、基于拉普拉斯变换LAP(Laplace Transform)的融合方法、基于离散小波变换DWT(Discrete Wavelet Transform)的融合方法以及基于多尺度变换(Multi-scale Decomposition)和脉冲耦合神经网络PCNN(PulseCoupled Neural Network)相结合的融合方法等。相比于彩色空间变换、PCA等方法易导致光谱失真，基于多尺度分解和PCNN相结合的方法能有效防止图像对比度的降低，保护图像细节特征，因而更适合用于多模态医学图像融合。

2002年，Do和Vetterli在继承小波多尺度分析思想的基础上提出了一种“真正”的二维图像稀疏表示方法-Contourlet变换，其基函数分布于多尺度、多方向上，仅用少量系数即可有效地捕获图像中的边缘轮廓及方向性纹理等细节信息，具有高度的方向敏感性和各向异性(Anisotropic)，很好地弥补了小波变换的不足，故本发明选用Contourlet作为多尺度变换工具。

Contourlet域分解系数在尺度间、方向间和空间邻域内存在极强的统计相关性，作为一种有效的统计建模技术，隐马尔科夫模型HMM(Hidden Markov Model)能够对图像多尺度变换所得系数的非高斯性(Non-Gaussian)、聚集性(Clustering)和持续性(Persistence)进行准确的描述。基于上下文的隐马尔可夫模型CHMM(Contextual HMM)，利用上下文信息(ContextInformation)来充分、有效地捕获Contourlet域高频方向子带系数在尺度间的持续性、尺度内的多方向选择性和空间邻域内的能量聚集特性。因此，本发明采用CHMM对医学图像高频方向子带进行统计建模。

为了提高融合图像的性能，融合规则的选取同样至关重要。通常，取平均、加权取平均或绝对值取大是比较常用的低频系数融合策略，本发明中选用基于区域绝对值和取大的融合规则，更好地保持了图像的对比度。而高频系数通常采用绝对值取大或基于区域特征取大的融合策略，但这些方案均未考虑子带系数不同尺度、不同方向之间潜在的相关性。为了充分提取待融合图像中的互补信息，避免有效细节的丢失，本发明提出改进的脉冲耦合神经网络M-PCNN(Modified PCN)并设计高频融合规则，并使用改进的拉普拉斯能量和SML(Sum-Modified-Laplacian)作为M-PCNN神经元的外部激励输入，CHMM训练参数计算边缘概率密度函数Edge PDF的归一化值作为其链接强度(Linking Strength)取值，获得了较优的融合结果，提高了融合图像的质量。

发明内容

本发明的目的是针对上述现有技术的不足，提出一种基于统计建模的Contourlet域多模态医学图像融合方法，以解决已有的多模态医学图像融合方法所得融合图像空间分辨率不高或光谱信息失真的问题，并充分整合不同模态医学图像的结构信息和功能信息，有效保护图像细节，增强图像对比度和边缘轮廓，改善其视觉效果，提高融合图像的质量。

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

一种基于统计建模的Contourlet域多模态医学图像融合方法，其特征在于包括以下步骤：

1)对待融合的两幅多模态医学图像进行IHS变换，得到相应的亮度、色调和饱和度分量；

2)对亮度分量执行Contourlet变换，分解得到不同尺度、不同方向的高、低频子带系数；

3)采用不同的融合规则分别确定融合图像的低频子带系数和高频子带系数；

3.1)对表征图像近似信息的低频子带，采用基于区域绝对值和取大的融合规则进行融合；

3.2)对表征图像细节特征的高频方向子带，采用基于CHMM和M-PCNN的融合规则进行融合；

a)构建高频方向子带系数的CHMM统计模型，利用优化的期望最大化EM算法分初始化和迭代训练两步来估计模型参数，得到每个系数的标准差和概率密度函数，并计算其边缘概率密度函数Edge PDF；

b)采用基于Edge PDF的归一化值作为链接强度，改进的拉普拉斯能量和SML作为外部激励输入的M-PCNN模型设计融合规则，确定各个高频子带融合系数；

4)对步骤3)所得高、低频子带融合系数，执行Contourlet逆变换重新获得新的亮度分量；

5)对新的亮度、色调和饱和度分量进行IHS逆变换获得最终的彩色融合图像。

作为一种优选方案，所述的基于统计建模的Contourlet域多模态医学图像融合方法，步骤2)所述对亮度分量执行Contourlet变换，分解得到不同尺度、不同方向的子带系数 (j=1,2...，J;k＝1，2，...，m_j；i，l＝1，...，N_j)，其中J表示分解尺度，k表示每尺度分解的方向子带，i,l表示位置索引，J表示最大分解尺度，m_j表示j尺度分解的方向子带最大数目，N_j表示相应子带的系数总个数，A、B分别代表SPECT和MRI图像；

作为一种优选方案，步骤3.1)所述基于区域绝对值和取大的融合规则，具体如下：

a)计算低频系数在3×3邻域内的绝对值之和；

$D_{j0}^X(i,l)=\underset{(i,l)\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left|C_{j0}^X(i,l)\right|,X=A,B$

其中，表示SPECT或MRI图像在最粗糙尺度j₀，位置(i,l)处的低频系数，Ω是3×3的窗口邻域；

b)确定融合图像在最粗糙尺度j₀，位置(i,l)处的低频子带系数：

$C_{j0}^F(i,l)=\left\{\begin{array}{c}{C}_{j0}^A(i,l),D_{j0}^A(i,l)\&GreaterEqual;D_{j0}^B(i,l)\\ C_{j0}^B(i,l),D_{j0}^A(i,l)<D_{j0}^B(i,l)\end{array}\right.$

作为一种优选方案，步骤3.2)中的步骤a)具体如下：

首先，定义轮廓波分解系数之间三种重要关系：对于当前Contourlet系数X，定义位于同一尺度同一方向子带内相邻位置上的8个Contourlet系数为其邻居系数(NX)，捕获空间相关性；定义相邻较粗糙尺度对应方向子带内相应空间位置上的Contourlet系数为其父系数(PX)，表征跨尺度相关性；定义位于同一尺度不同方向子带内相应空间位置上一组Contourlet系数为其表兄弟系数(CX)，体现不同方向间的相关性；

其次，根据新的设计方案计算相应的上下文变量取值，利用两状态、零均值的高斯混合模型GMM来刻画高频方向子带的非高斯分布特性，则每个Contourlet系数都与一个上下文变量和一个隐状态相关联，接着进行CHMM统计建模，由此可以计算其边缘概率密度函数：

$f_{C_{j,k,i,l}|V_{j,k,i,l}}\left(C_{j,k,i,l}\right|V_{j,k,i,l}=v)=\underset{m=0}{\overset{1}{\mathrm{\&Sigma;}}}{P}_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)\&CenterDot;g\left(C_{j,k,i,l}\right|0,{\mathrm{\&sigma;}}_{j,k,i,l}^2)$

$P_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)=\frac{P_{S_{j,k,i,l}}(S_{j,k,i,l}=m)\&CenterDot;P_{V_{j,k,i,l}|S_{j,k,i,l}}(V_{j,k,i,l}=v|S_{j,k,i,l}=m)}{\underset{m=0}{\overset{1}{\mathrm{\&Sigma;}}}{P}_{S_{j,k,i,l}}(S_{j,k,i,l}=m)\&CenterDot;P_{V_{j,k,i,l}|S_{j,k,i,l}}(V_{j,k,i,l}=v|S_{j,k,i,l}=m)}$

其中，(j,k,i,l)分别表示尺度、方向和空间位置索引，V_j,k,i,l为上下文变量，S_j,k,i,l为隐状态变量， $P_{S_{\text{j,k,i,l}}}(S_{j,k,i,l}=m)$ 是状态为m时的概率， $P_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)$ 是上下文变量取值v的条件下状态为m的概率， $g(C_{j,k,i,l,}0,{\mathrm{\&sigma;}}_{j,k,i,l,m}^2)=\frac{1}{\sqrt{{2\mathrm{\&pi;\&sigma;}}_{j,k,i,l,m}^2}}\&CenterDot;\exp \{-\frac{{(C_{j,k,i,l}-0)}^2}{2\&CenterDot;{\mathrm{\&sigma;}}_{j,k,i,l,m}^2}\}$ 表示高斯条件概率密度函数，均值为零，σ_j,k,i,l,m表示标准差；

最后，利用优化的期望最大化EM算法分初始化和迭代训练两步来估计模型参数。

作为一种优选方案，步骤3.2)中的步骤b)具体如下：

首先，初始化改进的脉冲耦合神经网络M-PCNN模型参数。所有神经元都处于熄火状态， $U_{i,l}^{j,k}\left(0\right)=0,Y_{i,l}^{j,k}\left(0\right)=0,{\mathrm{\&theta;}}_{i,l}^{j,k}\left(0\right)=1,T_{i,l}^{j,k}\left(0\right)=0,$ 设 $W=\left[\begin{array}{ccc}1/\sqrt{2}& 1& 1/\sqrt{2}\\ 1& 0& 1\\ 1/\sqrt{2}& 1& 1/\sqrt{2}\end{array}\right],$ α_L=1.0,α_θ=0.2,V_L=1.2,V_θ=20，迭代次数N_max=200，其中j为分解尺度，k为每尺度的方向子带，下标(i,l)为神经元位置索引；

再将边缘概率密度函数Edge PDF归一化至[0，1]区间，作为M-PCNN的链接强度取值；并按照下式计算改进的拉普拉斯能量和SML作为M-PCNN的外部激励输入；

${\&dtri;}_{\mathrm{ML}}^{j,k}(x,y)=|{2C}^{j,k}(x,y)-C^{j,k}(x-\mathrm{step},y)-C^{j,k}(x+\mathrm{step},y)|+|{2C}^{j,k}(x,y)-C^{j,k}(x,y-\mathrm{step})-C^{j,k}(x,y+\mathrm{step})|$

${\mathrm{SML}}_{i,l}^{j,k}=\underset{x=i-N}{\overset{i+N}{\mathrm{\&Sigma;}}}\underset{y=l-N}{\overset{l+N}{\mathrm{\&Sigma;}}}{\&dtri;}_{\mathrm{ML}}^{j,k}(x,y)$

其中，step表示系数间的可变间距，本发明中取值为1，N为区域窗口大小，取值为3。

然后，按照下列公式进行迭代运算，更新模型各个参数值，并统计每个神经元总的脉冲点火次数；

$\left\{\begin{array}{c}{F}_{i,l}^{j,k}\left(n\right)={\mathrm{SML}}_{i,l}^{j,k}\\ L_{i,l}^{j,k}\left(n\right)=\exp (-{\mathrm{\&alpha;}}_L)L_{i,l}^{j,k}(n-1)+V_L\underset{p,q}{\mathrm{\&Sigma;}}{W}_{i,l;p,q}{Y}_{p,q}(n-1)\\ U_{i,l}^{j,k}\left(n\right)=F_{i,l}^{j,k}\left(n\right)\&times;(1+\mathrm{\&beta;}\&CenterDot;L_{i,l}^{j,k}\left(n\right))\\ {\mathrm{\&theta;}}_{i,l}^{j,k}\left(n\right)=\exp (-{\mathrm{\&alpha;}}_{\mathrm{\&theta;}}){\mathrm{\&theta;}}_{i,l}^{j,k}(n-1)+V_{\mathrm{\&theta;}}{Y}_{i,l}^{j,k}\left(n\right)\\ Y_{i,l}^{j,k}\left(n\right)=\left\{\begin{array}{cc}1,& U_{i,l}^{j,k}\left(n\right)>{\mathrm{\&theta;}}_{i,l}^{j,k}\left(n\right)\\ 0,& \mathrm{otherwise}\end{array}\right.\end{array}\right.$

$T_{i,l}^{j,k}\left(n\right)=T_{i,l}^{j,k}(n-1)+Y_{i,l}^{j,k}\left(n\right)$

其中，下标(i,l)是神经元对应的图像像素标号，j,k是高频方向子带尺度、方向索引，p,q表示其它神经元与此神经元相连的范围，n为迭代次数；为神经元反馈输入，本实施例中采用能够较好描述图像的轮廓、方向性纹理等细节信息的SML作为M-PCNN的外部激励输入，为神经元链接输入，为神经元内部活动，为神经元输出，用于统计每个神经元总的脉冲点火次数，W_i,l;p,q是相应的突触联接权值矩阵，V_L是链接输入的放大系数，和V_θ分别是变阈值函数和阈值放大系数，α_L和α_θ分别表示链接输入和变阈值函数的时间常数，β为链接强度，其取值决定了神经元之间耦合关系的强弱，本实施例中选用高频方向子带系数的边缘概率密度函数Edge PDF的归一化值。

最后，根据如下公式，选取总的脉冲点火次数较大(即点火最频繁)的M-PCNN神经元所对应的系数作为融合高频方向子带系数：

$C_{j,k}^F(i,l)=\left\{\begin{array}{c}{C}_{j,k}^A(i,l),\mathrm{if}\left|T_{i,l}^{A;j,k}\left(N_{\max }\right)\right|\&GreaterEqual;\left|T_{i,l}^{B;j,k}\left(N_{\max }\right)\right|\\ C_{j,k}^B(i,l),\mathrm{if}\left|T_{i,l}^{A;j,k}\left(N_{\max }\right)\right|<\left|T_{i,l}^{B;j,k}\left(N_{\max }\right)\right|\end{array}\right.$

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

本发明相对比现有医学图像融合方法具有如下的优点：

1、本发明采用基于Contourlet域的图像融合方法，相比小波(Wavelet)变换仅能将图像分解为水平、垂直和对角三个高频子带，Contourlet变换在每一尺度都允许分解不同的方向数目，故其对方向信息的表现能力更强，能充分挖掘待融合图像的轮廓、纹理等细节特征，因而更适于处理图像的奇异性，得到信息量更丰富、清晰度更高、质量更好的融合图像。

2、本发明的医学图像融合方法采用Contourlet变换对图像进行多尺度多方向分解，相比于非采样Contourlet变换NSCT(Non-subsampled Contourlet Transform)，具有结构灵活、计算复杂度低的优点，因此，更能满足医学应用实时性的需求。

3、本发明的医学图像融合方法对低频近似子带系数，采用基于区域绝对值和取大的融合规则，能有效提高图像的亮度和对比度，改善整体视觉效果；对于高频方向子带系数，利用上下文隐马尔可夫模型来充分捕获Contourlet系数在尺度间、方向间和空间邻域内的统计相关性，进而提出采用基于边缘概率密度函数Edge PDF的归一化值作为链接强度，改进的拉普拉斯能量和SML作为外部激励输入的M-PCNN模型设计融合规则；可以最大限度地保护图像中的边缘轮廓和纹理信息，避免丢失细节，从而增加图像的可信度和可理解性。

附图说明

图1是本发明基于统计建模的Contourlet域多模态医学图像融合方法的流程图。

图2是本发明Contourlet域系数父子关系、最近邻关系和表兄弟关系示意图。

图3是本发明Contourlet域上下文隐马尔科夫模型C-CHMM的结构图。

图4是本发明改进的脉冲耦合神经网络M-PCNN的结构图。

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

图5(b)是本发明一个实施例的待融合SPECT图像。

图5(c)-(i)是本发明一个实施例的融合结果示意图。

图中：(c)基于IHS的融合图像；(d)基于PCA的融合图像；(e)基于DWT的融合图像；(f)基于IHS+Wavelet变换的融合图像；(g)基于IHS+Contourlet变换的融合图像；(h)基于IHS+NSCT的融合图像；(i)本发明方法的融合图像。

具体实施方式

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

步骤1，对待融合的两幅多模态医学图像进行IHS变换，得到相应的亮度、色调和饱和度分量；

步骤2，对亮度分量分别执行Contourlet变换，分解得到不同尺度、不同方向的高、低频子带系数(j=1,2...，J；k＝1，2，...，m_j；i，l＝１，...，N_j)，其中尺度分解LP选用“9-7”双正交滤波器，方向滤波器组DFB选用“pkva”，方向分解参数设为[2,2,3,3]，即进行4尺度分解，由粗糙到精细尺度所对应的方向子带数目依次为4,4,8,8；j表示分解尺度，k表示每尺度分解的方向子带，i,l表示位置索引，J表示最大分解尺度，m_j表示j尺度分解的方向子带最大数目，N_j表示相应子带的系数总个数，A、B分别代表SPECT和MRI图像；

步骤3，低频近似子带系数融合和高频方向子带系数融合：

1)对表征图像近似信息的低频子带，采用基于区域绝对值和取大的融合规则进行融合：

1.1)计算低频系数在3×3邻域内的绝对值之和；

$D_{j0}^X(i,l)=\underset{(i,l)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}\left|C_{j0}^X(i,l)\right|,X=A,B$

其中，表示SPECT或MRI图像在最粗糙尺度j₀，位置(i,l)处的低频系数，Ω是3×3的窗口邻域；

1.2)确定融合图像在最粗糙尺度j₀，位置(i,l)处的低频子带系数：

$C_{j0}^F(i,l)=\left\{\begin{array}{c}{C}_{j0}^A(i,l),D_{j0}^A(i,l)\&GreaterEqual;D_{j0}^B(i,l)\\ C_{j0}^B(i,l),D_{j0}^A(i,l)<D_{j0}^B(i,l)\end{array}\right.$

2)对表征图像细节特征的高频方向子带，采用基于CHMM和M-PCNN的融合规则进行融合；

2.1)构建高频方向子带系数的上下文隐马尔可夫统计模型，利用优化的期望最大化EM算法分初始化和迭代训练两步来估计CHMM参数，得到每个系数的标准差和概率密度函数，并计算其边缘概率密度函数Edge PDF，具体如下：

首先，定义轮廓波分解系数之间三种重要关系：对于当前Contourlet系数X，定义位于同一尺度同一方向子带内相邻位置上的8个Contourlet系数为其邻居系数(NX)，捕获空间相关性；定义相邻较粗糙尺度对应方向子带内相应空间位置上的Contourlet系数为其父系数(PX)，表征跨尺度相关性；定义位于同一尺度不同方向子带内相应空间位置上一组Contourlet系数为其表兄弟系数(CX)，体现不同方向间的相关性。将上述三种系数统称为广义邻居系数(Generalized Neighborhood Coefficients(GNC))，如图2所示；

其次，根据新的上下文设计方案，如图3所示，利用当前Contourlet系数的父节点、8个最近邻节点和两个最近表兄弟节点来计算其对应的V值，

$\mathrm{context}={\mathrm{\&omega;}}_0\&CenterDot;\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left|{\mathrm{NA}}_i\right|}^2+{\mathrm{\&omega;}}_1\&CenterDot;\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left|{\mathrm{NB}}_i\right|}^2+{\mathrm{\&omega;}}_2\&CenterDot;{\left|P\right|}^2+{\mathrm{\&omega;}}_3\&CenterDot;({\left|C_1\right|}^2+{\left|C_2\right|}^2)$

$V_{j,k,i,l}=\left\{\begin{array}{c}1,\mathrm{when}\mathrm{context}\&GreaterEqual;4\&CenterDot;{\mathrm{\&omega;}}_0\&CenterDot;E_n+4\&CenterDot;{\mathrm{\&omega;}}_1\&CenterDot;E_n\\ +{\mathrm{\&omega;}}_2\&CenterDot;E_p+{\mathrm{\&omega;}}_3\&CenterDot;(E_{C_1}+E_{C_2})\\ 0,\mathrm{otherwise}\end{array}\right.$

其中， $E_x=\frac{1}{N}\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(\mathrm{Coefficients}\left(i\right)\right)}^2,(E_x=E_n,E_p,E_{C_1},E_{C_2})$ 分别表示当前子带、较粗糙尺度同方向子带和同一尺度两个相邻方向子带系数的平均能量，N为对应子带系数总个数；ω₀、ω₁分别表示最近邻系数中的直接邻居NA(direct neighbors)和间接邻居NB(diagonal andcounter-diagonal neighbors)所对应的权值因子，ω₂、ω₃分别表示父系数和两个最近表兄弟系数所对应的权值因子。ω_i(i＝0,1,2,3)反映了GNC对应的相关性在上下文设计方案中所占比重。通过对多组经验值比较，本实施例中权值因子最终选取ω₀=1.0、ω₁=0.4、ω₂=1.0、ω₃=0.6。

再次，利用两状态、零均值的高斯混合模型GMM来刻画高频方向子带的非高斯分布特性，则每个Contourlet系数都与一个上下文变量和一个隐状态相关联，接着进行CHMM统计建模，由此可以计算其边缘概率密度函数：

$f_{C_{j,k,i,l}|V_{j,k,i,l}}\left(C_{j,k,i,l}\right|V_{j,k,i,l}=v)=\underset{m=0}{\overset{1}{\mathrm{\&Sigma;}}}{P}_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)\&CenterDot;g\left(C_{j,k,i,l}\right|0,{\mathrm{\&sigma;}}_{j,k,i,l}^2)$

$P_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)=\frac{P_{S_{j,k,i,l}}(S_{j,k,i,l}=m)\&CenterDot;P_{V_{j,k,i,l}|S_{j,k,i,l}}(V_{j,k,i,l}=v|S_{j,k,i,l}=m)}{\underset{m=0}{\overset{1}{\mathrm{\&Sigma;}}}{P}_{S_{j,k,i,l}}(S_{j,k,i,l}=m)\&CenterDot;P_{V_{j,k,i,l}|S_{j,k,i,l}}(V_{j,k,i,l}=v|S_{j,k,i,l}=m)}$

其中，(j,k,i,l)分别表示尺度、方向和空间位置索引，V_j,k,i,l为上下文变量，S_j,k,i,l为隐状态变量， $P_{S_{\text{j,k,i,l}}}(S_{j,k,i,l}=m)$ 是状态为m时的概率， $P_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)$ 是上下文变量取值v的条件下状态为m的概率， $g(C_{j,k,i,l,}0,{\mathrm{\&sigma;}}_{j,k,i,l,m}^2)=\frac{1}{\sqrt{{2\mathrm{\&pi;\&sigma;}}_{j,k,i,l,m}^2}}\&CenterDot;\exp \{-\frac{{(C_{j,k,i,l}-0)}^2}{2\&CenterDot;{\mathrm{\&sigma;}}_{j,k,i,l,m}^2}\}$ 表示高斯条件概率密度函数，均值为零，σ_j,k,i,l,m表示标准差；

最后，利用优化的期望最大化EM算法分初始化和迭代训练两步来估计模型参数，其中EM训练算法中由粗糙到精细各尺度对应的局部窗口大小依次设为5×5、7×7、9×9、11×11，迭代次数分别设为20次和5次。

2.2)采用基于Edge PDF的归一化值作为链接强度，改进的拉普拉斯能量和SML作为外部激励输入的M-PCNN模型设计融合规则，确定各个高频子带融合系数；

首先，初始化改进的脉冲耦合神经网络M-PCNN模型参数。所有神经元都处于熄火状态， $U_{i,l}^{j,k}\left(0\right)=0,Y_{i,l}^{j,k}\left(0\right)=0,{\mathrm{\&theta;}}_{i,l}^{j,k}\left(0\right)=1,T_{i,l}^{j,k}\left(0\right)=0,$ 设 $W=\left[\begin{array}{ccc}1/\sqrt{2}& 1& 1/\sqrt{2}\\ 1& 0& 1\\ 1/\sqrt{2}& 1& 1/\sqrt{2}\end{array}\right],$ α_L=1.0,α_θ=0.2,V_L=1.2,V_θ=20，迭代次数N_max=200，其中j为分解尺度，k为每尺度的方向子带，下标(i,l)为神经元位置索引；再将边缘概率密度函数Edge PDF归一化至[0，1]区间，作为M-PCNN的链接强度取值；并按照下式计算改进的拉普拉斯能量和SML作为M-PCNN的外部激励输入；

${\&dtri;}_{\mathrm{ML}}^{j,k}(x,y)=|{2C}^{j,k}(x,y)-C^{j,k}(x-\mathrm{step},y)-C^{j,k}(x+\mathrm{step},y)|+|{2C}^{j,k}(x,y)-C^{j,k}(x,y-\mathrm{step})-C^{j,k}(x,y+\mathrm{step})|$

${\mathrm{SML}}_{i,l}^{j,k}=\underset{x=i-N}{\overset{i+N}{\mathrm{\&Sigma;}}}\underset{y=l-N}{\overset{l+N}{\mathrm{\&Sigma;}}}{\&dtri;}_{\mathrm{ML}}^{j,k}(x,y)$

其中，step表示系数间的可变间距，本发明中取值为1，N为区域窗口大小，取值为3。

然后，M-PCNN中每个神经元都包括接受域、调制耦合器和脉冲产生器三部分，如图4所示。按照其描述公式进行迭代运算，更新各参数值，并统计每个神经元总的脉冲点火次数；

$\left\{\begin{array}{c}{F}_{i,l}^{j,k}\left(n\right)={\mathrm{SML}}_{i,l}^{j,k}\\ L_{i,l}^{j,k}\left(n\right)=\exp (-{\mathrm{\&alpha;}}_L)L_{i,l}^{j,k}(n-1)+V_L\underset{p,q}{\mathrm{\&Sigma;}}{W}_{i,l;p,q}{Y}_{p,q}(n-1)\\ U_{i,l}^{j,k}\left(n\right)=F_{i,l}^{j,k}\left(n\right)\&times;(1+\mathrm{\&beta;}\&CenterDot;L_{i,l}^{j,k}\left(n\right))\\ {\mathrm{\&theta;}}_{i,l}^{j,k}\left(n\right)=\exp (-{\mathrm{\&alpha;}}_{\mathrm{\&theta;}}){\mathrm{\&theta;}}_{i,l}^{j,k}(n-1)+V_{\mathrm{\&theta;}}{Y}_{i,l}^{j,k}\left(n\right)\\ Y_{i,l}^{j,k}\left(n\right)=\left\{\begin{array}{cc}1,& U_{i,l}^{j,k}\left(n\right)>{\mathrm{\&theta;}}_{i,l}^{j,k}\left(n\right)\\ 0,& \mathrm{otherwise}\end{array}\right.\end{array}\right.$

$T_{i,l}^{j,k}\left(n\right)=T_{i,l}^{j,k}(n-1)+Y_{i,l}^{j,k}\left(n\right)$

其中，下标(i,l)是神经元对应的图像像素标号，j,k是高频方向子带尺度、方向索引，p,q表示其它神经元与此神经元相连的范围，n为迭代次数；为神经元反馈输入，本实施例中采用能够较好描述图像的轮廓、方向性纹理等细节信息的SML作为M-PCNN的外部激励输入，为神经元链接输入，为神经元内部活动，为神经元输出，用于统计每个神经元总的脉冲点火次数，W_i,l;p,q是相应的突触联接权值矩阵，V_L是链接输入的放大系数，和V_θ分别是变阈值函数和阈值放大系数，α_L和α_θ分别表示链接输入和变阈值函数的时间常数，β为链接强度，其取值决定了神经元之间耦合关系的强弱，本实施例中选用高频方向子带系数的边缘概率密度函数Edge PDF的归一化值。M-PCNN模型在调制域进行耦合调制，计算内部活动输入到脉冲产生部分用于产生神经元的脉冲输出脉冲产生器中，若则对应神经元被激活，产生一个脉冲输出，称为一次点火，同时神经元的阈值通过反馈迅速提高，继续进行下一次迭代运算；当超过时，脉冲产生器关闭，神经元停止发放脉冲，阈值开始指数下降，当再次低于时，脉冲产生器再次被打开，对应神经元重新处于激活状态，如此迭代循环直至达到设定的最大迭代次数。

最后，根据如下公式，选取总的脉冲点火次数较大(即点火最频繁)的M-PCNN神经元所对应的系数作为融合高频方向子带系数：

$C_{j,k}^F(i,l)=\left\{\begin{array}{c}{C}_{j,k}^A(i,l),\mathrm{if}\left|T_{i,l}^{A;j,k}\left(N_{\max }\right)\right|\&GreaterEqual;\left|T_{i,l}^{B;j,k}\left(N_{\max }\right)\right|\\ C_{j,k}^B(i,l),\mathrm{if}\left|T_{i,l}^{A;j,k}\left(N_{\max }\right)\right|<\left|T_{i,l}^{B;j,k}\left(N_{\max }\right)\right|\end{array}\right.$

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

步骤4，对高、低频子带融合系数执行Contourlet逆变换重新获得新的亮度分量；

步骤5，对新的亮度、色调和饱和度分量进行IHS逆变换获得最终的彩色融合图像。

仿真实验

为了验证本发明的可行性和有效性，采用MRI和SPECT两幅医学图像，大小为256×256，如图5(a)和5(b)所示，根据本发明方法进行融合实验。

综上所述，通过图5的融合结果比较可以看出：本发明方法所得融合图像最大程度地忠于原始信息，更好地保持了待融合图像中的边缘轮廓、纹理等特征，有效避免了光谱信息失真，因而图像的对比度和清晰度更高，细节更突出，主观视觉效果最好，即融合结果更理想。

表1给出了采用各种融合方法所得融合结果的客观评价指标。其中，数据加粗表明对应的多模态医学图像融合方法所得评价指标取值最优。

表1基于各种融合方法的融合性能比较

表1中通过标准差、熵、清晰度、边缘强度、Qw、Qabf和光谱扭曲度来衡量融合图像的质量，进而验证本发明融合方法的可行性和有效性。上表中融合方法一为基于IHS的融合方法，融合方法二为基于主成分分析PCA的融合方法，融合方法三为基于离散Wavelet变换DWT融合方法，融合方法四为基于IHS和Wavelet变换的融合方法，融合方法五为基于IHS和Contourlet变换的融合方法，融合方法六为基于IHS和NSCT的融合方法，融合方法三、四、五、六均采用低频子带取平均、高频子带绝对值取大的融合规则。

由表1数据可知，本发明方法所获得的融合图像在标准差、熵、清晰度、边缘强度、光谱扭曲度等客观评价指标上要优于其它融合方法。标准差反映图像灰度相对于灰度均值的离散情况，其值越大，灰度级越分散，图像反差越大，可看出更多的信息；熵值反映了图像携带信息量的多少，熵值越大，说明包含的信息量越多，融合效果越好；清晰度反映图像对微小细节反差表达的能力，清晰度越高则图像融合效果越好；边缘强度用于衡量图像边缘细节的丰富程度，其值越大，则融合图像的边缘越清晰，效果越好；Qw指标取值越接近1，表明图像融合效果越好；而Qabf利用Sobel边缘检测算子来衡量有多少边缘信息从源图像转移到了融合图像，可综合评价源图像与融合图像之间的结构相似度，其值越接近1，说明图像融合效果越好。光谱扭曲度直接反映了融合影像的光谱失真程度，其值越小，表明融合效果越好。

Claims (7)

1.一种基于统计建模的Contourlet域多模态医学图像融合方法，其特征是：首先对待融合的多模态医学图像进行IHS变换，对得到的亮度分量执行Contourlet变换，并对高频子带系数建立上下文隐马尔可夫CHMM统计模型，然后分别采用区域绝对值和取大以及基于CHMM和改进的脉冲耦合神经网络M-PCNN的融合规则确定低、高频子带融合系数，再执行Contourlet逆变换重构新的亮度分量，最后进行IHS逆变换获得最终的彩色融合图像。

2.根据权利要求1所述的基于统计建模的Contourlet域多模态医学图像融合方法，其特征在于包括以下具体步骤：

1)对待融合的两幅多模态医学图像进行IHS变换，得到相应的亮度、色调和饱和度分量；

2)对亮度分量执行Contourlet变换，分解得到不同尺度、不同方向的高、低频子带系数；

3)采用不同的融合规则分别确定融合图像的低频子带系数和高频子带系数；

3.1)对表征图像近似信息的低频子带，采用基于区域绝对值和取大的融合规则进行融合；

3.2)对表征图像细节特征的高频方向子带，采用基于CHMM和M-PCNN的融合规则进行融合；

a)构建高频方向子带系数的CHMM统计模型，利用优化的期望最大化EM算法分初始化和迭代训练两步来估计模型参数，得到每个系数的标准差和概率密度函数，并计算其边缘概率密度函数Edge PDF；

b)采用基于Edge PDF的归一化值作为链接强度，改进的拉普拉斯能量和SML作为外部激励输入的M-PCNN模型设计融合规则，确定各个高频子带融合系数；

4)对步骤3)所得高、低频子带融合系数，执行Contourlet逆变换重新获得新的亮度分量；

5)对新的亮度、色调和饱和度分量进行IHS逆变换获得最终的彩色融合图像。

3.根据权利要求1或2所述的基于统计建模的Contourlet域多模态医学图像融合方法，其特征在于：步骤2)所述对亮度分量执行Contourlet变换，分解得到不同尺度、不同方向的子带系数(j=1,2...，J;k＝1，2，...，m_j；i，l＝1，...，N_j)，其中j表示分解尺度，k表示每尺度分解的方向子带，i,l表示位置索引，J表示最大分解尺度，m_j表示j尺度分解的方向子带最大数目，N_j表示相应子带的系数总个数，A、B分别代表SPECT和MRI图像。

4.根据权利要求1或2所述的基于统计建模的Contourlet域多模态医学图像融合方法，其特征在于：步骤3.1)所述基于区域绝对值和取大的融合规则，具体如下：

a)计算低频系数在3×3邻域内的绝对值之和；

$D_{j0}^X(i,l)=\underset{(i,l)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}\left|C_{j0}^X(i,l)\right|,X=A,B;$

其中，表示SPECT或MRI图像在最粗糙尺度j₀，位置(i,l)处的低频系数，Ω是3×3的窗口邻域；

b)确定融合图像在最粗糙尺度j₀，位置(i,l)处的低频子带系数；

$C_{j0}^F(i,l)=\left\{\begin{array}{c}{C}_{j0}^A(i,l),D_{j0}^A(i,l)\&GreaterEqual;D_{j0}^B(i,l)\\ C_{j0}^B(i,l),D_{j0}^A(i,l)<D_{j0}^B(i,l)\end{array}\right..$

5.根据权利要求1或2所述的基于统计建模的Contourlet域多模态医学图像融合方法，其特征在于：步骤3.2)中的步骤a)具体如下：

首先，定义轮廓波分解系数之间三种重要关系：对于当前Contourlet系数X，定义位于同一尺度同一方向子带内相邻位置上的8个Contourlet系数为其邻居系数(NX)，捕获空间相关性；定义相邻较粗糙尺度对应方向子带内相应空间位置上的Contourlet系数为其父系数(PX)，表征跨尺度相关性；定义位于同一尺度不同方向子带内相应空间位置上一组Contourlet系数为其表兄弟系数(CX)，体现不同方向间的相关性；

其次，根据新的设计方案计算相应的上下文变量取值，利用两状态、零均值的高斯混合模型GMM来刻画高频方向子带系数的非高斯分布特性，则每个Contourlet系数都与一个上下文变量和一个隐状态相关联，接着进行CHMM统计建模，由此可以计算其边缘概率密度函数Edge PDF；

$f_{C_{j,k,i,l}|V_{j,k,i,l}}\left(C_{j,k,i,l}\right|V_{j,k,i,l}=v)=\underset{m=0}{\overset{1}{\mathrm{\&Sigma;}}}{P}_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)\&CenterDot;g\left(C_{j,k,i,l}\right|0,{\mathrm{\&sigma;}}_{j,k,i,l}^2);$

$P_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)=\frac{P_{S_{j,k,i,l}}(S_{j,k,i,l}=m)\&CenterDot;P_{V_{j,k,i,l}|S_{j,k,i,l}}(V_{j,k,i,l}=v|S_{j,k,i,l}=m)}{\underset{m=0}{\overset{1}{\mathrm{\&Sigma;}}}{P}_{S_{j,k,i,l}}(S_{j,k,i,l}=m)\&CenterDot;P_{V_{j,k,i,l}|S_{j,k,i,l}}(V_{j,k,i,l}=v|S_{j,k,i,l}=m)};$

其中，(j,k,i,l)分别表示尺度、方向和空间位置索引，V_j,k,i,l为上下文变量，S_j,k,i,l为隐状态变量， $P_{S_{\text{j,k,i,l}}}(S_{j,k,i,l}=m)$ 是状态为m时的概率， $P_{S_{j,k,i,l}|V_{j,k,i,l}}(S_{j,k,i,l}=m|V_{j,k,i,l}=v)$ 是上下文变量取值v的条件下状态为m的概率， $g(C_{j,k,i,l,}0,{\mathrm{\&sigma;}}_{j,k,i,l,m}^2)=\frac{1}{\sqrt{{2\mathrm{\&pi;\&sigma;}}_{j,k,i,l,m}^2}}\&CenterDot;\exp \{-\frac{{(C_{j,k,i,l}-0)}^2}{2\&CenterDot;{\mathrm{\&sigma;}}_{j,k,i,l,m}^2}\}$ 表示高斯条件概率密度函数，均值为零，σ_j,k,i,l,m表示标准差；

最后，利用优化的期望最大化EM算法分初始化和迭代训练两步来估计模型参数。

6.根据权利要求1或2所述的基于统计建模的Contourlet域多模态医学图像融合方法，其特征在于：步骤3.2)中的步骤b)具体如下：

首先，初始化改进的脉冲耦合神经网络M-PCNN模型参数。所有神经元都处于熄火状态， $U_{i,l}^{j,k}\left(0\right)=0,Y_{i,l}^{j,k}\left(0\right)=0,{\mathrm{\&theta;}}_{i,l}^{j,k}\left(0\right)=1,T_{i,l}^{j,k}\left(0\right)=0,$ 设 $W=\left[\begin{array}{ccc}1/\sqrt{2}& 1& 1/\sqrt{2}\\ 1& 0& 1\\ 1/\sqrt{2}& 1& 1/\sqrt{2}\end{array}\right],$ α_L=1.0,α_θ=0.2,V_L=1.2,V_θ=20，迭代次数N_max=200，其中j为分解尺度，k为每尺度的方向子带，下标(i,l)为神经元位置索引；

再将边缘概率密度函数Edge PDF归一化至[0，1]区间，作为M-PCNN的链接强度取值；并按照下式计算改进的拉普拉斯能量和SML作为M-PCNN的外部激励输入；

${\&dtri;}_{\mathrm{ML}}^{j,k}(x,y)=|{2C}^{j,k}(x,y)-C^{j,k}(x-\mathrm{step},y)-C^{j,k}(x+\mathrm{step},y)|+|{2C}^{j,k}(x,y)-C^{j,k}(x,y-\mathrm{step})-C^{j,k}(x,y+\mathrm{step})|;$

${\mathrm{SML}}_{i,l}^{j,k}=\underset{x=i-N}{\overset{i+N}{\mathrm{\&Sigma;}}}\underset{y=l-N}{\overset{l+N}{\mathrm{\&Sigma;}}}{\&dtri;}_{\mathrm{ML}}^{j,k}(x,y);$

其中，step表示系数间的可变空间距离，本发明中取值为1，N为区域窗口大小，取值为3；

然后，按照下列公式进行迭代运算，更新参数值，并统计每个神经元总的脉冲点火次数；

$\left\{\begin{array}{c}{F}_{i,l}^{j,k}\left(n\right)={\mathrm{SML}}_{i,l}^{j,k}\\ L_{i,l}^{j,k}\left(n\right)=\exp (-{\mathrm{\&alpha;}}_L)L_{i,l}^{j,k}(n-1)+V_L\underset{p,q}{\mathrm{\&Sigma;}}{W}_{i,l;p,q}{Y}_{p,q}(n-1)\\ U_{i,l}^{j,k}\left(n\right)=F_{i,l}^{j,k}\left(n\right)\&times;(1+\mathrm{\&beta;}\&CenterDot;L_{i,l}^{j,k}\left(n\right))\\ {\mathrm{\&theta;}}_{i,l}^{j,k}\left(n\right)=\exp (-{\mathrm{\&alpha;}}_{\mathrm{\&theta;}}){\mathrm{\&theta;}}_{i,l}^{j,k}(n-1)+V_{\mathrm{\&theta;}}{Y}_{i,l}^{j,k}\left(n\right)\\ Y_{i,l}^{j,k}\left(n\right)=\left\{\begin{array}{cc}1,& U_{i,l}^{j,k}\left(n\right)>{\mathrm{\&theta;}}_{i,l}^{j,k}\left(n\right)\\ 0,& \mathrm{otherwise}\end{array}\right.\end{array}\right.$

$T_{i,l}^{j,k}\left(n\right)=T_{i,l}^{j,k}(n-1)+Y_{i,l}^{j,k}\left(n\right);$

其中，下标(i,l)是神经元对应的图像像素标号，j,k是高频方向子带尺度、方向索引，p,q表示其它神经元与此神经元相连的范围，n为迭代次数；为神经元反馈输入，本实施例中采用能够较好描述图像的轮廓、方向性纹理等细节信息的SML作为M-PCNN的外部激励输入，为神经元链接输入，为神经元内部活动，为神经元输出，用于统计每个神经元总的脉冲点火次数，W_i,l;p,q是相应的突触联接权值矩阵，V_L是链接输入的放大系数，和V_θ分别是变阈值函数和阈值放大系数，α_L和α_θ分别表示链接输入和变阈值函数的时间常数，β为链接强度，其取值决定了神经元之间耦合关系的强弱，本实施例中选用高频方向子带系数的边缘概率密度函数Edge PDF的归一化值；

最后，根据如下公式，选取总的脉冲点火次数较大(即点火最频繁)的M-PCNN神经元所对应的系数作为融合高频方向子带系数：

$C_{j,k}^F(i,l)=\left\{\begin{array}{c}{C}_{j,k}^A(i,l),\mathrm{if}\left|T_{i,l}^{A;j,k}\left(N_{\max }\right)\right|\&GreaterEqual;\left|T_{i,l}^{B;j,k}\left(N_{\max }\right)\right|\\ C_{j,k}^B(i,l),\mathrm{if}\left|T_{i,l}^{A;j,k}\left(N_{\max }\right)\right|<\left|T_{i,l}^{B;j,k}\left(N_{\max }\right)\right|\end{array}\right.;$

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

7.根据权利要求2或6所述的基于统计建模的Contourlet域多模态医学图像融合方法，其特征在于：上述采用M-PCNN模型的融合规则是依据若干次迭代后M-PCNN各个神经元产生的脉冲点火总次数决定高频方向子带融合系数的选取；脉冲点火次数反映了神经元所受外部刺激的大小，表征了Contourlet分解系数中所包含有效细节信息的多少，因此选择总的脉冲点火次数较大(即点火最频繁)的神经元所对应的高频方向子带系数可以包含待融合图像更多的细节特征，进而能获得更优的融合图像；M-PCNN模型在调制域进行耦合调制，计算内部活动输入到脉冲产生部分用于产生神经元的脉冲输出

在脉冲产生部分，若则对应神经元被激活，产生一个脉冲输出，称为一次点火，同时神经元的阈值通过反馈迅速提高，继续进行下一次迭代运算；当超过时，脉冲产生器关闭，神经元停止发放脉冲，阈值开始指数下降，当再次低于时，脉冲产生器再次被打开，对应的神经元重新处于激活状态，如此迭代循环直至达到设定的最大迭代次数。