CN111369488B - 用于多图像融合的二维多元信号经验模态快速分解方法 - Google Patents

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

Info

Publication number
CN111369488B
CN111369488B CN202010465445.2A CN202010465445A CN111369488B CN 111369488 B CN111369488 B CN 111369488B CN 202010465445 A CN202010465445 A CN 202010465445A CN 111369488 B CN111369488 B CN 111369488B
Authority
CN
China
Prior art keywords
signal
image
dimensional
mean
fused
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010465445.2A
Other languages
English (en)
Other versions
CN111369488A (zh
Inventor
方懿德
夏亦犁
裴文江
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jiangsu Jihui Mobile Communication Technology Research Institute Co ltd
Original Assignee
Jiangsu Jihui Mobile Communication Technology Research Institute Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jiangsu Jihui Mobile Communication Technology Research Institute Co ltd filed Critical Jiangsu Jihui Mobile Communication Technology Research Institute Co ltd
Priority to CN202010465445.2A priority Critical patent/CN111369488B/zh
Publication of CN111369488A publication Critical patent/CN111369488A/zh
Application granted granted Critical
Publication of CN111369488B publication Critical patent/CN111369488B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Quality & Reliability (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种用于多图像融合的二维多元信号经验模态快速分解方法，包括如下步骤：将需要融合的多幅图像组合成一个二维多元信号，作为MA‑BMEMD的输入信号，其中每一幅图像为一个二维信号；将输入的二维多元信号分解到多个空间尺度得到投影信号，采用基于均值估计的二维多元经验模态分解算法将每个二维信号在每个空间尺度的投影信号分解为若干本征模态函数；使用局部方差最大法将同一空间尺度下不同的本征模态函数进行融合，将所有空间尺度的融合后的信号叠加即可得到融合后的图像。本发明可以在保证BMEMD融合效果的前提下，减少图像空间信息的丢失，且较为显著地减少算法复杂度，降低运算成本。

Description

用于多图像融合的二维多元信号经验模态快速分解方法
技术领域
本发明属于多图像融合领域，具体涉及一种用于多图像融合的二维多元信号经验模态快速分解方法。
背景技术
图像融合技术时常被运用在两个场景下，一种情况是将多张焦点不同的散焦图片融合在一起形成一张清晰的对焦图片；另一种情况是将同一场景多次曝光的图片进行融合，过度曝光或者曝光不足都会使得图片信息缺失，融合之后会使得图像信息充足。
近些年提出来的经验模态分解（Empirical Mode Decomposition，EMD）是一种完全由数据驱动的信号分解算法。它将信号分解成不同振荡频率的成分，这些成分被称为本征模态函数（Intrinsic Mode Functions，IMFs）。在EMD的基础之上又有了许多扩展性的发展。多元经验模态分解（Multivariate empirical mode decomposition，MEMD）通过将多元信号同时映射到不同方向上，得到多元信号的局部均值和信号包络，以此实现对多元信号的分解。二维经验模态分解（Bidimensional Empirical Mode Decomposition，BEMD）将EMD运用到二维信号上，通过寻找二维曲面的局部极大值和极小值，在此基础上构造均值曲面，从而实现二维信号（图像）的分解。在BEMD基础上提出的一种基于均值估计的快速算法（Mean Approximation Bidimensional Empirical Mode Decomposition，MA-BEMD）可以实现二维信号的快速分解。但是使用该方法进行图像融合时无法保证待融合的每张图的本征模态函数都分解到相同空间尺度。二维多元经验模态分解（Bidimensional MultivariateEmpirical Mode Decomposition，BMEMD）将多元经验模态分解运用到图像融合上，实现了对不同图像在相同空间尺度上的融合，减少了图像融合时空间尺度的丢失。BMEMD在生成本征模态函数时需要不停地重复筛选过程，且每次筛选过程都需要对极大值和极小值分别进行插值，得到极小值和极大值包络，需要耗费较多地计算时间。
发明内容
发明目的：针对现有技术中存在的问题，本发明公开了一种用于多图像融合的二维多元信号经验模态快速分解方法，将一种基于均值估计的快速算法运用到二维多元经验模态分解上，快速地将多张图片分解到相同空间尺度，并在此基础上进行融合。
技术方案：本发明采用如下技术方案：一种用于多图像融合的二维多元信号经验模态快速分解方法，其特征在于，包括如下步骤：
步骤A、将需要融合的多幅图像组合成一个二维多元信号，作为MA-BMEMD的输入信号，其中每一幅图像为二维多元信号中的一个二维信号；
步骤B、将输入信号中的每个二维信号同时分解到若干个空间尺度中，得到投影信号，采用基于均值估计的二维多元经验模态分解算法将每个投影信号分解为若干本征模态函数；
步骤C、使用局部方差最大法或局部能量最大法将同一空间尺度下不同的本征模态函数进行融合，再将所有空间尺度的融合后的信号叠加即可得到融合后的图像；
其中，所述步骤B包括：
步骤B1、将输入的二维多元信号投影到若干个方向上，得到每个方向上的投影信号；
步骤B2、提取每个方向上的投影信号的极值点；
步骤B3、通过投影信号的极值点得到每个二维信号的极值点，结合均值估计快速算法，由极值点得到均值点；
步骤B4、使用三次样条插值对每个方向上的所有均值点进行插值，扩充成图像大小，得到每个方向上的均值曲面；
步骤B5、将所有方向上的均值曲面进行平均，得到二维多元信号的均值曲面；
步骤B6、输入信号减去均值曲面，即为分解出的一个本征模态函数，并判断该本征模态函数是否满足预先设定的停止标准，若不满足，则将该本征模态函数作为输入，若满足，则将均值曲面作为输入，重复步骤B1至B5，直到提取出所有本征模态函数。
优选的，所述步骤B3中，在得到的极值点的基础上构造德劳内三角，对于一个德劳 内三角，使用重心作为均值点
Figure 841241DEST_PATH_IMAGE001
，均值点坐标为
Figure 562073DEST_PATH_IMAGE002
，则：
Figure 727475DEST_PATH_IMAGE003
Figure 559165DEST_PATH_IMAGE004
Figure 860833DEST_PATH_IMAGE005
其中，德劳内三角的三个顶点分别为
Figure 752566DEST_PATH_IMAGE006
，
Figure 139684DEST_PATH_IMAGE007
，
Figure 775065DEST_PATH_IMAGE008
，顶点坐标分别为
Figure 196819DEST_PATH_IMAGE009
、
Figure 259453DEST_PATH_IMAGE010
和
Figure 133868DEST_PATH_IMAGE011
；
Figure 307361DEST_PATH_IMAGE012
为均值点的像素值，
Figure 583621DEST_PATH_IMAGE013
为德劳内三角的顶点的像素值，
Figure 817156DEST_PATH_IMAGE014
表示三个顶点的权重。
优选的，采用等权重划分方式时，
Figure 178868DEST_PATH_IMAGE015
；
采用基于欧式距离的权重划分方式时，
Figure 421630DEST_PATH_IMAGE016
Figure 817976DEST_PATH_IMAGE017
其中
Figure 222413DEST_PATH_IMAGE018
表示均值点到第p个顶点的欧式距离。
优选的，所述步骤B6中，停止标准为：两次相邻迭代运算的柯西标准差
Figure 805841DEST_PATH_IMAGE019
和
Figure 586715DEST_PATH_IMAGE020
之间的差异小于事先设定的阈值时，停止迭代，即
Figure 837568DEST_PATH_IMAGE021
其中，柯西标准差
Figure 412906DEST_PATH_IMAGE022
，M、N是图像的高 和宽；
Figure 483630DEST_PATH_IMAGE023
代表第x次筛选的结果；
Figure 802616DEST_PATH_IMAGE024
为事先设定的阈值。
优选的，所述步骤B5中，
Figure 907975DEST_PATH_IMAGE025
其中，
Figure 654214DEST_PATH_IMAGE026
为第k个方向上的均值曲面，A为二维多元信号的均值曲面。
优选的，所述步骤C中，对同一空间尺度的本征模态函数进行融合的步骤包括：
步骤C1、计算每幅图像的每个本征模态函数的局部方差或局部能量；
步骤C2、计算每幅图像的每个像素点对应的权重矩阵；
步骤C3、通过局部方差和权重矩阵计算该空间尺度下融合后的每个本征模态函数；
步骤C4、将所有空间尺度下融合后的本征模态函数进行叠加，即可得到融合后的图像。
有益效果：本发明具有如下有益效果：
（1）本发明结合了现有技术的优点，成功将二维多元信号（多张图片）同时分解到多个空间尺度，并在此基础上实现了对图像的融合，减少了图像空间信息的丢失；
（2）本发明在将图像分解得到多个本征模态函数的筛选过程中，使用了基于均值估计的快速算法，不再需要拟合极大值曲面和极小值曲面，通过均值估计的方法直接拟合均值曲面，较大地减少了算法的复杂度，提高了程序效率。
附图说明
图1为焦点位于图像左下部分的散焦图像；
图2为焦点位于图像中上部分的散焦图像；
图3为焦点位于图像右下部分的散焦图像；
图4为使用MA-BMEMD(
Figure 477814DEST_PATH_IMAGE027
)方式融合出来的图像；
图5为使用MA_BMEMD(
Figure 334911DEST_PATH_IMAGE028
)方式融合出来的图像；
图6为理想全焦点图像；
图7为使用BMEMD融合出来的图像；
图8为某次筛选过程中，使用均值估计算法计算均值点时产生的德劳内三角示意图；
图9为通过均值估计的方法由极值点获得均值点的示意图；
图10为现有的BMEMD技术融合多幅图像与本发明融合多幅图像定量评价指标示意图；
图11为本发明的流程示意图。
具体实施方式
下面结合附图对本发明作更进一步的说明。
本发明公开了一种用于多图像融合的二维多元信号经验模态快速分解方法，涉及二维多元信号的快速处理问题，可以快速将多幅图像融合，可应用于全焦点图像、适度曝光图像、全光谱图像等多种场景，本发明使用了一种通过均值估计得到均值曲面的近似算法，改善了经验模态分解时的效率。二维多元经验模态分解用于图像融合的主要思路是将不同的图像同时分解到不同的空间尺度，得到每张图像的二维本征模态函数（BIMFs），最后使用局部方差最大法或者局部能量最大法对相同空间尺度的BIMFs进行融合。
如图11所示，本发明包含如下步骤：
（1）将需要融合的多幅图像组合成一个二维多元信号，作为基于均值估计的二维多元经验模态分解（MA-BMEMD）算法的输入信号，其中每一幅图像为二维多元信号中的一个二维信号。
（2）MA-BMEMD算法将输入的二维多元信号中的每个二维信号，采用均值估计算法同时分解到多个空间尺度，即分解成多个本征模态函数（BIMFs），二维信号在空间尺度的BIMFs代表了该二维信号在该空间尺度的映射。
将二维多元信号分解成多个BIMFs，具体包括如下步骤：
（21）确定需要融合的图像的数目为n，投影的方向向量的个数为K。根据准蒙特卡 洛法，生成K个单位方向矢量，其中每个单位方向矢量对应的角度为
Figure 560356DEST_PATH_IMAGE029
。将二维多元信号I 作为输入投影到这
Figure 477497DEST_PATH_IMAGE030
个方向上，分别得到投影信号
Figure 257234DEST_PATH_IMAGE031
，其中，
Figure 183602DEST_PATH_IMAGE032
其中
Figure 263553DEST_PATH_IMAGE033
为二维多元信号I的第l个二维信号，
Figure 86016DEST_PATH_IMAGE034
即为归一化的单位方向矢量，
Figure 884207DEST_PATH_IMAGE035
为第l个二维信号在第k个方向上的投影信号。
归一化的单位方向矢量
Figure 83107DEST_PATH_IMAGE036
生成的具体方法如下：
首先生成非归一化的方向矢量：对于任一正整数c，其基b展开形式
Figure 283145DEST_PATH_IMAGE037
与c满足以下等式：
Figure 542088DEST_PATH_IMAGE038
式中，T表示c的基b展开形式
Figure 296417DEST_PATH_IMAGE039
的长度；那么，该正整数c的基b反转表示为
Figure 564587DEST_PATH_IMAGE040
， 如下：
Figure 353552DEST_PATH_IMAGE041
此时，非归一化的方向矢量构成的矩阵H表示为：
Figure 783396DEST_PATH_IMAGE042
式中，
Figure 290601DEST_PATH_IMAGE043
为互不相等的正整数。
然后进行归一化，对于给定的非归一化的方向矢量
Figure 96883DEST_PATH_IMAGE044
，可以表示为
Figure 740354DEST_PATH_IMAGE045
。该方向矢量对应的角度
Figure 606679DEST_PATH_IMAGE046
，其中：
Figure 335600DEST_PATH_IMAGE047
方向矢量
Figure 945573DEST_PATH_IMAGE044
对应的归一化的单位方向矢量为
Figure 443551DEST_PATH_IMAGE048
，其中：
Figure 480777DEST_PATH_IMAGE049
（22）提取投影信号
Figure 696994DEST_PATH_IMAGE050
的极值点，记录其坐标
Figure 845079DEST_PATH_IMAGE051
。
（23）根据极值点坐标，得到每个二维信号的极值点
Figure 463142DEST_PATH_IMAGE052
，其中
Figure 671270DEST_PATH_IMAGE053
。通过均值估计快速算法，由极值点
Figure 374783DEST_PATH_IMAGE052
得到均值点
Figure 326559DEST_PATH_IMAGE054
。
本发明提出了一种基于极值点
Figure 799128DEST_PATH_IMAGE052
得到近似估计均值点
Figure 178157DEST_PATH_IMAGE054
的快速算法，具 体如下：
使用德劳内三角化对若干极值点
Figure 380686DEST_PATH_IMAGE052
进行三角划分，得到多个极值点三角形即德 劳内三角，如图8所示。对于一个特定的德劳内三角，如图9所示，假设其顶点为
Figure 870573DEST_PATH_IMAGE006
，
Figure 463228DEST_PATH_IMAGE007
，
Figure 747579DEST_PATH_IMAGE008
，顶点坐标分别为
Figure 160106DEST_PATH_IMAGE009
、
Figure 453684DEST_PATH_IMAGE010
和
Figure 900846DEST_PATH_IMAGE055
，使 用其重心作为一个均值点
Figure 621677DEST_PATH_IMAGE054
，其坐标为
Figure 521500DEST_PATH_IMAGE056
，则：
Figure 353190DEST_PATH_IMAGE057
Figure 920438DEST_PATH_IMAGE058
其中，
Figure 812170DEST_PATH_IMAGE059
和
Figure 199289DEST_PATH_IMAGE060
均为向下取整的整数。
均值点
Figure 834670DEST_PATH_IMAGE054
的像素值大小为：
Figure 256424DEST_PATH_IMAGE061
其中，
Figure 319058DEST_PATH_IMAGE062
为均值点的像素值，
Figure 193473DEST_PATH_IMAGE063
为德劳内三角的顶点的像素值，
Figure 366965DEST_PATH_IMAGE064
表示三个顶点的权重。通常，
Figure 643226DEST_PATH_IMAGE065
有两种取值方式，一种是等权重
Figure 876761DEST_PATH_IMAGE066
方式，即
Figure 238472DEST_PATH_IMAGE067
另一种是基于欧式距离的权重划分
Figure 950076DEST_PATH_IMAGE068
方式：
Figure 80843DEST_PATH_IMAGE069
其中
Figure 485280DEST_PATH_IMAGE070
表示均值点与第p个顶点的欧式距离：
Figure 334287DEST_PATH_IMAGE071
本发明提出的一种基于极值点
Figure 115161DEST_PATH_IMAGE052
得到近似估计均值点
Figure 100435DEST_PATH_IMAGE054
的快速算法，其 算法复杂度推演如下：
由P个点构成的德劳内三角剖分共含有2P-2-y个三角形，其中y表示该P个点构成 的点集中位于边界上的点。因为图像像素点较多，可以近似忽略y，即产生2P-2个德劳内三 角形。假设共有P个均值点，使用
Figure 675773DEST_PATH_IMAGE066
方式得到均值点的计算复杂度为：
Figure 12076DEST_PATH_IMAGE072
Figure 331062DEST_PATH_IMAGE073
Figure 436421DEST_PATH_IMAGE074
Figure 182660DEST_PATH_IMAGE075
其中
Figure 6260DEST_PATH_IMAGE076
表示生成德劳内三角的计算复杂度，
Figure 863358DEST_PATH_IMAGE077
表示由德劳内三角得到均值点的计算复杂度，ADD表示一次加法运算，DIV表示一次除法运 算。
（24）基于德劳内三角划分，使用三次样条插值对均值点
Figure 823223DEST_PATH_IMAGE054
进行插值扩充至图 像大小，得到投影信号
Figure 5943DEST_PATH_IMAGE050
的均值曲面
Figure 51259DEST_PATH_IMAGE078
。
本发明基于德劳内三角划分，使用三次样条插值对均值点
Figure 712048DEST_PATH_IMAGE054
进行插值，得到 均值曲面
Figure 791999DEST_PATH_IMAGE078
的过程，其算法复杂度为：
Figure 614462DEST_PATH_IMAGE079
其中M、N分别表示二维信号每一维的长度，即图像的高和宽。
由此可以得到均值曲面估计的算法复杂度：
Figure 147074DEST_PATH_IMAGE080
Figure 611554DEST_PATH_IMAGE081
传统的BMEMD计算均值曲面时，需要先计算极大值曲面和极小值曲面，其算法复杂度计算如下：
Figure 546012DEST_PATH_IMAGE082
Figure 804955DEST_PATH_IMAGE083
由此可得到
Figure 824863DEST_PATH_IMAGE084
Figure 827454DEST_PATH_IMAGE085
其中，
Figure 616419DEST_PATH_IMAGE086
为计算极大值曲面的计算复杂度，
Figure 311842DEST_PATH_IMAGE087
为计算极小 值曲面的计算复杂度，
Figure 819047DEST_PATH_IMAGE088
为计算均值曲面的计算复杂度。
考虑到对于一张图像，由于
Figure 359750DEST_PATH_IMAGE089
，所以
Figure 268800DEST_PATH_IMAGE090
，换而言之，对于单次筛选过程，使用均值 估计求得均值曲面的方法比传统方法复杂度更低。
同时，在本发明所述的算法中，由极值点得到均值点后，只需要进行一次插值扩充就可以得到均值曲面；而在传统的BMEMD方法中，首先需要对极大值点进行差值扩充得到极大值曲面，再对极小值点进行差值扩充得到极小值曲面，需要进行两次插值扩充才能得到均值曲面。因此，本发明的计算效率高于传统的BMEMD方法。
（25）将K个方向上的均值曲面
Figure 135125DEST_PATH_IMAGE078
进行平均作为二维多元信号I的均值曲面A，即
Figure 129626DEST_PATH_IMAGE091
（26）计算
Figure 474019DEST_PATH_IMAGE092
，R即为分解出的一个BIMFs，判断R是否满足预先设定的 BIMFs停止标准。若不满足，则将R作为输入，重复步骤（21）-（25）；若满足，则将
Figure 237576DEST_PATH_IMAGE093
作为 输入，重复步骤（21）-（25），直到所有BIMFs被提取出来。其中，步骤（21）-（25）被称为一次筛 选过程。
本发明采用的停止准则如下：当两次相邻迭代运算的柯西标准差（Cauchy standard deviation, SD）
Figure 274802DEST_PATH_IMAGE094
和
Figure 491020DEST_PATH_IMAGE095
之间的差异小于事先设定的阈值时，停止迭代。柯 西标准差按如下公式计算：
Figure 639104DEST_PATH_IMAGE022
相邻两次迭代差异需满足：
Figure 257168DEST_PATH_IMAGE096
其中，
Figure 465295DEST_PATH_IMAGE097
代表第x次筛选的结果，M、N是图像的高和宽，
Figure 168809DEST_PATH_IMAGE098
为事先设 定的阈值。采用柯西标准差的停止准则，本发明所述的算法能够较完整的提取出原信号的 BIMFs。
（3）在每一个空间尺度，使用局部方差最大法将BIMFs进行融合，将融合后的信号作为融合后的图像在该空间尺度的BIMFs，将融合后的图像在所有空间尺度的BIMFs叠加即可得到融合后的图像。
本发明使用局部方差最大法对同一空间尺度的BIMFs进行融合，具体步骤如下：
（31）记第p张图像的第q个BIMFs在位置
Figure 855005DEST_PATH_IMAGE099
处的像素值大小为
Figure 593154DEST_PATH_IMAGE100
，对于每 个BIMFs计算其局部方差，局部方差定义如下：
Figure 972183DEST_PATH_IMAGE101
其中，z为正整数，用于衡量滑动窗口的大小；
Figure 162993DEST_PATH_IMAGE102
表示局部均值，具体计算公 式如下：
Figure 652880DEST_PATH_IMAGE103
（32）每个像素点对应的权重矩阵为：
Figure 979956DEST_PATH_IMAGE104
（33）将每幅图像的第q个BIMFs进行融合后，得到第q个空间尺度的BIMFs，计算公式为：
Figure 529886DEST_PATH_IMAGE105
（34）将融合后的m个尺度空间的BIMFs进行叠加，即可得到融合后的图像为
Figure 207992DEST_PATH_IMAGE106
其中m代表步骤（33）中融合得到的BIMFs的个数。
若使用局部能量最大法对同一空间尺度的BIMFs经行融合，具体步骤如下：
（31’）记第p张图像的第q个BIMFs在位置
Figure 501570DEST_PATH_IMAGE107
处的像素值大小为
Figure 948732DEST_PATH_IMAGE108
，对于每 个BIMFs计算其局部能量，局部能量定义如下：
Figure 403984DEST_PATH_IMAGE109
（32’）每个像素点
Figure 569386DEST_PATH_IMAGE110
对应的权重矩阵
Figure 401076DEST_PATH_IMAGE111
可以计算为：
Figure 702744DEST_PATH_IMAGE112
（33’）将每幅图像的第q个BIMFs进行融合后，得到第q个空间尺度的BIMFs，计算公式为：
Figure 328898DEST_PATH_IMAGE113
（34’）将融合后的m个尺度空间的BIMFs进行叠加，即可得到融合后的图像为
Figure 981596DEST_PATH_IMAGE114
其中m代表步骤（33’）中融合得到的BIMFs的个数。
本次实验尝试将多张散焦图像融合，通过与BMEMD算法比较融合结果，证明本发明 所述的方法在保证一定图片质量的基础上对计算效率的较大提升，体现该方法在工程上的 适用性。图1~图3为焦点分别位于左下、中上、右下的散焦图像。将本发明所述的方法应用于 图1~图3所述的图像融合时，图4和图5所示分别为使用MA-BMEMD（
Figure 616976DEST_PATH_IMAGE115
）方式和MA_ BMEMD（
Figure 773151DEST_PATH_IMAGE116
）方式融合出来的图像。图6为理想全焦点图像。图7为使用BMEMD方法融 合出来的图像。通过比较，可发现使用MA-BMEMD和BMEMD方法融合出来的图像和理想化的图 像有较强的一致性，但是MA-BMEMD的计算时间明显短于BMEMD。
对于多图像融合的结果，评价指标有平均信息熵（Entropy）、空间频率（SpatialFrequency，SF）、均方根误差（Root Mean Square Error，RMSE）、互相关系数（Correlation，Corr）、峰值信噪比（Peak Signal-to-Noise Ratio）和结构相似性（StructuralSimilarity Index，SSIM），其计算方式分别如下：
Figure 835785DEST_PATH_IMAGE117
其中，Entropy表示所求得的图像平均信息，
Figure 710200DEST_PATH_IMAGE118
表示像素值为i的点占图像总像素 的比率熵。
Figure 149272DEST_PATH_IMAGE119
其中，SF是所求得的图像空间频率，M、N分别是融合后图像的高和宽，
Figure 425533DEST_PATH_IMAGE120
是融合后 的图像。
Figure 659068DEST_PATH_IMAGE121
其中，RMSE是所求得的融合后的图像与理想图像之间的均方根误差，M、N分别是融 合后图像的高和宽，
Figure 20779DEST_PATH_IMAGE120
是融合后图像，
Figure 997962DEST_PATH_IMAGE122
是理想图像。
Figure 128729DEST_PATH_IMAGE123
其中，
Figure 533166DEST_PATH_IMAGE124
Figure 116594DEST_PATH_IMAGE125
Corr是所求得的融合后图像与理想图像之间的互相关系数，
Figure 897468DEST_PATH_IMAGE126
和
Figure 882742DEST_PATH_IMAGE127
分别表示融 合后图像和理想图像的均值，M、N分别是融合后图像的高和宽，
Figure 458079DEST_PATH_IMAGE120
是融合后图像，
Figure 794383DEST_PATH_IMAGE122
是理想图 像。
Figure 113369DEST_PATH_IMAGE128
其中，
Figure 218728DEST_PATH_IMAGE129
PSNR是所求得融合后图像的峰值信噪比，MSE表示当前图像和参考图像的均方误 差，M、N分别是融合后图像的高和宽，
Figure 964967DEST_PATH_IMAGE120
是融合后图像，
Figure 522987DEST_PATH_IMAGE122
是理想图像。
Figure 645664DEST_PATH_IMAGE130
Figure 605530DEST_PATH_IMAGE131
其中，
Figure 522670DEST_PATH_IMAGE132
Figure 833566DEST_PATH_IMAGE133
Figure 494354DEST_PATH_IMAGE134
SSIM是所求得融合后的图像结构相似性，通常取
Figure 308727DEST_PATH_IMAGE135
，
Figure 662348DEST_PATH_IMAGE136
，
Figure 929381DEST_PATH_IMAGE137
，一般的，
Figure 393860DEST_PATH_IMAGE138
，
Figure 62739DEST_PATH_IMAGE139
，
Figure 587261DEST_PATH_IMAGE140
，
Figure 607170DEST_PATH_IMAGE141
、
Figure 609761DEST_PATH_IMAGE142
分别 表示图像
Figure 398725DEST_PATH_IMAGE120
和
Figure 94149DEST_PATH_IMAGE122
的均值，
Figure 335774DEST_PATH_IMAGE143
、
Figure 142056DEST_PATH_IMAGE144
分别表示图像
Figure 51107DEST_PATH_IMAGE120
和
Figure 651852DEST_PATH_IMAGE122
的方差，
Figure 646353DEST_PATH_IMAGE145
表示图像
Figure 990747DEST_PATH_IMAGE120
和
Figure 754303DEST_PATH_IMAGE122
的协方 差。
如图10所示为现有的BMEMD技术融合多幅图像与本发明融合多幅图像定量评价指标示意图。对比可发现本方法所述算法和理想图像的有相近的信息熵和空间信息，对比已有的BMEMD算法，本发明提出的MA-BMEMD算法融合得到的图像虽然峰值信噪比、相关性及均方误差相近，同时该快速算法在极大地提升了运算效率，体现了该算法在工程上的适用性和优越性。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种用于多图像融合的二维多元信号经验模态快速分解方法，其特征在于，包括如下步骤：
步骤A、将需要融合的多幅图像组合成一个二维多元信号，作为基于均值估计的二维多元经验模态分解算法，即MA-BMEMD的输入信号，其中每一幅图像为二维多元信号中的一个二维信号；
步骤B、将输入信号中的每个二维信号同时分解到若干个空间尺度中，得到投影信号，采用基于均值估计的二维多元经验模态分解算法将每个投影信号分解为若干本征模态函数；
步骤C、使用局部方差最大法或局部能量最大法将同一空间尺度下不同的本征模态函数进行融合，再将所有空间尺度的融合后的信号叠加即可得到融合后的图像；
其中，所述步骤B包括：
步骤B1、将输入的二维多元信号投影到若干个方向上，得到每个方向上的投影信号；
步骤B2、提取每个方向上的投影信号的极值点；
步骤B3、通过投影信号的极值点得到每个二维信号的极值点，在极值点的基础上构造德劳内三角，在德劳内三角内结合均值估计快速算法，由极值点得到均值点，具体为：使用德劳内三角化对极值点进行三角划分，得到多个极值点三角形即德劳内三角，每个德劳内三角的重心作为一个均值点，对于每个德劳内三角，通过三个顶点的坐标与像素值计算对应的重心即均值点的坐标与像素值；
步骤B4、使用三次样条插值对每个方向上的所有均值点进行插值，扩充成图像大小，得到每个方向上的均值曲面；
步骤B5、将所有方向上的均值曲面进行平均，得到二维多元信号的均值曲面；
步骤B6、输入信号减去均值曲面，即为分解出的一个本征模态函数，并判断该本征模态函数是否满足预先设定的停止标准，若不满足，则将该本征模态函数作为输入，若满足，则将均值曲面作为输入，重复步骤B1至B5，直到提取出所有本征模态函数。
2.根据权利要求1所述的一种用于多图像融合的二维多元信号经验模态快速分解方 法，其特征在于，所述步骤B3中，在得到的极值点的基础上构造德劳内三角，对于一个德劳 内三角，使用重心作为均值点
Figure 565731DEST_PATH_IMAGE001
，均值点坐标为
Figure 833901DEST_PATH_IMAGE002
，则：
Figure 622866DEST_PATH_IMAGE003
Figure 52710DEST_PATH_IMAGE004
Figure 559915DEST_PATH_IMAGE005
其中，德劳内三角的三个顶点分别为
Figure 366197DEST_PATH_IMAGE006
，
Figure 9668DEST_PATH_IMAGE007
，
Figure 875992DEST_PATH_IMAGE008
，顶点坐标分别为
Figure 604914DEST_PATH_IMAGE009
、
Figure 214887DEST_PATH_IMAGE010
和
Figure 978444DEST_PATH_IMAGE011
；
Figure 750090DEST_PATH_IMAGE012
为均值点的像素值，
Figure 231887DEST_PATH_IMAGE013
为德劳内三角的顶点的像素值，
Figure 379972DEST_PATH_IMAGE014
表示三个顶点的权重。
3.根据权利要求2所述的一种用于多图像融合的二维多元信号经验模态快速分解方 法，其特征在于，采用等权重划分方式时，
Figure 998035DEST_PATH_IMAGE015
；
采用基于欧式距离的权重划分方式时，
Figure 675004DEST_PATH_IMAGE016
Figure 202737DEST_PATH_IMAGE017
其中
Figure 154512DEST_PATH_IMAGE018
表示均值点到第p个顶点的欧式距离。
4.根据权利要求1所述的一种用于多图像融合的二维多元信号经验模态快速分解方 法，其特征在于，所述步骤B6中，停止标准为：两次相邻迭代运算的柯西标准差
Figure 627082DEST_PATH_IMAGE019
和
Figure 6111DEST_PATH_IMAGE020
之间的差异小于事先设定的阈值时，停止迭代，即
Figure 196920DEST_PATH_IMAGE021
其中，柯西标准差
Figure 686808DEST_PATH_IMAGE022
，M、N是图像的高和宽；
Figure 279463DEST_PATH_IMAGE023
代表第x次筛选的结果；
Figure 829393DEST_PATH_IMAGE024
为事先设定的阈值。
5.根据权利要求1所述的一种用于多图像融合的二维多元信号经验模态快速分解方法，其特征在于，所述步骤B5中，
Figure 241920DEST_PATH_IMAGE025
其中，
Figure 535498DEST_PATH_IMAGE026
为第k个方向上的均值曲面，A为二维多元信号的均值曲面。
6.根据权利要求1所述的一种用于多图像融合的二维多元信号经验模态快速分解方法，其特征在于，所述步骤C中，对同一空间尺度的本征模态函数进行融合的步骤包括：
步骤C1、计算每幅图像的每个本征模态函数的局部方差或局部能量；
步骤C2、计算每幅图像的每个像素点对应的权重矩阵；
步骤C3、通过局部方差和权重矩阵计算该空间尺度下融合后的每个本征模态函数；
步骤C4、将所有空间尺度下融合后的本征模态函数进行叠加，即可得到融合后的图像。
CN202010465445.2A 2020-05-28 2020-05-28 用于多图像融合的二维多元信号经验模态快速分解方法 Active CN111369488B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010465445.2A CN111369488B (zh) 2020-05-28 2020-05-28 用于多图像融合的二维多元信号经验模态快速分解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010465445.2A CN111369488B (zh) 2020-05-28 2020-05-28 用于多图像融合的二维多元信号经验模态快速分解方法

Publications (2)

Publication Number Publication Date
CN111369488A CN111369488A (zh) 2020-07-03
CN111369488B true CN111369488B (zh) 2020-11-10

Family

ID=71212298

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010465445.2A Active CN111369488B (zh) 2020-05-28 2020-05-28 用于多图像融合的二维多元信号经验模态快速分解方法

Country Status (1)

Country Link
CN (1) CN111369488B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113139557B (zh) * 2021-05-10 2024-03-29 东南大学 一种基于二维多元经验模态分解的特征提取方法
CN113298890B (zh) * 2021-05-14 2022-07-15 之江实验室 非尺度混叠及边缘保留的图像多尺度分解方法及调色方法
CN113742802A (zh) * 2021-09-03 2021-12-03 国网经济技术研究院有限公司 用于工程图纸融合的二维多元信号经验模态分解快速方法
CN114372943A (zh) * 2021-12-30 2022-04-19 北京理工大学珠海学院 一种图像融合的方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110137980A (zh) * 2019-04-10 2019-08-16 国网辽宁省电力有限公司电力科学研究院 一种基于Hilbert-Hung和MEMD的电力系统低频振荡模式辨识方法
CN110287446A (zh) * 2019-05-29 2019-09-27 东南大学 一种快速二维多元经验模态分解算法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150293143A1 (en) * 2013-07-21 2015-10-15 National Chiao Tung University Feature Extraction from Human Gaiting Patterns using Principal Component Analysis and Multivariate Empirical Mode Decomposition
CN109767411B (zh) * 2018-12-27 2023-08-04 东南大学 一种用于多图像融合的二维多元经验模态分解算法
CN109659957B (zh) * 2019-01-10 2022-05-17 东北电力大学 基于apit-memd电力系统低频振荡模式辨识方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110137980A (zh) * 2019-04-10 2019-08-16 国网辽宁省电力有限公司电力科学研究院 一种基于Hilbert-Hung和MEMD的电力系统低频振荡模式辨识方法
CN110287446A (zh) * 2019-05-29 2019-09-27 东南大学 一种快速二维多元经验模态分解算法

Also Published As

Publication number Publication date
CN111369488A (zh) 2020-07-03

Similar Documents

Publication Publication Date Title
CN111369488B (zh) 用于多图像融合的二维多元信号经验模态快速分解方法
US9652890B2 (en) Methods and systems of generating an anatomically-constrained local model for performance capture
Shen et al. Layer depth denoising and completion for structured-light rgb-d cameras
US9639737B2 (en) Methods and systems of performing performance capture using an anatomically-constrained local model
Bai et al. Video snapcut: robust video object cutout using localized classifiers
US8929610B2 (en) Methods and apparatus for robust video stabilization
US20160321838A1 (en) System for processing a three-dimensional (3d) image and related methods using an icp algorithm
KR101195942B1 (ko) 카메라 보정 방법 및 이를 이용한 3차원 물체 재구성 방법
EP3367334B1 (en) Depth estimation method and depth estimation apparatus of multi-view images
Liu et al. Hybrid neural fusion for full-frame video stabilization
US20060193535A1 (en) Image matching method and image interpolation method using the same
CN106934827A (zh) 三维场景的重建方法和装置
US20060228002A1 (en) Simultaneous optical flow estimation and image segmentation
CN103702098A (zh) 一种时空域联合约束的三视点立体视频深度提取方法
Turkulainen et al. DN-Splatter: Depth and Normal Priors for Gaussian Splatting and Meshing
JP5561786B2 (ja) ３次元形状モデル高精度化方法およびプログラム
CN111460741B (zh) 一种基于数据驱动的流体模拟方法
Wang et al. Depth map super-resolution based on dual normal-depth regularization and graph Laplacian prior
CN107424122A (zh) 一种大位移下形变辅助的图像插补方法
US12073576B2 (en) Apparatus and method for generating depth map from multi-view image
Heimann et al. Joint Geometry and Attribute Upsampling of Point Clouds Using Frequency-Selective Models with Overlapped Support
Zhao et al. 3dfill: Reference-guided image inpainting by self-supervised 3d image alignment
Akyazi et al. Graph-based inpainting of disocclusion holes for zooming in 3d scenes
Chang et al. SC-AOF: A Sliding Camera and Asymmetric Optical-Flow-Based Blending Method for Image Stitching
Lee et al. Guess The Unseen: Dynamic 3D Scene Reconstruction from Partial 2D Glimpses

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: A Fast Empirical Mode Decomposition Method for 2D Multivariate Signals Used in Multi image Fusion

Effective date of registration: 20231220

Granted publication date: 20201110

Pledgee: Zijin Branch of Nanjing Bank Co.,Ltd.

Pledgor: Jiangsu Jihui Mobile Communication Technology Research Institute Co.,Ltd.

Registration number: Y2023980072754

PE01 Entry into force of the registration of the contract for pledge of patent right