用于多图像融合的二维多元信号经验模态分解快速算法的制作方法

本发明属于多图像融合领域，具体涉及一种用于多图像融合的二维多元信号经验模态分解快速算法。

背景技术：

图像融合技术时常被运用在两个场景下，一种情况是将多张焦点不同的散焦图片融合在一起形成一张清晰的对焦图片；另一种情况是将同一场景多次曝光的图片进行融合，过度曝光或者曝光不足都会使得图片信息缺失，融合之后会使得图像信息充足。

近些年提出来的经验模态分解（empiricalmodedecomposition，emd）是一种完全由数据驱动的信号分解算法。它将信号分解成不同振荡频率的成分，这些成分被称为本征模态函数（intrinsicmodefunctions，imfs）。在emd的基础之上又有了许多扩展性的发展。多元经验模态分解（multivariateempiricalmodedecomposition，memd）通过将多元信号同时映射到不同方向上，得到多元信号的局部均值和信号包络，以此实现对多元信号的分解。二维经验模态分解（bidimensionalempiricalmodedecomposition，bemd）将emd运用到二维信号上，通过寻找二维曲面的局部极大值和极小值，在此基础上构造均值曲面，从而实现二维信号（图像）的分解。在bemd基础上提出的一种基于均值估计的快速算法（meanapproximationbidimensionalempiricalmodedecomposition，ma-bemd）可以实现二维信号的快速分解。但是使用该方法进行图像融合时无法保证待融合的每张图的本征模态函数都分解到相同空间尺度。二维多元经验模态分解（bidimensionalmultivariateempiricalmodedecomposition，bmemd）将多元经验模态分解运用到图像融合上，实现了对不同图像在相同空间尺度上的融合，减少了图像融合时空间尺度的丢失。bmemd在生成本征模态函数时需要不停地重复筛选过程，且每次筛选过程都需要对极大值和极小值分别进行插值，得到极小值和极大值包络，需要耗费较多地计算时间。

技术实现要素：

发明目的：针对现有技术中存在的问题，本发明公开了一种用于多图像融合的二维多元信号经验模态分解快速算法，将一种基于均值估计的快速算法运用到二维多元经验模态分解上，快速地将多张图片分解到相同空间尺度，并在此基础上进行融合。

技术方案：本发明采用如下技术方案：一种用于多图像融合的二维多元信号经验模态分解快速算法，其特征在于，包括如下步骤：

步骤a、将需要融合的多幅图像组合成一个二维多元信号，作为ma-bmemd的输入信号，其中每一幅图像为二维多元信号中的一个二维信号；

步骤b、将输入信号中的每个二维信号同时分解到若干个空间尺度中，得到投影信号，采用基于均值估计的二维多元经验模态分解算法将每个投影信号分解为若干本征模态函数；

步骤c、使用局部方差最大法或局部能量最大法将同一空间尺度下不同的本征模态函数进行融合，再将所有空间尺度的融合后的信号叠加即可得到融合后的图像；

其中，所述步骤b包括：

步骤b1、将输入的二维多元信号投影到若干个方向上，得到每个方向上的投影信号；

步骤b2、提取每个方向上的投影信号的极值点；

步骤b3、通过投影信号的极值点得到每个二维信号的极值点，结合均值估计快速算法，由极值点得到均值点；

步骤b4、使用三次样条插值对每个方向上的所有均值点进行插值，扩充成图像大小，得到每个方向上的均值曲面；

步骤b5、将所有方向上的均值曲面进行平均，得到二维多元信号的均值曲面；

步骤b6、输入信号减去均值曲面，即为分解出的一个本征模态函数，并判断该本征模态函数是否满足预先设定的停止标准，若不满足，则将该本征模态函数作为输入，若满足，则将均值曲面作为输入，重复步骤b1至b5，直到提取出所有本征模态函数。

优选的，所述步骤b3中，在得到的极值点的基础上构造德劳内三角，对于一个德劳内三角，使用重心作为均值点，均值点坐标为，则：

其中，德劳内三角的三个顶点分别为，，，顶点坐标分别为、和；为均值点的像素值，为德劳内三角的顶点的像素值，表示三个顶点的权重。

优选的，采用等权重划分方式时，；

采用基于欧式距离的权重划分方式时，

其中表示均值点到第p个顶点的欧式距离。

优选的，所述步骤b6中，停止标准为：两次相邻迭代运算的柯西标准差和之间的差异小于事先设定的阈值时，停止迭代，即

其中，柯西标准差，m、n是图像的高和宽；代表第x次筛选的结果；为事先设定的阈值。

优选的，所述步骤b5中，

其中，为第k个方向上的均值曲面，a为二维多元信号的均值曲面。

优选的，所述步骤c中，对同一空间尺度的本征模态函数进行融合的步骤包括：

步骤c1、计算每幅图像的每个本征模态函数的局部方差或局部能量；

步骤c2、计算每幅图像的每个像素点对应的权重矩阵；

步骤c3、通过局部方差和权重矩阵计算该空间尺度下融合后的每个本征模态函数；

步骤c4、将所有空间尺度下融合后的本征模态函数进行叠加，即可得到融合后的图像。

有益效果：本发明具有如下有益效果：

（1）本发明结合了现有技术的优点，成功将二维多元信号（多张图片）同时分解到多个空间尺度，并在此基础上实现了对图像的融合，减少了图像空间信息的丢失；

（2）本发明在将图像分解得到多个本征模态函数的筛选过程中，使用了基于均值估计的快速算法，不再需要拟合极大值曲面和极小值曲面，通过均值估计的方法直接拟合均值曲面，较大地减少了算法的复杂度，提高了程序效率。

附图说明

图1为焦点位于图像左下部分的散焦图像；

图2为焦点位于图像中上部分的散焦图像；

图3为焦点位于图像右下部分的散焦图像；

图4为使用ma-bmemd()方式融合出来的图像；

图5为使用ma_bmemd()方式融合出来的图像；

图6为理想全焦点图像；

图7为使用bmemd融合出来的图像；

图8为某次筛选过程中，使用均值估计算法计算均值点时产生的德劳内三角示意图；

图9为通过均值估计的方法由极值点获得均值点的示意图；

图10为现有的bmemd技术融合多幅图像与本发明融合多幅图像定量评价指标示意图；

图11为本发明的流程示意图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

本发明公开了一种用于多图像融合的二维多元信号经验模态分解快速算法，涉及二维多元信号的快速处理问题，可以快速将多幅图像融合，可应用于全焦点图像、适度曝光图像、全光谱图像等多种场景，本发明使用了一种通过均值估计得到均值曲面的近似算法，改善了经验模态分解时的效率。二维多元经验模态分解用于图像融合的主要思路是将不同的图像同时分解到不同的空间尺度，得到每张图像的二维本征模态函数（bimfs），最后使用局部方差最大法或者局部能量最大法对相同空间尺度的bimfs进行融合。

如图11所示，本发明包含如下步骤：

（1）将需要融合的多幅图像组合成一个二维多元信号，作为基于均值估计的二维多元经验模态分解（ma-bmemd）算法的输入信号，其中每一幅图像为二维多元信号中的一个二维信号。

（2）ma-bmemd算法将输入的二维多元信号中的每个二维信号，采用均值估计算法同时分解到多个空间尺度，即分解成多个本征模态函数（bimfs），二维信号在空间尺度的bimfs代表了该二维信号在该空间尺度的映射。

将二维多元信号分解成多个bimfs，具体包括如下步骤：

（21）确定需要融合的图像的数目为n，投影的方向向量的个数为k。根据准蒙特卡洛法，生成k个单位方向矢量，其中每个单位方向矢量对应的角度为。将二维多元信号i作为输入投影到这个方向上，分别得到投影信号，其中，

其中为二维多元信号i的第l个二维信号，即为归一化的单位方向矢量，为第l个二维信号在第k个方向上的投影信号。

归一化的单位方向矢量生成的具体方法如下：

首先生成非归一化的方向矢量：对于任一正整数c，其基b展开形式与c满足以下等式：

式中，t表示c的基b展开形式的长度；那么，该正整数c的基b反转表示为，如下：

此时，非归一化的方向矢量构成的矩阵h表示为：

式中，为互不相等的正整数。

然后进行归一化，对于给定的非归一化的方向矢量，可以表示为。该方向矢量对应的角度，其中：

方向矢量对应的归一化的单位方向矢量为，其中：

（22）提取投影信号的极值点，记录其坐标。

（23）根据极值点坐标，得到每个二维信号的极值点，其中。通过均值估计快速算法，由极值点得到均值点。

本发明提出了一种基于极值点得到近似估计均值点的快速算法，具体如下：

使用德劳内三角化对若干极值点进行三角划分，得到多个极值点三角形即德劳内三角，如图8所示。对于一个特定的德劳内三角，如图9所示，假设其顶点为，，，顶点坐标分别为、和，使用其重心作为一个均值点，其坐标为，则：

其中，和均为向下取整的整数。

均值点的像素值大小为：

其中，为均值点的像素值，为德劳内三角的顶点的像素值，表示三个顶点的权重。通常，有两种取值方式，一种是等权重方式，即

另一种是基于欧式距离的权重划分方式：

其中表示均值点与第p个顶点的欧式距离：

本发明提出的一种基于极值点得到近似估计均值点的快速算法，其算法复杂度推演如下：

由p个点构成的德劳内三角剖分共含有2p-2-y个三角形，其中y表示该p个点构成的点集中位于边界上的点。因为图像像素点较多，可以近似忽略y，即产生2p-2个德劳内三角形。假设共有p个均值点，使用方式得到均值点的计算复杂度为：

其中表示生成德劳内三角的计算复杂度，表示由德劳内三角得到均值点的计算复杂度，add表示一次加法运算，div表示一次除法运算。

（24）基于德劳内三角划分，使用三次样条插值对均值点进行插值扩充至图像大小，得到投影信号的均值曲面。

本发明基于德劳内三角划分，使用三次样条插值对均值点进行插值，得到均值曲面的过程，其算法复杂度为：

其中m、n分别表示二维信号每一维的长度，即图像的高和宽。

由此可以得到均值曲面估计的算法复杂度：

传统的bmemd计算均值曲面时，需要先计算极大值曲面和极小值曲面，其算法复杂度计算如下：

由此可得到

考虑到对于一张图像，由于，所以，换而言之，对于单次筛选过程，使用均值估计求得均值曲面的方法比传统方法复杂度更低。

同时，在本发明所述的算法中，由极值点得到均值点后，只需要进行一次插值扩充就可以得到均值曲面；而在传统的bmemd方法中，首先需要对极大值点进行差值扩充得到极大值曲面，再对极小值点进行差值扩充得到极小值曲面，需要进行两次插值扩充才能得到均值曲面。因此，本发明的计算效率高于传统的bmemd方法。

（25）将k个方向上的均值曲面进行平均作为二维多元信号i的均值曲面a，即

（26）计算，r即为分解出的一个bimfs，判断r是否满足预先设定的bimfs停止标准。若不满足，则将r作为输入，重复步骤（21）-（25）；若满足，则将作为输入，重复步骤（21）-（25），直到所有bimfs被提取出来。其中，步骤（21）-（25）被称为一次筛选过程。

本发明采用的停止准则如下：当两次相邻迭代运算的柯西标准差（cauchystandarddeviation,sd）和之间的差异小于事先设定的阈值时，停止迭代。柯西标准差按如下公式计算：

相邻两次迭代差异需满足：

其中，代表第x次筛选的结果，m、n是图像的高和宽，为事先设定的阈值。采用柯西标准差的停止准则，本发明所述的算法能够较完整的提取出原信号的bimfs。

（3）在每一个空间尺度，使用局部方差最大法将bimfs进行融合，将融合后的信号作为融合后的图像在该空间尺度的bimfs，将融合后的图像在所有空间尺度的bimfs叠加即可得到融合后的图像。

本发明使用局部方差最大法对同一空间尺度的bimfs进行融合，具体步骤如下：

（31）记第p张图像的第q个bimfs在位置处的像素值大小为，对于每个bimfs计算其局部方差，局部方差定义如下：

其中，z为正整数，用于衡量滑动窗口的大小；表示局部均值，具体计算公式如下：

（32）每个像素点对应的权重矩阵为：

（33）将每幅图像的第q个bimfs进行融合后，得到第q个空间尺度的bimfs，计算公式为：

（34）将融合后的m个尺度空间的bimfs进行叠加，即可得到融合后的图像为

其中m代表步骤（33）中融合得到的bimfs的个数。

若使用局部能量最大法对同一空间的尺度的bimfs经行融合，具体步骤如下：

（31’）记第p张图像的第q个bimfs在位置处的像素值大小为，对于每个bimfs计算其局部能量，局部能量定义如下：

（32’）每个像素点对应的权重矩阵可以计算为：

（33’）将每幅图像的第q个bimfs进行融合后，得到第q个空间尺度的bimfs，计算公式为：

（34’）将融合后的m个尺度空间的bimfs进行叠加，即可得到融合后的图像为

其中m代表步骤（33’）中融合得到的bimfs的个数。

本次实验尝试将多张散焦图像融合，通过与bmemd算法比较融合结果，证明本发明所述的算法在保证一定图片质量的基础上对计算效率的较大提升，体现该算法在工程上的适用性。图1~图3为焦点分别位于左下、中上、右下的散焦图像。将本发明所述的算法应用于图1~图3所述的图像融合时，图4和图5所示分别为使用ma-bmemd（）方式和ma_bmemd（）方式融合出来的图像。图6为理想全焦点图像。图7为使用bmemd方法融合出来的图像。通过比较，可发现使用ma-bmemd和bmemd方法融合出来的图像和理想化的图像有较强的一致性，但是ma-bmemd的计算时间明显短于bmemd。

对于多图像融合的结果，评价指标有评价指标有平均信息熵（entropy）、空间频率（spatialfrequency，sf）、均方根误差（rootmeansquareerror，rmse）、互相关系数（correlation，corr）、峰值信噪比（peaksignal-to-noiseratio）、结构相似性（structuralsimilarityindex，ssim）其计算方式分别如下：

其中，entropy表示所求得的图像平均信息，表示像素值为i的点占图像总像素的比率熵。

其中，sf是所求得的图像空间频率，m、n分别是融合后图像的高和宽，是融合后的图像。

其中，rmse是所求得的融合后的图像与理想图像之间的均方根误差，m、n分别是融合后图像的高和宽，是融合后图像，是理想图像。

其中，

corr是所求得的融合后图像与理想图像之间的互相关系数，和分别表示融合后图像和理想图像的均值，m、n分别是融合后图像的高和宽，是融合后图像，是理想图像。

其中

psnr是所求得得融合后的图像峰值信噪比，mse表示当前图像和参考图像的均方误差，m、n分别是融合后图像的高和宽，是融合后图像，是理想图像。

其中

ssim是所求得融合后的图像结构相似性，通常取，，，一般的，，，，、分别表示图像和的均值，、分别表示图像和的方差，表示图像和的协方差。

如图10所示，现有的bmemd技术融合多幅图像与本发明融合多幅图像定量评价指标示意图。对比可发现本法所述算法和理想图像的有相近的信息熵和空间信息，对比已有的bmemd算法，本发明提出的ma-bmemd算法融合得到图像虽然峰值信噪比、相关性及均方误差相近，同时该快速算法在极大地提升了运算效率，体现了该算法在工程上的适用性和优越性。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。