CN110969669B - 基于互信息配准的可见光与红外相机联合标定方法 - Google Patents

基于互信息配准的可见光与红外相机联合标定方法 Download PDF

Info

Publication number
CN110969669B
CN110969669B CN201911153787.4A CN201911153787A CN110969669B CN 110969669 B CN110969669 B CN 110969669B CN 201911153787 A CN201911153787 A CN 201911153787A CN 110969669 B CN110969669 B CN 110969669B
Authority
CN
China
Prior art keywords
image
visible light
infrared
pixel
point
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
CN201911153787.4A
Other languages
English (en)
Other versions
CN110969669A (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 CN201911153787.4A priority Critical patent/CN110969669B/zh
Publication of CN110969669A publication Critical patent/CN110969669A/zh
Application granted granted Critical
Publication of CN110969669B publication Critical patent/CN110969669B/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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • 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/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10041—Panchromatic 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/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/20036—Morphological image processing
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30244—Camera pose

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了基于互信息配准的可见光与红外相机联合标定方法，属于图像处理和计算机视觉领域。通过提取并匹配特征点来对外参进行修正。为了缩小特征点的匹配范围，在特征点检测之前，本发明通过互信息来对红外图像和可见光图像进行配准。根据配准之后的结果进行匹配，可以有效利用可见光图像和红外图像之间的位置关系，从而更加有效的对红外相机和可见光相机进行联合自标定，操作简便，结果精确。

Description

基于互信息配准的可见光与红外相机联合标定方法
技术领域
本发明属于图像处理和计算机视觉领域，涉及基于互信息配准的可见光与红外相机联合标定方法。
背景技术
红外线(Infrared)是波长介于微波与可见光之间的电磁波，波长比红光要长。高于绝对零度(-273.15℃)的物质都可以产生红外线。红外图像由于其具有透过雾、雨等进行观察的能力而被广泛用于军事国防、资源勘探、气象预报、环境监测、医学诊治、海洋研究等不同领域。利用红外线可以隔着薄雾和烟雾拍摄景物，而且在夜间也可以进行红外摄影。红外相机成像的优点是在极端场景(低光、雨雪、浓雾等)也可以成像，缺点是分辨率低、图像细节较模糊。相比之下，可见光相机的优点是分辨率高、图像细节清晰，但是在极端场景下不能成像。因此，将红外相机和可见光相机结合起来具有重大的现实意义。
立体视觉是计算机视觉领域的重要主题。其目的是重建场景的3D几何信息。双目立体视觉是立体视觉的重要领域。在双目立体视觉中，左右摄像头用于模拟两只眼睛。通过计算双目图像之间的差异来计算深度图像。双目立体视觉具有效率高，准确度高，系统结构简单，成本低的优点。由于双目立体视觉需要匹配左右图像捕获点上的相同点，因此相机两个镜头的焦距和图像捕获中心，以及左右两个镜头之间的位置关系。为了得到以上数据，我们需要对相机进行标定。获取可见光相机和红外相机之间的位置关系称为联合标定。
在标定过程中获得了相机的两个镜头参数和相对位置参数，但这些参数不稳定。当温度、湿度等发生变化时，相机镜头的内部参数也会发生变化。另外，由于意外的相机碰撞，两个镜头之间的位置关系可能会改变。因此，每次使用摄像机时，都必须修改内部和外部参数，这就是自标定。在已知相机内部参数的情况下，通过分别提取红外图像特征和可见光图像特征来对红外镜头和可见光镜头的位置关系进行修正，即红外相机与可见光相机的联合自标定。
为了缩小特征点的匹配范围，在特征点检测之前，通过互信息来对红外图像和可见光图像进行配准。互信息(Mutual Information)是信息论里一种有用的信息度量，它可以看成是一个随机变量中包含的关于另一个随机变量的信息量，或者说是一个随机变量由于已知另一个随机变量而减少的不肯定性。
发明内容
本发明旨在解决由于温湿度、震动等因素造成红外相机和可见光相机位置关系的改变。从拍摄到的红外场景图像和可见光场景图像中提取并匹配特征点，并根据匹配的特征点对红外相机和可见光相机之间的位置关系进行修正，从而解决红外相机和可见光相机因温度和震动导致其外参发生变化的问题。
技术方案为，基于互信息配准的可见光与红外相机联合标定方法，包括如下步骤：
1)原图校正：将原图根据红外相机和可见光相机各自内参和原来的外参进行去畸变和双目校正；
2)根据红外图像与可见光图像的互信息确定红外图像在可见光图像的最佳对应位置，即对红外图像和可见光图像进行配准；
3)在配准后的红外图像和可见光图像上分别提取特征点；
4)对上一步提取的特征点进行匹配；
5)根据步骤2)和步骤3)的结果计算配准后的红外图像特征点对应于红外原图像的特征点；
6)判断特征点覆盖区域：将图像分成m*n个格子，如果特征点覆盖到所有格子，则进行下一步，否则继续拍摄图像，重复步骤1)～步骤5)；
7)修正标定结果：使用所有特征点的图像坐标计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加。
所述步骤1)原图校正，具体包括以下步骤：
1-1)计算图像的像素点对应的正规坐标系下的坐标。其中，正规坐标系是相机坐标系在平面Z＝1的投影；而相机坐标系是以相机的中心作为图像坐标系的原点，以图片方向为XY轴方向，以垂直于图像为Z轴方向的坐标系。像素坐标系以图片的左上角为原点，其x轴和y轴分别与图像坐标系的x轴和y轴平行。像素坐标系的单位是像素。像素坐标与正规坐标的关系如下：
u＝KX
Figure BDA0002284263860000031
其中，
Figure BDA0002284263860000032
表示图像的像素坐标；
Figure BDA0002284263860000033
表示相机的内参矩阵，fx和fy分别表示图像x方向和y方向的焦距，单位是像素，(cx,cy)表示相机的主点位置，即相机中心在图像上的对应位置；
Figure BDA0002284263860000034
是正规坐标系下的坐标。已知图像的像素坐标系以及相机的内参计算出像素点对应的正规坐标系，即X＝K-1u；
1-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)是实际带有畸变的正规坐标。
1-3)根据原来两相机的旋转关系将两图转回来：已知原来两个相机之间的旋转矩阵R和平移向量t，使得：
Xr＝RXl+t
其中，Xl表示红外相机的正规坐标，Xr表示可见光相机的正规坐标。将红外图像向R正方向旋转一半的角度，将可见光图像向R反方向旋转一半的角度；
1-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系。
所述步骤2)中对红外图像和可见光图像进行配准，具体包括以下步骤：
2-1)将红外图像和可见光图像按照同比例系数s缩小。这样做的目的是减小配准的计算量。
2-2)计算缩小后最佳位置：使用归一化互信息来描述红外图像和可见光图像重合区域的相似程度。
互信息是描述两个随机变量之间的相关熵,即任意变量包含另一个任意变量的信息多少。互信息的计算公式如下：
MI(A,B)＝H(A)+H(B)-H(A,B)
其中，H(A)和H(B)分别表示图像A和图像B的熵，H(A,B)表示图像A和图像B重叠区域的联合熵。熵(entropy)指的是体系的混乱的程度，它在控制论、概率论、数论、天体物理、生命科学等领域都有重要应用，在不同的学科中也有引申出的更为具体的定义。图像的熵表示图像包含的信息量，其计算公式如下
Figure BDA0002284263860000051
P(a)表示图像A上的灰度概率分布。a∈[0,255]
图像的联合熵的计算公式如下：
Figure BDA0002284263860000052
PAB(a,b)表示图像A和图像B重叠区域上的灰度联合概率分布，a∈[0,255]，b∈[0,255]。当图像A和图像B完全配准时，H(A,B)最小，即MI(A,B)最大。由此可见,互信息可以用来判定两幅图像间的相似性大小,当互信息取最大值时,相似性最大,即两幅图像实现配准。不过会随着重叠区域变化，联合熵也会变化，这时使用互信息会产生误配。为解决这个问题，可以使用归一化互信息来替代互信息：
Figure BDA0002284263860000053
记缩小的红外图像A0相对于缩小的可见光图像B0的位置是t0，A0(t0)和B0(t0)表示A0和B0在位置t0上的重叠区域，计算缩小后最佳位置t0 *：
Figure BDA0002284263860000054
2-3)计算原尺寸下的最佳位置：将缩小的红外图像A0和缩小的可见光图像B0缩放到原尺寸图像A和B，步骤2-2)计算出来的
Figure BDA0002284263860000055
缩放至
Figure BDA0002284263860000056
在t周围w×w范围内计算原尺寸下的最佳位置t*：
Figure BDA0002284263860000057
2-4)确定旋转角：对每一个侯选位置按照一个角度范围旋转多次(比如-10°～10°范围内分为200份，即从-10°位置开始每次转0.1°)，选取使NMI最大的对应位置和旋转角度。
所述步骤3)中提取特征点，具体包括以下步骤：
3-1)构建单尺度差分高斯金字塔(DoG)。差分高斯金字塔DoG是由相邻尺度空间做差得来的，常被用于尺度不变特征变换(Scale-invariant feature transform，SIFT)。一幅图像的尺度空间被定义为：高斯卷积核与该图像的卷积，它是高斯卷积核中的参数σ的函数。具体地，场景图像I(x,y的尺度空间为：
L(x,y,σ)＝G(x,y,σ)*I(x,y)
其中，
Figure BDA0002284263860000061
是高斯核函数，σ是尺度因子，σ的大小决定图像的平滑程度。大的σ值对应粗糙尺度(低分辨率)，小的σ值对应精细尺度(高分辨率)。*表示卷积操作。我们称L(x,y,σ)为图像I(x,y的尺度空间。
将不同尺度的尺度空间做差得到一层差分高斯金字塔(如图3所示)，此外还需乘一个归一化尺度因子λ，使得DoG图像的最大值是255。
D(x,y,σ)＝λ(L(x,y,kσ)-L(x,y,σ))
和SIFT不同的是，仅仅计算一种差分尺度特征。这么做的原因有两点：第一，计算多个尺度空间的计算量太大了，无法实现实时性；第二，使用多尺度空间得到的SIFT特征的准确度太低。
3-2)对于在上步得到的DoG中的每一个点，将其与邻域内的像素点作比较，判断其是否为局部极值点。若为局部极值点则为高斯特征点；
3-2-1)记上步得到的DoG为D；将D做膨胀操作，结果记为D1；将D1中每一个像素点与其8-邻域上的点作比较，如果改像素点是局部极大值，则将其加入候选点集P1里；
3-2-2)将D取反再做膨胀操作，结果记为D2。将D2中每一个像素点与其8-邻域上的点作比较，如果改像素点是局部极小值，则将其加入候选点集P2里。
3-2-3)将P1和P2取交集，得到P3＝P1∩P2；取P3中DoG灰度值大于15的点作为特征点集{P}；
3-3)由于仅仅根据高斯特征判断出的特征点里会出现噪点，所以我们需要对高斯特征点进行去噪。这里可以使用常用的滤波器来过滤噪点及边缘点。
所述步骤4)中特征点匹配，具体包括以下步骤：
4-1)将红外图像和可见光图像区域同时分为m×n个块。对于红外图每一个特征点
Figure BDA0002284263860000071
找到其在红外图对应的块
Figure BDA0002284263860000072
块
Figure BDA0002284263860000073
所对应的可见光图搜索范围记为
Figure BDA0002284263860000074
如图4所示。由于之前已经进行图像配准，
Figure BDA0002284263860000075
的范围就是块
Figure BDA0002284263860000076
在可见光图上的对应块。找到一个能够描述特征点相似程度的变量来评估
Figure BDA0002284263860000077
和
Figure BDA0002284263860000078
中任意一点的相似程度，如果相似程度最大值大于阈值t1，则视为粗匹配点
Figure BDA0002284263860000079
4-2)如果
Figure BDA00022842638600000710
和
Figure BDA00022842638600000711
中相似程度最大值sfirst和次大值ssecond满足：
F(sfirst,ssecond)≥t2
则保留该匹配，其中t2为阈值，F(sfirst，ssecond)用于描述sfirst和ssecond之间的关系。
按照该规则筛选之后，再按照步骤4-2)、4-3)的方法匹配
Figure BDA00022842638600000712
在红外图对应的特征点
Figure BDA00022842638600000713
如果满足
Figure BDA00022842638600000714
则保留该匹配
Figure BDA00022842638600000715
4-3)以红外图特征点
Figure BDA00022842638600000716
为基准，抛物线拟合优化对应可见光图的整数像素特征点
Figure BDA00022842638600000717
得到的对应可见光图的亚像素特征点
Figure BDA00022842638600000718
Figure BDA00022842638600000719
其中
Figure BDA00022842638600000720
为x方向上的亚像素偏移量，
Figure BDA00022842638600000721
为y方向上的亚像素偏移量。
4-4)以对应可见光图整数像素特征点
Figure BDA00022842638600000722
为基准，根据4-4)的方法计算出对应红外图的亚像素特征点
Figure BDA00022842638600000723
其中
Figure BDA00022842638600000724
为x方向上的亚像素偏移量，
Figure BDA00022842638600000725
为y方向上的亚像素偏移量。
4-5)最终的匹配点对为
Figure BDA00022842638600000726
所述步骤7)中修正标定结果，具体包括以下步骤：
7-1)使用随机抽样一致性(RANSAC)对点对做进一步筛选。
7-2)求解基础矩阵F和本质矩阵E：红外和可见光对应像素点对ul、ur和基础矩阵F的关系是：
Figure BDA00022842638600000727
可以将对应点坐标代入上式，构建齐次线性方程组求解F。
基础矩阵和本质矩阵的关系是：
Figure BDA0002284263860000081
其中，Kl、Kr分别是红外相机和可见光相机的内参矩阵。
7-3)从本质矩阵分解出旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：
E＝[t]×R
其中[t]×表示t的叉乘矩阵。
将E做奇异值分解，得
Figure BDA0002284263860000082
定义两个矩阵
Figure BDA0002284263860000083
和
Figure BDA0002284263860000084
ZW＝Σ
所以E可以写成以下两种形式
(1)E＝UZUTUWVT
令[t]×＝UZUT，R＝UWVT
(2)E＝-UZUTUWTVT
令[t]×＝-UzUT，R＝UWTVT
7-4)将分解出的旋转和平移关系叠加到原来的红外相机和可见光相机的位置关系里面；
记去畸变前的旋转矩阵为R0，平移向量为t0＝(tx,ty,tz)T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′x,t′y,t′z)T。则新的Rnew和tnew如下
Figure BDA0002284263860000085
Figure BDA0002284263860000086
此外，还需要将tnew乘一个系数，使得tnew在x方向上的分量
Figure BDA0002284263860000091
本发明的有益效果：本发明解决了由于温湿度、震动等因素造成红外相机和可见光相机位置关系的改变。具有速度快、结果精确、操作简单等优点。此外我们通过互信息来对红外图像和可见光图像进行配准。相比于普通的方法进一步缩小了特征点匹配范围。
附图说明
图1为整体流程图。
图2为校正流程图
图3表示高斯差分金字塔(DoG)。
图4是分块匹配的示意图。(a)左对应点所在的图像分块位置，(b)右对应点所在的图像分块位置。
具体实施方式
本发明旨在解决由于温湿度、震动等因素造成红外相机和可见光相机位置关系的改变。结合附图及实施例详细说明如下：
1)原图校正：将原图根据红外相机和可见光相机各自内参和原来的外参进行去畸变和双目校正。流程如图2所示。
1-1)计算图像的像素点对应的正规坐标系下的坐标。其中，正规坐标系是相机坐标系在平面Z＝1的投影；而相机坐标系是以相机的中心作为图像坐标系的原点，以图片方向为XY轴方向，以垂直于图像为Z轴方向的坐标系。像素坐标系以图片的左上角为原点，其x轴和y轴分别与图像坐标系的x轴和y轴平行。像素坐标系的单位是像素。像素坐标与正规坐标的关系如下：
u＝KX
Figure BDA0002284263860000092
其中，
Figure BDA0002284263860000101
表示图像的像素坐标；
Figure BDA0002284263860000102
表示相机的内参矩阵，fx和fy分别表示图像x方向和y方向的焦距，单位是像素，(cx,cy)表示相机的主点位置，即相机中心在图像上的对应位置；
Figure BDA0002284263860000103
是正规坐标系下的坐标。已知图像的像素坐标系以及相机的内参计算出像素点对应的正规坐标系，即X＝K-1u；
1-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)是实际带有畸变的正规坐标。
1-3)根据原来两相机的旋转关系将两图转回来：已知原来两个相机之间的旋转矩阵R和平移向量t，使得：
Xr＝RXl+t
其中，Xl表示红外相机的正规坐标，Xr表示可见光相机的正规坐标。将红外图像向R正方向旋转一半的角度，将可见光图像向R反方向旋转一半的角度；
1-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系。
2)根据红外图像与可见光图像的互信息来确定红外图像在可见光图像的最佳对应位置，即对红外图像和可见光图像进行配准。
2-1)将红外图像和可见光图像按照同比例系数s缩小。这样做的目的是减小配准的计算量。
2-2)计算缩小后最佳位置：使用归一化互信息来描述红外图像和可见光图像重合区域的相似程度。
互信息是描述两个随机变量之间的相关熵,即任意变量包含另一个任意变量的信息多少。互信息的计算公式如下：
MI(A,B)＝H(A)+H(B)-H(A,B)
其中，H(A)和H(B)分别表示图像A和图像B的熵，H(A,B)表示图像A和图像B重叠区域的联合熵。熵(entropy)指的是体系的混乱的程度，它在控制论、概率论、数论、天体物理、生命科学等领域都有重要应用，在不同的学科中也有引申出的更为具体的定义。图像的熵表示图像包含的信息量，其计算公式如下
Figure BDA0002284263860000111
P(a)表示图像A上的灰度概率分布。a∈[0,255]
图像的联合熵的计算公式如下：
Figure BDA0002284263860000112
PAB(a,b)表示图像A和图像B重叠区域上的灰度联合概率分布，a∈[0,255]，b∈[0,255]。当图像A和图像B完全配准时，H(A,B)最小，即MI(A,B)最大。由此可见,互信息可以用来判定两幅图像间的相似性大小,当互信息取最大值时,相似性最大,即两幅图像实现配准。不过会随着重叠区域变化，联合熵也会变化，这时使用互信息会产生误配。为解决这个问题，可以使用归一化互信息来替代互信息：
Figure BDA0002284263860000121
记缩小的红外图像A0相对于缩小的可见光图像B0的位置是t0，A0(t0)和B0(t0)表示A0和B0在位置t0上的重叠区域，计算缩小后最佳位置t0 *：
Figure BDA0002284263860000122
2-3)计算原尺寸下的最佳位置：将缩小的红外图像A0和缩小的可见光图像B0缩放到原尺寸图像A和B，步骤2-2)计算出来的
Figure BDA0002284263860000123
缩放至
Figure BDA0002284263860000124
在t周围w×w范围内计算原尺寸下的最佳位置t*：
Figure BDA0002284263860000125
2-4)确定旋转角：对每一个侯选位置按照一个角度范围旋转多次(比如-10°～10°范围内分为200份，即从-10°位置开始每次转0.1°)，选取使NMI最大的对应位置和旋转角度。
3)在配准后的红外图像和可见光图像上分别提取特征点。
3-1)构建单尺度差分高斯金字塔(DoG)。差分高斯金字塔差分高斯金字塔是由相邻尺度空间做差得来的，常被用于尺度不变特征变换(Scale-invariant featuretransform，SIFT)。一幅图像的尺度空间被定义为：高斯卷积核与该图像的卷积，它是高斯卷积核中的参数σ的函数。具体地，场景图像I(x,y)的尺度空间为：
L(x,y,σ)＝G(x,y,σ)*I(x,y)
其中，
Figure BDA0002284263860000126
是高斯核函数，σ是尺度因子，σ的大小决定图像的平滑程度。大的σ值对应粗糙尺度(低分辨率)，小的σ值对应精细尺度(高分辨率)。*表示卷积操作。我们称L(x,y,σ)为图像I(x,y)的尺度空间。
将不同尺度的尺度空间做差得到一层差分高斯金字塔(如图3所示)，此外还需乘一个归一化尺度因子λ，使得DoG图像的最大值是255。
D(x,y,σ)＝λ(L(x,y,kσ)-L(x,y,σ))
和SIFT不同的是，我们仅仅计算一种差分尺度特征。这么做的原因有两点：第一，计算多个尺度空间的计算量太大了，无法实现实时性；第二，使用多尺度空间得到的SIFT特征的准确度太低。
3-2)对于在上步得到的DoG中的每一个点，将其与邻域内的像素点作比较，判断其是否为局部极值点。
3-3)由于仅仅根据高斯特征判断出的特征点里会出现噪点，所以我们需要对高斯特征点进行去噪。这里可以使用常用的滤波器来过滤噪点及边缘点。
3-2-1)记上步得到的DoG为D。将D做膨胀操作，结果记为D1。将D1中每一个像素点与其8-邻域上的点作比较，如果改像素点是局部极大值，则将其加入候选点集P1里。
3-2-2)将D取反再做膨胀操作，结果记为D2。将D2中每一个像素点与其8-邻域上的点作比较，如果改像素点是局部极小值，则将其加入候选点集P2里。
3-2-3)将P1和P2取交集，得到P3＝P1∩P2。取P3中DoG灰度值大于15的点作为特征点集{P}。
4)对上一步提取的特征点进行匹配。
4-1)将红外图像和可见光图像区域同时分为m×n个块。对于红外图每一个特征点
Figure BDA0002284263860000131
找到其在红外图对应的块
Figure BDA0002284263860000132
块
Figure BDA0002284263860000133
所对应的可见光图搜索范围记为
Figure BDA0002284263860000134
如图4所示。由于之前已经进行图像配准，
Figure BDA0002284263860000135
的范围就是块
Figure BDA0002284263860000136
在可见光图上的对应块。找到一个能够描述特征点相似程度的变量来评估
Figure BDA0002284263860000137
和
Figure BDA0002284263860000138
中任意一点的相似程度，如果相似程度最大值大于阈值t1，则视为粗匹配点
Figure BDA0002284263860000139
4-2)如果
Figure BDA00022842638600001310
和
Figure BDA00022842638600001311
中相似程度最大值sfirst和次大值ssecond满足：
F(sfirst,ssecond)≥t2
则保留该匹配，其中t2为阈值，F(sfirst,ssecond)用于描述sfirst和ssecond之间的关系。
按照该规则筛选之后，再按照步骤4-2)、4-3)的方法匹配
Figure BDA0002284263860000141
在红外图对应的特征点
Figure BDA0002284263860000142
如果满足
Figure BDA0002284263860000143
则保留该匹配
Figure BDA0002284263860000144
4-3)以红外图特征点
Figure BDA0002284263860000145
为基准，抛物线拟合优化对应可见光图的整数像素特征点
Figure BDA0002284263860000146
得到的对应可见光图的亚像素特征点
Figure BDA0002284263860000147
Figure BDA0002284263860000148
其中
Figure BDA0002284263860000149
为x方向上的亚像素偏移量，
Figure BDA00022842638600001410
为y方向上的亚像素偏移量。
4-4)以对应可见光图整数像素特征点
Figure BDA00022842638600001411
为基准，根据4-4)的方法计算出对应红外图的亚像素特征点
Figure BDA00022842638600001412
其中
Figure BDA00022842638600001413
为x方向上的亚像素偏移量，
Figure BDA00022842638600001414
为y方向上的亚像素偏移量。
4-5)最终的匹配点对为
Figure BDA00022842638600001415
5)根据步骤2)和步骤3)的结果计算配准后的红外图像特征点对应于红外原图像的特征点。
6)判断特征点覆盖区域：将图像分成m*n个格子，如果特征点覆盖到所有格子，则进行下一步，否则继续拍摄图像，重复步骤1)～步骤5)。
7)修正标定结果：使用所有特征点的图像坐标来计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加。
7-1)使用随机抽样一致性(RANSAC)对点对做进一步筛选。
7-2)求解基础矩阵F和本质矩阵E：红外和可见光对应像素点对ul、ur和基础矩阵F的关系是：
Figure BDA00022842638600001416
可以将对应点坐标代入上式，构建齐次线性方程组求解F。
基础矩阵和本质矩阵的关系是：
Figure BDA0002284263860000151
其中，Kl、Kr分别是红外相机和可见光相机的内参矩阵。
7-3)从本质矩阵分解出旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：
E＝[t]×R
其中[t]×表示t的叉乘矩阵。
将E做奇异值分解，得
Figure BDA0002284263860000152
定义两个矩阵
Figure BDA0002284263860000153
和
Figure BDA0002284263860000154
ZW＝Σ
所以E可以写成以下两种形式
(1)E＝UZUTUWVT
令[t]×＝UZUT，R＝UWVT
(2)E＝-UZUTUWTVT
令[t]×＝-UZUT，R＝UWTVT
7-4)将分解出的旋转和平移关系叠加到原来的红外相机和可见光相机的位置关系里面；
记去畸变前的旋转矩阵为R0，平移向量为t0＝(tx,ty,tz)T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′x,t′y,t′z)T。则新的Rnew和tnew如下
Figure BDA0002284263860000155
Figure BDA0002284263860000156
此外，还需要将tnew乘一个系数，使得tnew在x方向上的分量
Figure BDA0002284263860000157

Claims (5)

1.基于互信息配准的可见光与红外相机联合标定方法，其特征在于，包括如下步骤：
1)原图校正：将原图根据红外相机和可见光相机各自内参和原来的外参进行去畸变和双目校正；
2)根据红外图像与可见光图像的互信息确定红外图像在可见光图像的最佳对应位置，即对红外图像和可见光图像进行配准；
3)在配准后的红外图像和可见光图像上分别提取特征点；
4)对上一步提取的特征点进行匹配；包括以下步骤：
4-1)将红外图像和可见光图像区域同时分为m×n个块；对于红外图每一个特征点
Figure FDA0003262804880000011
找到其在红外图对应的块
Figure FDA0003262804880000012
块
Figure FDA0003262804880000013
所对应的可见光图搜索范围记为
Figure FDA0003262804880000014
经图像配准后，
Figure FDA0003262804880000015
的范围就是块
Figure FDA0003262804880000016
在可见光图上的对应块；找到一个能够描述特征点相似程度的变量评估
Figure FDA0003262804880000017
和
Figure FDA0003262804880000018
中任意一点的相似程度，如果相似程度最大值大于阈值t1，则视为粗匹配点
Figure FDA0003262804880000019
4-2)如果
Figure FDA00032628048800000110
和
Figure FDA00032628048800000111
中相似程度最大值sfirst和次大值ssecond满足：
F(sfirst,ssecond)≥t2
则保留该匹配，其中t2为阈值，F(sfirst,ssecond)用于描述sfirst和ssecond之间的关系；
按照该规则筛选之后，再按照步骤4-1)、4-2)的方法匹配pi r在红外图对应的特征点
Figure FDA00032628048800000112
如果满足
Figure FDA00032628048800000113
则保留该匹配
Figure FDA00032628048800000114
4-3)以红外图特征点
Figure FDA00032628048800000115
为基准，抛物线拟合优化对应可见光图的整数像素特征点
Figure FDA00032628048800000116
得到的对应可见光图的亚像素特征点
Figure FDA00032628048800000117
Figure FDA00032628048800000118
其中
Figure FDA00032628048800000119
为x方向上的亚像素偏移量，
Figure FDA00032628048800000120
为y方向上的亚像素偏移量；
4-4)以对应可见光图整数像素特征点
Figure FDA00032628048800000121
为基准，根据4-3)的方法计算出对应红外图的亚像素特征点
Figure FDA00032628048800000122
其中
Figure FDA00032628048800000123
为x方向上的亚像素偏移量，
Figure FDA00032628048800000124
为y方向上的亚像素偏移量；
4-5)最终的匹配点对为
Figure FDA0003262804880000021
5)根据步骤2)和步骤3)的结果计算配准后的红外图像特征点对应于红外原图像的特征点；
6)判断特征点覆盖区域：将图像分成m*n个格子，如果特征点覆盖到所有格子，则进行下一步，否则继续拍摄图像，重复步骤1)～步骤5)；
7)修正标定结果：使用所有特征点的图像坐标计算校正之后的两相机之间的位置关系，然后与原来的外参相叠加。
2.根据权利要求1所述的基于互信息配准的可见光与红外相机联合标定方法，其特征在于，步骤1)原图校正，包括以下步骤：
1-1)计算图像的像素点对应的正规坐标系下的坐标；其中，正规坐标系是相机坐标系在平面Z＝1的投影；而相机坐标系是以相机的中心作为图像坐标系的原点，以图片方向为XY轴方向，以垂直于图像为Z轴方向的坐标系；像素坐标系以图片的左上角为原点，其x轴和y轴分别与图像坐标系的x轴和y轴平行；像素坐标系的单位是像素；像素坐标与正规坐标的关系如下：
u＝KX
Figure FDA0003262804880000022
其中，
Figure FDA0003262804880000023
表示图像的像素坐标；
Figure FDA0003262804880000024
表示相机的内参矩阵，fx和fy分别表示图像x方向和y方向的焦距，单位是像素，(cx,cy)表示相机的主点位置，即相机中心在图像上的对应位置；
Figure FDA0003262804880000025
是正规坐标系下的坐标；已知图像的像素坐标系以及相机的内参计算出像素点对应的正规坐标系，即X＝K-1u；
1-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)是实际带有畸变的正规坐标；
1-3)根据原来两相机的旋转关系将两图转回来：已知原来两个相机之间的旋转矩阵R和平移向量t，使得：
Xr＝RXl+t
其中，Xl表示红外相机的正规坐标，Xr表示可见光相机的正规坐标；将红外图像向R正方向旋转一半的角度，将可见光图像向R反方向旋转一半的角度；
1-4)根据公式u＝KX将去畸旋转后的图像还原至像素坐标系。
3.根据权利要求1所述的基于互信息配准的可见光与红外相机联合标定方法，其特征在于，步骤2)中对红外图像和可见光图像进行配准，包括以下步骤：
2-1)将红外图像和可见光图像按照同比例系数s缩小；
2-2)计算缩小后最佳位置：使用归一化互信息描述红外图像和可见光图像重合区域的相似程度；
归一化互信息：
Figure FDA0003262804880000041
其中，互信息是描述两个随机变量之间的相关熵,即任意变量包含另一个任意变量的信息多少；H(A)和H(B)分别表示图像A和图像B的熵，H(A,B)表示图像A和图像B重叠区域的联合熵；图像的熵表示图像包含的信息量，其计算公式如下
Figure FDA0003262804880000042
P(a)表示图像A上的灰度概率分布；a∈[0,255]
图像的联合熵的计算公式如下：
Figure FDA0003262804880000043
PAB(a,b)表示图像A和图像B重叠区域上的灰度联合概率分布，a∈[0,255]，b∈[0,255]；
记缩小的红外图像A0相对于缩小的可见光图像B0的位置是t0，A0(t0)和B0(t0)表示A0和B0在位置t0上的重叠区域，计算缩小后最佳位置t0 *：
Figure FDA0003262804880000044
2-3)计算原尺寸下的最佳位置：将缩小的红外图像A0和缩小的可见光图像B0缩放到原尺寸图像A和B，步骤2-2)计算出来的
Figure FDA0003262804880000045
缩放至
Figure FDA0003262804880000046
在t周围w×w范围内计算原尺寸下的最佳位置t*：
Figure FDA0003262804880000047
2-4)确定旋转角：对每一个侯选位置按照一个角度范围旋转N次，N是大于零的整数，选取使NMI最大的对应位置和旋转角度。
4.根据权利要求1所述的基于互信息配准的可见光与红外相机联合标定方法，其特征在于，步骤3)中提取特征点，包括以下步骤：
3-1)构建单尺度差分高斯金字塔DoG；差分高斯金字塔DoG是由相邻尺度空间做差得来的，场景图像I(x,y)的尺度空间为：
L(x,y,σ)＝G(x,y,σ)*I(x,y)
其中，
Figure FDA0003262804880000051
是高斯核函数，σ是尺度因子，σ的大小决定图像的平滑程度；*表示卷积操作；
计算一层差分高斯金字塔：
D(x,y,σ)＝λ(L(x,y,kσ)-L(x,y,σ))
其中，λ是归一化因子，使得DoG图像的最大值是255；
3-2)对于在3-1)中得到的DoG中的每一个点，将其与邻域内的像素点作比较，判断其是否为局部极值点；若为局部极值点则为高斯特征点；
3-2-1)记上步得到的DoG为D；将D做膨胀操作，结果记为D1；将D1中每一个像素点与其8-邻域上的点作比较，如果改像素点是局部极大值，则将其加入候选点集P1里；
3-2-2)将D取反再做膨胀操作，结果记为D2；将D2中每一个像素点与其8-邻域上的点作比较，如果改像素点是局部极小值，则将其加入候选点集P2里；
3-2-3)将P1和P2取交集，得到P3＝P1∩P2；取P3中DoG灰度值大于15的点作为特征点集{P}；
3-3)对高斯特征点去噪；用滤波器过滤噪点及边缘点。
5.根据权利要求1所述的基于互信息配准的可见光与红外相机联合标定方法，其特征在于，步骤7)中修正标定结果，包括以下步骤：
7-1)使用随机抽样一致性方法对点对做进一步筛选；
7-2)求解基础矩阵F和本质矩阵E：红外和可见光对应像素点对ul、ur和基础矩阵F的关系是：
Figure FDA0003262804880000052
将对应点坐标代入上式，构建齐次线性方程组求解F；
基础矩阵和本质矩阵的关系是：
Figure FDA0003262804880000053
其中，Kl、Kr分别是红外相机和可见光相机的内参矩阵；
7-3)从本质矩阵分解出旋转和平移关系：本质矩阵E与旋转R和平移t的关系如下：
E＝[t]×R
其中[t]×表示t的叉乘矩阵；
将E做奇异值分解，得
Figure FDA0003262804880000061
定义两个矩阵
Figure FDA0003262804880000062
和
Figure FDA0003262804880000063
ZW＝Σ
所以E写成以下两种形式
(1)E＝UZUTUWVT
令[t]×＝UZUT，R＝UWVT
(2)E＝-UZUTUWTVT
令[t]×＝-UZUT，R＝UWTVT
7-4)将分解出的旋转和平移关系叠加到原来的红外相机和可见光相机的位置关系里面；
记去畸变前的旋转矩阵为R0，平移向量为t0＝(tx,ty,tz)T；上一步计算出的旋转矩阵为R，平移向量为t＝(t′x,t′y,t′z)T；则新的Rnew和tnew如下：
Figure FDA0003262804880000064
Figure FDA0003262804880000065
此外，要将tnew乘一个系数，使得tnew在x方向上的分量
Figure FDA0003262804880000066
CN201911153787.4A 2019-11-22 2019-11-22 基于互信息配准的可见光与红外相机联合标定方法 Active CN110969669B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911153787.4A CN110969669B (zh) 2019-11-22 2019-11-22 基于互信息配准的可见光与红外相机联合标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911153787.4A CN110969669B (zh) 2019-11-22 2019-11-22 基于互信息配准的可见光与红外相机联合标定方法

Publications (2)

Publication Number Publication Date
CN110969669A CN110969669A (zh) 2020-04-07
CN110969669B true CN110969669B (zh) 2021-12-03

Family

ID=70031243

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911153787.4A Active CN110969669B (zh) 2019-11-22 2019-11-22 基于互信息配准的可见光与红外相机联合标定方法

Country Status (1)

Country Link
CN (1) CN110969669B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113537204A (zh) * 2020-04-20 2021-10-22 富华科精密工业（深圳）有限公司 基于红外特征、机器学习的小火焰检测方法及计算机装置
CN111627072B (zh) * 2020-04-30 2023-10-24 贝壳技术有限公司 一种对多传感器进行标定的方法、装置和存储介质
CN111915683A (zh) * 2020-07-27 2020-11-10 湖南大学 图像位置标定方法、智能设备及存储介质
CN112634374B (zh) * 2020-12-18 2023-07-14 杭州海康威视数字技术股份有限公司 双目相机的立体标定方法、装置、系统及双目相机
CN113744349A (zh) * 2021-08-31 2021-12-03 湖南航天远望科技有限公司 一种红外光谱图像测量对准方法、装置及介质
CN116704048B (zh) * 2023-08-09 2023-11-17 四川元祉智慧科技有限公司 一种双光配准方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108828606A (zh) * 2018-03-22 2018-11-16 中国科学院西安光学精密机械研究所 一种基于激光雷达和双目可见光相机联合测量方法
WO2019081305A1 (en) * 2017-10-24 2019-05-02 Schreder S.A. METHOD AND SYSTEM FOR CONTROLLING A LUMINAIRE, AND LUMINAIRE COMPRISING SUCH A CONTROL SYSTEM
CN110349221A (zh) * 2019-07-16 2019-10-18 北京航空航天大学 一种三维激光雷达与双目可见光传感器的融合标定方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7961134B2 (en) * 2009-03-18 2011-06-14 The United States Of America As Represented By The Secretary Of The Army Metric and self-calibration for an automatic, surveillance-based change detection system operating on noisy imagery

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019081305A1 (en) * 2017-10-24 2019-05-02 Schreder S.A. METHOD AND SYSTEM FOR CONTROLLING A LUMINAIRE, AND LUMINAIRE COMPRISING SUCH A CONTROL SYSTEM
CN108828606A (zh) * 2018-03-22 2018-11-16 中国科学院西安光学精密机械研究所 一种基于激光雷达和双目可见光相机联合测量方法
CN110349221A (zh) * 2019-07-16 2019-10-18 北京航空航天大学 一种三维激光雷达与双目可见光传感器的融合标定方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Automatic calibration and registration of lidar and stereo camera without calibration objects;Vijay John 等;《2015 IEEE International Conference on Vehicular Electronics and Safety (ICVES)》;20160218;全文 *
大空间建筑火源的精确定位方法;毛颖 等;《计算机应用与软件》;20160229;第33卷(第2期);全文 *

Also Published As

Publication number Publication date
CN110969669A (zh) 2020-04-07

Similar Documents

Publication Publication Date Title
CN110969669B (zh) 基于互信息配准的可见光与红外相机联合标定方法
CN110969668B (zh) 一种长焦双目相机的立体标定算法
CN110969670B (zh) 基于显著特征的多光谱相机动态立体标定方法
US11398053B2 (en) Multispectral camera external parameter self-calibration algorithm based on edge features
CN110992409B (zh) 基于傅里叶变换配准的多光谱立体相机动态配准方法
CN110135455B (zh) 影像匹配方法、装置及计算机可读存储介质
CN111080709B (zh) 基于轨迹特征配准的多光谱立体相机自标定算法
CN110956661B (zh) 基于双向单应矩阵的可见光与红外相机动态位姿计算方法
CN111080529A (zh) 一种加强鲁棒性的无人机航拍图像拼接方法
CN110880191B (zh) 基于直方图均衡化的红外立体相机动态外参计算方法
CN110910456B (zh) 基于Harris角点互信息匹配的立体相机动态标定方法
CN109118544B (zh) 基于透视变换的合成孔径成像方法
CN111369605B (zh) 一种基于边缘特征的红外与可见光图像的配准方法和系统
CN112016478B (zh) 一种基于多光谱图像融合的复杂场景识别方法及系统
CN111860651B (zh) 一种基于单目视觉的移动机器人半稠密地图构建方法
CN112254656A (zh) 一种基于结构表面点特征的立体视觉三维位移测量方法
CN110136048B (zh) 一种图像配准方法及系统、存储介质及终端
CN116958419A (zh) 一种基于波前编码的双目立体视觉三维重建系统和方法
CN113793266A (zh) 一种多目机器视觉图像拼接方法、系统及存储介质
CN111127353A (zh) 一种基于块配准和匹配的高动态图像去鬼影方法
CN110910457B (zh) 基于角点特征的多光谱立体相机外参计算方法
CN108447084B (zh) 基于orb特征的立体匹配补偿方法
CN111833384B (zh) 一种可见光和红外图像快速配准方法及装置
Flores et al. Generating a full spherical view by modeling the relation between two fisheye images
Mingquan et al. Optimized detection method of automotive floor solder joints based on digital-analog information

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