CN110956661A - 基于双向单应矩阵的可见光与红外相机动态位姿计算方法 - Google Patents

Abstract

本发明属于图像处理和计算机视觉领域，公开了基于双向单应矩阵的可见光与红外相机动态位姿计算方法。通过提取并匹配红外图像和可见光图像上的直线特征来计算红外图像和可见光图像之间的变换矩阵。在红外图像上检测Harris角点，然后根据变换矩阵计算检测出的Harris角点在可见光图像上的初始对应，并在初始对应周围检测Harris角点作为最终对应。之后在可见光图像上做同样的操作。提高了红外特征点与可见光特征点的匹配数量和匹配效率，从而更有效地对红外相机和可见光相机进行联合自标定。

Description

基于双向单应矩阵的可见光与红外相机动态位姿计算方法

技术领域

本发明属于图像处理和计算机视觉领域，涉及基于双向单应矩阵的可见光与红外相机动态位姿计算方法。

背景技术

红外线(Infrared)是波长介于微波与可见光之间的电磁波，波长比红光要长。高于绝对零度(-273.15℃)的物质都可以产生红外线。红外图像由于其具有透过雾、雨等进行观察的能力而被广泛用于军事国防、资源勘探、气象预报、环境监测、医学诊治、海洋研究等不同领域。利用红外线可以隔着薄雾和烟雾拍摄景物，而且在夜间也可以进行红外摄影。红外相机成像的优点是在极端场景(低光、雨雪、浓雾等)也可以成像，缺点是分辨率低、图像细节较模糊。相比之下，可见光相机的优点是分辨率高、图像细节清晰，但是在极端场景下不能成像。因此，将红外相机和可见光相机结合起来具有重大的现实意义。

立体视觉是计算机视觉领域的重要主题。其目的是重建场景的3D几何信息。双目立体视觉是立体视觉的重要领域。在双目立体视觉中，左右摄像头用于模拟两只眼睛。通过计算双目图像之间的差异来计算深度图像。双目立体视觉具有效率高，准确度高，系统结构简单，成本低的优点。由于双目立体视觉需要匹配左右图像捕获点上的相同点，因此相机两个镜头的焦距和图像捕获中心，以及左右两个镜头之间的位置关系。为了得到以上数据，需要对相机进行标定。获取可见光相机和红外相机之间的位置关系称为联合标定。

在标定过程中获得了相机的两个镜头参数和相对位置参数，但这些参数不稳定。当温度、湿度等发生变化时，相机镜头的内部参数也会发生变化。另外，由于意外的相机碰撞，两个镜头之间的位置关系可能会改变。因此，每次使用摄像机时，都必须修改内部和外部参数，这就是自标定。在已知相机内部参数的情况下，通过分别提取红外图像特征和可见光图像特征来对红外镜头和可见光镜头的位置关系进行修正，即红外相机与可见光相机的联合自标定。

为了缩小特征点的匹配范围，在特征点检测之前，对红外图像和可见光图像进行配准。由于直线特征在不同模态下具有鲁棒性，通过直线特征的提取与匹配计算红外图像和可见光图像的配准关系。然后分别提取红外图像和可见光图像上的Harris角点，然后根据配准关系确定角点在另一台相机下的的对应点。

发明内容

本发明旨在解决由于温湿度、震动等因素造成红外相机和可见光相机位置关系的改变。首先利用直线的匹配关系对可见光图像和红外图像进行配准。然后根据配准关系确定特征点的对应点，从而对原有的标定结果进行修正。

包括下列步骤：

1)原图校正：使用红外相机和可见光相机同时拍摄一组有运动物体的场景连续帧；将原图根据红外相机和可见光相机各自内参和原来的外参进行去畸变和双目校正；流程如图2所示。

2)提取匹配直线，并根据配对的直线计算单应矩阵H^*

3)分别在红外图像和可见光图像检测Harris角点

4)根据单应矩阵寻找上一步检测出的Harris角点的匹配点。

5)判断特征点覆盖区域：将图像分成m*n个格子，如果特征点覆盖到所有格子，则进行下一步，否则继续拍摄图像，重新进行步骤1)～步骤4)。

6)修正标定结果：使用所有特征点的图像坐标来计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加。

步骤1)中原图校正，包括如下步骤：

1-1)计算图像的像素点对应的正规坐标系下的坐标；其像素坐标系以图片的左上角为原点，其x轴和y轴分别与图像坐标系的x轴和y轴平行；像素坐标系的单位是像素；以相机的光心作为图像坐标系的原点，且将光心到图像平面的距离缩放到1；像素坐标与正规坐标的关系如下：

u＝KX

其中，

表示图像的像素坐标；

表示相机的内参矩阵，f_x和f_y分别表示图像x方向和y方向的焦距，单位是像素，(c_x,c_y)表示相机的主点位置，即相机中心在图像上的对应位置；

是正规坐标系下的坐标；已知图像的像素坐标系以及相机的内参计算出像素点对应的正规坐标系，即X＝K^-1u；

1-2)去除图像畸变：由于镜头生产工艺的限制，实际情况下的镜头会存在一些失真现象导致非线性的畸变；因此纯线性模型不能完全准确地描述成像几何关系；非线性畸变可大致分为径向畸变和切向畸变；

图像径向畸变是图像像素点以畸变中心为中心点，沿着径向产生的位置偏差，从而导致图像中所成的像发生形变；径向畸变的表述如下：

x_d＝x(1+k₁r²+k₂r⁴+k₃r⁶)

y_d＝y(1+k₁r²+k₂r⁴+k₃r⁶)

其中，r²＝x²+y²，k₁、k₂、k₃为径向畸变参数；

图像切向畸变是由于摄像机制造上的缺陷使得透镜本身与图像平面不平行而产生的，定量描述为：

x_d＝x+(2p₁xy+p₂(r²+2x²))

y_d＝y+(p₁(r²+2y²)+2p₂xy)

其中，p₁、p₂为切向畸变系数；

畸变前后的坐标关系如下：

x_d＝x(1+k₁r²+k₂r⁴+k₃r⁶)+(2p₁xy+p₂(r²+2x²))

y_d＝y(1+k₁r²+k₂r⁴+k₃r⁶)+(p₁(r²+2y²)+2p₂xy)

其中，(x,y)是理想状态下的正规坐标，(x_d,y_d)是实际带有畸变的正规坐标；

1-3)根据原来两相机的旋转关系将两图转回来：已知原来两个相机之间的旋转矩阵R和平移向量t，使得：

X_r＝RX_l+t

其中，X_l表示红外相机的正规坐标，X_r表示可见光相机的正规坐标；将可见光图像旋转R反方向一半的角度，将红外图像旋转R正方向一半的角度；

1-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系。

所述步骤2)，具体包括以下步骤：

2-1)将红外图像和可见光图像分别做直方图均衡化并提取Canny边缘点；

2-2)使用类RANSAC的方法提取出直线：首先随机选出两个像素点并连线，然后计算该连线的得分：

其中，Ω是边缘点的集合d(g,x,y)表示点(x,y)到直线g的距离，τ是阈值；每次选取若干点对连线，取得分最大的作为检测出来的直线；

为每条直线分配四个属性：

第一个，L_ls＝|x_end-x_start|/D_start→end：其中x_start和x_end分别表示直线的起点和终点，D_start→end表示起点到终点的距离；

第二个，sp：如果L_ls＞0.7，则sp表示直线与x轴夹角，否则sp表示直线与y轴夹角；

第三个，os：如果L_ls＞0.7，则os表示直线相对于y轴的偏移，否则os表示直线相对于x轴的偏移；

第四个，pec_edge＝N_edge/N_total：表示直线附近的边缘点占比，其中N_edge表示直线附近边缘点个数，N_total表示直线附近像素点的总个数；

2-3)直线配对：直线(L_ls,sp,pec_edge)和

的配对得分如下所示：

其中K(·)表示Epanechnikov核函数，σ表示核的大小；

2-4)计算H^*＝min_H(error(l′,Hl))，其中l表示红外直线，l′表示可见光直线。

所述步骤3)，具体包括以下步骤：

3-1)构建梯度矩阵M；在图像上，角点表现为不同边缘之间的交点；而且不论从什么角度来观测，它都是不同边缘之间的交点，不会因为视角的变化而变化；此外，角点邻域内的点的梯度会有大幅变化；特征点应满足：当移动窗口时，特征点所在窗口与其周围各个方向窗口的亮度分布差别很大；将窗口移动[u,v]时，灰度变化如下所示：

将上式泰勒展开，得：

其中，ω(x,y)表示(x,y)点对应的权值，权值是常数，或者是高斯核的对应系数；I_x和I_y分别表示图像在x方向和y方向上的梯度，使用高斯函数计算梯度，也可以两个3*3得卷积核(比如Prewitt算子或者Sobel算子)求得，矩阵M表示为：

3-2)计算矩阵M的两个特征值λ₁和λ₂，λ₁和λ₂所对应的特征向量分别代表着灰度变化最快和最慢的两个方向；λ₁和λ₂的大小关系和对应点的属性存在一下的对应关系：

当λ₁和λ₂的值都很小时，该点落在平滑区域内；

当λ₁>>λ₂或者λ₂>>λ₁时，该点落在图像的边缘上；

当λ₁和λ₂的值都很大，且处于同一大小水平时，可以认为该点属于角点；

3-3)使用一个值R来描述该点的角点响应，然后通过阈值σ1和σ2来判断该点是不是一个角点；角点响应值R＝det(M)-k*trace(M)²，其中det(M)表示矩阵M对应的行列式的值，trace(M)表示矩阵M的迹，即：

det(M)＝λ₁*λ₂

trace(M)＝λ₁+λ₂

其中k是常数，一般取0.04～0.06

当|R|<σ1时，该区域是平面；

当R<0时，该区域是直线；

当R＞σ2时，该区域是角点。

步骤4)中寻找上一步检测出的Harris角点的匹配点，包括以下步骤：

4-1)记步骤3)在红外图像检测出来的角点集合为

对集合

里的每个角点

根据步骤2)计算出的单应矩阵计算

在可见光图像上的初始对应点，在初始对应点周围w*w窗口内检测Harris角点，记角点集为

如果

为空，则取下一个角点，重复本步骤，否则计算

中各个点和

的相似度，如果最大的相似度大于t1，则将其对应点

作为

在可见光图像上的初始对应点并加入到匹配点集中，选取下一点重复该步骤直到遍历所有红外角点；

4-2)记步骤3)在可见光图像检测出来的角点集合为

对集合

里的每个角点

根据步骤2)计算出的单应矩阵计算

在红外图像上的初始对应点，在初始对应点周围w*w窗口内检测Harris角点，记角点集为

如果

为空，则取下一个角点，重复本步骤，否则计算

中各个点和

的相似度，如果最大的相似度大于t1，则将其对应点

作为

在可见光图像上的初始对应点并加入到匹配点集中，选取下一点重复该步骤直到遍历所有可见光角点；

4-3)将匹配点集去重；

4-4)以红外图特征点

为基准，抛物线拟合优化对应可见光图的整数像素特征点

得到的对应可见光图的亚像素特征点

其中

为x方向上的亚像素偏移量，

为y方向上的亚像素偏移量；

4-5)以对应可见光图整数像素特征点

为基准，根据4-4)的方法计算出对应红外图的亚像素特征点

其中

为x方向上的亚像素偏移量，

为y方向上的亚像素偏移量；

4-6)得到最终的匹配点对为

步骤6)中修正标定结果，包括以下步骤：

6-1)使用随机抽样一致性对点对做进一步筛选；

6-2)求解基础矩阵F和本质矩阵E：红外和可见光对应像素点对u_l、u_r和基础矩阵F的关系是：

将对应点坐标代入上式，构建齐次线性方程组求解F；

基础矩阵和本质矩阵的关系是：

其中，K_l、K_r分别是红外和可见光相机的内参矩阵。

6-3)从本质矩阵分解出旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：

E＝[t]_×R

其中[t]_×表示t的叉乘矩阵。

将E做奇异值分解，得

定义两个矩阵

和

ZW＝Σ

所以E可以写成以下两种形式

(1)E＝UZU^TUWV^T

令[t]_×＝UZU^T，R＝UWV^T

(2)E＝-UZU^TUW^TV^T

令[t]_×＝-UZU^T，R＝UW^TV^T

6-4)将分解出的旋转和平移关系叠加到原来的外参里面；

记去畸变前的旋转矩阵为R₀，平移向量为t₀＝(t_x,t_y,t_z)^T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′_x,t′_y,t′_z)^T；则新的R_new和t_new如下：

此外，还需要将t_new乘一个系数，使得t_new在x方向上的分量

本发明的有益效果是：

本发明解决了由于温湿度、震动等因素造成红外双目相机位置关系的改变。具有速度快、结果精确、操作简单等优点。使用直线这一跨模态鲁棒的特征进行图像配准，减小了匹配范围。根据配准的单应矩阵确定Harris角点对应点的初始值，并根据该初始值周围窗口进行Harris角点检测，这样既提高了特征点匹配效率，又提高了特征点对的个数。

附图说明

图1为整体流程图。

图2为校正流程示意图。

图3为Harris角点条件判断示意图。

具体实施方式

1)原图校正：将原图根据红外相机和可见光相机各自内参和原来的外参进行去畸变和双目校正。流程如图2所示。

1-1)计算图像的像素点对应的正规坐标系下的坐标。其中，正规坐标系是相机坐标系在平面Z＝1的投影；而相机坐标系是以相机的中心作为图像坐标系的原点，以图片方向为XY轴方向，以垂直于图像为Z轴方向的坐标系。像素坐标系以图片的左上角为原点，其x轴和y轴分别与图像坐标系的x轴和y轴平行。像素坐标系的单位是像素。像素坐标与正规坐标的关系如下：

u＝KX

其中，

表示图像的像素坐标；

表示相机的内参矩阵，f_x和f_y分别表示图像x方向和y方向的焦距，单位是像素，(c_x,c_y)表示相机的主点位置，即相机中心在图像上的对应位置；

是正规坐标系下的坐标。已知图像的像素坐标系以及相机的内参计算出像素点对应的正规坐标系，即X＝K^-1u；

1-2)去除图像畸变：由于镜头生产工艺的限制，实际情况下的镜头会存在一些失真现象导致非线性的畸变。因此纯线性模型不能完全准确地描述成像几何关系。非线性畸变可大致分为径向畸变和切向畸变。

图像径向畸变是图像像素点以畸变中心为中心点，沿着径向产生的位置偏差，从而导致图像中所成的像发生形变。径向畸变的大致表述如下：

x_d＝x(1+k₁r²+k₂r⁴+k₃r⁶)

y_d＝y(1+k₁r²+k₂r⁴+k₃r⁶)

其中，r²＝x²+y²，k₁、k₂、k₃为径向畸变参数。

图像切向畸变是由于摄像机制造上的缺陷使得透镜本身与图像平面不平行而产生的，可定量描述为：

x_d＝x+(2p₁xy+p₂(r²+2x²))

y_d＝y+(p₁(r²+2y²)+2p₂xy)

其中，p₁、p₂为切向畸变系数。

综上，畸变前后的坐标关系如下：

x_d＝x(1+k₁r²+k₂r⁴+k₃r⁶)+(2p₁xy+p₂(r²+2x²))

y_d＝y(1+k₁r²+k₂r⁴+k₃r⁶)+(p₁(r²+2y²)+2p₂xy)

其中，(x,y)是理想状态下的正规坐标，(x_d,y_d)是实际带有畸变的正规坐标。

1-3)根据原来两相机的旋转关系将两图转回来：已知原来两个相机之间的旋转矩阵R和平移向量t，使得：

X_r＝RX_l+t

其中，X_l表示红外相机的正规坐标，X_r表示可见光相机的正规坐标。将红外图像向R正方向旋转一半的角度，将可见光图像向R反方向旋转一半的角度；

1-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系。

2)提取匹配直线，并根据配对的直线计算单应矩阵H^*

2-1)将红外图像和可见光图像分别做直方图均衡化并提取Canny边缘点。

2-2)使用类RANSAC的方法提取出直线：首先随机选出两个像素点并连线，然后计算该连线的得分：

其中Ω是边缘点的集合d(g,x,y)表示点(x,y)到直线g的距离，τ是阈值。每次选取若干点对连线，取得分最大的作为检测出来的直线。

为每条直线分配四个属性：

第一个，L_ls＝|x_end-x_start|/D_start→end：其中x_start和x_end分别表示直线的起点和终点，D_start→end表示起点到终点的距离；

第二个，sp：如果L_ls＞0.7，则sp表示直线与x轴夹角，否则sp表示直线与y轴夹角；

第三个，os：如果L_ls＞0.7，则os表示直线相对于y轴的偏移，否则os表示直线相对于x轴的偏移；

第四个，pec_edge＝N_edge/N_total：表示直线附近的边缘点占比，其中N_edge表示直线附近边缘点个数，N_total表示直线附近像素点的总个数；

2-3)直线配对：直线(L_ls,sp,pec_edge)和

的配对得分如下所示：

其中K(·)表示Epanechnikov核函数，σ表示核的大小。

2-4)计算H^*＝min_H(error(l′,Hl))，其中l表示红外直线，l’表示可见光直线。

3)分别在红外图像和可见光图像检测Harris角点

3-1)构建梯度矩阵M。在图像上，角点表现为不同边缘之间的交点。而且不论从什么角度来观测，它都是不同边缘之间的交点，不会因为视角的变化而变化。此外，角点邻域内的点的梯度会有大幅变化。特征点应满足：当移动窗口时，特征点所在窗口与其周围各个方向窗口的亮度分布差别很大。将窗口移动[u,v]时，灰度变化如下所示：

将上式泰勒展开，得：

其中，ω(x,y)表示(x,y)点对应的权值，权值可以是常数，也可以是高斯核的对应系数。I_x和I_y分别表示图像在x方向和y方向上的梯度，这里可以使用高斯函数来计算梯度，也可以两个3*3得卷积核(比如Prewitt算子或者Sobel算子)求得，矩阵M可表示为：

3-2)计算矩阵M的两个特征值λ₁和λ₂，λ₁和λ₂所对应的特征向量分别代表着灰度变化最快和最慢的两个方向。λ₁和λ₂的大小关系和对应点的属性存在一下的对应关系，如图3所示：

(1)当λ₁和λ₂的值都很小时，该点落在平滑区域内。

(2)当λ₁>>λ₂或者λ₂>>λ₁时，该点落在图像的边缘上。

(3)当λ₁和λ₂的值都很大，且处于同一大小水平时，可以认为该点属于角点。

3-3)使用一个值R来描述该点的角点响应，然后通过阈值σ1和σ2来判断该点是不是一个角点；角点响应值R＝det(M)-k*trace(M)²，其中det(M)表示矩阵M对应的行列式的值，trace(M)表示矩阵M的迹，即：

det(M)＝λ₁*λ₂

trace(M)＝λ₁+λ₂

其中k是常数，一般取0.04～0.06

当|R|<σ1时，该区域是平面；

当R<0时，该区域是直线；

当R＞σ2时，该区域是角点。

4)根据单应矩阵寻找上一步检测出的Harris角点的匹配点。

4-1)记步骤3)在红外图像检测出来的角点集合为

对集合

里的每个角点

根据步骤2)计算出的单应矩阵计算

在可见光图像上的初始对应点，在初始对应点周围w*w窗口内检测Harris角点，记角点集为

如果

为空，则取下一个角点，重复本步骤，否则计算

中各个点和

的相似度，如果最大的相似度大于t1，则将其对应点

作为

在可见光图像上的初始对应点并加入到匹配点集中，选取下一点重复该步骤直到遍历所有红外角点；

4-2)记步骤3)在可见光图像检测出来的角点集合为

对集合

里的每个角点

根据步骤2)计算出的单应矩阵计算

在红外图像上的初始对应点，在初始对应点周围w*w窗口内检测Harris角点，记角点集为

如果

为空，则取下一个角点，重复本步骤，否则计算

中各个点和

的相似度，如果最大的相似度大于t1，则将其对应点

作为

在可见光图像上的初始对应点并加入到匹配点集中，选取下一点重复该步骤直到遍历所有可见光角点；

4-3)将匹配点集去重；

4-4)以红外图特征点

为基准，抛物线拟合优化对应可见光图的整数像素特征点

得到的对应可见光图的亚像素特征点

其中

为x方向上的亚像素偏移量，

为y方向上的亚像素偏移量；

4-5)以对应可见光图整数像素特征点

为基准，根据4-4)的方法计算出对应红外图的亚像素特征点

其中

为x方向上的亚像素偏移量，

为y方向上的亚像素偏移量；

4-6)得到最终的匹配点对为

5)判断特征点覆盖区域：将图像分成m*n个格子，如果特征点覆盖到所有格子，则进行下一步，否则继续拍摄图像，重新进行步骤1)～步骤4)。

6)修正标定结果：使用所有特征点的图像坐标来计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加。

6-1)使用随机抽样一致性(RANSAC)对点对做进一步筛选。

6-2)求解基础矩阵F和本质矩阵E：红外和可见光对应像素点对u_l、u_r和基础矩阵F的关系是：

可以将对应点坐标代入上式，构建齐次线性方程组求解F。

基础矩阵和本质矩阵的关系是：

其中，K_l、K_r分别是红外和可见光相机的内参矩阵。

6-3)从本质矩阵分解出旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：

E＝[t]_×R

其中[t]_×表示t的叉乘矩阵。

将E做奇异值分解，得

定义两个矩阵

和

ZW＝Σ

所以E可以写成以下两种形式

(1)E＝UZU^TUWV^T

令[t]_×＝UZU^T，R＝UWV^T

(2)E＝-UZU^TUW^TV^T

令[t]_×＝-UZU^T，R＝UW^TV^T

6-4)将分解出的旋转和平移关系叠加到原来的外参里面；

记去畸变前的旋转矩阵为R₀，平移向量为t₀＝(t_x,t_y,t_z)^T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′_x,t′_y,t′_z)^T；则新的R_new和t_new如下：

此外，还需要将t_new乘一个系数，使得t_new在x方向上的分量

Claims (5)

1.基于双向单应矩阵的可见光与红外相机动态位姿计算方法，其特征在于，包括如下步骤：

1)原图校正：使用红外相机和可见光相机同时拍摄一组有运动物体的场景连续帧；将原图根据红外相机和可见光相机各自内参和原来的外参进行去畸变和双目校正；

2)提取匹配直线，并根据配对的直线计算单应矩阵H^*；

2-1)将红外图像和可见光图像分别做直方图均衡化并提取Canny边缘点；

2-2)使用类RANSAC的方法提取出直线：首先随机选出两个像素点并连线，然后计算该连线的得分：

其中，Ω是边缘点的集合d(g，x，y)表示点(x，y)到直线g的距离，τ是阈值；每次选取若干点对连线，取得分最大的作为检测出来的直线；

为每条直线分配四个属性：

第一个，L_ls＝|x_end-x_start|/D_start→end：其中x_start和x_end分别表示直线的起点和终点，D_start→end表示起点到终点的距离；

第二个，sp：如果L_ls＞0.7，则sp表示直线与x轴夹角，否则sp表示直线与y轴夹角；

第三个，os：如果L_ls＞0.7，则os表示直线相对于y轴的偏移，否则os表示直线相对于x轴的偏移；

第四个，pec_edge＝N_edge/N_total：表示直线附近的边缘点占比，其中N_edge表示直线附近边缘点个数，N_total表示直线附近像素点的总个数；

2-3)直线配对：直线(L_ls，sp，pec_edge)和

的配对得分如下所示：

其中K(·)表示Epanechnikov核函数，σ表示核的大小；

2-4)计算H^*＝min_H(error(l′，Hl))，其中l表示红外直线，l′表示可见光直线；

3)分别在红外图像和可见光图像检测Harris角点；

4)根据步骤2)得到的单应矩阵H^*，寻找步骤3)检测出的Harris角点的匹配点，即特征点；

5)判断特征点覆盖区域：将图像分成m*n个格子，如果特征点覆盖到所有格子，则进行下一步，否则继续拍摄图像，重新进行步骤1)～步骤4)；

6)修正标定结果：使用所有特征点的图像坐标来计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加。

2.根据权利要求1所述的基于双向单应矩阵的可见光与红外相机动态位姿计算方法，其特征在于，步骤1)中原图校正，包括如下步骤：

1-1)计算图像的像素点对应的正规坐标系下的坐标；其像素坐标系以图片的左上角为原点，其x轴和y轴分别与图像坐标系的x轴和y轴平行；像素坐标系的单位是像素；以相机的光心作为图像坐标系的原点，且将光心到图像平面的距离缩放到1；像素坐标与正规坐标的关系如下：

u＝KX

其中，

表示图像的像素坐标；

表示相机的内参矩阵，f_x和f_y分别表示图像x方向和y方向的焦距，单位是像素，(c_x，c_y)表示相机的主点位置，即相机中心在图像上的对应位置；

是正规坐标系下的坐标；已知图像的像素坐标系以及相机的内参计算出像素点对应的正规坐标系，即X＝K^-1u；

1-2)去除图像畸变：图像径向畸变是图像像素点以畸变中心为中心点，沿着径向产生的位置偏差，从而导致图像中所成的像发生形变；径向畸变的表述如下：

x_d＝x(1+k₁r²+k₂r⁴+k₃r⁶)

y_d＝y(1+k₁r²+k₂r⁴+k₃r⁶)

其中，r²＝x²+y²，k₁、k₂、k₃为径向畸变参数；

图像切向畸变是由于摄像机制造上的缺陷使得透镜本身与图像平面不平行而产生的，定量描述为：

x_d＝x+(2p₁xy+p₂(r²+2x²))

y_d＝y+(p₁(r²+2y²)+2p₂xy)

其中，p₁、p₂为切向畸变系数；

畸变前后的坐标关系如下：

x_d＝x(1+k₁r²+k₂r⁴+k₃r⁶)+(2p₁xy+p₂(r²+2x²))

y_d＝y(1+k₁r²+k₂r⁴+k₃r⁶)+(p₁(r²+2y²)+2p₂xy)

其中，(x，y)是理想状态下的正规坐标，(x_d，y_d)是实际带有畸变的正规坐标；

1-3)根据原来两相机的旋转关系将两图转回来：已知原来两个相机之间的旋转矩阵R和平移向量t，使得：

X_r＝RX_l+t

其中，X_l表示红外相机的正规坐标，X_r表示可见光相机的正规坐标；将可见光图像旋转R反方向一半的角度，将红外图像旋转R正方向一半的角度；

1-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系。

3.根据权利要求1所述的基于双向单应矩阵的可见光与红外相机动态位姿计算方法，其特征在于，步骤3)中检测Harris角点，包括以下步骤：

3-1)构建梯度矩阵M；将窗口移动[u，v]时，灰度变化如下所示：

将上式泰勒展开，得：

其中，ω(x，y)表示(x，y)点对应的权值，权值是常数，或者是高斯核的对应系数；I_x和I_y分别表示图像在x方向和y方向上的梯度，矩阵M表示为：

3-2)计算矩阵M的两个特征值λ₁和λ₂，λ₁和λ₂所对应的特征向量分别代表着灰度变化最快和最慢的两个方向；

3-3)使用一个值R来描述该点的角点响应，然后通过阈值σ1和σ2来判断该点是不是一个角点；角点响应值R＝det(M)-k*trace(M)²，其中det(M)表示矩阵M对应的行列式的值，trace(M)表示矩阵M的迹，即：

det(M)＝λ₁*λ₂

trace(M)＝λ₁+λ₂

其中k是常数，一般取0.04～0.06

当|R|＜σ1时，该区域是平面；

当R＜0时，该区域是直线；

当R＞σ2时，该区域是角点。

4.根据权利要求1所述的基于双向单应矩阵的可见光与红外相机动态位姿计算方法，其特征在于，步骤4)中寻找上一步检测出的Harris角点的匹配点，包括以下步骤：

4-1)记步骤3)在红外图像检测出来的角点集合为

对集合

里的每个角点

根据步骤2)计算出的单应矩阵计算

在可见光图像上的初始对应点，在初始对应点周围w*w窗口内检测Harris角点，记角点集为

如果

为空，则取下一个角点，重复本步骤，否则计算

中各个点和

的相似度，如果最大的相似度大于t1，则将其对应点

作为

在可见光图像上的初始对应点并加入到匹配点集中，选取下一点重复该步骤直到遍历所有红外角点；

4-2)记步骤3)在可见光图像检测出来的角点集合为

对集合

里的每个角点

根据步骤2)计算出的单应矩阵计算

在红外图像上的初始对应点，在初始对应点周围w*w窗口内检测Harris角点，记角点集为

如果

为空，则取下一个角点，重复本步骤，否则计算

中各个点和

的相似度，如果最大的相似度大于t1，则将其对应点

作为

在可见光图像上的初始对应点并加入到匹配点集中，选取下一点重复该步骤直到遍历所有可见光角点；

4-3)将匹配点集去重；

4-4)以红外图特征点

为基准，抛物线拟合优化对应可见光图的整数像素特征点

得到的对应可见光图的亚像素特征点

其中

为x方向上的亚像素偏移量，

为y方向上的亚像素偏移量；

4-5)以对应可见光图整数像素特征点

为基准，根据4-4)的方法计算出对应红外图的亚像素特征点

其中

为x方向上的亚像素偏移量，

为y方向上的亚像素偏移量；

4-6)得到最终的匹配点对为

5.根据权利要求1所述的基于双向单应矩阵的可见光与红外相机动态位姿计算方法，其特征在于，步骤6)中修正标定结果，包括以下步骤：

6-1)使用随机抽样一致性对点对做进一步筛选；

6-2)求解基础矩阵F和本质矩阵E：红外和可见光对应像素点对u₁、u_r和基础矩阵F的关系是：

将对应点坐标代入上式，构建齐次线性方程组求解F；

基础矩阵和本质矩阵的关系是：

其中，K₁、K_r分别是红外和可见光相机的内参矩阵；

6-3)从本质矩阵分解出旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：

E＝[t]_×R

其中[t]_×表示t的叉乘矩阵；

将E做奇异值分解，得

定义两个矩阵

和

ZW＝∑

所以E可以写成以下两种形式

(1)E＝UZU^TUWV^T

令[t]_×＝UZU^T，R＝UWV^T

(2)E＝-UZU^TUW^TV^T

令[t]_×＝-UZU^T，R＝UW^TV^T

6-4)将分解出的旋转和平移关系叠加到原来的外参里面；

记去畸变前的旋转矩阵为R₀，平移向量为t₀＝(t_x，t_y，t_z)^T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′_x，t′_y，t′_z)^T；则新的R_new和t_new如下：

此外，还需要将t_new乘一个系数，使得t_new在x方向上的分量

Images