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

本发明涉及一种基于弹性配准和网格优化的邻接图像拼接方法，属于图像智能化处理技术领域。

背景技术：

目前，随着无人机技术的发展，无人机在各类场合下的应用越来越广泛。无人机遥感是一种低空遥感，具有图像获取快速、定位准确、操作简单等优点，与航天航空遥感相比，空间分辨率更高，成本更低。由于无人机航高、焦距和视角的限制，单张无人机图像很难反映整个待测区域的情况，为了扩大视场，需要依靠图像拼接技术将多幅航拍图像融合成具有宽视角的、同时保持地面分辨率的全景图像。

图像拼接是将多幅图像拼接成具有更宽视场的全景图的过程。图像拼接包括三个阶段:特征提取与匹配、图像配准和图像合成。传统的图像拼接方法是利用全局单应性对两幅图像进行对齐，如仿射或投影变换，代表方法是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)

其中，ea为对齐项、els为局部相似项、egs为全局相似项、el为直线保持项，λls和λl为对应的权重系数，v为所有图像的网格顶点坐标集，n为所有图像的数量；j为邻接图像之间的邻接关系；m^ij为网格匹配点集合；

对齐项ea(v)中，为网格顶点或网格对应点所在网格的四个顶点的双线性插值；

局部相似项els(v)和全局相似项egs(v)中，ei为ii的所有边的集合，和表示原始图像的某条边及其变形后的边；为边经历的相似变换，和为相似变换中的对应元素，表示为顶点变量的线性组合，si、θi为最佳尺度和最佳旋转角度，为边的权重函数；

直线保持项el(v)中，li表示图像ii中的直线集合，lu、lv、lk分别为直线上的起点、终点和中间采样点，u为直线局部坐标系下的1维坐标，为直线起点、终点和中间采样点所在网格的四个顶点的双线性插值。

作为优选方案，所述si、θi最佳尺度和最佳旋转角度获取的步骤如下：

从多平面单应性矩阵h中估算焦距的初值并分别形成ii、ij的内参矩阵ki、kj，通过以下公式获得ii和ij之间的3d旋转rij的初始估计：

其中，ii、ij为相邻接的两张图像，r代表3d旋转矩阵的参数；

初始化之后，对所有的ki和3d旋转rij初始值执行光束平差法得到得每幅邻接图像ii精细化的焦距fi和3d旋转ri；

每幅邻接图像的最佳尺度计算公式如下：

si＝f1/fi

其中，f1为参考图像的焦距；

使用lsd检测图像的直线，通过弹性配准可得到两幅邻接图像ii和ij之间线的对应关系，每对线段对唯一地确定相对旋转角度，再按照ransac算法投票筛选得到每幅邻接图像的最佳旋转角度θi。

作为优选方案，所述m^ij网格匹配点集合获取步骤如下：

根据特征匹配内点集和多平面单应性变换矩阵h，通过下述两式计算ij上的特征目标点qi在ii上的对应投影点q′i以及投影点q′i和特征源点pi的偏差ri：

其中，和为投影点q′i在x方向和y方向上的分量，和为投影点误差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

其中，jλ是变形能量函数，jd是对齐项，js是平滑项，λ是权重系数；表示投影点q′i在x方向上经历的变形，为投影点误差ri在x方向上的分量。i是第i个特征匹配内点，n代表特征匹配内点的数量；

同理，对于y方向上的变形h(x,y)，能量函数如下：

jλ＝jd+λjs

其中，jλ是变形能量函数，jd是对齐项，js是平滑项，λ是权重系数；表示投影点q′i在y方向上经历的变形，为投影点误差ri在y方向上的分量。

根据薄板样条理论，通过最小化jλ分别得到变形函数g(x，y)和h(x，y)唯一的解析解：

其中，di＝‖x-q′i‖表示任意像素坐标x＝(x，y)^t与第i个投影点q′i的距离，n为投影点的数量；

上述公式具有2(n+3)个待求解系数α＝(α1，α2，α3)^t、β＝(β1，β2,β3)^t，通过求解下述矩阵方程得到：

其中，dij＝‖q′j-q′i‖,c是权重系数，i为单位矩阵；q＝(q′1，…，q′n)^t表示齐次投影点，n为投影点的数量，和为投影点误差在x方向和y方向上的分量；

在求得变形函数之后，对于ii上的网格顶点通过加上形变后再进行多平面单应性变换，得到其在ij上的对应点

如果网格顶点vⁱ及其对应点vⁱ′落在ii和ij的重叠区域中，则收集vⁱ形成网格匹配点集合m^ij，m^ij中元素作为网格顶点，为网格顶点在ij上投影点。

作为优选方案，所述特征匹配内点集获取步骤如下：

1-1)对邻接关系j中一对邻接图像ii和ij，通过sift算法检测每幅邻接图像的特征，通过2nn算法进行匹配获取每幅邻接图像初始匹配点集；

1-2)利用ransac算法从初始匹配点集中估算出一个单应性矩阵对应的内点集；

1-3)从初始匹配点集中提取单应性矩阵对应的内点集，将剩余的匹配点再次进行ransac算法估算出另一个单应性矩阵对应的内点集；

1-4)重复步骤1-3)多次，直至剩余的匹配点数量小于设定阈值；

1-5)将每一步估算提取的单应性矩阵对应的内点集合并成特征匹配内点集

其中，pi为ii上的第i个特征源点，qi为ij上的第i个特征目标点，n为特征匹配内点的数量。

作为优选方案，所述多平面单应性矩阵h获取步骤如下：

通过特征匹配内点集用直接线性变换方法计算出邻接图像之间的多平面单应性变换矩阵h。

作为优选方案，v所有图像的网格顶点坐标集获取步骤如下：

输入n张邻接图像i1，i2，...，in，图像数量n≥2；获取邻接图像之间的邻接关系j；

通过邻接图像和邻接图像之间的邻接关系构建全局匹配图

以work_scale尺度对邻接图像进行降采样，对降采样图像以分割像素大小划分网格，获取每张图片的网格顶点坐标其中vⁱ表示第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像素大小划分网格，获取每张图片的网格顶点坐标其中vⁱ表示第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)将每一步估算提取的单应性矩阵对应的内点集合并成特征匹配内点集

其中，pi为ii上的第i个特征源点，qi为ij上的第i个特征目标点，n为特征匹配内点的数量。

通过最终得到的特征匹配内点集用直接线性变换方法计算出邻接图像之间的多平面单应性变换矩阵h。

如图2所示，步骤(3)基于薄板样条的弹性模型进行邻接图像配准的流程，由于多平面单应性仍然难以将ii和ij的所有像素位置对齐，因此设计该步骤的目的就是利用基于薄板样条的弹性模型将两幅图像更精确地配准，具体步骤如下：

(3-1)根据步骤(2)求得的特征匹配内点集和多平面单应性变换矩阵h，通过下述两式计算ij上的特征目标点qi在ii上的对应投影点q′i以及投影点q′i和特征源点pi的偏差ri：

其中，和为投影点q′i在x方向和y方向上的分量，和为投影点误差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

其中，jλ是变形能量函数，jd是对齐项，js是平滑项，λ是权重系数，用以控制平滑程度。表示投影点q′i在x方向上经历的变形，为投影点误差ri在x方向上的分量。i是第i个特征匹配内点，n代表特征匹配内点的数量。

同理，对于y方向上的变形h(x，y)，能量函数如下：

jλ＝jd+λjs

其中，jλ是变形能量函数，jd是对齐项，js是平滑项，λ是权重系数，用以控制平滑程度。表示投影点q′i在y方向上经历的变形，为投影点误差ri在y方向上的分量。

(3-3)根据薄板样条理论，通过最小化jλ分别得到变形函数g(x，y)和h(x，y)唯一的解析解：

其中，di＝‖x-q′i‖表示任意像素坐标x＝(x，y)^t与第i个投影点q′i的距离，n为投影点的数量。

上述公式具有2(n+3)个待求解系数α＝(α1,α2,α3)^t、β＝(β1,β2,β3)^t，可以通过求解下述矩阵方程得到：

其中，dij＝‖q′j-q′i‖,c是权重系数，i为单位矩阵；q＝(q′1，…，q′n)^t表示齐次投影点，n为投影点的数量，和为投影点误差在x方向和y方向上的分量。

(3-4)在求得变形函数之后，对于ii上的网格顶点通过加上形变后再进行多平面单应性变换，得到其在ij上的对应点

如果网格顶点vⁱ及其对应点vⁱ′落在ii和ij的重叠区域中，则收集vⁱ形成网格顶点作为网格匹配点集合m^ij中元素。网格顶点在ij上对应点为

参考图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的初始估计：

其中，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)

其中，ea为对齐项、els为局部相似项、egs为全局相似项、el为直线保持项，λls和λl为对应的权重系数，v为所有图片的网格顶点坐标集，n为所有图片的数量；

对齐项ea(v)中，为网格顶点或网格顶点对应点或所在网格的四个顶点的双线性插值；

局部相似项els(v)和全局相似项egs(v)中，ei为ii的所有边的集合，和表示原始图片的某条边及其变形后的边。为边经历的相似变换，和为相似变换中的对应元素，可表示为顶点变量的线性组合，si、θi为由步骤s4得到的最佳尺度和旋转，为边的权重函数；

直线保持项el(v)中，li表示图像ii中的直线集合，lu、lv、lk分别为直线上的起点、终点和中间采样点，u为直线局部坐标系下的1维坐标，为直线采样点所在网格的四个顶点的双线性插值。在直线采样过程中，基于网格优化算法，选取直线的长度大于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)为网格优化得到的变形网格顶点。可见有益效果如下：在重叠区域能将图像精确对齐，去除拼接断裂、重影现象；图像在全局范围内保持自然的形状。

本发明实验结果表明，相比于传统方法，能有效消除拼接重影、不对齐现象，具有一定的视差容忍度，并且能够减少多图拼接产生的失真，保持图像形状，具有自然观感。

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