CN112862683B - 一种基于弹性配准和网格优化的邻接图像拼接方法 - Google Patents

一种基于弹性配准和网格优化的邻接图像拼接方法 Download PDF

Info

Publication number
CN112862683B
CN112862683B CN202110174293.5A CN202110174293A CN112862683B CN 112862683 B CN112862683 B CN 112862683B CN 202110174293 A CN202110174293 A CN 202110174293A CN 112862683 B CN112862683 B CN 112862683B
Authority
CN
China
Prior art keywords
image
point
scale
adjacent
grid
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
CN202110174293.5A
Other languages
English (en)
Other versions
CN112862683A (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.)
Tongji University
Original Assignee
Tongji University
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 Tongji University filed Critical Tongji University
Priority to CN202110174293.5A priority Critical patent/CN112862683B/zh
Publication of CN112862683A publication Critical patent/CN112862683A/zh
Application granted granted Critical
Publication of CN112862683B publication Critical patent/CN112862683B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/40Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4038Image mosaicing, e.g. composing plane images from plane sub-images
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/20Image enhancement or restoration using local operators
    • G06T5/30Erosion or dilatation, e.g. thinning
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing

Landscapes

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

Abstract

本发明公开了一种基于弹性配准和网格优化的邻接图像拼接方法,首先使用SIFT算法进行特征提取与匹配,通过序列RANSAC算法获取特征匹配内点,计算多平面单应性;接着对无人机图像划分网格,使用基于弹性模型的方法对相邻的无人机图像进行配准;然后针对网格顶点坐标集构造四个约束项来建立能量函数,通过最小化能量函数进行网格优化,得到变形网格顶点;最后通过三角形纹理映射、最佳缝合线和多通道融合算法的处理步骤后得到高分辨率的无人机图像拼接结果。本发明实验结果表明,相比于传统方法,能有效消除拼接重影、不对齐现象,具有一定的视差容忍度,并且能够减少多图拼接产生的失真,保持图像形状,具有自然观感。

Description

一种基于弹性配准和网格优化的邻接图像拼接方法
技术领域
本发明涉及一种基于弹性配准和网格优化的邻接图像拼接方法,属于图像智能化处理技术领域。
背景技术
目前,随着无人机技术的发展,无人机在各类场合下的应用越来越广泛。无人机遥感是一种低空遥感,具有图像获取快速、定位准确、操作简单等优点,与航天航空遥感相比,空间分辨率更高,成本更低。由于无人机航高、焦距和视角的限制,单张无人机图像很难反映整个待测区域的情况,为了扩大视场,需要依靠图像拼接技术将多幅航拍图像融合成具有宽视角的、同时保持地面分辨率的全景图像。
图像拼接是将多幅图像拼接成具有更宽视场的全景图的过程。图像拼接包括三个阶段:特征提取与匹配、图像配准和图像合成。传统的图像拼接方法是利用全局单应性对两幅图像进行对齐,如仿射或投影变换,代表方法是AutoStitch。该方法假设拍摄场景处于同一平面,或者图像是围绕摄像机投影中心旋转拍摄,即输入图像之间是包含很少或没有视差的。在此假设下,全局单应性取得了很好的效果,但当违反上述成像假设时,很可能会出现误对准伪影。当这种情况发生时,这些方法尝试使用后处理图像混合来隐藏错位的区域,但当视差存在或较大时,该类方法仍会失效。
针对传统全局拼接算法的上述缺陷,现有技术采用了一种局部变换模型,如平滑变化仿射拼接(SVA)、尽可能投影的移动直接线性变换(APAP)算法、鲁棒弹性变换(REW)算法等。这些方法基于网格变形,利用局部单应性对图像进行配准,能够处理一定的视差场景。在此基础上,对网格施加不同的约束可以达到不同的效果,如保形半投影变换拼接(SPHP)算法、尽可能自然的图像拼接(AANAP)算法、具有全局相似性的自然图像拼接(NISwGSP)算法,这些算法对网格施加不同的限制,如多平面单应性、线性度等,达到了减少投影失真的效果,增加图像的自然拼接度。但是,对于无人机图像拼接,由于其分辨率较高、地形复杂等,更容易产生拼接不齐、地形畸变等问题,直接运用现有技术仍不能达到预期的效果,因此,本领域技术人员急需要设计出新的图像拼接方法。
发明内容
目的:为了克服现有技术中存在的不足,本发明提供一种基于弹性配准和网格优化的邻接图像拼接方法,使用基于薄板样条的弹性模型进行邻接图像的配准,使用网格优化方法进行全局形状的保持,解决了无人机图像拼接过程中存在的不对齐伪影和形状失真问题。
技术方案:为解决上述技术问题,本发明采用的技术方案为:
一种基于弹性配准和网格优化的邻接图像拼接方法,包括如下步骤:
将原始邻接图像降采样至seam_scale尺度,得到seam_scale尺度图像,将seam_scale尺度和work_scale尺度求比值,得到seam_work_aspect,seam_scale尺度小于work_scale尺度;
将变形网格顶点乘以seam_work_aspect得到seam_scale尺度图像新的网格顶点坐标,根据seam_scale尺度图像新的网络顶点坐标利用三角形仿射变换方法对seam_scale尺度图像进行纹理映射得到变换图像;
对变换图像执行基于图割法的最佳缝合线算法,得到缝合线两侧图像的缝合线掩码;
将原始邻接图像采样至compose_scale尺度,得到compose_scale尺度图像,将compose_scale尺度和work_scale尺度求比值,得到compose_work_aspect,compose_scale尺度大于work_scale尺度;
将变形网格顶点乘以compose_work_aspect,得到compose_scale尺度图像新的网格顶点坐标,根据compose_scale尺度图像新的网格顶点坐标利用三角形仿射变换方法对compose_scale尺度图像进行纹理映射得到高分辨率的变换图像;
将缝合线掩码膨胀放大至compose_scale尺度,在compose_scale尺度上针对高分辨率的变换图像和膨胀放大后的缝合线掩码执行多通道融合算法,得到拼接结果图。
作为优选方案,所述变形网格顶点获取步骤如下:
对构造的网格优化能量函数E(V),使用稀疏线性求解器进行优化求解,得到每幅图像相对于参考图像的变形网格顶点坐标,再通过归一化得到每幅图像的变形网格顶点;
所述网格优化能量函数E(V)计算公式如下:
E(V)=Ea(V)+λlsEls(V)+Egs(V)+λlEl(V)
Figure BDA0002940422410000031
Figure BDA0002940422410000032
Figure BDA0002940422410000033
Figure BDA0002940422410000034
其中,Ea为对齐项、Els为局部相似项、Egs为全局相似项、El为直线保持项,λls和λl为对应的权重系数,V为所有图像的网格顶点坐标集,N为所有图像的数量;J为邻接图像之间的邻接关系;Mij为网格匹配点集合;
对齐项Ea(V)中,
Figure BDA0002940422410000035
为网格顶点
Figure BDA0002940422410000036
或网格对应点
Figure BDA0002940422410000037
所在网格的四个顶点的双线性插值;
局部相似项Els(V)和全局相似项Egs(V)中,Ei为Ii的所有边的集合,
Figure BDA0002940422410000038
Figure BDA0002940422410000039
表示原始图像的某条边及其变形后的边;
Figure BDA00029404224100000310
为边
Figure BDA00029404224100000311
经历的相似变换,
Figure BDA00029404224100000312
Figure BDA00029404224100000313
为相似变换
Figure BDA00029404224100000314
中的对应元素,表示为顶点变量的线性组合,si、θi为最佳尺度和最佳旋转角度,
Figure BDA00029404224100000315
为边
Figure BDA00029404224100000316
的权重函数;
直线保持项El(V)中,Li表示图像Ii中的直线集合,lu、lv、lk分别为直线
Figure BDA00029404224100000317
上的起点、终点和中间采样点,u为直线局部坐标系下的1维坐标,
Figure BDA00029404224100000318
为直线起点、终点和中间采样点所在网格的四个顶点的双线性插值。
作为优选方案,所述si、θi最佳尺度和最佳旋转角度获取的步骤如下:
从多平面单应性矩阵H中估算焦距的初值并分别形成Ii、Ij的内参矩阵Ki、Kj,通过以下公式获得Ii和Ij之间的3D旋转Rij的初始估计:
Figure BDA0002940422410000041
其中,Ii、Ij为相邻接的两张图像,R代表3D旋转矩阵的参数;
初始化之后,对所有的Ki和3D旋转Rij初始值执行光束平差法得到得每幅邻接图像Ii精细化的焦距fi和3D旋转Ri
每幅邻接图像的最佳尺度计算公式如下:
si=f1/fi
其中,f1为参考图像的焦距;
使用LSD检测图像的直线,通过弹性配准可得到两幅邻接图像Ii和Ij之间线的对应关系,每对线段对唯一地确定相对旋转角度,再按照RANSAC算法投票筛选得到每幅邻接图像的最佳旋转角度θi
作为优选方案,所述Mij网格匹配点集合获取步骤如下:
根据特征匹配内点集
Figure BDA0002940422410000042
和多平面单应性变换矩阵H,通过下述两式计算Ij上的特征目标点qi在Ii上的对应投影点q′i以及投影点q′i和特征源点pi的偏差ri
Figure BDA0002940422410000043
Figure BDA0002940422410000044
其中,
Figure BDA0002940422410000045
Figure BDA0002940422410000046
为投影点q′i在X方向和Y方向上的分量,
Figure BDA0002940422410000047
Figure BDA0002940422410000048
为投影点误差ri在X方向和Y方向上的分量;i是第i个特征匹配内点;
构造用于计算最优变形的能量函数,使得能够通过投影点误差拟合图像平面Ii上的任意像素坐标x=(x,y)T在X方向上的变形g(x,y)和Y方向上变形h(x,y),x、y代表X方向的坐标和Y方向的坐标。对于X方向上的变形g(x,y),能量函数如下:
Jλ=JD+λJS
Figure BDA0002940422410000049
Figure BDA00029404224100000410
其中,Jλ是变形能量函数,JD是对齐项,JS是平滑项,λ是权重系数;
Figure BDA00029404224100000411
表示投影点q′i在X方向上经历的变形,
Figure BDA00029404224100000412
为投影点误差ri在X方向上的分量。i是第i个特征匹配内点,n代表特征匹配内点的数量;
同理,对于Y方向上的变形h(x,y),能量函数如下:
Jλ=JD+λJS
Figure BDA0002940422410000051
Figure BDA0002940422410000052
其中,Jλ是变形能量函数,JD是对齐项,JS是平滑项,λ是权重系数;
Figure BDA0002940422410000053
表示投影点q′i在Y方向上经历的变形,
Figure BDA0002940422410000054
为投影点误差ri在Y方向上的分量。
根据薄板样条理论,通过最小化Jλ分别得到变形函数g(x,y)和h(x,y)唯一的解析解:
Figure BDA0002940422410000055
其中,di=‖x-q′i‖表示任意像素坐标x=(x,y)T与第i个投影点q′i的距离,n为投影点的数量;
上述公式具有2(n+3)个待求解系数
Figure BDA0002940422410000056
α=(α1,α2,α3)T、β=(β1,β23)T,通过求解下述矩阵方程得到:
Figure BDA0002940422410000057
其中,
Figure BDA0002940422410000058
dij=‖q′j-q′i‖,C是权重系数,I为单位矩阵;Q=(q′1,…,q′n)T表示齐次投影点,n为投影点的数量,
Figure BDA0002940422410000059
Figure BDA00029404224100000510
为投影点误差在X方向和Y方向上的分量;
在求得变形函数之后,对于Ii上的网格顶点
Figure BDA00029404224100000511
通过加上形变后再进行多平面单应性变换,得到其在Ij上的对应点
Figure BDA00029404224100000512
Figure BDA00029404224100000513
如果网格顶点vi及其对应点vi′落在Ii和Ij的重叠区域中,则收集vi形成网格匹配点集合Mij,Mij中元素
Figure BDA00029404224100000514
作为网格顶点,
Figure BDA00029404224100000515
为网格顶点
Figure BDA00029404224100000516
在Ij上投影点。
作为优选方案,所述特征匹配内点集
Figure BDA00029404224100000517
获取步骤如下:
1-1)对邻接关系J中一对邻接图像Ii和Ij,通过SIFT算法检测每幅邻接图像的特征,通过2NN算法进行匹配获取每幅邻接图像初始匹配点集;
1-2)利用RANSAC算法从初始匹配点集中估算出一个单应性矩阵对应的内点集;
1-3)从初始匹配点集中提取单应性矩阵对应的内点集,将剩余的匹配点再次进行RANSAC算法估算出另一个单应性矩阵对应的内点集;
1-4)重复步骤1-3)多次,直至剩余的匹配点数量小于设定阈值;
1-5)将每一步估算提取的单应性矩阵对应的内点集合并成特征匹配内点集
Figure BDA0002940422410000061
其中,pi为Ii上的第i个特征源点,qi为Ij上的第i个特征目标点,n为特征匹配内点的数量。
作为优选方案,所述多平面单应性矩阵H获取步骤如下:
通过特征匹配内点集
Figure BDA0002940422410000062
用直接线性变换方法计算出邻接图像之间的多平面单应性变换矩阵H。
作为优选方案,V所有图像的网格顶点坐标集获取步骤如下:
输入N张邻接图像I1,I2,...,IN,图像数量N≥2;获取邻接图像之间的邻接关系J;
通过邻接图像和邻接图像之间的邻接关系构建全局匹配图
以work_scale尺度对邻接图像进行降采样,对降采样图像以分割像素大小划分网格,获取每张图片的网格顶点坐标
Figure BDA0002940422410000063
其中vi表示第i张图片的一个顶点坐标,m为顶点个数;
所有图片的网格顶点坐标集为V=(V1,…,VN)。
作为优选方案,所述参考图像采用一组邻接图像的第一张图像。
作为优选方案,所述邻接图像通过无人机拍摄获得。
有益效果:本发明提供的一种基于弹性配准和网格优化的邻接图像拼接方法,通过基于薄板样条的弹性模型进行邻接图像的配准,提升重叠区域的对齐精度;通过构造对齐项、相似项、直线保持项等约束项建立网格优化框架,保持图像的原始形状。实验结果表明,该方法具有鲁棒、高效的特点,适用于无人机图像拼接场景。
附图说明
图1是本发明的整体算法框架的示意图。
图2是本发明的弹性配准流程的示意图。
图3是本发明的对比实验结果的示意图。
具体实施方式
下面结合具体实施例对本发明作更进一步的说明。
如图1所示,一种基于弹性配准和网格优化的邻接图像拼接方法,主要包括以下步骤:
(1)输入图像和邻接图;
(2)特征检测与匹配;
(3)弹性模型配准生成网格匹配点;
(4)全局相似性估计;
(5)构造约束项网格优化;
(6)纹理映射合成图像。
S1.获取一组邻接图像及图像之间的邻接关系,选取参考图像。对邻接图像进行降采样,并划分网格得到网格顶点坐标集;
S2.对每幅邻接图像,使用SIFT算法提取特征,对特征进行匹配,获得初始匹配点集,对初始匹配点集使用序列RANSAC算法滤除错误匹配点,得到特征匹配内点集,利用特征匹配内点集计算邻接图像之间的多平面单应性变换矩阵;
S3.利用基于薄板样条的弹性模型进行邻接图像的配准,将邻接图像精确对齐,生成网格匹配点集合;
S4.利用S2求得的多平面单应性变换矩阵和S3求得的网格匹配点估计每幅邻接图像的焦距和3D旋转,选取每幅图像相对于参考图像的最佳尺度和旋转;
S5.针对网格顶点坐标集构建网格顶点能量函数,该能量函数包含四个约束项:对齐项、局部相似项、全局相似项、直线保持项,通过稀疏线性求解器进行网格优化求解,得到相对于参考图像的变形网格顶点;
S6.利用S5求得的变形网格顶点,通过三角形纹理映射、最佳缝合线和多通道融合的处理步骤后得到最终拼接结果。
实施例1:
步骤(1)中,输入N张用无人机拍摄的邻接图像I1,I2,...,IN,图像数量N≥2;获取邻接图像之间的邻接关系J,可由无人机航线规划方法获取,通过邻接图像和邻接图像之间的邻接关系构建全局匹配图;在不失一般性的情况下,使用I1作为参考图像。由于无人机原始图像分辨率较高,以800×600像素大小对邻接图像进行降采样,对降采样图像以40×40像素大小划分网格,获取每张图片的网格顶点坐标
Figure BDA0002940422410000081
其中vi表示第i张图片的一个顶点坐标,m为顶点个数。所有图片的网格顶点坐标集为V=(V1,…,VN)。记降采样图像相对于原始图像的缩小尺度为work_scale。在接下来的(1)-(5)步骤过程中,均在此work_scale尺度上对图像进行操作。
步骤(2)中,对邻接关系J中一对邻接图像Ii和Ij,通过SIFT算法检测每幅邻接图像的特征,通过2NN算法进行匹配获取每幅邻接图像初始匹配点集。由于经典RANSAC算法在滤除错误匹配点时,只计算出邻接图像中一个平面的单应性,牺牲了大量正确的特征匹配内点,因此该处使用序列RANSAC算法获取特征匹配内点集,具体方法如下:
(2-1)利用RANSAC算法从初始匹配点集中估算出一个单应性矩阵对应的内点集;
(2-2)从初始匹配点集中提取单应性矩阵对应的内点集,将剩余的匹配点再次进行RANSAC算法估算出另一个单应性矩阵对应的内点集;
(2-3)重复步骤(2-2)多次,直至剩余的匹配点数量小于设定阈值,设定阈值设置为40;
(2-4)将每一步估算提取的单应性矩阵对应的内点集合并成特征匹配内点集
Figure BDA0002940422410000082
其中,pi为Ii上的第i个特征源点,qi为Ij上的第i个特征目标点,n为特征匹配内点的数量。
通过最终得到的特征匹配内点集
Figure BDA0002940422410000083
用直接线性变换方法计算出邻接图像之间的多平面单应性变换矩阵H。
如图2所示,步骤(3)基于薄板样条的弹性模型进行邻接图像配准的流程,由于多平面单应性仍然难以将Ii和Ij的所有像素位置对齐,因此设计该步骤的目的就是利用基于薄板样条的弹性模型将两幅图像更精确地配准,具体步骤如下:
(3-1)根据步骤(2)求得的特征匹配内点集
Figure BDA0002940422410000084
和多平面单应性变换矩阵H,通过下述两式计算Ij上的特征目标点qi在Ii上的对应投影点q′i以及投影点q′i和特征源点pi的偏差ri
Figure BDA0002940422410000091
Figure BDA0002940422410000092
其中,
Figure BDA0002940422410000093
Figure BDA0002940422410000094
为投影点q′i在X方向和Y方向上的分量,
Figure BDA0002940422410000095
Figure BDA0002940422410000096
为投影点误差ri在X方向和Y方向上的分量。i是第i个特征匹配内点。
(3-2)构造用于计算最优变形的能量函数,使得能够通过投影点误差拟合图像平面Ii上的任意像素坐标x=(x,y)T在X方向上的变形g(x,y)和Y方向上变形h(x,y),x、y代表X方向的坐标和Y方向的坐标。对于X方向上的变形g(x,y),能量函数如下:
Jλ=JD+λJS
Figure BDA0002940422410000097
Figure BDA0002940422410000098
其中,Jλ是变形能量函数,JD是对齐项,JS是平滑项,λ是权重系数,用以控制平滑程度。
Figure BDA0002940422410000099
表示投影点q′i在X方向上经历的变形,
Figure BDA00029404224100000910
为投影点误差ri在X方向上的分量。i是第i个特征匹配内点,n代表特征匹配内点的数量。
同理,对于Y方向上的变形h(x,y),能量函数如下:
Jλ=JD+λJS
Figure BDA00029404224100000911
Figure BDA00029404224100000912
其中,Jλ是变形能量函数,JD是对齐项,JS是平滑项,λ是权重系数,用以控制平滑程度。
Figure BDA00029404224100000913
表示投影点q′i在Y方向上经历的变形,
Figure BDA00029404224100000914
为投影点误差ri在Y方向上的分量。
(3-3)根据薄板样条理论,通过最小化Jλ分别得到变形函数g(x,y)和h(x,y)唯一的解析解:
Figure BDA0002940422410000101
其中,di=‖x-q′i‖表示任意像素坐标x=(x,y)T与第i个投影点q′i的距离,n为投影点的数量。
上述公式具有2(n+3)个待求解系数
Figure BDA0002940422410000102
α=(α123)T、β=(β123)T,可以通过求解下述矩阵方程得到:
Figure BDA0002940422410000103
其中,
Figure BDA0002940422410000104
dij=‖q′j-q′i‖,C是权重系数,I为单位矩阵;Q=(q′1,…,q′n)T表示齐次投影点,n为投影点的数量,
Figure BDA0002940422410000105
Figure BDA0002940422410000106
为投影点误差在X方向和Y方向上的分量。
(3-4)在求得变形函数之后,对于Ii上的网格顶点
Figure BDA0002940422410000107
通过加上形变后再进行多平面单应性变换,得到其在Ij上的对应点
Figure BDA0002940422410000108
Figure BDA0002940422410000109
如果网格顶点vi及其对应点vi′落在Ii和Ij的重叠区域中,则收集vi形成网格顶点
Figure BDA00029404224100001010
作为网格匹配点集合Mij中元素。网格顶点
Figure BDA00029404224100001011
在Ij上对应点为
Figure BDA00029404224100001012
参考图2,(a)为步骤(2-4)求得的Ii与Ij的特征匹配内点集,(b)为步骤(3-1)计算得到的Ij在Ii上的投影点误差,(c)为步骤(3-3)根据投影点误差拟合的网格形变量,(d)为步骤(3-4)计算得到的均匀分布的网格匹配点。
步骤(4)中,由于适当的尺度和旋转能够保持图像形状,因此设计该步骤来选取每幅图像相对于参考图像的最佳尺度和旋转。具体步骤如下:
(4-1)估计焦距和3D旋转。从步骤(2)计算出的多平面单应性矩阵H中估算焦距的初值并分别形成Ii、Ij的内参矩阵Ki、Kj,通过以下公式获得Ii和Ij之间的3D旋转Rij的初始估计:
Figure BDA00029404224100001013
其中,R代表3D旋转矩阵的参数。该式可通过SVD分解求解。相比于AutoStitch的初始化,该式使用更好的初始化单应性变换和更均匀分布的网格匹配点。在初始化之后,对所有的Ki和3D旋转Rij初始值执行光束平差法得到得每幅邻接图像Ii精细化的焦距fi和3D旋转Ri
(4-2)选择尺度si。每幅邻接图像的尺度可设置为:
si=f1/fi (9)
其中,f1为参考图像的焦距。
(4-3)选择旋转θi。使用LSD检测图像的直线,通过弹性配准可得到两幅邻接图像Ii和Ij之间线的对应关系,每对线段对唯一地确定相对旋转角度,再按照RANSAC算法投票筛选得到每幅邻接图像的最佳旋转角度θi
步骤(5)中所述网络顶点能量函数如式(10)-(14):
E(V)=Ea(V)+λlsEls(V)+Egs(V)+λlEl(V) (10)
Figure BDA0002940422410000111
Figure BDA0002940422410000112
Figure BDA0002940422410000113
Figure BDA0002940422410000114
其中,Ea为对齐项、Els为局部相似项、Egs为全局相似项、El为直线保持项,λls和λl为对应的权重系数,V为所有图片的网格顶点坐标集,N为所有图片的数量;
对齐项Ea(V)中,
Figure BDA0002940422410000115
为网格顶点或网格顶点对应点
Figure BDA0002940422410000116
Figure BDA0002940422410000117
所在网格的四个顶点的双线性插值;
局部相似项Els(V)和全局相似项Egs(V)中,Ei为Ii的所有边的集合,
Figure BDA0002940422410000118
Figure BDA0002940422410000119
表示原始图片的某条边及其变形后的边。
Figure BDA00029404224100001110
为边
Figure BDA00029404224100001111
经历的相似变换,
Figure BDA0002940422410000121
Figure BDA0002940422410000122
为相似变换
Figure BDA0002940422410000123
中的对应元素,可表示为顶点变量的线性组合,si、θi为由步骤S4得到的最佳尺度和旋转,
Figure BDA0002940422410000124
为边
Figure BDA0002940422410000125
的权重函数;
直线保持项El(V)中,Li表示图像Ii中的直线集合,lu、lv、lk分别为直线
Figure BDA0002940422410000126
上的起点、终点和中间采样点,u为直线局部坐标系下的1维坐标,
Figure BDA0002940422410000127
为直线采样点所在网格的四个顶点的双线性插值。在直线采样过程中,基于网格优化算法,选取直线的长度大于60像素。限于LSD提取到的直线可能较为片段化,本发明还提供了一种交互式方式,可以由用户对需要保护的直线进行手动提取;
对构造的网格优化能量函数E(V),使用稀疏线性求解器进行优化求解,得到每幅图像相对于参考图像的变形网格顶点坐标,再通过归一化得到每幅图像的变形网格顶点。
以上(1)-(5)步骤的图像操作均在800×600像素的work_scale尺度上进行,在得到该尺度下的变形网格顶点坐标后,还需进行步骤(6)中的以下处理步骤:
(6-1)将原始邻接图像降采样至seam_scale尺度得到seam_scale尺度图像,seam_scale尺度小于work_scale尺度,得到seam_scale和work_scale之比seam_work_aspect,将步骤(5)得到的变形网格顶点结果乘以seam_work_aspect得到seam_scale尺度图像新的网格顶点坐标,根据seam_scale尺度图像新的网络顶点坐标利用三角形仿射变换方法对seam_scale尺度图像进行纹理映射得到变换图像;对变换图像执行基于图割法的最佳缝合线算法,得到缝合线两侧图像的缝合线掩码;
(6-2)将原始邻接图像采样至compose_scale尺度得到compose_scale尺度图像,compose_scale尺度大于work_scale尺度,得到compose_scale与work_scale之比compose_work_aspect,将步骤(5)得到的变形网格顶点结果乘以compose_work_aspect,得到compose_scale尺度图像新的网格顶点坐标,根据compose_scale尺度图像新的网格顶点坐标利用三角形仿射变换方法对compose_scale尺度图像进行纹理映射得到高分辨率的变换图像;将(6-1)得到的缝合线掩码膨胀放大至compose_scale尺度,在compose_scale尺度上针对高分辨率的变换图像和膨胀放大后的缝合线掩码执行多通道融合算法,得到最终完整清晰、高分辨率的无人机拼接结果图。
一种基于弹性配准和网格优化的邻接图像拼接方法,其一组对比实验结果如图3所示。其中,图(a)为AutoStitch方法得到的拼接结果,图(b)为本发明得到的拼接结果,图(c)为网格优化得到的变形网格顶点。可见有益效果如下:在重叠区域能将图像精确对齐,去除拼接断裂、重影现象;图像在全局范围内保持自然的形状。
本发明实验结果表明,相比于传统方法,能有效消除拼接重影、不对齐现象,具有一定的视差容忍度,并且能够减少多图拼接产生的失真,保持图像形状,具有自然观感。
以上所述仅是本发明的优选实施方式,应当指出:对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。

Claims (3)

1.一种基于弹性配准和网格优化的邻接图像拼接方法,其特征在于:包括如下步骤:
将原始邻接图像降采样至seam_scale尺度,得到seam_scale尺度图像,将seam_scale尺度和work_scale尺度求比值,得到seam_work_aspect,seam_scale尺度小于work_scale尺度;
将变形网格顶点乘以seam_work_aspect得到seam_scale尺度图像新的网格顶点坐标,根据seam_scale尺度图像新的网络顶点坐标利用三角形仿射变换方法对seam_scale尺度图像进行纹理映射得到变换图像;
对变换图像执行基于图割法的最佳缝合线算法,得到缝合线两侧图像的缝合线掩码;
将原始邻接图像采样至compose_scale尺度,得到compose_scale尺度图像,将compose_scale尺度和work_scale尺度求比值,得到compose_work_aspect,compose_scale尺度大于work_scale尺度;
将变形网格顶点乘以compose_work_aspect,得到compose_scale尺度图像新的网格顶点坐标,根据compose_scale尺度图像新的网格顶点坐标利用三角形仿射变换方法对compose_scale尺度图像进行纹理映射得到高分辨率的变换图像;
将缝合线掩码膨胀放大至compose_scale尺度,在compose_scale尺度上针对高分辨率的变换图像和膨胀放大后的缝合线掩码执行多通道融合算法,得到拼接结果图;
所述变形网格顶点获取步骤如下:
对构造的网格优化能量函数E(V),使用稀疏线性求解器进行优化求解,得到每幅图像相对于参考图像的变形网格顶点坐标,再通过归一化得到每幅图像的变形网格顶点;
所述网格优化能量函数E(V)计算公式如下:
E(V)=Ea(V)+λlsEls(V)+Egs(V)+λlEl(V)
Figure FDA0003911207880000011
Figure FDA0003911207880000021
Figure FDA0003911207880000022
Figure FDA0003911207880000023
其中,Ea为对齐项、Els为局部相似项、Egs为全局相似项、El为直线保持项,λls和λl为对应的权重系数,V为所有图像的网格顶点坐标集,N为所有图像的数量;J为邻接图像之间的邻接关系;Mij为网格匹配点集合;
对齐项Ea(V)中,
Figure FDA0003911207880000024
为网格顶点
Figure FDA0003911207880000025
或网格对应点
Figure FDA0003911207880000026
所在网格的四个顶点的双线性插值;
局部相似项Els(V)和全局相似项Egs(V)中,Ei为Ii的所有边的集合,
Figure FDA0003911207880000027
Figure FDA0003911207880000028
表示原始图像的某条边及其变形后的边;
Figure FDA0003911207880000029
为边
Figure FDA00039112078800000210
经历的相似变换,
Figure FDA00039112078800000211
Figure FDA00039112078800000212
Figure FDA00039112078800000213
为相似变换
Figure FDA00039112078800000214
中的对应元素,表示为顶点变量的线性组合,si、θi为最佳尺度和最佳旋转角度,
Figure FDA00039112078800000215
为边
Figure FDA00039112078800000216
的权重函数;
直线保持项El(V)中,Li表示图像Ii中的直线集合,lu、lv、lk分别为直线
Figure FDA00039112078800000217
上的起点、终点和中间采样点,u为直线局部坐标系下的1维坐标,
Figure FDA00039112078800000218
为直线起点、终点和中间采样点所在网格的四个顶点的双线性插值;
所述si、θi最佳尺度和最佳旋转角度获取的步骤如下:
从多平面单应性矩阵H中估算焦距的初值并分别形成Ii、Ij的内参矩阵Ki、Kj,通过以下公式获得Ii和Ij之间的3D旋转Rij的初始估计:
Figure FDA00039112078800000219
其中,Ii、Ij为相邻接的两张图像,R代表3D旋转矩阵的参数;
初始化之后,对所有的Ki和3D旋转Rij初始值执行光束平差法得到得每幅邻接图像Ii精细化的焦距fi和3D旋转Ri
每幅邻接图像的最佳尺度计算公式如下:
si=f1/fi
其中,f1为参考图像的焦距;
使用LSD检测图像的直线,通过弹性配准可得到两幅邻接图像Ii和Ij之间线的对应关系,每对线段对唯一地确定相对旋转角度,再按照RANSAC算法投票筛选得到每幅邻接图像的最佳旋转角度θi
所述Mij网格匹配点集合获取步骤如下:
根据特征匹配内点集
Figure FDA0003911207880000031
和多平面单应性变换矩阵H,通过下述两式计算Ij上的特征目标点qi在Ii上的对应投影点q′i以及投影点q′i和特征源点pi的偏差ri
Figure FDA0003911207880000032
Figure FDA0003911207880000033
其中,
Figure FDA0003911207880000034
Figure FDA0003911207880000035
为投影点q′i在X方向和Y方向上的分量,
Figure FDA0003911207880000036
Figure FDA0003911207880000037
为投影点误差ri在X方向和Y方向上的分量;i是第i个特征匹配内点;
构造用于计算最优变形的能量函数,使得能够通过投影点误差拟合图像平面Ii上的任意像素坐标x=(x,y)T在X方向上的变形g(x,y)和Y方向上变形h(x,y),x、y代表X方向的坐标和Y方向的坐标;对于X方向上的变形g(x,y),能量函数如下:
Jλ=JD+λJS
Figure FDA0003911207880000038
Figure FDA0003911207880000039
其中,Jλ是变形能量函数,JD是对齐项,JS是平滑项,λ是权重系数;
Figure FDA00039112078800000310
表示投影点q′i在X方向上经历的变形,
Figure FDA00039112078800000311
为投影点误差ri在X方向上的分量;i是第i个特征匹配内点,n代表特征匹配内点的数量;
同理,对于Y方向上的变形h(x,y),能量函数如下:
Jλ=JD+λJS
Figure FDA00039112078800000312
Figure FDA0003911207880000041
其中,Jλ是变形能量函数,JD是对齐项,JS是平滑项,λ是权重系数;
Figure FDA0003911207880000042
表示投影点q′i在Y方向上经历的变形,
Figure FDA0003911207880000043
为投影点误差ri在Y方向上的分量;
根据薄板样条理论,通过最小化Jλ分别得到变形函数g(x,y)和h(x,y)唯一的解析解:
Figure FDA0003911207880000044
其中,di=‖x-q′i‖表示任意像素坐标x=(x,y)T与第i个投影点q′i的距离,n为投影点的数量;
上述公式具有2(n+3)个待求解系数
Figure FDA0003911207880000045
α=(α123)T、β=(β123)T,通过求解下述矩阵方程得到:
Figure FDA0003911207880000046
其中,
Figure FDA0003911207880000047
dij=‖q′j-q′i‖,C是权重系数,I为单位矩阵;Q=(q′1,…,q′n)T表示齐次投影点,n为投影点的数量,
Figure FDA0003911207880000048
Figure FDA0003911207880000049
为投影点误差在X方向和Y方向上的分量;
在求得变形函数之后,对于Ii上的网格顶点
Figure FDA00039112078800000410
通过加上形变后再进行多平面单应性变换,得到其在Ij上的对应点
Figure FDA00039112078800000411
Figure FDA00039112078800000412
如果网格顶点vi及其对应点vi′落在Ii和Ij的重叠区域中,则收集vi形成网格匹配点集合Mij,Mij中元素
Figure FDA00039112078800000413
作为网格顶点,
Figure FDA00039112078800000414
为网格顶点
Figure FDA00039112078800000415
在Ij上投影点;
所述特征匹配内点集
Figure FDA00039112078800000416
获取步骤如下:
1-1)对邻接关系J中一对邻接图像Ii和Ij,通过SIFT算法检测每幅邻接图像的特征,通过2NN算法进行匹配获取每幅邻接图像初始匹配点集;
1-2)利用RANSAC算法从初始匹配点集中估算出一个单应性矩阵对应的内点集;
1-3)从初始匹配点集中提取单应性矩阵对应的内点集,将剩余的匹配点再次进行RANSAC算法估算出另一个单应性矩阵对应的内点集;
1-4)重复步骤1-3)多次,直至剩余的匹配点数量小于设定阈值;
1-5)将每一步估算提取的单应性矩阵对应的内点集合并成特征匹配内点集
Figure FDA0003911207880000051
其中,pi为Ii上的第i个特征源点,qi为Ij上的第i个特征目标点,n为特征匹配内点的数量;
所述多平面单应性矩阵H获取步骤如下:
通过特征匹配内点集
Figure FDA0003911207880000052
用直接线性变换方法计算出邻接图像之间的多平面单应性变换矩阵H;
所有图像的网格顶点坐标集获取步骤如下:
输入N张邻接图像I1,I2,...,IN,图像数量N≥2;获取邻接图像之间的邻接关系J;
通过邻接图像和邻接图像之间的邻接关系构建全局匹配图;
以work_scale尺度对邻接图像进行降采样,对降采样图像以分割像素大小划分网格,获取每张图片的网格顶点坐标
Figure FDA0003911207880000053
其中vi表示第i张图片的一个顶点坐标,m为顶点个数;
所有图片的网格顶点坐标集为V=(V1,…,VN)。
2.根据权利要求1所述的一种基于弹性配准和网格优化的邻接图像拼接方法,其特征在于:所述参考图像采用一组邻接图像的第一张图像。
3.根据权利要求1所述的一种基于弹性配准和网格优化的邻接图像拼接方法,其特征在于:所述邻接图像通过无人机拍摄获得。
CN202110174293.5A 2021-02-07 2021-02-07 一种基于弹性配准和网格优化的邻接图像拼接方法 Active CN112862683B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110174293.5A CN112862683B (zh) 2021-02-07 2021-02-07 一种基于弹性配准和网格优化的邻接图像拼接方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110174293.5A CN112862683B (zh) 2021-02-07 2021-02-07 一种基于弹性配准和网格优化的邻接图像拼接方法

Publications (2)

Publication Number Publication Date
CN112862683A CN112862683A (zh) 2021-05-28
CN112862683B true CN112862683B (zh) 2022-12-06

Family

ID=75989340

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110174293.5A Active CN112862683B (zh) 2021-02-07 2021-02-07 一种基于弹性配准和网格优化的邻接图像拼接方法

Country Status (1)

Country Link
CN (1) CN112862683B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113506216B (zh) * 2021-06-24 2024-03-12 煤炭科学研究总院 一种全景图像拼接的快速缝合线寻优方法
CN114387153B (zh) * 2021-12-13 2023-07-04 复旦大学 一种用于插管机器人的视野拓展方法
CN114913064B (zh) * 2022-03-15 2024-07-02 天津理工大学 基于结构保持和多对多匹配的大视差图像拼接方法及装置
CN115393196B (zh) * 2022-10-25 2023-03-24 之江实验室 一种无人机面阵摆扫的红外多序列影像无缝拼接方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107730528A (zh) * 2017-10-28 2018-02-23 天津大学 一种基于grabcut算法的交互式图像分割与融合方法
CN109961398A (zh) * 2019-02-18 2019-07-02 鲁能新能源(集团)有限公司 风机叶片图像分割与网格优化拼接方法
CN110136090A (zh) * 2019-04-11 2019-08-16 中国地质大学(武汉) 以局部保持配准的鲁棒弹性模型无人机图像拼接方法
CN110428367A (zh) * 2019-07-26 2019-11-08 北京小龙潜行科技有限公司 一种图像拼接方法及装置
CN110781903A (zh) * 2019-10-12 2020-02-11 中国地质大学(武汉) 基于网格优化和全局相似性约束的无人机图像拼接方法
CN111915484A (zh) * 2020-07-06 2020-11-10 天津大学 基于密集匹配与自适应融合的参考图引导超分辨率方法
CN112308775A (zh) * 2020-09-23 2021-02-02 中国石油大学(华东) 水下图像拼接方法及装置

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107734268A (zh) * 2017-09-18 2018-02-23 北京航空航天大学 一种结构保持的宽基线视频拼接方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107730528A (zh) * 2017-10-28 2018-02-23 天津大学 一种基于grabcut算法的交互式图像分割与融合方法
CN109961398A (zh) * 2019-02-18 2019-07-02 鲁能新能源(集团)有限公司 风机叶片图像分割与网格优化拼接方法
CN110136090A (zh) * 2019-04-11 2019-08-16 中国地质大学(武汉) 以局部保持配准的鲁棒弹性模型无人机图像拼接方法
CN110428367A (zh) * 2019-07-26 2019-11-08 北京小龙潜行科技有限公司 一种图像拼接方法及装置
CN110781903A (zh) * 2019-10-12 2020-02-11 中国地质大学(武汉) 基于网格优化和全局相似性约束的无人机图像拼接方法
CN111915484A (zh) * 2020-07-06 2020-11-10 天津大学 基于密集匹配与自适应融合的参考图引导超分辨率方法
CN112308775A (zh) * 2020-09-23 2021-02-02 中国石油大学(华东) 水下图像拼接方法及装置

Also Published As

Publication number Publication date
CN112862683A (zh) 2021-05-28

Similar Documents

Publication Publication Date Title
CN112862683B (zh) 一种基于弹性配准和网格优化的邻接图像拼接方法
CN110211043B (zh) 一种用于全景图像拼接的基于网格优化的配准方法
CN110648398B (zh) 基于无人机航摄数据的正射影像实时生成方法及系统
CN104732482B (zh) 一种基于控制点的多分辨率图像拼接方法
KR101175097B1 (ko) 파노라마 영상 생성 방법
US20240169674A1 (en) Indoor scene virtual roaming method based on reflection decomposition
CN106447601B (zh) 一种基于投影-相似变换的无人机遥感影像拼接方法
CN110111250B (zh) 一种鲁棒的自动全景无人机图像拼接方法及装置
CN103985133B (zh) 基于图割能量优化的影像间最优拼接线寻找方法及系统
CN107767339B (zh) 一种双目立体图像拼接方法
CN104463859B (zh) 一种基于跟踪指定点的实时视频拼接方法
JP2017503290A (ja) 無特徴抽出の高密度sfm三次元再構成法
CN110781903B (zh) 基于网格优化和全局相似性约束的无人机图像拼接方法
CN110717936B (zh) 一种基于相机姿态估计的图像拼接方法
CN110246161B (zh) 一种360度全景图像无缝拼接的方法
Li et al. A study on automatic UAV image mosaic method for paroxysmal disaster
CN110415304B (zh) 一种视觉标定方法及系统
CN115456870A (zh) 基于外参估计的多图像拼接方法
Wan et al. Drone image stitching using local mesh-based bundle adjustment and shape-preserving transform
CN110580715A (zh) 一种基于照度约束和格网变形的图像对齐方法
CN107067368B (zh) 基于影像变形的街景影像拼接方法及系统
CN106595602B (zh) 基于同名直线特征的相对定向方法
CN107330856A (zh) 一种基于投影变换和薄板样条的全景成像方法
CN111047513A (zh) 一种用于柱面全景拼接的鲁棒性图像对齐方法及装置
EP2879090B1 (en) Aligning ground based images and aerial imagery

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