CN110910456A - 基于Harris角点互信息匹配的立体相机动态标定算法 - Google Patents

基于Harris角点互信息匹配的立体相机动态标定算法 Download PDF

Info

Publication number
CN110910456A
CN110910456A CN201911152551.9A CN201911152551A CN110910456A CN 110910456 A CN110910456 A CN 110910456A CN 201911152551 A CN201911152551 A CN 201911152551A CN 110910456 A CN110910456 A CN 110910456A
Authority
CN
China
Prior art keywords
image
point
corner
camera
matching
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.)
Granted
Application number
CN201911152551.9A
Other languages
English (en)
Other versions
CN110910456B (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.)
Dalian University of Technology
Original Assignee
Dalian University of Technology
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 Dalian University of Technology filed Critical Dalian University of Technology
Priority to CN201911152551.9A priority Critical patent/CN110910456B/zh
Publication of CN110910456A publication Critical patent/CN110910456A/zh
Application granted granted Critical
Publication of CN110910456B publication Critical patent/CN110910456B/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
    • G06T7/00—Image analysis
    • G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
    • G06T7/85—Stereo camera calibration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/80—Geometric correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • 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/10—Image acquisition modality
    • G06T2207/10004—Still image; Photographic image
    • G06T2207/10012—Stereo images
    • 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/10—Image acquisition modality
    • G06T2207/10048—Infrared image
    • 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/20112—Image segmentation details
    • G06T2207/20164—Salient point detection; Corner detection

Landscapes

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

Abstract

本发明属于图像处理和计算机视觉领域，涉及基于Harris角点互信息匹配的立体相机动态标定算法。所述方法包括下列步骤：第一步：Harris角点检测。第二步：基于互信息的角点匹配。第三步：原图校正：输入上一步得到的匹配的左右角点以及红外双目相机各自内参和原来的外参，计算左右两图的角点经过去畸校正后的坐标。第四步：判断角点覆盖区域：将图像分成m*n个格子，如果角点覆盖到所有格子，则进行下一步，否则继续拍摄图像，提取角点。本发明解决了由于温湿度、震动等因素造成红外双目相机位置关系的改变，具有速度快、结果精确、操作简单等优点。

Description

基于Harris角点互信息匹配的立体相机动态标定算法
技术领域
本发明属于图像处理和计算机视觉领域，涉及基于Harris角点互信息匹配的立体相机动态标定算法。
背景技术
立体视觉是计算机视觉领域的重要主题。其目的是重建场景的3D几何信息。在双目立体视觉中，左右摄像头用于模拟两只眼睛。通过计算双目图像之间的差异来计算深度图像。双目立体视觉具有效率高，准确度高，系统结构简单，成本低的优点。由于双目立体视觉需要匹配左右图像捕获点上的相同点，因此相机两个镜头的焦距和图像捕获中心，以及左右两个镜头之间的位置关系。为了得到以上数据，我们需要对相机进行标定。
红外线(Infrared)是波长介于微波与可见光之间的电磁波，波长比红光要长。高于绝对零度(-273.15℃)的物质都可以产生红外线。红外图像由于其具有透过雾、雨等进行观察的能力而被广泛用于军事国防、资源勘探、气象预报、环境监测、医学诊治、海洋研究等不同领域。利用红外线可以隔着薄雾和烟雾拍摄景物，而且在夜间也可以进行红外摄影。将红外双目相机进行标定和校正，可以在低光、浓雾及雨雪等恶劣环境下估计视差及深度，从而实现全天候的三维视觉感知。
在标定过程中获得了相机的两个镜头参数和相对位置参数，但这些参数不稳定。当温度，湿度等发生变化时，相机镜头的内部参数也会发生变化。另外，由于意外的相机碰撞，两个镜头之间的位置关系可能会改变。因此，每次使用摄像机时，都必须修改内部和外部参数，这就是自标定。在已知相机内部参数的情况下，通过提取红外图像的角点来对两个红外镜头的位置关系进行修正，即红外双目相机的自标定。
发明内容
本发明旨在解决由于温湿度、震动等因素造成红外双目相机位置关系的改变。通过提取左右红外相机各自的角点并进行匹配，并根据这些角点对原有的标定结果进行修正。
基于Harris角点互信息匹配的立体相机动态标定算法，包括下列步骤：
第一步：Harris角点检测：使用红外双目相机拍摄场景图像，并在红外图像上检测Harris角点以待匹配。
第一步中Harris角点检测，具体包括以下步骤：
1-1)使用左右相机拍摄图像，获取左右图，分别在左右图上进行角点检测。
1-2)为图像上的每一个像素点构建梯度矩阵M。
在图像上，角点表现为不同边缘之间的交点。而且不论从什么角度来观测，它都是不同边缘之间的交点，不会因为视角的变化而变化。此外，角点邻域内的点的梯度会有大幅变化。角点应满足：当移动窗口时，角点所在窗口与其周围各个方向窗口的亮度分布差别很大。将窗口移动[u,v]时，灰度变化如下所示：
Figure BDA0002283946090000021
将上式泰勒展开，得：
Figure BDA0002283946090000022
其中，(x,y)表示窗口内的一点，ω(x,y)表示(x,y)点对应的权值，权值可以是常数，也可以是高斯核的对应系数。Ix和Iy分别表示图像(x,y)点在x方向和y方向上的梯度，矩阵M可表示为：
Figure BDA0002283946090000023
矩阵M计算方法如下：
计算图像I在x方向和y方向上的梯度图像：
Figure BDA0002283946090000024
Figure BDA0002283946090000025
其中
Figure BDA0002283946090000026
表示卷积。
Figure BDA0002283946090000027
1-3)根据每一个像素点的矩阵M来判断该像素点是否为角点。
计算矩阵M的两个特征值λ1和λ2，λ1和λ2所对应的特征向量分别代表着灰度变化最快和最慢的两个方向。λ1和λ2的大小关系和对应点的属性存在以下的对应关系，如图2所示：
(1)当λ1和λ2的值都很小时，该点落在平滑区域内。
(2)当λ1>>λ2或者λ2>>λ1时，该点落在图像的边缘上。
(3)当λ1和λ2的值都很大，且处于同一大小水平时，可以认为该点属于角点。
使用一个值R来描述该点的角点响应，然后通过阈值σ1和σ2来判断该点是不是一个角点。角点响应值R＝det(M)-k*trace(M)2，其中det(M)表示矩阵M对应的行列式的值，trace(M)表示矩阵M的迹，即：
det(M)＝λ1*λ2
trace(M)＝λ1+λ2
其中k是常数，一般取0.04～0.06
当|R|<σ1时，该区域是平面。
当R<0时，该区域是直线。
当R＞σ2时，该区域是角点。
1-4)记左图的Harris角点集为
Figure BDA0002283946090000031
右图的Harris角点集为
Figure BDA0002283946090000032
第二步：基于互信息的角点匹配：输入左右灰度图像以及上一步得到的左右图像的角点集
Figure BDA0002283946090000033
和
Figure BDA0002283946090000034
根据匹配窗口分别计算左右相机带匹配点对的互信息并根据其计算结果进行匹配，得到匹配关系{(Pl,Pr)}。
第二步中角点匹配，具体包括以下步骤：
2-1)将左图像和右图像都分为m×n个块。对于左图每一个角点
Figure BDA0002283946090000035
进行步骤2-2)～2-6)。
2-2)找到
Figure BDA0002283946090000036
在左图对应的块
Figure BDA0002283946090000037
(如图4(a)所示)。块
Figure BDA0002283946090000038
在右图相同位置的块为
Figure BDA0002283946090000039
与块
Figure BDA00022839460900000310
具有相同横坐标和纵坐标的块集合
Figure BDA00022839460900000311
(如图4(b)所示)，其角点集记为
Figure BDA00022839460900000312
我们使用像素的互信息来评估像素点之间的相似程度。互信息是一种对影像明暗变化不敏感的相关性测度，它通过两张影像各自的熵H以及两者的联合熵来定义，熵代表影像的信息量，图像的熵越大代表包含的像素灰度越丰富，灰度分布越均匀。像素点的互信息的计算公式如下：
Figure BDA0002283946090000041
其中
Figure BDA0002283946090000042
Figure BDA0002283946090000043
PI(i)表示图像I的灰度概率密度分布，
Figure BDA0002283946090000044
表示图像I1和I2的灰度联合概率分布。
Figure BDA0002283946090000045
表示卷积运算，g(i)和g(i,k)表示高斯核，n表示对应像素点个数。两个像素点p1和p2的相似程度s可以用互信息表示，即
Figure BDA0002283946090000046
其中I1(p1)和I2(p2)分别表示p1点在I1的灰度值以及p2点在I2的灰度值。
计算
Figure BDA0002283946090000047
和
Figure BDA0002283946090000048
中任意一点的互信息作为相似程度，如果相似程度大于阈值t1，则视为粗匹配点，其集合记为
Figure BDA0002283946090000049
否则舍弃该点，选取下一个角点重新进行步骤2-2)。
2-3)如果
Figure BDA00022839460900000410
和
Figure BDA00022839460900000411
中相似程度最大值sfirst和次大值ssecond满足：
F(sfirst,ssecond)≥t2
则保留该匹配，取
Figure BDA00022839460900000412
中相似程度最大的点
Figure BDA00022839460900000413
作为匹配点，其中t2为阈值,F(sfirst,ssecond)用于描述sfirst和ssecond之间的关系。如不满足，则舍弃该点，选取下一个角点重新进行步骤2-2)。
按照该规则筛选之后，再按照步骤2-2)～2-3)匹配
Figure BDA00022839460900000414
在左图对应的角点
Figure BDA00022839460900000415
如果满足
Figure BDA00022839460900000416
则保留该匹配
Figure BDA00022839460900000417
如不满足，则舍弃该点，选取下一个角点重新进行步骤2-2)。
2-4)以左图角点
Figure BDA00022839460900000418
为基准，抛物线拟合优化对应右图的整数像素角点
Figure BDA00022839460900000419
得到的对应右图的亚像素角点
Figure BDA00022839460900000420
其中
Figure BDA00022839460900000421
为x方向上的亚像素偏移量，
Figure BDA00022839460900000422
为y方向上的亚像素偏移量。
2-5)以对应右图整数像素角点
Figure BDA00022839460900000423
为基准，根据2-4)的方法计算出对应左图的亚像素角点
Figure BDA00022839460900000424
其中
Figure BDA00022839460900000425
为x方向上的亚像素偏移量，
Figure BDA00022839460900000426
为y方向上的亚像素偏移量。
2-6)得到最终的匹配点对为
Figure BDA00022839460900000427
选取下一个左图角点重新进行步骤2-2)～2-6)。
第三步：原图校正：输入上一步得到的匹配的左右角点以及红外双目相机各自内参和原来的外参，计算左右两图的角点经过去畸校正后的坐标。流程如图3所示。
步骤3)中原图校正，具体包括以下步骤：
3-1)计算匹配的左右角点
Figure BDA00022839460900000428
对应的正规坐标系下的坐标。
像素坐标系以图片的左上角为原点，其x轴和y轴分别与图像坐标系的x轴和y轴平行。像素坐标系的单位是像素，像素是图像显示的基本且不可分割的单位。正规坐标系以相机的光心作为图像坐标系的原点，且将光心到图像平面的距离缩放到1。像素坐标与正规坐标的关系如下：
u＝KX
Figure BDA0002283946090000051
其中，
Figure BDA0002283946090000052
表示图像的像素坐标；
Figure BDA0002283946090000053
表示相机的内参矩阵，fx和fy分别表示图像x方向和y方向的焦距(单位是像素)，(cx,cy)表示相机住店的位置；
Figure BDA0002283946090000054
是正规坐标系下的坐标。已知图像的像素坐标系以及相机的内参可以计算出像素点对应的正规坐标系，即
X＝K-1u
对于每一对左右相机匹配角点
Figure BDA0002283946090000055
它们的正规坐标系为：
Figure BDA0002283946090000056
Figure BDA0002283946090000057
其中，
Figure BDA0002283946090000058
和
Figure BDA0002283946090000059
分别是
Figure BDA00022839460900000510
和
Figure BDA00022839460900000511
的像素坐标，
Figure BDA00022839460900000512
和
Figure BDA00022839460900000513
分别是
Figure BDA00022839460900000514
和
Figure BDA00022839460900000515
的正规坐标，Kl和Kr分别是左相机和右相机的内参矩阵。
3-2)去除图像畸变：根据左右图像角点的正规坐标和左右相机各自的畸变系数来计算出左右图像角点去畸变后的正规坐标。
由于镜头生产工艺的限制，实际情况下的镜头会存在一些失真现象导致非线性的畸变，可大致分为径向畸变和切向畸变。
图像径向畸变是图像像素点以畸变中心为中心点，沿着径向产生的位置偏差，从而导致图像中所成的像发生形变。径向畸变的大致表述如下：
xd＝x(1+k1r2+k2r4+k3r6)
yd＝y(1+k1r2+k2r4+k3r6)
其中，r2＝x2+y2，k1、k2、k3为径向畸变参数。
切向畸变是由于摄像机制造上的缺陷使得透镜本身与图像平面不平行而产生的，可定量描述为：
xd＝x+(2p1xy+p2(r2+2x2))
yd＝y+(p1(r2+2y2)+2p2xy)
其中，p1、p2为切向畸变系数。
综上，畸变前后的坐标关系如下：
xd＝x(1+k1r2+k2r4+k3r6)+(2p1xy+p2(r2+2x2))
yd＝y(1+k1r2+k2r4+k3r6)+(p1(r2+2y2)+2p2xy)
其中，(x,y)是理想状态下的正规坐标，(xd,yd)是实际带有畸变的正规坐标。我们以(xd,yd)作为(x,y)的初值，迭代计算若干次(例如20次)得到实际的(x,y)。
3-3)根据原来两相机的旋转关系将左右两图旋转：已知原来左右相机之间的旋转矩阵R和平移向量t，使得
Xr＝RXl+t
其中，Xl表示左相机的正规坐标，Xr表示右相机的正规坐标。将左图旋转R正方向一半的角度，将右图旋转R反方向一半的角度。
对于上一步得到的每一对去畸变之后的左右角点
Figure BDA0002283946090000061
的正规坐标
Figure BDA0002283946090000062
Figure BDA0002283946090000063
3-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系。根据上一步更新的左右角点
Figure BDA0002283946090000064
的正规坐标
Figure BDA0002283946090000065
计算去畸校正后的图像坐标
Figure BDA0002283946090000066
Figure BDA0002283946090000067
第四步：判断角点覆盖区域：将图像分成m*n个格子，如果角点覆盖到所有格子，则进行下一步，否则继续拍摄图像，提取角点。
第五步：修正标定结果：使用所有角点的图像坐标来计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加。
步骤5)中修正外参结果，具体包括以下步骤：
5-1)根据左右两图的角点对坐标以及左右相机的内参矩阵来求解基础矩阵F和本质矩阵E：左右对应像素点对ul、ur和基础矩阵F的关系是：
Figure BDA0002283946090000071
使用随机抽样一致性(RANSAC)对点对做进一步筛选，之后将对应点坐标代入上式，构建齐次线性方程组求解F。
基础矩阵和本质矩阵的关系是：
Figure BDA0002283946090000072
其中，Kl、Kr分别是左右相机的内参矩阵。
5-2)从本质矩阵分解出校正之后左右相机旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：
E＝[t]×R
其中[t]×表示t的叉乘矩阵。
将E做奇异值分解，得
Figure BDA0002283946090000073
定义两个矩阵
Figure BDA0002283946090000074
和
Figure BDA0002283946090000075
ZW＝Σ
所以E可以写成以下两种形式
(1)E＝UZUTUWVT
令[t]×＝UZUT，R=UWVT
(2)E＝-UZUTUWTVT
令[t]×＝-UZUT，R=UWTVT
得到四对R和t，选取具有三维意义的解。
5-3)将分解出的旋转和平移关系叠加到原来的外参里面。
记去畸变前的旋转矩阵为R0，平移向量为t0＝(tx,ty,tz)T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′x,t′y,t′z)T。则新的Rnew和tnew如下
Figure BDA0002283946090000076
Figure BDA0002283946090000077
还需要将tnew乘一个系数，使得tnew在x方向上的分量
Figure BDA0002283946090000078
本发明的有益效果是：
本发明解决了由于温湿度、震动等因素造成红外双目相机位置关系的改变。具有速度快、结果精确、操作简单等优点。
附图说明
图1为整体流程示意图。
图2为Harris角点条件判断示意图。
图3为双目校正流程示意图。
图4(a)为分块匹配的左图示意图。
图4(b)为分块匹配的右图示意图。
具体实施方式
本发明旨在解决由于温湿度、震动等因素造成红外双目相机位置关系的改变。结合附图及实施例详细说明如下：
基于Harris角点互信息匹配的立体相机动态标定算法，包括下列步骤：
第一步：Harris角点检测：使用红外双目相机拍摄场景图像，并在红外图像上检测Harris角点以待匹配；
第一步中Harris角点检测，具体包括以下步骤：
1-3)使用左右相机拍摄图像，获取左右图，分别在左右图上进行角点检测；
1-4)为图像上的每一个像素点构建梯度矩阵M；
矩阵M计算方法如下：
计算图像I在x方向和y方向上的梯度图像：
Figure BDA0002283946090000081
Figure BDA0002283946090000082
其中
Figure BDA0002283946090000083
表示卷积；
Figure BDA0002283946090000084
1-3)根据每一个像素点的矩阵M来判断该像素点是否为角点；
使用一个值R来描述该点的角点响应，然后通过阈值σ1和σ2来判断该点是角点，角点响应值R＝det(M)-k*trace(M)2，其中det(M)表示矩阵M对应的行列式的值，trace(M)表示矩阵M的迹，即：
det(M)＝λ1*λ2
trace(M)＝λ1+λ2
其中k是常数，一般取0.04～0.06
当|R|<σ1时，该区域是平面；
当R<0时，该区域是直线；
当R＞σ2时，该区域是角点；
1-4)记左图的Harris角点集为
Figure BDA0002283946090000091
右图的Harris角点集为
Figure BDA0002283946090000092
第二步：基于互信息的角点匹配：输入左右灰度图像以及上一步得到的左右图像的角点集
Figure BDA0002283946090000093
和
Figure BDA0002283946090000094
根据匹配窗口分别计算左右相机带匹配点对的互信息并根据其计算结果进行匹配，得到匹配关系{(Pl,Pr)}；
2-1)将左图像和右图像都分为m×n个块；对于左图每一个角点
Figure BDA0002283946090000095
进行步骤2-2)～2-6)；
2-2)找到
Figure BDA0002283946090000096
在左图对应的块
Figure BDA0002283946090000097
块
Figure BDA0002283946090000098
在右图相同位置的块为
Figure BDA0002283946090000099
与块
Figure BDA00022839460900000910
具有相同横坐标和纵坐标的块集合
Figure BDA00022839460900000911
其角点集记为
Figure BDA00022839460900000912
计算
Figure BDA00022839460900000913
和
Figure BDA00022839460900000914
中任意一点的互信息作为相似程度，如果相似程度大于阈值t1，则视为粗匹配点，其集合记为
Figure BDA00022839460900000915
否则舍弃该点，选取下一个角点重新进行步骤2-2)；
2-3)如果
Figure BDA00022839460900000916
和
Figure BDA00022839460900000917
中相似程度最大值sfirst和次大值ssecond满足：
F(sfirst,ssecond)≥t2
则保留该匹配，取
Figure BDA00022839460900000918
中相似程度最大的点
Figure BDA00022839460900000919
作为匹配点，其中t2为阈值,F(sfirst,ssecond)用于描述sfirst和ssecond之间的关系；如不满足，则舍弃该点，选取下一个角点重新进行步骤2-2)；
按照该规则筛选之后，再按照步骤2-2)～2-3)匹配
Figure BDA00022839460900000920
在左图对应的角点
Figure BDA00022839460900000921
如果满足
Figure BDA00022839460900000922
则保留该匹配
Figure BDA00022839460900000923
如不满足，则舍弃该点，选取下一个角点重新进行步骤2-2)；
2-4)以左图角点
Figure BDA00022839460900000924
为基准，抛物线拟合优化对应右图的整数像素角点
Figure BDA00022839460900000925
得到的对应右图的亚像素角点
Figure BDA00022839460900000926
其中
Figure BDA00022839460900000927
为x方向上的亚像素偏移量，
Figure BDA00022839460900000928
为y方向上的亚像素偏移量；
2-5)以对应右图整数像素角点
Figure BDA0002283946090000101
为基准，根据2-4)的方法计算出对应左图的亚像素角点
Figure BDA0002283946090000102
其中
Figure BDA0002283946090000103
为x方向上的亚像素偏移量，
Figure BDA0002283946090000104
为y方向上的亚像素偏移量；
2-6)得到最终的匹配点对为
Figure BDA0002283946090000105
选取下一个左图角点重新进行步骤2-2)～2-6)。
第三步：原图校正：
输入上一步得到的匹配的左右角点以及红外双目相机各自内参和原来的外参，计算左右两图的角点经过去畸校正后的坐标；
第三步中原图校正，具体包括以下步骤：
3-1)计算匹配的左右角点
Figure BDA0002283946090000106
对应的正规坐标系下的坐标；
对于每一对左右相机匹配角点
Figure BDA0002283946090000107
它们的正规坐标系为：
Figure BDA0002283946090000108
Figure BDA0002283946090000109
其中，
Figure BDA00022839460900001010
和
Figure BDA00022839460900001011
分别是
Figure BDA00022839460900001012
和
Figure BDA00022839460900001013
的像素坐标，
Figure BDA00022839460900001014
和
Figure BDA00022839460900001015
分别是
Figure BDA00022839460900001016
和
Figure BDA00022839460900001017
的正规坐标，Kl和Kr分别是左相机和右相机的内参矩阵；
3-2)去除图像畸变：根据左右图像角点的正规坐标和左右相机各自的畸变系数来计算出左右图像角点去畸变后的正规坐标；
以(xd,yd)作为(x,y)的初值，迭代计算若干次得到实际的(x,y)；
3-3)根据原来两相机的旋转关系将左右两图旋转：已知原来左右相机之间的旋转矩阵R和平移向量t，使得
Xr＝RXl+t
其中，Xl表示左相机的正规坐标，Xr表示右相机的正规坐标；将左图旋转R正方向一半的角度，将右图旋转R反方向一半的角度；
对于上一步得到的每一对去畸变之后的左右角点
Figure BDA00022839460900001018
的正规坐标
Figure BDA00022839460900001019
Figure BDA00022839460900001020
3-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系；根据上一步更新的左右角点
Figure BDA00022839460900001021
的正规坐标
Figure BDA00022839460900001022
计算去畸校正后的图像坐标
Figure BDA00022839460900001023
Figure BDA00022839460900001024
第四步：判断角点覆盖区域：将图像分成m*n个格子，如果角点覆盖到所有格子，则进行下一步，否则继续拍摄图像，提取角点；
第五步：修正标定结果：使用所有角点的图像坐标来计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加；
步骤5)中修正外参结果，具体包括以下步骤：
5-1)使用随机抽样一致性(RANSAC)对点对做进一步筛选，之后将对应点坐标代入上式，构建齐次线性方程组求解F；
基础矩阵和本质矩阵的关系是：
Figure BDA0002283946090000111
其中，Kl、Kr分别是左右相机的内参矩阵；
5-2)从本质矩阵分解出校正之后左右相机旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：
E＝[t]×R
其中[t]×表示t的叉乘矩阵；
将E做奇异值分解，得
Figure BDA0002283946090000112
定义两个矩阵
Figure BDA0002283946090000113
和
Figure BDA0002283946090000114
ZW＝Σ
所以E可以写成以下两种形式
(1)E＝UZUTUWVT
令[t]×＝UZUT，R＝UWVT
(2)E＝-UZUTUWTVT
令[t]×＝-UZUT，R=UWTVT
得到四对R和t，选取具有三维意义的解；
5-3)将分解出的旋转和平移关系叠加到原来的外参里面；
去畸变前的旋转矩阵R和平移向量t
Figure BDA0002283946090000115
t＝[-335.5808 1.5591 -0.4805]T
上一步计算出的旋转矩阵为R′和平移向量为t′
Figure BDA0002283946090000121
t′＝[-1.0000 -0.0021 -0.0042]T
新的Rnew和tnew
Figure BDA0002283946090000122
tnew＝[-335.5808 -1.4520 -0.4218]T。

Claims (3)

1.基于Harris角点互信息匹配的立体相机动态标定算法，其特征在于，包括下列步骤：
第一步：Harris角点检测：使用红外双目相机拍摄场景图像，并在红外图像上检测Harris角点以待匹配；
第一步中Harris角点检测，具体包括以下步骤：
1-1)使用左右相机拍摄图像，获取左右图，分别在左右图上进行角点检测；
1-2)为图像上的每一个像素点构建梯度矩阵M；
矩阵M计算方法如下：
计算图像I在x方向和y方向上的梯度图像：
Figure FDA0002283946080000011
Figure FDA0002283946080000012
其中
Figure FDA0002283946080000013
表示卷积；
Figure FDA0002283946080000014
1-3)根据每一个像素点的矩阵M来判断该像素点是否为角点；
使用一个值R来描述该点的角点响应，然后通过阈值σ1和σ2来判断该点是角点，角点响应值R＝det(M)-k*trace(M)2，其中det(M)表示矩阵M对应的行列式的值，trace(M)表示矩阵M的迹，即：
det(M)＝λ1*λ2
trace(M)＝λ1+λ2
其中k是常数，一般取0.04～0.06
当|R|<σ1时，该区域是平面；
当R<0时，该区域是直线；
当R＞σ2时，该区域是角点；
1-4)记左图的Harris角点集为
Figure FDA0002283946080000015
右图的Harris角点集为
Figure FDA0002283946080000016
第二步：基于互信息的角点匹配：输入左右灰度图像以及上一步得到的左右图像的角点集
Figure FDA0002283946080000017
和
Figure FDA0002283946080000018
根据匹配窗口分别计算左右相机带匹配点对的互信息并根据其计算结果进行匹配，得到匹配关系{(Pl,Pr)}；
第三步：原图校正：
输入上一步得到的匹配的左右角点以及红外双目相机各自内参和原来的外参，计算左右两图的角点经过去畸校正后的坐标；
第四步：判断角点覆盖区域：将图像分成m*n个格子，如果角点覆盖到所有格子，则进行下一步，否则继续拍摄图像，提取角点；
第五步：修正标定结果：使用所有角点的图像坐标来计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加；
步骤5)中修正外参结果，具体包括以下步骤：
5-1)使用随机抽样一致性(RANSAC)对点对做进一步筛选，之后将对应点坐标代入上式，构建齐次线性方程组求解F；
基础矩阵和本质矩阵的关系是：
Figure FDA0002283946080000021
其中，K1、Kr分别是左右相机的内参矩阵；
5-2)从本质矩阵分解出校正之后左右相机旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：
E＝[t]×R
其中[t]×表示t的叉乘矩阵；
将E做奇异值分解，得
Figure FDA0002283946080000022
定义两个矩阵
Figure FDA0002283946080000023
和
Figure FDA0002283946080000024
ZW＝∑
所以E可以写成以下两种形式
(1)E＝UZUTUWVT
令[t]×＝UZUT，R＝UWVT
(2)E＝-UZUTUWTVT
令[t]×＝-UZUT，R＝UWTVT
得到四对R和t，选取具有三维意义的解；
5-3)将分解出的旋转和平移关系叠加到原来的外参里面；
记去畸变前的旋转矩阵为R0，平移向量为t0＝(tx，ty，tz)T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′x，t′y，t′z)T，则新的Rnew和tnew如下
Figure FDA0002283946080000025
Figure FDA0002283946080000026
还需要将tnew乘一个系数，使得tnew在x方向上的分量
Figure FDA0002283946080000027
2.如权利要求1所述的基于Harris角点互信息匹配的立体相机动态标定算法，其特征在于，第二步中角点匹配，具体包括以下步骤：
2-1)将左图像和右图像都分为m×n个块；对于左图每一个角点
Figure FDA0002283946080000031
进行步骤2-2)～2-6)；
2-2)找到
Figure FDA0002283946080000032
在左图对应的块
Figure FDA0002283946080000033
块
Figure FDA0002283946080000034
在右图相同位置的块为
Figure FDA0002283946080000035
与块
Figure FDA00022839460800000341
具有相同横坐标和纵坐标的块集合
Figure FDA0002283946080000036
其角点集记为
Figure FDA0002283946080000037
计算
Figure FDA0002283946080000038
和
Figure FDA0002283946080000039
中任意一点的互信息作为相似程度，如果相似程度大于阈值t1，则视为粗匹配点，其集合记为
Figure FDA00022839460800000310
否则舍弃该点，选取下一个角点重新进行步骤2-2)；
2-3)如果
Figure FDA00022839460800000311
和
Figure FDA00022839460800000312
中相似程度最大值sfirst和次大值ssecond满足：
F(sfirst,ssecond)≥t2
则保留该匹配，取
Figure FDA00022839460800000313
中相似程度最大的点
Figure FDA00022839460800000314
作为匹配点，其中t2为阈值,F(sfirst,ssecond)用于描述sfirst和ssecond之间的关系；如不满足，则舍弃该点，选取下一个角点重新进行步骤2-2)；
按照该规则筛选之后，再按照步骤2-2)～2-3)匹配
Figure FDA00022839460800000315
在左图对应的角点
Figure FDA00022839460800000316
如果满足
Figure FDA00022839460800000317
则保留该匹配
Figure FDA00022839460800000318
如不满足，则舍弃该点，选取下一个角点重新进行步骤2-2)；
2-4)以左图角点
Figure FDA00022839460800000319
为基准，抛物线拟合优化对应右图的整数像素角点
Figure FDA00022839460800000320
得到的对应右图的亚像素角点
Figure FDA00022839460800000321
其中
Figure FDA00022839460800000322
为x方向上的亚像素偏移量，
Figure FDA00022839460800000323
为y方向上的亚像素偏移量；
2-5)以对应右图整数像素角点
Figure FDA00022839460800000324
为基准，根据2-4)的方法计算出对应左图的亚像素角点
Figure FDA00022839460800000325
其中
Figure FDA00022839460800000326
为x方向上的亚像素偏移量，
Figure FDA00022839460800000327
为y方向上的亚像素偏移量；
2-6)得到最终的匹配点对为
Figure FDA00022839460800000328
选取下一个左图角点重新进行步骤2-2)～2-6)。
3.如权利要求2所述的基于Harris角点互信息匹配的立体相机动态标定算法，其特征在于，第三步中原图校正，具体包括以下步骤：
3-1)计算匹配的左右角点
Figure FDA00022839460800000329
对应的正规坐标系下的坐标；
对于每一对左右相机匹配角点
Figure FDA00022839460800000330
它们的正规坐标系为：
Figure FDA00022839460800000331
Figure FDA00022839460800000332
其中，
Figure FDA00022839460800000333
和
Figure FDA00022839460800000334
分别是Pi l和
Figure FDA00022839460800000336
的像素坐标，
Figure FDA00022839460800000337
和
Figure FDA00022839460800000338
分别是Pi l和
Figure FDA00022839460800000340
的正规坐标，Kl和Kr分别是左相机和右相机的内参矩阵；
3-2)去除图像畸变：根据左右图像角点的正规坐标和左右相机各自的畸变系数来计算出左右图像角点去畸变后的正规坐标；
以(xd,yd)作为(x,y)的初值，迭代计算若干次得到实际的(x,y)；
3-3)根据原来两相机的旋转关系将左右两图旋转：已知原来左右相机之间的旋转矩阵R和平移向量t，使得
Xr＝RXl+t
其中，Xl表示左相机的正规坐标，Xr表示右相机的正规坐标；将左图旋转R正方向一半的角度，将右图旋转R反方向一半的角度；
对于上一步得到的每一对去畸变之后的左右角点
Figure FDA0002283946080000041
的正规坐标
Figure FDA0002283946080000042
Figure FDA0002283946080000043
3-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系；根据上一步更新的左右角点
Figure FDA0002283946080000044
的正规坐标
Figure FDA0002283946080000045
计算去畸校正后的图像坐标
Figure FDA0002283946080000046
Figure FDA0002283946080000047
CN201911152551.9A 2019-11-22 2019-11-22 基于Harris角点互信息匹配的立体相机动态标定方法 Active CN110910456B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911152551.9A CN110910456B (zh) 2019-11-22 2019-11-22 基于Harris角点互信息匹配的立体相机动态标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911152551.9A CN110910456B (zh) 2019-11-22 2019-11-22 基于Harris角点互信息匹配的立体相机动态标定方法

Publications (2)

Publication Number Publication Date
CN110910456A true CN110910456A (zh) 2020-03-24
CN110910456B CN110910456B (zh) 2020-09-29

Family

ID=69818903

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911152551.9A Active CN110910456B (zh) 2019-11-22 2019-11-22 基于Harris角点互信息匹配的立体相机动态标定方法

Country Status (1)

Country Link
CN (1) CN110910456B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113284189A (zh) * 2021-05-12 2021-08-20 深圳市格灵精睿视觉有限公司 畸变参数标定方法、装置、设备及存储介质
CN113409399A (zh) * 2021-06-10 2021-09-17 武汉库柏特科技有限公司 一种双相机联合标定方法、系统及装置
CN113450416A (zh) * 2020-06-15 2021-09-28 天津工业大学 一种应用于三目相机立体标定的tcsc方法
CN113766209A (zh) * 2020-05-29 2021-12-07 上海汉时信息科技有限公司 相机偏移量处理方法及装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR100755450B1 (ko) * 2006-07-04 2007-09-04 중앙대학교 산학협력단 평면 호모그래피를 이용한 3차원 재구성 장치 및 방법
US20100259624A1 (en) * 2007-10-24 2010-10-14 Kai Li Method and apparatus for calibrating video camera
CN102509304A (zh) * 2011-11-24 2012-06-20 江南大学 基于智能优化的摄像机标定方法
EP2660776A1 (en) * 2012-05-01 2013-11-06 Universität Bern Image distortion correction and robust phantom detection
CN109064516A (zh) * 2018-06-28 2018-12-21 北京航空航天大学 一种基于绝对二次曲线像的相机自标定方法
CN110456330A (zh) * 2019-08-27 2019-11-15 中国人民解放军国防科技大学 一种相机与激光雷达之间外参无目标自动标定方法及系统

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR100755450B1 (ko) * 2006-07-04 2007-09-04 중앙대학교 산학협력단 평면 호모그래피를 이용한 3차원 재구성 장치 및 방법
US20100259624A1 (en) * 2007-10-24 2010-10-14 Kai Li Method and apparatus for calibrating video camera
CN102509304A (zh) * 2011-11-24 2012-06-20 江南大学 基于智能优化的摄像机标定方法
EP2660776A1 (en) * 2012-05-01 2013-11-06 Universität Bern Image distortion correction and robust phantom detection
CN109064516A (zh) * 2018-06-28 2018-12-21 北京航空航天大学 一种基于绝对二次曲线像的相机自标定方法
CN110456330A (zh) * 2019-08-27 2019-11-15 中国人民解放军国防科技大学 一种相机与激光雷达之间外参无目标自动标定方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
PANDEY 等: "Automatic Targetless Extrinsic Calibration of a 3D Lidar and Camera by Maximizing Mutual Information", 《 TWENTY-SIXTH AAAI CONFERENCE ON ARTIFICIAL INTELLIGENCE. AAAI PRESS》 *
韩松 等: "基于自适应互信息的红外/深度双摄像机时空配准", 《华东理工大学学报(自然科学版)》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113766209A (zh) * 2020-05-29 2021-12-07 上海汉时信息科技有限公司 相机偏移量处理方法及装置
CN113766209B (zh) * 2020-05-29 2024-04-30 上海汉时信息科技有限公司 相机偏移量处理方法及装置
CN113450416A (zh) * 2020-06-15 2021-09-28 天津工业大学 一种应用于三目相机立体标定的tcsc方法
CN113450416B (zh) * 2020-06-15 2024-03-15 天津工业大学 一种应用于三目相机立体标定的tcsc方法
CN113284189A (zh) * 2021-05-12 2021-08-20 深圳市格灵精睿视觉有限公司 畸变参数标定方法、装置、设备及存储介质
CN113409399A (zh) * 2021-06-10 2021-09-17 武汉库柏特科技有限公司 一种双相机联合标定方法、系统及装置

Also Published As

Publication number Publication date
CN110910456B (zh) 2020-09-29

Similar Documents

Publication Publication Date Title
CN110910456B (zh) 基于Harris角点互信息匹配的立体相机动态标定方法
CN110969668B (zh) 一种长焦双目相机的立体标定算法
WO2021098083A1 (zh) 基于显著特征的多光谱相机动态立体标定算法
US11398053B2 (en) Multispectral camera external parameter self-calibration algorithm based on edge features
CN110956661B (zh) 基于双向单应矩阵的可见光与红外相机动态位姿计算方法
CN110135455B (zh) 影像匹配方法、装置及计算机可读存储介质
CN110969669B (zh) 基于互信息配准的可见光与红外相机联合标定方法
CN111080709B (zh) 基于轨迹特征配准的多光谱立体相机自标定算法
CN110880191B (zh) 基于直方图均衡化的红外立体相机动态外参计算方法
CN107146200B (zh) 一种基于图像拼接质量评价的无人机遥感图像拼接方法
CN110992409B (zh) 基于傅里叶变换配准的多光谱立体相机动态配准方法
CN107767339B (zh) 一种双目立体图像拼接方法
WO2021017588A1 (zh) 一种基于傅立叶频谱提取的图像融合方法
CN107560592B (zh) 一种用于光电跟踪仪联动目标的精确测距方法
CN113744315B (zh) 一种基于双目视觉的半直接视觉里程计
CN112016478B (zh) 一种基于多光谱图像融合的复杂场景识别方法及系统
CN111899345B (zh) 一种基于2d视觉图像的三维重建方法
CN110910457B (zh) 基于角点特征的多光谱立体相机外参计算方法
CN111047636A (zh) 基于主动红外双目视觉的避障系统和避障方法
Hsu et al. Object detection using structure-preserving wavelet pyramid reflection removal network
CN116433822B (zh) 一种神经辐射场训练方法、装置、设备及介质
CN113962904B (zh) 一种高光谱图像滤波降噪的方法
CN113670268B (zh) 基于双目视觉的无人机和电力杆塔距离测量方法
CN113808070A (zh) 一种双目数字散斑图像相关的视差测量方法
CN113763261A (zh) 一种海雾气象条件下的远小目标实时检测方法

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