CN102129676B - 一种基于二维经验模态分解的显微图像融合方法 - Google Patents

Abstract

本发明涉及一种基于二维经验模态分解的显微图像融合方法，该方法采用二维经验模态分解方法对采集的序列显微源图像进行多尺度分解，获得源图像的多级尺度的高频分量，按照局部显著性准则进行融合处理，并对源图像的低频分量采用主成分分析方法进行融合处理，最后反向重构获取融合图像。本发明采用二维经验模态分解方法对采集的序列显微图像进行多尺度分解，分解过程是自适应的；按照基于区域极大值的局部显著性准则进行高频融合处理，考虑了相邻系数的相关性，可以提取出各幅源图像的聚焦清晰的细节信息；采用主成分分析方法进行低频融合处理，能够很好的利用源图像像素的相关信息，提高融合图像的目视解译效果，使得融合后图像的质量得到提高。

Description

一种基于二维经验模态分解的显微图像融合方法

技术领域

本发明涉及图像融合的技术，特别涉及一种基于二维经验模态分解的显微图像融合方法，用以融合显微成像过程中聚焦位置不同的多幅图像，由局部清晰的序列图像提取得到整幅清晰的图像，特别适合红外光源条件下的显微成像。

背景技术

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

至今为止，人们已研究出了多种方法用于融合多聚焦图像，可主要分为空域融合方法和多分辨分析方法两大类。

各类空域融合方法计算速度快，但在处理离焦程度较大的显微图像时，容易产生块效应，效果不理想。

以小波变换为代表的多分辨分析方法具有良好的平滑性能，目前得到广泛的应用和认可，但这种方法有两个问题：一是所抽取的特征只能反映图像在单一函数经过平移和伸缩所构成函数上的能量；二是需要对小波基和一些参数进行选取。

1998年美国学者Norden E.Huang等人提出了一种用于分析非线性非平稳数据的经验模态分解(Empirical Mode Decomposition，EMD)方法，它是基于数据时域局部特征、自适应的时频分析工具。二维经验模态分解(Bidimensional Empirical ModeDecomposition，BEMD)是一维EMD分解思想与方法在二维平面上的推广，可用于图像数据的分析和处理，通过将原始图像自适应地分解为有限数量的子图像，可以将图像从高频到低频的局部窄带的各个细节清晰地分解出来，残差表示图像趋势。提取来的各级固有模态函数(Intrinsic Mode Function，IMF)都具有当前图像中局部最高的空间振荡频率，也就是当前图像的纹理特征。这种新的图像多尺度分析方法，其与传统多尺度分析技术区别主要在于：首先，极值点距离被引入进行局部尺度的确定从而使分解具有自适应和完全数据驱动的特性；其次，每一成分的提取都使用迭代计算方法，并使用某种准则确定迭代的终止。由于这些原因，BEMD在自适应的提取图像符合视觉感知的成分上有其独特的优势。但作为一种新兴的技术，BEMD方法在理论上还不尽完善，其应用更处于起步阶段，且很少有关于图像融合的讨论和具体应用的尝试，需深入研究探索。

主成分分析(PCA)是统计学中分析数据的一种有效的方法，其主要目的是降维，并且在降维后保存了原数据中的主要信息，从而使数据更易于处理。但是，在非多尺度分解的框架下，基于PCA分解确定加权系数的图像融合方法对显微图像聚焦细节信息的提取能力有限。但若能将该方法融入多尺度分解算法，必将取得更好的融合图像，值得在这方面继续加以研究。

综上，目前的现有融合技术应用于序列显微图像融合，特别是红外光源条件下的显微图像的融合时，还存在一些不足：各类空域融合方法存在块效应问题；小波变换方法只提取水平、垂直和对角线三个方向细节信息且需要对小波基等参数进行选择，并且参数选择的好坏对融合结果有很大的影响。

发明内容

本发明的目的在于，针对现有技术的不足，本发明立足提出一种显微图像融合方法，以期达到既不产生块效应现象，又无需人为进行参数的选择，并且能够更好地提取显微图像的聚焦信息的目标，实现全二维的自适应显微图像融合，从而提供一种基于二维经验模态分解的显微图像融合方法，其基于全新的多尺度分解结构，具有完全数据驱动的自适应性，具有更强的细节获取能力，并且反映显微图像独特的视觉意义，运用主成分分析方法和局部显著性准则进行分解后各成分的处理，使得融合后图像的质量得到提高。

对于普通显微镜，本发明可以有效扩大显微镜的景深，特别地，可以应用于使用红外光源的显微观察领域，如对具有趋光性的微小生物体的显微镜观察，还可以用于生物学、病理学、药物化学、材料检测等领域。

因为很多需要显微镜观察的微小生物体具有趋光性，此时显微镜的可见光光源就会对实验产生无法预料的影响而得不到预期的结果。使用红外光源代替普通光源，可以有效的避免光源对实验结果的影响。植物、细胞生长普遍具有趋光性，很多植物、细胞和微生物等的生长实验必须在避光的条件下进行，采用可见光进行实时观察，会对植物、细胞和微生物等的行为和效应产生影响，从而影响实验结果。大量的文献报道：生命体对红外波段的敏感性低，采用红外光源进行观察，对植物、细胞和微生物等的行为和效应产生的影响非常小，是一种适宜进行生命实验的观测手段。但是，使用红外光源代替普通光源，由于光信号较弱，照度较低，显微成像的图像质量受到一定程度的影响，亮度、对比度较低，需要更好的显微图像融合方法。

在显微光学成像过程中，显微镜物镜放大倍数越高，景深越小，只有那些在聚焦平面或其附近的结构才是可见的。而生物医学的发展要求显微镜成像既要有更高的分辨率又要有足够的景深，这是传统光学硬件的矛盾。解决上述问题的一个有效方法就是采集覆盖显微样本全部纵向信息的序列图像，利用数字图像处理技术，通过一定的融合规则对序列图像进行合成，从而重建出一幅每一景深部位均十分清晰的图像。

为实现上述发明的目的，本发明提供了一种基于二维经验模态分解的显微图像融合方法，该方法采用二维经验模态分解方法对采集的序列显微源图像进行多尺度分解，获得源图像的多级尺度的高频分量，按照局部显著性准则进行融合处理，并对源图像的低频分量采用主成分分析方法进行融合处理，最后反向重构获取融合图像，该方法包括如下步骤：

(1)对采集的序列显微源图像X₁，X₂，...，X_n分别进行二维经验模态分解处理，二维经验模态分解BEMD的过程反映了提取局部最高频、次高频的过程，得到每幅源图像的n级固有模态函数分量IMF和一个残差分量；

(2)对不同源图像对应的各级固有模态函数分量IMF中的像素，采用基于区域极大值的局部显著性选择准则进行融合处理，使得将融合后的固有模态函数分量具备所有源图像的清晰聚焦的细节信息；

(3)利用主成分分析PCA方法，分别计算出不同源图像对应的残差分量的自适应融合权重，按权重进行残差分量的融合处理；

(4)将融合后的各级固有模态函数分量和残差分量反向重构获取融合图像。

作为上述技术方案的一种改进，所述的步骤(1)中对每一幅源图像进行图像高频到低频的自然尺度分离，首先，分解出来的第1级固有模态函数IMF₁是图像所含有的最高频率分量，该分量的各处频率都对应着图像在各处的局部最高频，源图像减去第1级固有模态函数得到第1级残差分量；对第1级残差分量再进行分解，得到第2级固有模态函数和第2级残差分量；依此类推，得到n级固有模态函数和第n级残差；

所述的二维经验模态分解BEMD的处理过程包括如下步骤：

(1-1)：为了避免二维经验模态分解产生边界效应，采用局部镜像延拓对原始图像进行四周边界处理：I_原图＝F_镜像(I_原图)；

(1-2)：初始化：I＝I_原图，I_残差＝I，j＝0，j表示IMF的分解级数；

(1-3)：对所处理的残差图像曲面I_残差求取曲面局部极值点，包括所有局部极大值和极小值，初始时，I_残差就是源图像曲面I_原图；

(1-4)：对各极大值点和各极小值点分别进行曲面拟合，经插值后得到极大值点对应的上包络曲面E_u和极小值点对应的下包络曲面E_l；

(1-5)：将两曲面数据求平均得到均值包络曲面数据E_m：E_m＝(E_u+E_l)/2；

(1-6)：计算筛分终止条件标准偏差SD：

(1-7)：提取细节，用残差图像曲面减去均值包络曲面求得差值：I_残差＝I_残差-E_m；

(1-8)：重复上述步骤(1-3)～步骤(1-7)，直到满足给定的终止条件：(a)IMF的极值点和过零点数目必须相等或至多只相差一点；(b)在每一像素点，由极大值点定义的上包络线和由极小值点定义的下包络线的平均值为零；

(1-9)：计算残差，用图像I减去第j层固有模态函数(即I_残差)得到第j层残差分量并赋值给I：I＝I-I_残差；

(1-10)：对残差分量I重复步骤(1-3)～步骤(1-9)，直到满足残差不含IMF分量或已达到所需要的运算级数，依次得到图像的n级固有模态函数和第n级残差分量。

作为上述技术方案进一步的改进，所述的步骤(1-3)中，图像的局部极大值点为灰度值比周围3×3区域8个相邻像素点灰度值都高的点，图像的极小值点为灰度值比周围8个相邻像素点灰度值都低的点。

作为上述技术方案进一步的改进，所述的步骤(1-4)基于Delaunay三角剖分的插值方法对各极大值点和各极小值点分别进行曲面拟合；计算上包络曲面E_u时，将各极大值点分割成若干简单的小三角形区域，三角形的顶点即是图像的局部极大值点，在每个小三角形中立方插值构造插值曲面，再进一步把这些曲面片拼接起来，构造出一张大的插值曲面，即上包络曲面E_u；同样，用此方法再构造出极小值点对应的下包络曲面E_l。

作为上述技术方案进一步的改进，所述的步骤(1-8)采用简化的终止条件：当步骤(1-6)中计算得到的SD小于0.3时，就认为满足了收敛条件，令j＝j+1，此时得到的I_残差即为第j层二维固有模态函数IMF_j。

作为上述技术方案的另一种改进，所述的步骤(2)包括：

(2-1)首先，对于显微源图像X₁，为其各级IMF系数定义一个衡量其显著性的变量S：

$S_j(X_1,p)=\underset{q\∈Q}{\max }\left(\underset{q\∈Q}{\mathrm{\Σ}}\right|{\mathrm{IMF}}_j(X_1,q)\left|\right);$

其中，j表示IMF系数的级数，j＝1表示该分量为第一级IMF分量，反映显微图像X₁的局部最高频信息，j＝2表示该分量为第二级IMF分量，反映显微图像X1的局部次高频信息，以此类推；p＝(m，n)表示IMF系数的空间位置；Q表示以p为中心的一个3×3的方形窗口，q为窗口内的任意一点；IMF_j表示图像X₁位于q点位置的第j级IMF系数值；

对于图像X_i中对应的IMF系数同样定义S_j(X_i，p)；

(2-2)然后，在不同源图像的各级IMF系数中选择对应的显著性变量S值较大的IMF系数作为合成图像中对应位置的IMF系数；如果用C_j ^*(X_i，p)表示图像X_i相应位置上的决策表的值，则可表示为：

$C_j^{*}(X_1,p)=\left\{\begin{array}{cc}1,& S_j(X_1,p)\≥S_j(X_i,q),i=\mathrm{2,3},...n\\ 0,& \mathrm{other}\end{array}\right.$

$C_j^{*}(X_2,p)=\left\{\begin{array}{cc}1,& S_j(X_2,p)\≥S_j(X_i,p),i=1,3,...n\\ 0,& \mathrm{other}\end{array}\right.;$

…

$C_j^{*}(X_n,p)=\left\{\begin{array}{cc}1,& S_j(X_n,p)\≥S_j(X_i,p),i=1,2,...n-1\\ 0,& \mathrm{other}\end{array}\right.$

对得到的决策表进行基于多数表决原则的一致性验证，令修正后决策表的值为C_j(X_i，p)，则：

$C_j(X_1,p)=\left\{\begin{array}{cc}1,& \underset{q\∈Q}{\mathrm{\Σ}}{C}_j^{*}(X_1,q)\≥5\\ 0,& \mathrm{other}\end{array}\right.$

$C_j(X_2,p)=\left\{\begin{array}{cc}1,& \underset{q\∈Q}{\mathrm{\Σ}}{C}_j^{*}(X_2,q)\≥5\\ 0,& \mathrm{other}\end{array}\right.;$

…

$C_j(X_n,p)=\left\{\begin{array}{cc}1,& \underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{C}_j(X_i,p)=0\\ 0,& \mathrm{other}\end{array}\right.$

(2-3)得到决策表中各点的值后，就可计算融合图像F的IMF系数：

${\mathrm{IMF}}_j(F,p)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_j(X_i,p)\·{\mathrm{IMF}}_j(X_i,p).$

作为上述技术方案的再一改进，所述的步骤(3)对n幅显微源图像，把每幅图像看作一维向量记做x_i，i＝1，2，...，n，残差分量的融合处理包括：

(3-1)由源图像构造数据矩阵X：X＝(x₁，x₂，…，x_n)^T；

(3-2)计算数据矩阵X的协方差矩阵C：

σ_i，j ²为图像的方差，

x_i为第i幅源图像的平均灰度值；

$X={(x_1,x_2,...,x_n)}^T=\left(\begin{array}{ccccc}{x}_{11}& ...& x_{1j}& ...& x_{1m}\\ ...& ...& ...& ...& ...\\ x_{i1}& ...& x_{\mathrm{ij}}& ...& x_{\mathrm{im}}\\ ...& ...& ...& ...& ...\\ x_{n1}& ...& x_{\mathrm{nj}}& ...& x_{\mathrm{nm}}\end{array}\right);$ $C=\left(\begin{array}{ccccc}{\mathrm{\σ}}_{11}& ...& {\mathrm{\σ}}_{1j}& ...& {\mathrm{\σ}}_{1m}\\ ...& ...& ...& ...& ...\\ {\mathrm{\σ}}_{i1}& ...& {\mathrm{\σ}}_{\mathrm{ij}}& ...& {\mathrm{\σ}}_{\mathrm{im}}\\ ...& ...& ...& ...& ...\\ {\mathrm{\σ}}_{n1}& ...& {\mathrm{\σ}}_{\mathrm{nj}}& ...& {\mathrm{\σ}}_{\mathrm{nm}}\end{array}\right);$

(3-3)计算协方差矩阵的C的特征值λ及相应的特征向量ξ：

由特征值方程|λI-C|＝0，求出特征值λ_i和对应的特征向量ξ_i(i＝1，2，...，n)；

(3-4)确定加权系数ω_i：

(3-5)计算最终融合残差分量I_残差：

作为上述技术方案的又一改进，所述的步骤(4)反向重构获取融合图像I_融合结果：

本发明的方法首先对采集的序列显微图像进行二维经验模态分解处理，之后对不同图像对应的代表图像高频细节信息的各级固有模态函数分量IMFs中的像素按照局部显著性准则进行选取并进行融合处理，使得将融合后的固有模态函数分量具备所有源图像的清晰聚焦的细节信息。再利用主成分分析(PCA)方法确定不同图像对应的残差分量的自适应融合权重，按权重进行残差分量的融合处理，以提高融合图像的目视解译效果。最后将融合后的各级固有模态函数分量和残差分量反向重构获取融合图像。

本发明的方法包括如下具体步骤：

1、对采集的序列显微图像进行二维经验模态分解处理，BEMD分解的过程反映了提取局部最高频、次高频等的过程，得到n级IMF分量和一个残差分量；

2、对不同源图像对应的各级固有模态函数分量IMFs中的像素，本发明采用基于区域极大值的局部显著性选择准则进行融合处理，使得将融合后的固有模态函数分量具备所有源图像的清晰聚焦的细节信息。这是因为，在一幅显微图像的二维经验模态分解结果中，各级固有模态函数分量代表图像的高频细节信息，而且绝对值较大的IMF系数对应于图像中对比度变化较大的边缘等特征，表明此幅显微图像在该区域聚焦清晰。

基于区域极大值的局部显著性选择规则对于代表高频率域的各级IMF分量的融合非常适合。这是因为，二维经验模态分解后得到的IMF分量的系数在零值左右波动，绝对值越大表示该处灰度变化越剧烈，越有可能包含图像边缘等细节信息，即表明聚焦清晰程度越高。不同聚焦程度显微图像的各级IMF分量存在显著差异，IMF分量融合规则的选择对于最终的融合图像的质量是至关重要的。对于在红外光源条件下所成的显微图像而言，由于光信号较弱，照度较低，图像质量没有普通光源好，亮度、对比度较低。本发明考虑到显微图像特别是红外光源显微图像的特点，在处理各级IMF分量时，不仅考虑到当前像素点位置的IMF系数，还考虑了与它相邻位置的IMF系数。这种基于区域极大值的局部显著性选择规则体现了图像像素与它相邻像素的相关性，能更好的避免融合后图像的块效应现象。另外，为了使IMF系数在融合时对一个区域内的点采取相同的选择方案，本发明还对得到的决策表进行了基于多数表决原则的一致性验证。

3、利用主成分分析(PCA)方法，分别计算出不同源图像对应的残差分量的自适应融合权重，按权重进行残差分量的融合处理。这种方法很好的结合了源图像像素相关信息，计算的融合系数权使融合图像灰度分布方差最大化，因此能使融合图像的反差更好，细节更清晰，可以提高融合图像的目视解译效果。

4、将融合后的固有模态函数分量和残差分量反向重构获取融合图像。

本发明首先对多幅显微图像进行BEMD分解，获得每幅图像的多尺度信息；对各级固有模态函数按照局部显著性准则进行融合处理，保留了所有源图像的清晰聚焦的细节信息；利用PCA方法确定残差分量的自适应融合权重并进行融合，使得源图像像素的相关信息得到充分利用；将融合后的各级分量反向重构获取融合图像。本发明采用自适应的BEMD分解，完全由数据驱动，不需要预设任何滤波器，且是全二维处理，与只提取三个方向细节信息的小波方法相比，能获得更为丰富的显微聚焦信息，使得融合后图像的质量得到提高。

本发明的优点在于，本发明提出的基于二维经验模态分解的显微图像融合方法，采用二维经验模态分解方法对采集的序列显微图像进行多尺度分解，分解过程是自适应的，完全由数据驱动，不需要预设任何滤波器或小波函数；求取各级IMF分量类似高频滤波过程，获得各级尺度的高频部分，按照基于区域极大值的局部显著性准则进行融合处理，考虑了相邻系数的相关性，可以提取出各幅源图像的聚焦清晰的细节信息，而且BEMD方法是一种全二维的处理方法，与只提取水平、垂直和对角线三个方向的细节信息的小波方法相比，能获得更为丰富的显微聚焦信息；求取残差分量类似低频滤波过程，对应源图像的低频部分，采用主成分分析方法进行融合处理，能够很好的利用源图像像素的相关信息，提高融合图像的目视解译效果。本发明使得融合后图像的质量得到提高，对于多聚焦红外显微图像融合领域的应用具有重要意义和实用价值。

附图说明

图1是本发明的基于二维经验模态分解的显微图像融合方法流程示意图；

图2是BEMD分解流程框图；

图3是BEMD多尺度分解示意图；

图4是BEMD对Cameraman标准图像三级分解结果示例；

其中，(a)为待处理Cameraman标准图像，(b)为第1级IMF子图像，(c)为第2级IMF子图像，(d)为第3级IMF子图像，(e)为残差子图像。

图5是局部显著性融合规则示意图；

图6是本发明方法使用一种藻类显微序列图像作为融合源图像的融合结果示例；

其中，(a)、(b)和(c)为三幅源图像，(d)为本发明方法的融合结果图像。

具体实施方式

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

如图1所示，本发明首先对源图像X₁，X₂，...，X_n分别进行二维经验模态分解，得到每幅源图像的n级IMF分量和一个残差分量；对各级IMF分量按照局部显著性准则进行融合处理；对残差分量利用PCA方法进行融合处理；将融合后的各级IMF分量和残差分量进行逆向重构得到最终融合图像。

本发明的具体实施如下：

1、对采集的多幅显微图像分别进行二维经验模态分解处理

BEMD分解的过程反映了提取局部最高频、次高频等的过程，得到n级IMF分量和一个第n层残差分量。对一幅图像进行BEMD分解的实现过程如下(如图2所示)：

步骤1：为了避免二维经验模态分解产生边界效应，采用局部镜像延拓对原始图像进行四周边界处理：I_原图＝F_镜像(I_原图)。

步骤2：初始化：I＝I_原图，I_残差＝I，j＝0(j表示IMF的分解级数)。

步骤3：对所处理的残差图像曲面I_残差(初始时就是源图像曲面I_原图)求取曲面局部极值点，包括所有局部极大值和极小值。

图像的局部极大值点就是灰度值比周围3×3区域8个相邻像素点灰度值都高的点，灰度值比周围8个相邻像素点灰度值都低的点就是图像的极小值点，其他都不是极值点。逐行扫描I_残差的每一个像素点，将像素点I_残差(i，j)的灰度值与其邻域的8个像素的灰度值相比较，确定出图像的极大值点和极小值点。

步骤4：对各极大值点和各极小值点分别进行曲面拟合，经插值后得到极大值点对应的上包络曲面E_u和极小值点对应的下包络曲面E_l。

在BEMD的实现过程中，极值点找出来后，它们是散乱分布的，需要把它们组织起来，在空间上进行曲面拟合。包络曲面的拟合是实现BEMD的关键步骤，本发明采用目前常用的基于Delaunay三角剖分的插值方法。计算上包络曲面E_u时，将各极大值点分割成若干简单的小三角形区域，三角形的顶点即是图像的局部极大值点，在每个小三角形中立方插值构造插值曲面，再进一步把这些曲面片拼接起来，构造出一张大的插值曲面，即上包络曲面E_u。同样，用此方法再构造出极小值点对应的下包络曲面E_l。

步骤5：将两曲面数据求平均得到均值包络曲面数据E_m：E_m＝(E_u+E_l)/2。

步骤6：计算筛分终止条件：

这里的SD(标准偏差，standard deviation)就是均值包络曲面数据的绝对值的最大值除以I_残差的绝对值的最大值得到的结果。

步骤7：提取细节，用残差图像曲面减去均值包络曲面求得差值：I_残差＝I_残差-E_m。

步骤8：重复步骤3～步骤7，直到满足给定的终止条件。理论上，各级IMF需要满足下列两个条件：(1)IMF的极值点和过零点数目必须相等或至多只相差一点；(2)在每一像素点，由极大值点定义的上包络线和由极小值点定义的下包络线的平均值为零。为了减少计算量，使收敛速度较快些，本发明选择了简化的终止条件，即步骤6中计算得到的SD小于0.3时，就认为满足了收敛条件，令j＝j+1，此时得到的I_残差即为第j层二维固有模态函数IMF_j。

步骤9：计算残差，用图像I减去第j层固有模态函数(即I_残差)得到第j层残差分量并赋给I：I＝I-I_残差；

步骤10：对残差分量(即I)重复步骤3～步骤9，直到满足给定的终止条件，即残差不含IMF分量或已达到所需要的运算级数，依次得到图像的n级固有模态函数和第n级残差分量。

BEMD的分解过程实现了图像高频到低频的自然尺度分离过程。首先分解出来的第1级固有模态函数IMF₁是图像所含有的最高频率分量，该分量的各处频率都对应着图像在各处的局部最高频。原图像减去第1级固有模态函数得到第1级残差分量；对第1级残差分量再进行分解，得到第2级固有模态函数和第2级残差分量；依此类推，得到n级固有模态函数和第n级残差。每一级固有模态函数都由上一级残差分解得到，如图3所示BEMD多尺度分解示意图。对Cameraman标准图像做3层BEMD分解，得到的3级IMF和最后一级残差如图4所示。BEMD没有确定的基，它的“基”是根据信号而自适应产生的，这使得它具有良好的时频局部性。

2、对多幅源图像的各级固有模态函数IMF系数进行融合处理

在一幅显微图像的二维经验模态分解中，绝对值较大的固有模态函数IMF系数对应于图像中对比度变化较大的边缘等特征，表明此幅显微图像在该区域聚焦清晰。因此，对不同源图像对应的各级固有模态函数IMFs中的像素进行融合处理，本发明采用基于区域极大值的局部显著性选择准则。

设有多幅显微源图像X₁，X₂，...，X_n，对于显微图像X₁，可为其各级IMF系数定义一个衡量其显著性的变量S：

$S_j(X_1,p)=\underset{q\∈Q}{\max }\left(\underset{q\∈Q}{\mathrm{\Σ}}\right|{\mathrm{IMF}}_j(X_1,q)\left|\right);$

其中，j表示IMF系数的级数，j＝1表示该分量为第一级IMF分量，反映显微图像X₁的局部最高频信息，j＝2表示该分量为第二级IMF分量，反映显微图像X₁的局部次高频信息，以此类推；p＝(m，n)表示IMF系数的空间位置；Q表示以p为中心的一个3×3的方形窗口，q为窗口内的任意一点；IMF_j表示图像X₁位于q点位置的第j级IMF系数值。

对于图像X_i中对应的IMF系数同样可定义S_j(X_i，p)。为了在最后的融合图像中保留不同源图像中最显著的特征，本发明在不同源图像的IMF系数中选择值较大的IMF系数作为合成图像中对应位置的IMF系数。如果用C_j ^*(X_i，p)表示图像X_i相应位置上的决策表的值，上述思想用数学公式就可表示为：

$C_j^{*}(X_1,p)=\left\{\begin{array}{cc}1,& S_j(X_1,p)\≥S_j(X_i,q),i=\mathrm{2,3},...n\\ 0,& \mathrm{other}\end{array}\right.$

$C_j^{*}(X_2,p)=\left\{\begin{array}{cc}1,& S_j(X_2,p)\≥S_j(X_i,p),i=1,3,...n\\ 0,& \mathrm{other}\end{array}\right.;$

…

$C_j^{*}(X_n,p)=\left\{\begin{array}{cc}1,& S_j(X_n,p)\≥S_j(X_i,p),i=1,2,...n-1\\ 0,& \mathrm{other}\end{array}\right.$

为了在选择IMF系数时对一个区域内的点采取相同的选择方案，本发明对得到的决策表还进一步进行一致性验证。在这里，本发明采用多数表决原则。令修正后决策表的值为C_j(X_i，p)，则：

$C_j(X_1,p)=\left\{\begin{array}{cc}1,& \underset{q\∈Q}{\mathrm{\Σ}}{C}_j^{*}(X_1,q)\≥5\\ 0,& \mathrm{other}\end{array}\right.$

$C_j(X_2,p)=\left\{\begin{array}{cc}1,& \underset{q\∈Q}{\mathrm{\Σ}}{C}_j^{*}(X_2,q)\≥5\\ 0,& \mathrm{other}\end{array}\right.$

…

$C_j(X_n,p)=\left\{\begin{array}{cc}1,& \underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{C}_j(X_i,p)=0\\ 0,& \mathrm{other}\end{array}\right.$

得到决策表中各点的值后，就可计算融合图像F的IMF系数：

${\mathrm{IMF}}_j(F,p)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_j(X_i,p)\·{\mathrm{IMF}}_j(X_i,p)$

采用基于区域极大值的局部显著性准则进行IMF分量融合的思想如图5所示。

3、对多幅源图像的残差分量进行融合处理

采用主成分分析(PCA)方法，分别计算出不同源图像对应的残差分量的自适应融合权重，按权重进行残差分量的融合处理。

设有n幅显微源图像，把每幅图像看作一维向量记做x_i，i＝1，2，...，n，残差分量的融合处理实现过程如下：

步骤1：由源图像构造数据矩阵X：X＝(x₁，x₂，…，x_n)^T。

步骤2：计算数据矩阵X的协方差矩阵C：

σ_i，j ²为图像的方差，

x_i为第i幅源图像的平均灰度值。

$X={(x_1,x_2,...,x_n)}^T=\left(\begin{array}{ccccc}{x}_{11}& ...& x_{1j}& ...& x_{1m}\\ ...& ...& ...& ...& ...\\ x_{i1}& ...& x_{\mathrm{ij}}& ...& x_{\mathrm{im}}\\ ...& ...& ...& ...& ...\\ x_{n1}& ...& x_{\mathrm{nj}}& ...& x_{\mathrm{nm}}\end{array}\right);$ $C=\left(\begin{array}{ccccc}{\mathrm{\σ}}_{11}& ...& {\mathrm{\σ}}_{1j}& ...& {\mathrm{\σ}}_{1m}\\ ...& ...& ...& ...& ...\\ {\mathrm{\σ}}_{i1}& ...& {\mathrm{\σ}}_{\mathrm{ij}}& ...& {\mathrm{\σ}}_{\mathrm{im}}\\ ...& ...& ...& ...& ...\\ {\mathrm{\σ}}_{n1}& ...& {\mathrm{\σ}}_{\mathrm{nj}}& ...& {\mathrm{\σ}}_{\mathrm{nm}}\end{array}\right)$

步骤3：计算协方差矩阵的C的特征值λ及相应的特征向量ξ：

由特征值方程|λI-C|＝0，求出特征值λ_i和对应的特征向量ξ_i(i＝1，2，...，n)。

步骤4：确定加权系数ω_i：

步骤5：计算最终融合残差分量I_残差：

4、各融合后分量反向重构获取融合结果

将融合后的残差分量I_残差和融合后的各级固有模态函数IMF_j反向重构获取融合结果图像I_融合结果：

为了验证本发明方法的有效性，使用一组藻类显微源图像进行图像融合处理。图6中(a)、(b)与(c)为三幅待处理多聚焦藻类显微源图像，每幅图像中都有部分细胞处于聚焦位置，成清晰像，也有部分细胞处于离焦位置，成模糊像。(d)为使用本发明方法进行图像融合处理之后得到的各部位聚焦都很清晰的融合结果图像。

目前，在图像融合领域，小波变换方法是最受欢迎且使用最广的方法。为了证明本发明方法性能上的优越，使用相同的多聚焦显微源图像(图6中所示藻类显微源图像)，分别使用主分量分析方法、小波变换方法和本发明方法进行图像融合处理，得到的客观评价结果对照如表1所示。其中小波变换方法使用了bior4.4小波基，尺度系数的融合采用常见的平均融合规则，小波系数的融合采用常见的区域绝对值最大融合规则。这里选用了信息熵和标准差这两种常用的评价指标来进行融合结果客观评价。图像的熵值是衡量图像信息丰富程度的一个重要指标，若融合图像的熵越大，表示所含的信息越丰富，融合质量越好。标准差反映了图像灰度相对于灰度平均值的离散情况，若标准差大，则图像灰度级分布分散，图像的反差大，对比度大，可以看出更多的信息。它们的定义如下：

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

$\mathrm{STD}=\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(F(i,j)-\widehat{\mathrm{\μ}})}^2/(M\×N)}$

表1.本发明方法与其他方法的融合客观评价值对比(源图像为图6所示显微图像)

从表1中可以看出，本发明方法使得融合后图像的质量得到提高，在两个客观评价指标上都优于主分量分析方法和相同分解级数的小波变换方法。另外，随着分解级数的增加，尽管需要花费更多的运算时间，本发明方法可以获得更好的融合结果。而且，本发明方法由于使用二维经验模态分解方法作为多聚焦显微图像的多尺度分解手段，可以达到自适应提取图像细节信息的目的，从而避免小波方法需要选取合适的小波基等麻烦。同时BEMD是一种全二维的更为满足人的视觉特性的多尺度分解，相比小波方法只提取水平、垂直、对角线三个方向的细节信息，本发明方法使用BEMD来分解红外显微图像，能获得更为丰富的细节信息，进而得到质量更高的融合结果图像。

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

Claims (7)

1.一种基于二维经验模态分解的显微图像融合方法，该方法采用二维经验模态分解方法对采集的序列显微源图像进行多尺度分解，获得源图像的多级尺度的高频分量，按照局部显著性准则进行融合处理，并对源图像的低频分量采用主成分分析方法进行融合处理，最后反向重构获取融合图像，该方法包括如下步骤：

（1）对采集的序列显微源图像X₁，X₂，…，X_n分别进行二维经验模态分解处理，二维经验模态分解BEMD的过程反映了提取局部最高频、次高频的过程，得到每幅源图像的n级固有模态函数分量IMF和一个残差分量；

（2）对不同源图像对应的各级固有模态函数分量IMF中的像素，采用基于区域极大值的局部显著性选择准则进行融合处理，使得将融合后的固有模态函数分量具备所有源图像的清晰聚焦的细节信息，所述融合处理包括：

（2－1）首先，对于显微源图像X₁，为其各级IMF系数定义一个衡量其显著性的变量S：

$S_j(X_1,p)=\underset{q\∈Q}{\max }\left(\underset{q\∈Q}{\mathrm{\Σ}}\right|{\mathrm{IMF}}_j(X_1,q)\left|\right);$

其中，j表示IMF系数的级数，j=1表示该分量为第一级IMF分量，反映显微图像X₁的局部最高频信息，j=2表示该分量为第二级IMF分量，反映显微图像X₁的局部次高频信息，以此类推；p=(m,n)表示IMF系数的空间位置；Q表示以p为中心的一个3×3的方形窗口，q为窗口内的任意一点；IMF_j表示图像X₁位于q点位置的第j级IMF系数值；

对于图像X_i中对应的IMF系数同样定义S_j(X_i,p)；

（2－2）然后，在不同源图像的各级IMF系数中选择对应的显著性变量S值较大的IMF系数作为合成图像中对应位置的IMF系数；如果用

表示图像X_i相应位置上的决策表的值，则可表示为：

$C_j^{*}(X_1,p)=\left\{\begin{array}{cc}1,& S_j(X_1,p)\≥S_j(X_i,p),i=\mathrm{2,3},...n\\ 0,& \mathrm{other}\end{array}\right.$

$C_j^{*}(X_2,p)=\left\{\begin{array}{cc}1,& S_j(X_2,p)\≥S_j(X_i,p),i=\mathrm{1,3},...n\\ 0,& \mathrm{other}\end{array}\right.;$

...

$C_j^{*}(X_n,p)=\left\{\begin{array}{cc}1,& S_j(X_n,p)\≥S_j(X_i,p),i=\mathrm{1,2},...n-1\\ 0,& \mathrm{other}\end{array}\right.$

对得到的决策表进行基于多数表决原则的一致性验证，令修正后决策表的值为C_j(X_i,p)，则：

$C_j(X_1,p)=\left\{\begin{array}{cc}1,& \underset{q\∈Q}{\mathrm{\Σ}}{C}_j^{*}(X_1,q)\≥5\\ 0,& \mathrm{other}\end{array}\right.$

$C_j(X_2,p)=\left\{\begin{array}{cc}1,& \underset{q\∈Q}{\mathrm{\Σ}}{C}_j^{*}(X_2,q)\≥5\\ 0,& \mathrm{other}\end{array}\right.;$

...

$C_j(X_n,p)=\left\{\begin{array}{cc}1,& \underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{C}_j(X_i,p)=0\\ 0,& \mathrm{other}\end{array}\right.$

（2－3）得到决策表中各点的值后，就可计算融合图像F的IMF系数：

${\mathrm{IMF}}_j(F,p)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_j(X_i,p)\·{\mathrm{IMF}}_j(X_i,p);$

（3）利用主成分分析PCA方法，分别计算出不同源图像对应的残差分量的自适应融合权重，按权重进行残差分量的融合处理；

（4）将融合后的各级固有模态函数分量和残差分量反向重构获取融合图像。

2.根据权利要求1所述的基于二维经验模态分解的显微图像融合方法，其特征在于，所述的步骤（1）中对每一幅源图像进行图像高频到低频的自然尺度分离，首先，分解出来的第1级固有模态函数IMF₁是图像所含有的最高频率分量，该分量的各处频率都对应着图像在各处的局部最高频，源图像减去第1级固有模态函数得到第1级残差分量；对第1级残差分量再进行分解，得到第2级固有模态函数和第2级残差分量；依此类推，得到n级固有模态函数和第n级残差；

所述的二维经验模态分解BEMD的处理过程包括如下步骤：

（1－1）：为了避免二维经验模态分解产生边界效应，采用局部镜像延拓对原始图像进行四周边界处理：I_原图=F_镜像（I_原图）；

（1－2）：初始化：I=I_原图，I_残差=I，j=0，j表示IMF的分解级数；

（1－3）：对所处理的残差图像曲面I_残差求取曲面局部极值点，包括所有局部极大值和极小值，初始时，I_残差就是源图像曲面I_原图；

（1－4）：对各极大值点和各极小值点分别进行曲面拟合，经插值后得到极大值点对应的上包络曲面E_u和极小值点对应的下包络曲面E_l；

（1－5）：将两曲面数据求平均得到均值包络曲面数据E_m：E_m=（E_u+E_l）/2；

（1－6）：计算筛分终止条件标准偏差SD：

（1－7）：提取细节，用残差图像曲面减去均值包络曲面求得差值：I_残差=I_残差-E_m；

（1－8）：重复上述步骤（1－3）～步骤（1－7），直到满足给定的终止条件：（a）IMF的极值点和过零点数目必须相等或至多只相差一点；（b）在每一像素点，由极大值点定义的上包络线和由极小值点定义的下包络线的平均值为零；

（1－9）：计算残差，用图像I减去第j层固有模态函数即I_残差得到第j层残差分量并赋值给I：I=I-I_残差；

（1－10）：对残差分量I重复步骤（1－3）～步骤（1－9），直到满足残差不含IMF分量或已达到所需要的运算级数，依次得到图像的n级固有模态函数和第n级残差分量。

3.根据权利要求2所述的基于二维经验模态分解的显微图像融合方法，其特征在于，所述的步骤（1－3）中，图像的局部极大值点为灰度值比周围3×3区域8个相邻像素点灰度值都高的点，图像的极小值点为灰度值比周围8个相邻像素点灰度值都低的点。

4.根据权利要求2所述的基于二维经验模态分解的显微图像融合方法，其特征在于，所述的步骤（1－4）基于Delaunay三角剖分的插值方法对各极大值点和各极小值点分别进行曲面拟合；计算上包络曲面E_u时，将各极大值点分割成若干简单的小三角形区域，三角形的顶点即是图像的局部极大值点，在每个小三角形中立方插值构造插值曲面，再进一步把这些曲面片拼接起来，构造出一张大的插值曲面，即上包络曲面E_u；同样，用此方法再构造出极小值点对应的下包络曲面E_l。

5.根据权利要求2所述的基于二维经验模态分解的显微图像融合方法，其特征在于，所述的步骤（1－8）采用简化的终止条件：当步骤（1－6）中计算得到的SD小于0.3时，就认为满足了收敛条件，令j＝j+1，此时得到的I_残差即为第j层二维固有模态函数IMF_j。

6.根据权利要求1所述的基于二维经验模态分解的显微图像融合方法，其特征在于，所述的步骤（3）对n幅显微源图像，把每幅图像看作一维向量记做x_i，i=1,2，...,n，残差分量的融合处理包括：

（3－1）由源图像构造数据矩阵X：X=(x₁,x₂,…,x_n)^T；

（3－2）计算数据矩阵X的协方差矩阵C：

为图像的方差，

为第i幅源图像的平均灰度值；

$X={(x_1,x_2,...,x_n)}^T=\left(\begin{array}{ccccc}{x}_{11}& ...& x_{1j}& ...& x_{1m}\\ ...& ...& ...& ...& ...\\ x_{i1}& ...& x_{\mathrm{ij}}& ...& x_{\mathrm{im}}\\ ...& ...& ...& ...& ...\\ x_{n1}& ...& x_{\mathrm{nj}}& ...& x_{\mathrm{nm}}\end{array}\right);$ $C=\left(\begin{array}{ccccc}{\mathrm{\σ}}_{11}& ...& {\mathrm{\σ}}_{1j}& ...& {\mathrm{\σ}}_{1m}\\ ...& ...& ...& ...& ...\\ {\mathrm{\σ}}_{i1}& ...& {\mathrm{\σ}}_{\mathrm{ij}}& ...& {\mathrm{\σ}}_{\mathrm{im}}\\ ...& ...& ...& ...& ...\\ {\mathrm{\σ}}_{n1}& ...& {\mathrm{\σ}}_{\mathrm{nj}}& ...& {\mathrm{\σ}}_{\mathrm{nm}}\end{array}\right);$

（3－3）计算协方差矩阵的C的特征值λ及相应的特征向量ξ：

由特征值方程|λI-C|=0，求出特征值λ_i和对应的特征向量ξ_i（i=1,2,…,n）；

（3－4）确定加权系数ω_i：

（3－5）计算最终融合残差分量I_残差：

7.根据权利要求5所述的基于二维经验模态分解的显微图像融合方法，其特征在于，所述的步骤（4）反向重构获取融合图像I_融合结果：

Images