本发明属于图像信息处理领域，涉及图像拼接方法，尤其涉及基于空间三角面片拟合的透视图像拼接方法。
背景技术：
全景图像因其远大于人眼的视场角，在监控、侦查、影视制作以及虚拟现实等领域具有巨大应用价值。拼接式全景成像兼具大视场角和高分辨率的优点，正逐渐成为主流的全景图像采集方式。计算机硬件能力的不断提高和并行计算技术的快速发展也为多路图像甚至视频数据的快速处理提供了有利的现实条件。在计算机视觉和图像处理理论与技术的发展过程中，关于图像拼接技术的创新受到越来越广泛的关注。
图像拼接的目的是以尽可能小的计算代价输出高质量的全景图像。在多幅图像的拍摄过程中，相机的拍摄角度和位置的变化会导致图像之间出现结构偏差，即同一物体在不同图像中呈现出形状或位置上的差异。在场景深度未知的情况下，得到各路图像之间的精确投影关系是困难的。
在过去几十年中，一系列基于不同数学理论的图像拼接方法被提出。其中，传统的全局方法通常使用全局变换将所有输入图像投影到统一的观察平面，目前主流的图像拼接软件通常基于这一方案实现。若相机的光心保持不动或被拍摄场景是近似平面的，全局变换能够得到令人满意的结果；但若图像间存在明显的视差，全局变换的结果就会导致所得全景图像中出现纹理的模糊、重影或断裂等问题。
为消除视差对拼接质量的影响，研究者们提出了一系列局部自适应的图像变换方案，取代单一的全局变换，以达到更高的配准精度。但现有局部变换模型在匹配信息不足或纹理不丰富的区域可能产生不稳定的结果，且需要消耗大量的计算资源，无法满足实际应用对算法效率的要求。
技术实现要素：
本发明提出一种基于空间三角面片拟合的透视图像拼接方法。将被拍摄场景近似为一系列空间面片的组合，并在相机坐标系的单位投影面上逐面片地求解3参数的平面投影变换。与图像坐标系下8参数的单应变换相比，所提方法在应用于图像配准时具有更高的稳定性，因此能够达到更好的拼接效果。与现有的基于匹配数据加权的局部变换模型相比，所提方法具有更低的计算复杂度，因此能够达到更高的计算效率。
具体地，本发明采用的技术方案是：
一种基于空间三角面片拟合的透视图像拼接方法，其特征在于包括以下步骤：
S1.对有重叠区域的两幅输入图像进行特征提取与匹配，得到图像间的匹配特征点；
S2.在S1中得到的匹配特征点基础上，通过最小化关于相机参数的代价函数，得到一组近似的相机参数；
S3.利用S2中得到的相机参数，通过两步RANSAC，剔除S1中得到匹配点中的误匹配点；
S4.在S3中得到的正确匹配特征点基础上，对相机坐标系下的单位投影平面进行三角划分；
S5.在S2中得到的相机参数与S4中得到的三角划分基础上，求解每一个空间三角面片的参数；
S6.在S5中得到的三角面片参数基础上，逐面片地渲染得到精确配准的重投影图像；
S7.在S6中得到的重投影图像基础上，通过图像融合得到最终的全景图像。
优选地，所述步骤S2的具体方法包括：
用于拍摄图像I与I′的相机分别为C与C′；根据特征匹配结果进行相机参数估计的核心理论依据是极线约束：
x′TFx＝0
其中，x与x′为图像I与I′中的一对匹配点，F为基本矩阵。除此之外，增加关于全局投影的二级约束：
x′～Hx
其中，～表示正比关系，H为单应矩阵。则对于一对匹配点{x,x′}，需最小化的代价函数为：
E(x,x′)＝h(re,σe)+λh(rp,σp)
其中，系数λ用于控制rp的影响，取为λ＝0.01，函数h(·)为具有误匹配抑制能力的Huber误差函数：
参数σe与σp分别表示极线约束与全局投影约束下正常误差的上限，残差二者的度量单位均为像素，y＝K-1x，y′＝K′-1x′，其中K与K′分别为两相机C与C′的内参数矩阵，Q为本征矩阵，且Q＝[t]×R，其中R与t分别为C与C′的相机坐标系间的旋转矩阵与平移向量；在re的表达式中，记(a,b,c)T＝Qy，则标准化参数其中f′x与f′y分别为相机C′在x，y方向上的等效焦距；在rp的表达式中，系数w＝r3Ty，其中r3T为R的第三行；
对于所有匹配点{xi,x′i}，以及它们对应的残差re,i与rp,i，i＝1,…,n，总的代价函数为：
Eall是关于残差re与rp的二次函数，因此Eall的最小化通过LM迭代实现，所得结果是一组真实相机参数的近似值；其中K，K′为内参数，R，t为外参数。
优选地，所述步骤S3的具体方法包括：
两步RANSAC，包括
第一步以图像间的极线约束作为RANSAC的目标函数，即最小化
其中，yi＝K-1xi，y′i＝K′-1x′i；
第二步以全局平面变换模型作为目标函数，即最小化
取两次筛选后所得匹配数据的交集用于建立图像间的局部配准模型。
优选地，所述步骤S4的具体方法包括：
给定剔除误匹配之后的匹配点{xi,x′i|i＝1,…,n}，其中n表示正确匹配点的数量，以yi＝K-1xi，i＝1,…,n为顶点，对单位投影平面进行delaunay三角划分，所得平面三角形代表空间三角面片在C的相机坐标系内单位投影平面上的投影；为保持空间三角面片模型的一致性，对于C′的相机坐标系下的顶点y′i＝K′-1x′i，i＝1,…,n，适用与yi相一致的三角划分。
优选地，所述步骤S5的具体方法包括：
图像I与I′内一对匹配点x与x′是同一空间点P在不同相机成像平面上的投影。x与x′之间的变换关系可表示为
x′～Hx
其中
H＝K′(R+tmT)K-1
在已估计得到相机参数的情况下，描述点P所在三角面片的参数m是H各组成部分中唯一的未知量，其自由度为3；记y＝K-1x，y′＝K′-1x′，则投影关系转化为
y′～(R+tmT)y
根据上式，向量y′与(R+tmT)y平行，则
y′×(R+tmT)y＝03×1
根据点积和向量积对加法的分配率，有
-y′×(Ry)＝(y′×t)(mTy)
mTy是一个标量，上式为一个关于mTy的过约束线性方程组，最小二乘意义上求解可得
上式确定了关于m的一个线性约束，该约束可表示为
两相机单位投影平面上的一对三角形ΔABC与ΔA′B′C′是同一个空间三角面片的投影；根据三对匹配的顶点{yA,y′A}，{yB,y′B}，{yC,y′C}确立三个关于m的线性约束；记
其中，i∈{A,B,C}。则三对匹配点确定的线性约束可写为
Am＝13×1
其中，A＝(aAT,aBT,aCT)T。上式的约束个数与未知数个数均为3，可求得该空间三角面片的参数mΔ。
优选地，所述步骤S6的具体方法包括：
判断单位投影平面上一点P是否在ΔABC内部，包括
点P的空间位置可表示为对于下式
u(yB-yA)+v(yC-yA)＝yP-yA
其中，yP，yA，yB，yC分别为点P，A，B，C在单位投影面上的坐标；求解上式得到位置参数u与v，则点P是ΔABC内点的充要条件是0≤u，0≤v且u+v≤1；对于每个空间三角形，首先计算该三角形三个顶点在图像平面的投影点，并据此构建图像平面上当前空间三角形投影的矩形包围盒；然后标记包围盒内位于空间三角形投影内部的像素点；最后计算这些像素点上的映射并将其添加到用于图像重投影的映射表中。
根据所得映射表对输入图像进行重投影，得到精确配准后的重投影图像。
优选地，所述步骤S7的具体方法包括：
在观察平面上对重投影图像进行基于图像金字塔的多尺度融合，得到最终的全景图像。
本发明所述方法提出一种基于空间三角面片拟合的透视图像拼接方法。该方法通过对被拍摄场景结构的面片化近似处理，将被拍摄场景近似为一系列空间面片的组合，并在相机坐标系的单位投影面上逐面片地求解3参数的平面投影变换。在配准精度方面，达到了比图像坐标系下8参数单应变换更高的稳定性，因此能够达到更好的拼接效果。在计算效率方面，达到了比现有的基于匹配数据加权的局部变换模型更低的计算复杂度，能够满足更多实际应用的需求。
附图说明：
图1为根据本发明实施方式的总体流程图。
图2为根据本发明实施方式的输入图像。
图3为根据本发明实施方式的根据极线约束剔除误匹配点的结果图像。
图4为根据本发明实施方式的两步RANSAC后的误匹配点剔除的结果图像。
图5为根据本发明实施方式的平面三角划分的结果图像。
图6为根据本发明实施方式的平面三角形内点判断示意图。
具体实施方式
下面将结合附图和实例对本发明的具体实施方式作进一步的详细说明。
本发明所述方法通过对被拍摄场景结构的面片化近似处理，提出一种基于空间三角面片拟合的透视图像拼接方法。首先，利用图像间的点匹配结果，通过最小化关于相机参数的代价函数，得到一组近似的相机参数；然后，通过两步RANSAC，剔除误匹配点；然后，对相机坐标系下的单位投影平面进行三角划分，并求解每一个空间面片的参数；最后，逐面片地渲染得到精确配准的重投影图像，并在观察平面上融合得到最终的全景图像。总体流程图如图1所示。
本发明提供的基于空间三角面片拟合的透视图像拼接方法通过以下步骤实现：
S1.对有重叠区域的两幅输入图像进行特征提取与匹配，得到图像间的匹配特征点；
对于如图2所示输入图像I与I′，使用SIFT特征算法进行特征点的检测与匹配，得到图像间的一系列匹配点{xi,x′i|i＝1,…,n}，其中n为所得匹配点对的数量。
S2.在S1中得到的匹配特征点基础上，通过最小化关于相机参数的代价函数，得到一组近似的相机参数；
记用于拍摄图像I与I′的相机分别为C与C′。根据特征匹配结果进行相机参数估计的核心理论依据是极线约束：
x′TFx＝0
其中，x与x′为图像I与I′中的一对匹配点，F为基本矩阵。除此之外，增加关于全局投影的二级约束：
x′～Hx
其中，～表示正比关系，H为单应矩阵。则对于一对匹配点{x,x′}，需最小化的代价函数为：
E(x,x′)＝h(re,σe)+λh(rp,σp)
其中，系数λ用于控制rp的影响，取为λ＝0.01，函数h(·)为具有误匹配抑制能力的Huber误差函数：
参数σe与σp分别表示极线约束与全局投影约束下正常误差的上限，残差二者的度量单位均为像素，y＝K-1x，y′＝K′-1x′，其中K与K′分别为两相机C与C′的内参数矩阵，Q为本征矩阵，且Q＝[t]×R，其中R与t分别为C与C′的相机坐标系间的旋转矩阵与平移向量。在re的表达式中，记(a,b,c)T＝Qy，则标准化参数其中f′x与f′y分别为相机C′在x，y方向上的等效焦距。在rp的表达式中，系数w＝r3Ty，其中r3T为R的第三行。
对于所有匹配点{xi,x′i}，以及它们对应的残差re,i与rp,i，i＝1,…,n，总的代价函数为：
注意到Eall是关于残差re与rp的二次函数，因此Eall的最小化可通过LM迭代实现，所得结果是一组真实相机参数(内参数K，K′以及外参数R，t)的近似值。
S3.利用S2中得到的相机参数，通过两步RANSAC，剔除S1中得到匹配点中的误匹配点；
空间三角面片拟合的基础是图像间正确的点匹配关系，而误匹配的出现是不可避免的。因此，在图像配准之前需要首先进行误匹配的剔除。这一过程通过两步RANSAC实现：
第一步以图像间的极线约束作为RANSAC的目标函数，即最小化
其中，yi＝K-1xi，y′i＝K′-1x′i。如图2所示，远离其对应极线的匹配点被认为是错误的。
第二步以全局平面变换模型作为目标函数，即最小化
取两次筛选后所得匹配数据的交集用于建立图像间的局部配准模型。对于图2中的输入图像，使用上述两步RANSAC方法进行误匹配剔除的结果如图4所示。
S4.在S3中得到的正确匹配特征点基础上，对相机坐标系下的单位投影平面进行三角划分；
现实场景可近似地看作由大量不同方向和位置的空间面片拼合而成，本发明采用三角面片对场景进行拟合。
给定剔除误匹配之后的匹配点{xi,x′i|i＝1,…,n}(此时n表示正确匹配点的数量)，以yi＝K-1xi，i＝1,…,n为顶点，对单位投影平面进行delaunay三角划分，所得平面三角形代表空间三角面片在C的相机坐标系内单位投影平面上的投影。为保持空间三角面片模型的一致性，对于C′的相机坐标系下的顶点y′i＝K′-1x′i，i＝1,…,n，适用与yi相一致的三角划分。将单位投影面上的三角划分投影到全景图像平面的效果如图5所示。
S5.在S2中得到的相机参数与S4中得到的三角划分基础上，求解每一个空间三角面片的参数；
图像I与I′内一对匹配点x与x′是同一空间点P在不同相机成像平面上的投影。x与x′之间的变换关系可表示为
x′～Hx
其中
H＝K′(R+tmT)K-1
在已估计得到相机参数的情况下，描述点P所在三角面片的参数m是H各组成部分中唯一的未知量，其自由度为3。进一步地，记y＝K-1x，y′＝K′-1x′，则投影关系转化为
y′～(R+tmT)y
根据上式，向量y′与(R+tmT)y平行，则
y′×(R+tmT)y＝03×1
根据点积和向量积对加法的分配率，有
-y′×(Ry)＝(y′×t)(mTy)
注意到mTy是一个标量，上式为一个关于mTy的过约束线性方程组，最小二乘意义上求解可得
上式确定了关于m的一个线性约束，该约束可表示为
两相机单位投影平面上的一对三角形ΔABC与ΔA′B′C′是同一个空间三角面片的投影。则根据三对匹配的顶点{yA,y′A}，{yB,y′B}，{yC,y′C}能够确立三个关于m的线性约束。具体地，记
其中，i∈{A,B,C}。则三对匹配点确定的线性约束可写为
Am＝13×1
其中，A＝(aAT,aBT,aCT)T。上式的约束个数与未知数个数均为3，可求得唯一解，即该空间三角面片的参数mΔ。
S6.在S5中得到的三角面片参数基础上，逐面片地渲染得到精确配准的重投影图像；
在已知面片参数的基础上，为由输入图像渲染得到精确配准的重投影图像，还需建立图像点与各空间面片间的归属关系，即判断单位投影平面上一点P是否在ΔABC内部。如图6所示，点P的空间位置可表示为即
u(yB-yA)+v(yC-yA)＝yP-yA
其中，yP，yA，yB，yC分别为点P，A，B，C在单位投影面上的坐标。求解上式得到位置参数u与v，则点P是ΔABC内点的充要条件是0≤u，0≤v且u+v≤1。具体实现时，对于每个空间三角形，首先计算该三角形三个顶点在图像平面的投影点，并据此构建图像平面上当前空间三角形投影的矩形包围盒；然后标记包围盒内位于空间三角形投影内部的像素点；最后计算这些像素点上的映射并将其添加到用于图像重投影的映射表中。
根据所得映射表对输入图像进行重投影，得到精确配准后的重投影图像。
S7.在S6中得到的重投影图像基础上，通过图像融合得到最终的全景图像。
在观察平面上对重投影图像进行基于图像金字塔的多尺度融合，得到最终的全景图像。
以上所述，仅为本发明示例性的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。