CN112927133B - 一种基于一体化标定参数的图像空间投影拼接方法 - Google Patents

一种基于一体化标定参数的图像空间投影拼接方法 Download PDF

Info

Publication number
CN112927133B
CN112927133B CN202110169365.7A CN202110169365A CN112927133B CN 112927133 B CN112927133 B CN 112927133B CN 202110169365 A CN202110169365 A CN 202110169365A CN 112927133 B CN112927133 B CN 112927133B
Authority
CN
China
Prior art keywords
axis
coordinate system
rotation
point
image
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
CN202110169365.7A
Other languages
English (en)
Other versions
CN112927133A (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.)
Hunan Qiaokang Intelligent Technology Co ltd
Original Assignee
Hunan Qiaokang Intelligent Technology Co ltd
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 Hunan Qiaokang Intelligent Technology Co ltd filed Critical Hunan Qiaokang Intelligent Technology Co ltd
Priority to CN202110169365.7A priority Critical patent/CN112927133B/zh
Publication of CN112927133A publication Critical patent/CN112927133A/zh
Application granted granted Critical
Publication of CN112927133B publication Critical patent/CN112927133B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/40Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4038Image mosaicing, e.g. composing plane images from plane sub-images
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06KGRAPHICAL DATA READING; PRESENTATION OF DATA; RECORD CARRIERS; HANDLING RECORD CARRIERS
    • G06K7/00Methods or arrangements for sensing record carriers, e.g. for reading patterns
    • G06K7/10Methods or arrangements for sensing record carriers, e.g. for reading patterns by electromagnetic radiation, e.g. optical sensing; by corpuscular radiation
    • G06K7/14Methods or arrangements for sensing record carriers, e.g. for reading patterns by electromagnetic radiation, e.g. optical sensing; by corpuscular radiation using light without selection of wavelength, e.g. sensing reflected white light
    • G06K7/1404Methods for optical code recognition
    • G06K7/1408Methods for optical code recognition the method being specifically adapted for the type of code
    • G06K7/14172D bar codes
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06KGRAPHICAL DATA READING; PRESENTATION OF DATA; RECORD CARRIERS; HANDLING RECORD CARRIERS
    • G06K7/00Methods or arrangements for sensing record carriers, e.g. for reading patterns
    • G06K7/10Methods or arrangements for sensing record carriers, e.g. for reading patterns by electromagnetic radiation, e.g. optical sensing; by corpuscular radiation
    • G06K7/14Methods or arrangements for sensing record carriers, e.g. for reading patterns by electromagnetic radiation, e.g. optical sensing; by corpuscular radiation using light without selection of wavelength, e.g. sensing reflected white light
    • G06K7/1404Methods for optical code recognition
    • G06K7/1439Methods for optical code recognition including a method step for retrieval of the optical code
    • G06K7/1452Methods for optical code recognition including a method step for retrieval of the optical code detecting bar code edges
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/60Rotation of whole images or parts thereof
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/80Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
    • G06T7/85Stereo camera calibration
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00Indexing scheme for image data processing or generation, in general
    • G06T2200/32Indexing scheme for image data processing or generation, in general involving image mosaicing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10028Range image; Depth image; 3D point clouds

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Health & Medical Sciences (AREA)
  • Electromagnetism (AREA)
  • General Health & Medical Sciences (AREA)
  • Toxicology (AREA)
  • Artificial Intelligence (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

本发明公开了一种基于一体化标定参数的图像空间投影拼接方法,通过三维标定板,标定相机内参数;根据2D‑LiDAR单次扫描标定板板的点云数据确定2D‑LiDAR和相机之间的旋转矩阵和平移矩阵,根据旋转后相机图像中关键点坐标与采集末端零位时图像关键点坐标对比,结合相机与2D‑LiDAR的外参数,实现旋转云台的旋转参数的标定;使得标定结果处于旋转云台一体化坐标系中,完成旋转图像的空间投影及拼接。本发明提高了标定的效率和准确性,解决了传统的多传感器与运动机构需分步标定的问题;将所有标定结果,包括相机图像与激光点云,在基于旋转云台一体化坐标系内进行表示,并将其运动学方程用于旋转图像的空间投影及拼接方法。

Description

一种基于一体化标定参数的图像空间投影拼接方法
技术领域
本发明属于图像处理领域,尤其涉及一种基于一体化标定参数的图像空间投影拼接方法。
背景技术
获取室外真实场景的三维模型是计算机视觉、视觉导航、机器人等领域的一个重要研究课题,具有广阔的应用场景:如机器人导航,无人驾驶,移动道路测量车和城市或地图3D数字化等。对于大型场景的空间测量,单一传感器存在不可避免的局限性,2D-LiDAR与相机联合测量技术结合了两种类型传感器的优点,能同时获取纹理丰富的图像和几何结构信息。
在利用2D-LiDAR与相机联合测量技术进行空间测量前,需要对其几何安装参数进行标定。
目前,大多数相关工作主要是利用特定的标定板来计算两传感器间精准的相对位姿(旋转矩阵和平移矢量),该标定板一般是三角板或由三角板组成的V 型板等等,通过特定的结构提取激光点云和相机图像中的角点,从而完成标定。但是,激光点云存在稀疏性,会影响点云数据和图像数据中角点的准确性,从而影响标定的精确性。
此外,对建筑物,如桥梁、外墙面等进行外观检测时,由于摄像头经常需要调整角度和转向,导致拍摄得到的图片会出现倾斜、扭曲等,导致需要将图像拼接形成建筑的整体图形时,需要根据各图像边缘的相似度寻找到邻近的图片,然后将图片边缘、扭曲、拉伸或缩减等将邻近图片拼接在一起。其缺点在于,处理流程繁琐,耗费时间长,且当建筑有大量相似的部分时,很容易使得图片拼接错误,而且缺少部分图片时,即无法经图片拼接。
名词解释:
2D-LiDAR:二维激光雷达。
ChArUco标定板:棋盘格与ArUco二维码相结合的标定板。
发明内容
为解决上述问题,本发明公开了一种基于一体化标定参数的图像空间投影拼接方法,本发明采用特殊的U型三维标定板,仅需采集一次标定数据,即实现多传感器几何参数的自动标定,提出采集末端的旋转云台一体化坐标系,确定了旋转云台的运动学方程,形成了图像与2D激光点云的有效融合,以实现多自由度旋转图像的快速精确拼接。
为实现上述目的,本发明的技术方案为:
一种基于一体化标定参数的图像空间投影拼接方法,包括如下步骤:
步骤一、建立一体化标定***;
所述一体化标定***包括标定板和多自由度采集末端,多自由度采集末端包括旋转云台,旋转云台包括2D-LiDAR和至少两个相机;旋转云台固定连接在Y 轴旋转轴,Y轴旋转轴轴接连接有Z轴旋转轴;所述标定板包括ChArUco标定板, ChArUco标定板一侧固定有矩形侧板,另一侧固定有等腰直角三角形侧板,等腰直角三角形侧板的一个直角边为ChArUco标定板的侧边;矩形侧板外侧固定有第一边侧板;等腰直角三角形侧板的斜面上固定有第二边侧板;矩形侧板和等腰直角三角形侧板所在平面均与ChArUco标定板所在平面垂直设置;第一边侧板所在平面与ChArUco标定板坐在平面平行,第二边侧板所在平面与等腰直角三角形侧板所在平面垂直设置;
步骤二、获取2D-LiDAR单次扫描标定板的点云数据,提取标定板的角点信息,根据角点信息确定多自由度采集末端在各个旋转自由度时的零点位置;
步骤三、驱动Y轴旋转轴和Z轴旋转轴预设轨迹旋转,并采集相机和 2D-LiDAR数据形成完整的标定数据集;
步骤四、
提取完整的标定数据集中每张图像中的ArUco特征点和ChArUco角点坐标序列,所述ArUco特征点包括二维码的ID号和二维码四个顶点的平面坐标;根据采集末端不同旋转角度时相机图像中相应二维码四个顶点的平面坐标信息,标定相机内参;
步骤五、在获取2D-LiDAR扫描标定板的点云数据,提取标定板角点,计算各标定板角点在相应相机图像中的平面坐标信息,结合各个标定板角点的空间和平面坐标信息,确定2D-LiDAR和相机之间的旋转矩阵和平移矩阵;
步骤六:获取2D-LiDAR绕多自由度采集末端的Z轴旋转轴和Y轴旋转轴旋转得到的点云数据和相机图像,进行多自由度采集末端的旋转云台的Y轴和Z 轴旋转参数的联合标定;
步骤七:将Z轴旋转轴和Y轴旋转轴标定结果在2D-LiDAR坐标系内表示;
步骤八:多自由度采集末端采集目标物的图像,并将采集到图像的像素点均转化到2D-LiDAR坐标系内,完成目标物图像的拼接。
进一步的改进,所述步骤二包括如下步骤:
相机朝向ChArUco标定板,2D-LiDAR朝向与矩形侧板所在平面平行的方向,采集2D-LiDAR扫描标定板的激光点云;对激光点云进行稀疏激光均匀化采样,激光滤波,得到均匀连续的激光点云;将均匀连续的激光点云转换成图像,对图像进行形态学标定板角点取激光点云清晰连续的骨架,获得获取图像中激光点云的轮廓,当轮廓中包含的点的个数大于或等于预设值时,则为有效有效激光轮廓,否则为无效激光轮廓,重新采集2D-LiDAR扫描标定板的激光点云;
根据有效激光轮廓,进行多边形拟合后,提取有效激光轮廓中关于标定板的多边形轮廓线段,通过外扩矩阵确定标定板角点位置:
Figure GDA0003541226060000031
式1中,K1与K2分别表示多边形轮廓上任一线段的起始点和终止点的坐标;
K3和K4分别表示线段K1K2经过外扩后对应的起始点和终止点的坐标;K3和K4作为更新后的标定板角点位置;
然后根据给定的阈值s,计算外扩矩阵四个顶点的坐标,得到获取包含4个顶点的最小外接矩形;s表示此外接矩形高度的二分之一。(最小外接矩阵的二分;
通过给定的阈值s,计算中矩形的4个顶点在2D-LiDAR坐标系中的坐标值,并获取包含4个顶点的最小外接矩形,以判断激光点云在哪个矩形内,将所有激光点云分配到对应的矩形内,然后将矩形内的激光点拟合得到直线,根据两两直线的交点得到有效激光轮廓的角点,分别为p0、p1、p2、p3、p4、p5,其中p1、 p2、p3、p4即标定板角点;
其中,d=3δ,s=2δ,δ表示2D-LiDAR的测量误差;
驱动多自由度采集末端的Y轴旋转轴旋转,采集序列激光剖面,并提取处于ChArUco标定板上p2p3连线与2D-LiDAR坐标系中X轴方向的夹角,取夹角最小的位置为Y轴的零点,记录Y轴为零点时对应的Y轴机械角度θ,并驱动Y轴旋转至Y轴的零点;然后驱动Z轴旋转轴旋转,采集序列激光剖面,并提取线段 p2p3的长度,取p2p3长度最小时的位置为Z轴的零点,Z轴的零点对应的Z轴机械角度β,将{θ,β}角度作为多自由度采集末端的零点位置,即完成多自由度采集末端的零点位置标定。
进一步的改进,所述步骤三包括如下步骤:
将多自由度采集末端的Y轴旋转轴与Z轴旋转轴均置于零点位置,然后分两步进行标定采集相机和2D-LiDAR数据形成完整的标定数据集:
保持Z轴旋转轴处于Z轴的零点位置不动,驱动多自由度采集末端的Y轴旋转轴转动,每次转动1±0.5°;每次Y轴旋转轴转动后,多自由度采集末端上的两个相机拍摄图像并保存,记录2D-LiDAR一帧点云数据,同时记录多自由度采集末端在Y轴自由度和Z轴自由度的旋转角度,作为一组标定数据;在每个旋转位置保存对应的一组标定数据,作为Y轴转动的标定数据集;
然后保持Y轴旋转轴处于Y轴的零点位置不动,驱动多自由度采集末端的Z 轴旋转轴转动,每次转动1±0.5°,将转动过程中每组相机图像、激光点云和多自由度采集末端的姿态数据保存下来,作为Z轴转动的标定数据集;
Y轴转动的标定数据和Z轴转动的标定数据集形成完整的标定数据集。
进一步的改进,所述步骤四包括如下步骤:
提取Y轴转动的标定数据中每张图像中的ArUco特征点,所述ArUco特征点包括二维码的ID号和二维码四个顶点的平面坐标;同时提取出每张图像的 ChArUco角点坐标序列;基于Y轴旋转轴旋转的多个位置,图像坐标系上的角点序列和相对应的世界坐标序列,根据棋盘格标定方法标定出相机的内参数矩阵K 和畸变参数矩阵D,完成相机内参的标定:
Figure GDA0003541226060000041
D=[k1 k2 k3 m1 m2] (式2)
式2中,fx和fy分别表示x与y方向上像素与实际距离的比例;u0,v0分别为图像中心在像素平面即图像坐标系的u、v轴坐标,u为图像宽度方向为方向,v 为图像高度方向;k1,k2,k3表示镜头的径向畸变参数;m1,m2表示镜头的切向畸变参数。
进一步的改进,所述步骤五包括如下步骤:
根据Y轴转动的标定数据集,获取激光点云数据,提取有效激光轮廓中的标定板角点,即得到标定板角点p1、p2、p3和p4在2D-LiDAR坐标系中的空间坐标{Pi L|i=1,2,3,4};
根据得到的激光轮廓的角点坐标,通过标定板结构的几何约束,设ChArUco 标定板处于标定板世界坐标系Z=0平面,于是通过{Pi L|i=1,2,3,4}以及标定板的几何约束得到所述标定板角点在标定板世界坐标系的坐标{Pi W|i=1,2,3,4},即有:
Figure GDA0003541226060000051
Figure GDA0003541226060000052
Figure GDA0003541226060000053
Figure GDA0003541226060000054
式3中,d为二维码边缘与ChArUco标定板边缘的距离,w为ChArUco标定板的边长;
以相机图像中提取的二维码四个顶点作为控制点,结合标定的相机内参数及畸变参数,通过相机标定方法标定出相机外参,即相机坐标系与标定板世界坐标系之间的旋转矩阵RW2C与平移矩阵TW2C,然后确定标定板角点在图像坐标系中的实际坐标:
Figure GDA0003541226060000055
式4中,WXi,WYi,WZ分别为标定板角点在标定板世界坐标系中的X、Y、 Z分量;
Figure GDA0003541226060000062
分别表示标定板角点在相机坐标系中的X、Y、Z分量;
Figure GDA0003541226060000063
分别表示标定板角点在图像中的u轴和v轴坐标值;
得到标定板角点在图像坐标系中的理论坐标信息:
Figure GDA0003541226060000064
式5中,RL2C和TL2C分别表示2D-LiDAR与相机间的旋转矩阵和平移矩阵;LXi,LYi,LZi分别为标定板角点的X、Y、Z分量;CXi,CYi,CZi分别表示标定板角点在相机坐标系中对应的理论坐标值;Iui,Ivi分别表示标定板角点在图像中的理论像素坐标值;
至此,定义损失函数为:
Figure GDA0003541226060000065
式6中,Δu和Δv分别表示像素点在u和v方向的偏差;
采用最小二乘法,迭代求解得到2D-LiDAR与相机间旋转矩阵和平移矩阵。
进一步的改进,所述步骤六包括如下步骤:
根据完整的标定数据集,进行旋转云台绕Z轴和Y轴的旋转参数标定,标定结果包括代表旋转轴方向的单位向量和绕轴旋转点,旋转云台绕Z轴和Y轴的标定结果最后都统一到同一个坐标系2D-LiDAR坐标系下:
首先,设绕Y轴旋转轴旋转的旋转向量就是2D-LiDAR坐标系中Y轴(0,1,0),通过由粗到细的搜索方法求解绕得到Y轴旋转轴的旋转点,同样通过由粗到细的搜索方法得到Z轴旋转轴的旋转点;利用相机图像中的ArUco特征点作为外部控制点,将旋转时相机图像的ArUco特征点,与旋转到零点角度时相机图像的ArUco 特征点之间的偏差作为误差函数,利用最小二乘法,迭代求解误差函数最小时的旋转轴方向与旋转点;求解得到旋转云台的Y轴、Z轴的旋转参数,作为迭代求解时Y、Z轴的旋转轴向量和旋转点的初始值。
进一步的改进,所述由粗到细的搜索方法包括如下步骤:
一)获取Z轴旋转轴旋转不同角度的一系列激光点云数据,读取激光的单次点云数据,记录Z轴旋转轴第i次旋转的角度θi;对第i次转台绕Z轴旋转后采集到的原始点云数据,进行稀疏激光均匀化上采样和滤波处理,得到有效的点云数据,提取得到标定板角点p2、p3;最后得到标定板角点p2、p3与Z轴旋转轴旋转角度的序列{LSPi ji}(j=2,3),其中,i表示第i次旋转,j表示p2或p3点,LSPi j则表示第i次旋转时p2或p3点在2D-LiDAR坐标系中的坐标值;
二)将2D-LiDAR绕多自由度的旋转简化为空间上一点绕两向量旋转,存在如下关系:
Figure GDA0003541226060000071
Figure GDA0003541226060000072
where:
Figure GDA0003541226060000073
Sz=(XAZ,YAZ,ZAZ)T
r0Lnx·Lny·(1-cosθ)+cosθ,r1Lnx·Lny·(1-cosθ)+Lnz·sinθ
r2Lnz·Lnx·(1-cosθ)-Lny·sinθ,r3Lnx·Lny·(1-cosθ)-Lnz·sinθ
r4=(Lny)2·(1-cosθ)+cosθ,r5Lny·Lnz·(1-cosθ)+Lnx·sinθ
r6Lnz·Lnx·(1-cosθ)+Lny·sinθ,r7Lny·Lnz·(1-cosθ)-Lnx·sinθ
r8=(Lnz)2·(1-cosθ)+cosθ
(式7)
式7中,θ为Z轴旋转轴旋转的角度,
Figure GDA0003541226060000074
表示Z轴旋转轴旋转的方向向量, Sz表示Z轴旋转轴的旋转点,
Figure GDA0003541226060000075
则表示旋转云台绕
Figure GDA0003541226060000076
方向旋转θ角度对应的旋转矩阵,P1表示2D-LiDAR坐标系内任意一点,P1
Figure GDA0003541226060000077
轴旋转θ角度对应点为P0;r0、r1、r2...r8分别表示旋转矩阵
Figure GDA0003541226060000078
的各元素XAZ,YAZ,ZAZ,分别表示绕Z轴旋转点SZ在2D-LiDAR坐标系的X、Y、Z值,LnxLnyLnz分别表示旋转轴向量
Figure GDA0003541226060000079
在2D-LiDAR坐标系的X、Y、Z轴分量;
确定点P0绕Y轴旋转轴
Figure GDA00035412260600000710
轴旋转α后的几何关系:
Figure GDA00035412260600000711
式8中,α表示旋转云台绕Y轴旋转的角度,
Figure GDA0003541226060000082
表示绕Y轴旋转的方向向量, SY则表示绕Y轴旋转点,
Figure GDA0003541226060000083
表示旋转云台绕
Figure GDA0003541226060000084
方向旋转α角度对应的旋转矩阵;P2则表示P0在绕
Figure GDA0003541226060000085
旋转α角度后的对应点;
将绕Z轴旋转轴旋转后的标定板角点p2、p3还原到旋转零角度的公式即为:
Figure GDA0003541226060000086
公式9中,θ0表示在Z轴旋转轴旋转标定为零角度时,传感器记录的实际角度,而θi则表示在第i次Z轴旋转轴旋转时所记录的角度值;
LSPi j(X),LSPi j(Y),LSPi j(Z)分别表示Z轴旋转轴第i次旋转后,标定板第j 个角点在2D-LiDAR坐标系中的坐标值,
Figure GDA0003541226060000087
分别表示标定板角点在还原至绕Z轴旋转零角度后,在2D-LiDAR坐标系中的坐标值;
将所有采集末端绕Z轴旋转采集的激光点云数据均还原到绕Z轴旋转零角度,即形成点云数据;
三)设旋Z轴旋转轴向量为
Figure GDA0003541226060000088
计算Z轴旋转轴的旋转点,根据由粗到细的搜索策略,先利用循环迭代的方法进行粗略的求解:
A.初始时,设旋Z轴旋转轴的旋转点为SZ(0,0,0);
B.Z轴旋转轴旋转采集到的所有激光点云数据,均还原到Z轴旋转轴旋转角度为零时的状态;
C.根据还原后的点云数据,计算得到包含还原后的点云数据的最小外接矩阵,记最小外接矩阵的四个顶点为{Lqk}(k=1,2,3,4);
D.将四个顶点连线,形成四条线段{lk}(k=1,2,3,4),随后给定阈值d和s,得到lk线段各自的外扩矩形;
E.依次判断激光点云落在哪个外扩矩形内,并形成每个外扩矩形的激光点集;根据每个矩形内点集拟合直线,并计算矩形内所有点到对应直线的距离累加值;
F.计算四个外扩矩形的顶点到拟合直线距离累加值之和,如果小于上一次迭代的累加值之和,则根据阈值Δy将旋转点SZ的Y分量加上Δy;Δy=1-3mm;
G.若迭代次数超过设定阈值次数则结束,否则回到步骤B运行循环;
迭代过程停止之后,得到粗略求解的采集末端绕Z轴旋转点,记为
Figure GDA0003541226060000091
其中
Figure GDA0003541226060000092
表示
Figure GDA0003541226060000093
点在2D-LiDAR坐标系中的X 分量,
Figure GDA0003541226060000094
表示
Figure GDA0003541226060000095
点在2D-LiDAR坐标系中的Y分量,
Figure GDA0003541226060000096
表示
Figure GDA0003541226060000097
点在2D-LiDAR坐标系中的Z分量;
Figure GDA0003541226060000098
作为细化搜索的迭代初值,设定一个绕Z轴旋转轴旋转点Y分量的迭代阈值dy,dy=0.1*Δy,同时将细化搜索时旋转点Y值搜索范围控制在
Figure GDA0003541226060000099
之内;于是按照步骤A至步骤F的迭代过程,细化搜索得到旋转云台绕Z轴旋转轴的旋转点,最终结果得到的旋转点记为
Figure GDA00035412260600000910
其中
Figure GDA00035412260600000911
表示
Figure GDA00035412260600000912
点在2D-LiDAR坐标系中的X分量,
Figure GDA00035412260600000913
表示
Figure GDA00035412260600000914
点在2D-LiDAR坐标系中的Y分量,
Figure GDA00035412260600000915
表示
Figure GDA00035412260600000916
点在2D-LiDAR坐标系中的Z分量;
根据步骤一)至步骤三),同样得到旋转云台绕Y轴的旋转点,记为
Figure GDA00035412260600000917
其中
Figure GDA00035412260600000918
表示
Figure GDA00035412260600000919
点在2D-LiDAR坐标系中的X分量,
Figure GDA00035412260600000920
表示
Figure GDA00035412260600000921
点在2D-LiDAR坐标系中的Y分量,
Figure GDA00035412260600000922
表示
Figure GDA00035412260600000923
点在2D-LiDAR坐标系中的Z分量;完成对Z轴旋转轴和Y轴旋转轴的标零。
进一步的改进,包括如下步骤:
a.在对Z轴旋转轴和Y轴旋转轴标零之后,记录下标零时Z轴旋转轴和Y轴旋转轴的传感器角度值,记为θ0,α0,运用图像处理方法提取此时相机图像中 ArUco特征点的像素坐标,记为:
Figure GDA0003541226060000101
其中j表示旋转轴在零角度时,相机图像中第j个ArUco二维码,k则表示第j个ArUco二维码的第k个顶点,LC表示左相机图像,RC则表示右相机图像;
Figure GDA0003541226060000102
表示在旋转零角度时,左相机图像中,第j个ArUco二维码的第k个顶点在图像坐标系中的u轴坐标,
Figure GDA0003541226060000103
表示左相机图像中,第j个ArUco二维码的第k个顶点在图像坐标系中的v轴坐标;
Figure GDA0003541226060000104
分别表示右相机图像中相应点的u轴和v轴坐标;
根据Y轴转动的标定数据和Z轴转动的标定数据集,当第i次Z轴旋转轴和 Y轴旋转轴旋转角度分别αii,提取相机图像ArUco特征点:
Figure GDA0003541226060000105
i表示第i次旋转。
根据左相机的内参与外参,将左相机的图像坐标转化为标定板世界坐标:
Figure GDA0003541226060000106
式11中,i表示第i次旋转,j表示图像中的第j个ArUco二维码,k则表示二维码的第k个顶点;RW2C与TW2C分别表示左相机的旋转矩阵和平移矩阵, KLC为左相机内参数矩阵;
Figure GDA0003541226060000107
为左相机图像ArUco特征点的像素坐标,
Figure GDA0003541226060000108
分别为特征点在左相机坐标系中X、Y、Z轴的坐标,
Figure GDA0003541226060000109
分别为特征点在标定板世界坐标系中X、Y、Z轴对应坐标值;
因图像中ArUco二维码的四个顶点世界坐标均已知,故求得
Figure GDA00035412260600001010
根据2D-LiDAR与相机间关系标定结果,得到ArUco二维码各顶点在2D-LiDAR 坐标系中的空间坐标表示:
Figure GDA0003541226060000111
式12中,RL2LC与TL2LC分别表示2D-LiDAR与左相机之间的旋转矩阵和平移矩阵,
Figure GDA0003541226060000112
表示第i次旋转时,左相机图像的第j个ArUco二维码的第k个顶点在2D-LiDAR坐标系中的对应坐标值;
同理根据右相机的内参与外参,将右相机的图像坐标转化为2D-LiDAR坐标;在2D-LiDAR坐标系下,旋转云台绕Y轴的旋转轴为
Figure GDA0003541226060000113
绕轴旋转点为 [XAY,YAY,ZAY]T,将在2D-LiDAR坐标系下的ArUco特征点还原到绕Y轴旋转轴零角度时的坐标,坐标表示为:
Figure GDA0003541226060000114
式13中,α0表示Y轴旋转轴在零角度时,传感器记录的旋转角度值,αi则为第i次旋转时传感器记录的旋转角度值;
Figure GDA0003541226060000115
分别对应着在 2D-LiDAR坐标系下相应特征点还原至Y轴旋转轴在零角度时的坐标;
同样的,在2D-LiDAR坐标系下,旋转云台绕Z轴的旋转轴为
Figure GDA0003541226060000116
绕轴旋转点为[XAZ,YAZ,ZAZ]T,则ArUco二维码各顶点还原到绕Z轴旋转零角度,坐标表示为:
Figure GDA0003541226060000121
式14中,θ0表示旋转云台绕Z轴旋转零角度时,传感器记录的旋转角度值,θi则为第i次旋转时传感器记录的旋转角度值;
Figure GDA0003541226060000122
对应着相应特征点还原至2D-LiDAR坐标系下绕Z轴旋转轴在零角度时的坐标;利用相机内参数得到ArUco特征点还原至Z轴旋转轴和Y轴旋转轴均在零角度时的像素坐标,称之为理论像素坐标;具体推导公式如下:
Figure GDA0003541226060000123
式15中,
Figure GDA0003541226060000124
即为ArUco特征点在Z轴旋转轴和Y轴旋转轴均在零角度时的理论像素坐标,
Figure GDA0003541226060000125
为对应的理论2D-LiDAR坐标值;
关注理论像素坐标,与Z轴旋转轴和Y轴旋转轴均在零角度时所提取的ArUco 特征点像素坐标,建立误差函数如下:
Figure GDA0003541226060000126
Figure GDA0003541226060000127
Figure GDA0003541226060000131
式16中,ΔLu表示左相机图像在图像坐标系u方向的像素差,ΔLv表示左相机图像在图像坐标系v方向的像素差,ΔRu和ΔRv则表示右相机图像的相应差值,Δu和Δv则反映两个相机在u和v方向的整体像素偏差;根据式16,利用最小二乘法迭代求解旋转云台绕Z轴和Y轴的旋转参数,将Z轴旋转轴和Y轴旋转轴均在零点时的数据作为初始值,即:
Figure GDA0003541226060000132
Figure GDA0003541226060000133
最终解得旋转云台绕Z轴和Y轴的旋转参数,包括旋转轴
Figure GDA0003541226060000134
及其对应的绕轴旋转点SY,SZ
将所有在2D-LiDAR坐标系下的点按照如下公式还原到Z轴旋转轴和Y轴旋转轴均在零角度时的状态:
Figure GDA0003541226060000135
式17中,LX,LY,LZ表示在2D-LiDAR坐标系下表示的任意一点的X、Y、Z 轴坐标,α和θ则表示此时传感器记录的旋转云台绕Y轴和Z轴的旋转角度,LX0LY0LZ0则表示将所述任意一点还原至旋转平台绕Y、Z轴旋转零角度时的坐标值;
b.将旋转云台的旋转参数结果转换到旋转云台一体化坐标系内,其中Y轴旋转轴通过基座滑动连接有直线滑轨;
滑轨坐标系Osld与基座坐标系Obase之间旋转和平移均已知:
Figure GDA0003541226060000141
式18中,Mbase2sld(d)表示滑轨坐标系与基座坐标系之间的平移矩阵, Xsld,Ysld,Zsld表示滑轨坐标系中点的X、Y、Z轴的坐标值,而Xbase,Ybase,Zbase代表滑轨坐标系中的点在基座坐标系中的对应位置坐标;
基座坐标系Obase与云台偏转坐标系Oyaw之间旋转和平移均已知:
Figure GDA0003541226060000142
式19中,Ryaw2base和tyaw2base分别表示基座坐标系与云台偏转坐标系之间的旋转和平移矩阵,Xbase,Ybase,Zbase分别表示基座坐标系中点的X、Y、Z轴坐标值,而Xyaw,Yyaw,Zyaw则代表基座坐标系中的点在绕Z轴旋转的偏转坐标系中的对应位置坐标;
云台偏转坐标系Oyaw与旋转云台一体化坐标系Omot之间旋转和平移也已知:
Figure GDA0003541226060000143
式20中,Rmot2yaw和tmot2yaw分别表示偏转坐标系与旋转云台一体化坐标系之间的旋转和平移矩阵,Xyaw,Yyaw,Zyaw表示云台偏转坐标系中点的X、Y、Z轴坐标值,
Figure GDA0003541226060000144
则表示云台偏转坐标系的点在旋转云台Y、Z轴旋转角度均为零角度时,在一体化运动坐标中的对应坐标值;
还原旋转后的旋转云台一体化坐标系中点表示为:
Figure GDA0003541226060000151
式21中,Xmot,Ymot,Zmot表示在旋转云台一体化坐标系中点的X、Y、Z轴坐标值;β为旋转云台绕旋转云台一体化坐标系Z轴旋转的传感器记录角度,β0表示Z轴旋转零角度时传感器记录角度,α与α0分别表示旋转云台一体化坐标系Y 轴旋转的相应角度;
Figure GDA0003541226060000152
Figure GDA0003541226060000153
分别表示旋转云台一体化坐标系Y轴和Z轴的旋转矩阵;
在旋转云台一体化坐标系中
Figure GDA0003541226060000154
与所标定旋转云台三个方向的旋转轴
Figure GDA0003541226060000155
完全重合,其中,
Figure GDA0003541226060000156
Figure GDA0003541226060000157
分别为基于2D-LiDAR 坐标系的y轴和z轴旋转方向单位向量,且
Figure GDA0003541226060000158
的交点与旋转云台一体化坐标系原点Omot重合,计算
Figure GDA0003541226060000159
最小二乘交点:SAYAZ=[XAYAZ,YAYAZ,ZAYAZ]T,有如下推导:
Figure GDA00035412260600001510
Figure GDA00035412260600001511
式22中,[XLS,YLS,ZLS]分别表示2D-LiDAR坐标系内任一点的X、Y、Z轴坐标值,在旋转云台一体化坐标系中的对应的X、Y、Z轴坐标为[Xmot,Ymot,Zmot]; RLS2mot即为2D-LiDAR坐标系何旋转云台一体化坐标系之间的旋转矩阵,tLS2mot为平移矩阵;XAYAZ、YAYAZ和ZAYAZ分别表示
Figure GDA00035412260600001512
最小二乘交点SAYAZ在 2D-LiDAR坐标系中的X、Y、Z轴坐标;
综合式18-22,即有2D-LiDAR坐标系到滑轨坐标系的转换公式,表述如下:
Figure GDA0003541226060000161
Figure GDA0003541226060000162
Xsld、Ysld、Zsld分别表示2D-LiDAR坐标系中任意一点[XLS,YLS,ZLS],在还原至旋转云台旋转零角度时,对应滑轨坐标系中点的X、Y、Z轴坐标;
将相机图像坐标中任意点以及2D-LiDAR坐标系中任意点,至滑轨坐标系中的转换为结果坐标。
进一步的改进,所述步骤八包括如下步骤:
Ⅰ.将左右相机的视场统一为2D-LiDAR坐标系中的虚拟单目相机模式;
以HL,WL分别表示左相机图像在纵向和横向的像素点个数,HR,WR分别表示右相机图像在横向和纵向的像素点个数,左右相机t0~t4点在相机图像平面坐标系中的坐标分别为:
Figure GDA0003541226060000163
Figure GDA0003541226060000164
其中,
Figure GDA0003541226060000165
分别表示左相机图像中的四个角点,以及这些点在图像坐标系中u、v轴坐标位置,
Figure GDA0003541226060000166
则表示左相机图像中心点在图像坐标系中u、v 轴坐标位置;同样的
Figure GDA0003541226060000167
Figure GDA0003541226060000168
则表示右相机图像中相应点在图像坐标系中u、v轴坐标位置;
将图像点转换到相机坐标系中,随后转换到2D-LiDAR坐标系中,按照如下公式求得左右相机图像四个顶点与中心点在2D-LiDAR坐标系内坐标位置:
Figure GDA0003541226060000171
式25中,
Figure GDA0003541226060000172
分别表示前述左相机
Figure GDA0003541226060000173
在图像平面中的u、v轴像素坐标,
Figure GDA0003541226060000174
则分别表示右相机
Figure GDA0003541226060000175
在图像平面中的u、x轴像素坐标;
Figure GDA0003541226060000176
分别表示左、右相机内参的逆矩阵;
Figure GDA0003541226060000177
分别表示左、右相机与2D-LiDAR之间旋转矩阵的逆矩阵;TLS2LC,TLS2RC分别表示左、右相机与 2D-LiDAR之间的平移矩阵;
Figure GDA0003541226060000178
表示左相机
Figure GDA0003541226060000179
点在2D-LiDAR 坐标系中对应点的X、Y、Z轴坐标,而
Figure GDA00035412260600001710
则表示右相机
Figure GDA00035412260600001711
点在2D-LiDAR坐标系中对应点的X、Y、Z轴坐标;
左相机坐标系原点OLC和右相机坐标系原点ORC,在2D-LiDAR坐标系中对应点LCOLSRCOLS,其转换坐标计算如下:
Figure GDA00035412260600001712
Figure GDA00035412260600001713
式26中,
Figure GDA00035412260600001714
分别为左相机坐标系原点在2D-LiDAR坐标系中的X、Y、Z轴坐标,
Figure GDA00035412260600001715
为右相机坐标系原点在2D-LiDAR坐标系中X、Y、Z轴坐标;根据
Figure GDA00035412260600001716
Figure GDA00035412260600001717
与左、右相机坐标系原点在2D-LiDAR坐标系中相连所形成射线,射线与2D-LiDAR坐标中ZLS=h 平面的交点;ZLS表示2D-LiDAR坐标系的Z轴,h表示相机拍摄平面与2D-LiDAR 坐标系原点的垂直距离;ZLS=h平面即虚拟单目相机平面;
射线与ZLS=h平面的交点记为:
Figure GDA0003541226060000181
左相机的交点
Figure GDA0003541226060000182
的X、 Y与Z轴坐标计算如下:
Figure GDA0003541226060000183
同理得到右相机的相关交点
Figure GDA0003541226060000184
计算得到射线与ZLS=h平面的交点在2D-LiDAR坐标系中具***置后,得到 OE,T5,T6,T7,T8在2D-LiDAR坐标系位置;
OE表示左相机图像中心点与左相机坐标系原点OLC连线的直线与右相机图像中心点与右相机坐标系原点ORC连线的直线,在2D-LiDAR坐标系中的最小二乘交点;T5
Figure GDA0003541226060000185
Figure GDA0003541226060000186
连线在2D-LiDAR坐标系YLS=0平面的交点;T6
Figure GDA0003541226060000187
Figure GDA0003541226060000188
的连线中点在2D-LiDAR坐标系YLS=0平面的交点;T7
Figure GDA0003541226060000189
Figure GDA00035412260600001810
的连线在 2D-LiDAR坐标系XLS=0平面的交点;T8
Figure GDA00035412260600001811
Figure GDA00035412260600001812
的连线在2D-LiDAR坐标系 XLS=0平面的交点;
OE在2D-LiDAR坐标系中位置,通过两条直线的最小二乘交点确定。所述两条直线分别过左右相机的原点与T0点,具体计算方法:
Figure GDA00035412260600001813
Figure GDA00035412260600001814
式28中,X,Y,Z分别表示需要计算的T0在2D-LiDAR坐标X、Y、Z轴的坐标值;完成基于旋转云台一体化坐标系的虚拟单目相机视场的关键点计算;
Ⅱ.旋转图像的空间投影及拼接
由于相机在拍摄时,相平面与被拍摄平面并非平行,需要将图像投影到被拍摄物体的平面,随后将左右相机的图像进行拼接,以实现在三维模型表面纹理贴图等目的;
2.1.单相机旋转图像的空间投影变换
相机图像从各自的相机坐标系,投影虚拟单目相机平面:
计算左、右相机相平面四个角点在2D-LiDAR坐标系的对应坐标:
Figure GDA0003541226060000191
Figure GDA0003541226060000192
以及左、右相机坐标系原点在 2D-LiDAR坐标系中的对应空间坐标:LCOLSRCOLS
此时传感器记录的旋转云台绕Y轴旋转角度α和绕Z轴旋转角度θ,根据式 17,将左、右相机相平面角点和左、右相机坐标系原点在2D-LiDAR坐标系的位置还原到旋转云台在旋转零角度状态时的对应坐标,分别记为
Figure GDA0003541226060000193
Figure GDA0003541226060000194
LCOLS-0RCOLS-0
得到左、右相机相平面角点和坐标系原点的连线,在还原至旋转零角度后的 2D-LiDAR坐标系Z=h平面的交点位置;对于左相机交点位置分别为:
Figure GDA0003541226060000195
计算公式如下:
Figure GDA0003541226060000196
式29中,
Figure GDA0003541226060000197
表示
Figure GDA0003541226060000198
点坐标值的 XYZ分量,
Figure GDA0003541226060000199
则分别表示左相机相平面四个角点还原至旋转零角度后的点
Figure GDA00035412260600001910
的XYZ分量,LCOLS-0(X),LCOLS-0(Y),LCOLS-0(Z)则分别表示左相机坐标系原点在还原至旋转零角度后的LCOLS-0XYZ分量;
左相机相平面四个角点和左相机坐标系原点在还原至旋转零角度后的点在 2D-LiDAR坐标系中的最大和最小X、Y坐标值,并作为透视变换后图像尺寸的变换依据:
Figure GDA0003541226060000201
Figure GDA0003541226060000202
式30中,
Figure GDA0003541226060000203
表示在
Figure GDA0003541226060000204
四个点中,X分量最大值,
Figure GDA0003541226060000205
表示X分量的最小值,
Figure GDA0003541226060000206
Figure GDA0003541226060000207
分别为Y分量的最大值和最小值;
同样得到右相机相平面四个角点和右相机坐标系原点在还原至旋转零角度后的点在2D-LiDAR坐标系中的最大和最小X、Y坐标值;
2.2.设定透视变换后的图像尺寸与原图像尺寸相同,则原图像四个顶点对应的透视变换后图像坐标,对于左相机按如下公式进行计算:
Figure GDA0003541226060000208
式31中,
Figure GDA0003541226060000209
Figure GDA00035412260600002010
分别表示原左相机图像四个角点在透视变换后图像中的像素坐标;
根据四个控制点在原始图像和透视变换后图像中的坐标位置,计算透视变换矩阵,即得到左相机图像的透视变换结果;同理可得右相机图像的透视变换结果;
2.3.旋转图像的拼接
将左右相机的图像都变换到同一图像上,再对两个变换后的图像进行权重相加形成一张拼接后的图,于是有:
Figure GDA00035412260600002011
Figure GDA00035412260600002012
式32中,
Figure GDA0003541226060000211
Figure GDA0003541226060000212
分别表示在拼接图像四个角点与Z=h平面交点在X 轴方向的最大值和最小值,
Figure GDA0003541226060000213
Figure GDA0003541226060000214
表示拼接图像四个角点与Z=h平面交点在Y轴方向的最大值和最小值;
设定拼接后图像纵向的像素个数为H=HL,则横向的像素个数为:
Figure GDA0003541226060000215
式33中,HL表示左相机原始图像在纵向的像素点数,W为新拼接图像在横向的像素点数;
得到左相机图像四个顶点在新拼接图像中的像素坐标:
Figure GDA0003541226060000216
在式34中,
Figure GDA0003541226060000217
Figure GDA0003541226060000218
分别表示左相机图像的四个角点在新拼接图像中的像素坐标;
同样得到右相机图像的四个角点在新拼接图像中的像素坐标;
将左相机图像和右相机图像通过加权相加的方法拼接得到新拼接图像。
进一步的改进,将多组旋转云台在不同旋转角度位置拍摄的图像投影到 ZLS=h平面,得到对应的新拼接图像;在世界坐标系内旋转新拼接图像至旋转云台处于零点时的位置,实现旋转图像在3D模型中的还原贴图。
附图说明
图1为本发明实施例中采集末端结构示意图;
图2为本发明实施例中采集末端及其各运动关节坐标系示意图;
图3为本发明实施例中标定板示意图;
图4为本发明实施例中旋转矩形求解标定板角点示意图;
图5为本发明实施例中ChArUco板定板图;
图6为本发明实施例中相机图像ArUco检测结果图;
图7为本发明实施例中相机图像ChArUco关键角点检测结果图;
图8为本发明实施例中标定板几何约束示意图;
图9为本发明实施例中相机与2D-LiDAR关系标定结果校验图;
图10为本发明实施例中空间点绕Y轴和Z轴旋转简图;
图11为本发明实施例中2D-LiDAR坐标系中的虚拟单目相机关键点示意图;
图12为本发明实施例中左相机图像透视变换效果图;
图13为本发明实施例中右相机图像透视变换效果图;
图14为为本发明实施例中左相机图像透视变换效果图和右相机图像透视变换效果图合并的效果图;
图15为本发明实施例中激光点云标定板角点绕Z轴还原结果图;
图16为本发明实施例中激光点云标定板角点绕Y轴还原结果图;
图17为本发明实施例中多旋转图像的投影与拼接应用图;
图18为本发明实施例中一体化自动标定流程图。
其中,1相机,22D-LiDAR,3标定板,4Z轴旋转轴,5Y轴旋转轴,6基座, 7直线滑轨,31ChArUco标定板,32矩形侧板,33等腰直角三角形侧板,34第一边侧板,35第二边侧板。
具体实施方式
以下结合附图及实施例对本发明做进一步说明。
图18是本发明实施例提供的一种2D-LiDAR和相机的多自由度采集末端的联合标定方法的流程图。所述标定方法包括:
步骤一:获取2D-LiDAR单次扫描标定板的点云数据,提取标定板结构角点信息,根据角点信息确定多自由度采集末端在各个旋转自由度时的零点位置;
步骤二:根据后续各项标定需求,驱动采集末端按照特定轨迹运动,并采集相机和2D-LiDAR数据形成完整的标定数据集;
步骤三:根据相机获取的图像,提取出图像中ChArUco关键角点坐标信息,根据采集末端不同旋转角度时相机图像中相应二维码特征点的平面坐标信息,计算相机内参;
步骤四:在获取激光扫描U型标定板的点云数据后,提取U型标定板角点,根据标定板特定结构计算各角点在相应相机图像中的平面坐标信息,结合各个角点的空间和平面坐标信息,确定2D-LiDAR和相机之间的旋转矩阵和平移矩阵;
步骤五:获取2D-LiDAR绕采集末端Y轴和Z轴旋转得到的点云数据和相机图像,采用由粗到细的搜索策略,实现旋转云台的Y轴和Z轴旋转参数的联合标定;
步骤六:采集末端标定结果在其一体化的运动坐标系内表示;旋转云台一体化坐标系内关键点计算及运用;
步骤七:旋转后,左右相机图像的投影变换;投影变换后多图像的拼接。
下面进行详细说明:
(1)采集末端零点位置标定:
图3中的U型标定板中间为ChArUco标定板(如图5所示),一侧为两个互相垂直且尺寸已知的矩形板,另一侧为一块直角边与ChArUco标定板边长相同的等腰直角三角形板,有一块矩形板与其斜边相接。
将如图3所示U型标定板放置在竖直支架上,并将处于机械零点状态的采集末端设备放置在支架大概中心位置之下,采集2D-LiDAR扫描标定板的激光点云。
由于2D-LiDAR采样的角度以及激光与标定板轮廓之间的角度变化,得到的点云数据有的区域比较稀疏,影响后续提取激光轮廓关键点,同时有效的激光点也减少了,需要对激光点云进行稀疏激光均匀化上采样,激光滤波,得到均匀连续的激光点云。将均匀连续的激光点云转换成图像,对图像进行膨胀、二值化和细化等形态学处理,提取点云清晰连续的骨架。获取图像中点云的轮廓,根据轮廓中包含的点的个数判断该轮廓是否为有效激光轮廓。
根据得到的有效激光轮廓,通过直线拟合提取轮廓中的角点。在直线拟合过程中,提出了可旋转矩形的思想如图4。其中,d和s是我们给定的距离值,计算外扩矩阵:
Figure GDA0003541226060000231
公式1中,K1与K2表示经过前述图像处理后,计算得到的初始U型标定板轮廓关键点(如图3所示),K3和K4则表示经过外扩后U型标定板轮廓关键点的更新位置。
随后通过给定的阈值s,计算图4中矩形的4个顶点,并获取包含这4个顶点的最小外接矩阵。通过求解出的矩形方程,来判断激光点云在哪个矩形内,由此将所有激光点云以此分配带其所属矩形,随后将矩形内的激光点拟合得到直线,根据两两直线的交点确定相应的角点,分别为p0,p1,p2,p3,p4,p5
驱动采集末端设备首先绕其俯仰Y轴(如图2中的Ymot所示)旋转某一正负角度范围,采集序列激光剖面,并提取线段p2p3(如图3中所示)的角度,取线段角度最小的点为Y轴的零点,记录该帧点云对应的Y轴机械角度θ,并驱动Y轴旋转至该角度;然后绕偏转Z轴(如图2中Zmot所示)旋转某一正负角度范围,采集序列激光剖面,并提取线段p2p3的长度,取长度最小的点为Z轴的零点,并记录该帧点云对应的Z轴机械角度β,将{θ,β}角度作为采集末端的零点位置,即完成采集末端的零点位置标定。
(2)标定数据采集
根据(1)中所述零点标定方法,将采集末端的俯仰自由度(Y轴)与偏转自由度(Z轴)均置于零位置后,分两步进行标定所需要数据的采集。
保持偏转自由度(Z轴)处于零位置不动,驱动采集末端绕俯仰自由度(Y 轴)转动,每次转动约1°。每次绕Y轴转动后,控制采集末端上的左右相机拍摄图像并保存,记录2D-LiDAR一帧点云数据,同时记录采集末端在俯仰自由度和偏转自由度的旋转角度,以此作为一组标定数据。在每个旋转位置按照此方式保存对应的一组标定数据,以此为绕Y轴转动的标定数据集。
类似的,保持俯仰自由度(Y轴)处于零位置不动,驱动采集末端绕偏转自由度(Z轴)转动,每次转动约1°。将过程中每组相机图像、激光点云和采集末端姿态数据保存下来,作为绕Z轴转动的标定数据集。
(3)相机内参标定
本发明采用的是棋盘格与ArUco二维码相结合的ChArUco标定板(如图5)。
进行相机内参标定前,需要对要标定相机所拍摄的图像进行一系列处理:获取(2)所述绕Y轴转动的标定数据集中,由需标定相机所拍摄的标定板图像;提取每张图像中ArUco特征点(如图6所示),包括其ID号和二维码四个顶点的平面坐标;进一步的,根据之前所述ArUco特征点,提取出图像的ChArUco 角点坐标序列(如图7中绿色圆圈所示结果);基于采集末端在绕Y轴旋转的多个位置时,图像坐标系上的角点序列和相对应的世界坐标序列,根据棋盘格标定方法标定出相机的内参数矩阵K和畸变参数矩阵D:具体
Figure GDA0003541226060000251
D=[k1 k2 k3 p1 p2]
公式2
公式2中,fx,fy分别表示x与y方向上像素与实际距离的比例,也叫做x 轴和y轴上的归一化焦距;u0,v0分别为图像中心在像素平面的位置;k1,k2,k3则表示镜头的径向畸变参数;p1,p2表示镜头的切向畸变参数。
(4)2D-LiDAR与相机的外参标定
根据(2)所述采集末端绕Y轴旋转所采集到的数据集,获取激光点云数据,按照(1)所述方法提取激光轮廓中的标定板角点,即可得到所述角点(图3中所示p1、p2、p3和p4)在2D-LiDAR坐标系中的空间坐标{Pi L|i=1,2,3,4}。
根据得到的所述激光轮廓的角点坐标,通过标定板结构的几何约束(图8 所示)。假定标定板处于U型标定板世界坐标系Z=0平面,于是通过前述 {Pi L|i=1,2,3,4}以及标定板的几何约束,即可得到所述标定板角点在世界坐标系的坐标{Pi W|i=1,2,3,4},即有:
Figure GDA0003541226060000252
Figure GDA0003541226060000253
Figure GDA0003541226060000254
Figure GDA0003541226060000261
公式3中,d为标定二维码边缘与标定板边缘的距离,w为标定板的边长。
通过前述的相机图像中提取的ChArUco关键角点作为控制点,结合之前标定的相机内参数及其畸变参数,通过相机标定方法[1]标定出相机外参,即相机坐标系与U型标定板世界坐标系之间的旋转矩阵RW2C与平移矩阵TW2C,然后确定所述U型标定板角点在图像坐标系中的实际坐标:
Figure GDA0003541226060000262
公式4中,WXi,WYi,WZ对应着公式3中所观测到的标定板角点坐标X、Y、Z 分量;
Figure GDA0003541226060000263
表示观测点在相机坐标系中坐标值;
Figure GDA0003541226060000264
表示观测点在图像中的像素坐标值。
同样的,可以推导上述标定板角点在图像坐标系中的理论坐标信息:
Figure GDA0003541226060000265
公式5中,RL2C和TL2C分别表示激光与相机间旋转矩阵和平移矩阵;LXi,LYi,LZi对应着公式3中标定板角点的点云数据{Pi L|i=1,2,3,4}的X、Y、Z分量;CXi,CYi,CZi表示标定板角点在相机坐标系中对应的理论坐标值;Iui,Ivi表示其在图像中的理论像素坐标值。
至此,定义损失函数为:
Figure GDA0003541226060000266
采用最小二乘法,迭代求解激光与相机间旋转矩阵和平移矩阵。根据求解结果,将激光点云数据转换到对应相机的图像坐标系上,可以直观看到校验结果,如图9所示,图中绿色与蓝色的点线即为激光点云数据。
(5)旋转轴参数标定
根据(2)中采集到的绕Y轴和Z轴旋转标定数据集,分别用来进行旋转云台的Y轴和Z轴旋转参数标定。需要特别指出的是,标定结果包括代表旋转轴方向的单位向量和绕轴旋转点,Y轴和Z轴的旋转参数标定结果最后都统一到同一个坐标系下。
首先,根据旋转云台的结构,初步假设绕Y轴旋转的旋转向量就是2D-LiDAR 坐标系中Y轴(0,1,0),通过由粗到细的搜索方法求解绕Y轴的旋转点,对于 Z轴也是如此;随后,利用相机图像中的ArUco特征点作为外部控制点,将旋转时相机图像的ArUco特征点,与旋转零角度时相机图像的ArUco特征点之间的偏差作为误差函数,利用最小二乘法,迭代求解误差函数最小时的旋转轴与旋转点。而前述由粗到细搜索方法求解旋转云台的Y轴、Z轴的旋转参数,则作为迭代求解时Y、Z轴的旋转轴向量和旋转点的初始值。
①采用由粗到细的搜索方法计算旋转轴参数初值
获取转台绕Z轴旋转不同角度的一系列激光点云数据,读取激光的单次点云数据,记录转台第i次旋转的角度θi。按照(1)中所述方法与步骤,对第i次转台绕Z轴旋转后采集到的原始点云数据,进行稀疏激光均匀化上采样和滤波处理,得到有效合理的点云数据,提取激光轮廓的关键点(如图3所示,仅需要 p2p3)。最后得到关于标定板结构的关键角点与采集末端绕Z轴旋转角度的序列 {LSPi ji}(j=2,3),这里i表示第i次旋转,j表示p2或p3点,LSPi j则表示第 i次旋转时p2或p3点从激光点云获取的坐标(2D-LiDAR坐标系下)。
接下来,将对旋转云台的旋转运动带来空间坐标的变化进行说明。将所述 2D-LiDAR绕多自由度的旋转简化抽象为空间上一点绕两向量旋转,所述旋转过程的几何模型如图10所示,存在如下关系:
Figure GDA0003541226060000281
Figure GDA0003541226060000282
where:
Figure GDA0003541226060000283
Sz=(XAZ,YAZ,ZAZ)T
r0Lnx·Lny·(1-cosθ)+cosθ,r1Lnx·Lny·(1-cosθ)+Lnz·sinθ
r2Lnz·Lnx·(1-cosθ)-Lny·sinθ,r3Lnx·Lny·(1-cosθ)-Lnz·sinθ
r4=(Lny)2·(1-cosθ)+cosθ,r5Lny·Lnz·(1-cosθ)+Lnx·sinθ
r6Lnz·Lnx·(1-cosθ)+Lny·sinθ,r7Lny·Lnz·(1-cosθ)-Lnx·sinθ
r8=(Lnz)2·(1-cosθ)+cosθ
公式7
公式7中,θ为旋转云台绕Z轴旋转的角度,
Figure GDA0003541226060000284
表示绕Z轴旋转的方向向量,SZ表示绕Z轴旋转点,
Figure GDA0003541226060000285
则表示旋转矩阵。P1表示任意一点,其绕
Figure GDA0003541226060000286
轴旋转θ后对应点为P0
类似的,可确定点P0绕采集末端
Figure GDA0003541226060000287
轴旋转α后的几何关系:
Figure GDA0003541226060000288
公式8中,α表示旋转云台绕Y轴旋转的角度,
Figure GDA0003541226060000289
表示绕Y轴旋转的方向向量,SY则表示绕Y轴旋转点,
Figure GDA00035412260600002810
表示其对应的旋转矩阵。P2则表示P0在绕
Figure GDA00035412260600002811
旋转α后的对应点。
于是将绕Z轴旋转后的U型标定板角点(如图3所示,仅需要关注p2p3)还原到旋转零角度的公式即为:
Figure GDA00035412260600002812
公式9中,θ0表示在旋转云台绕Z轴旋转标定为零角度时,传感器记录的实际角度,而θi则表示在第i次绕Z轴旋转时所记录的角度值;
LSPi j(X),LSPi j(Y),LSPi j(Z)表示旋转云台第i次绕Z轴旋转后,U型标定板第j个关键角点在2D-LiDAR坐标系中的坐标值,
Figure GDA0003541226060000291
则表示这些关键角点在还原至绕Z轴旋转零角度后的坐标值。
依照上述方法,将所有采集末端绕Z轴旋转采集的激光点云数据均还原到绕Z轴旋转零角度,即形成如图15所示的点云数据。
假设旋转云台绕Z轴的旋转轴向量为
Figure GDA0003541226060000292
只需要计算绕Z轴旋转点即可。接下来,根据由粗到细的搜索策略,先利用循环迭代的方法进行粗略的求解:
1)初始时,假设旋转云台绕Z轴的旋转点为SZ(0,0,0);
2)按照之前所述公式,将采集末端绕Z轴旋转采集到的所有激光点云数据,均还原到绕Z轴旋转角度为零时的状态;
3)根据还原后的点云数据,计算得到包含这些点云数据的最小外接矩阵,记此矩阵的四个顶点为{Lqk}(k=1,2,3,4),如图15;
4)将上个步骤中四个顶点连线,形成四条线段{lk}(k=1,2,3,4),随后给定阈值d和s,按照(1)中计算标定板角点相同的方法,获得如图 4所示lk线段各自的外扩矩形;
5)依次判断激光点云落在哪个外扩矩形内,并形成每个外扩矩形的激光点集。根据每个矩形内点集拟合直线,并计算矩形内所有点到其对应直线的距离累加值;
6)计算四个外扩矩形的点到拟合直线距离累加值之和,如果其小于上一次迭代的累加值之和,则根据阈值Δy将旋转点SZ的Y分量加上Δy
7)若迭代次数超过一定阈值,则结束,否则回到2)继续。
上述迭代过程停止之后,可以得到粗略求解的采集末端绕Z轴旋转点,记为
Figure GDA0003541226060000293
进一步的,将
Figure GDA0003541226060000294
作为细化搜索的迭代初值,设定一个绕Z轴旋转点Y分量的迭代阈值dy,dy远小于上述迭代中的Δy。同时将细化搜索时旋转点Y值搜索范围控制在
Figure GDA0003541226060000301
之内。于是按照上述粗略搜索的迭代过程,进一步细化搜索采集末端绕Z轴的旋转点,最终结果为
Figure GDA0003541226060000302
对于旋转云台绕Y轴的旋转轴向量与旋转点,其计算方法与Z轴类似。同样获取标定板的关键角点与绕Y轴旋转的序列{LSPi ji}(j=1,2,3,4)。类似的,i表示第i次旋转,j表示图3中所示的p1、p2、p3或p4点,LSPi j则表示第i次旋转时U型标定板角点从激光点云获取的坐标(2D-LiDAR坐标系下)。如图16所示,设定其初始旋转轴向量和旋转点分别为:
Figure GDA0003541226060000303
SY(0,0,0),这里
Figure GDA0003541226060000304
SY含义与公式8中相同。
采用与绕Z轴标定计算时同样的公式,计算激光点云还原到绕Y轴零角度时的坐标,如图16所示,随后进行粗略搜索得到的绕Y轴旋转点,记为
Figure GDA0003541226060000305
进一步细化搜索,最终旋转云台绕Y轴旋转点的结果,记为
Figure GDA0003541226060000306
②基于图像控制点的旋转参数优化
根据①所述,在对Y,Z轴标零之后,记录下此时的旋转云台绕Y轴和Z轴的传感器角度值,记为θ0,α0,运用图像处理方法提取此时相机图像中ArUco 特征点(ArUco的四个顶点)的像素坐标,记为:
Figure GDA0003541226060000307
其中j表示旋转轴在零角度时,相机中第j个ArUco二维码,k则表示此二维码的第k个顶点,LC表示左相机图像,RC则表示右相机图像。
根据(2)所述采集的标定数据集,当第i次绕Y轴和Z轴角度分别αii, 提取相机图像ArUco特征点:
Figure GDA0003541226060000308
与前述类似,旋转零角度时特征点表示类似,这里i表示第i次旋转。
以左相机为例,根据左相机的内参与外参,可以将图像坐标转化为世界坐标:
Figure GDA0003541226060000311
公式11中,i表示第i次旋转,j表示图像中的第j个ArUco二维码,k则表示二维码的第k个顶点;RW2C与TW2C表示左相机的外参数(旋转矩阵和平移矩阵),KLC为左相机内参数矩阵;
Figure GDA0003541226060000312
为左相机图像ArUco特征点的像素坐标,
Figure GDA0003541226060000313
是这些特征点在左相机坐标系中的坐标表示,
Figure GDA0003541226060000314
则为特征点在U型标定板世界坐标系中的对应坐标值。
因图像中ArUco二维码的四个顶点世界坐标均已知,故可从上述公式求得
Figure GDA0003541226060000315
进一步的,根据(4)所述2D-LiDAR与相机间关系标定结果,可以推导出上述ArUco二维码各顶点在2D-LiDAR坐标系中的空间坐标表示:
Figure GDA0003541226060000316
公式12中,RL2LC与TL2LC分别表示2D-LiDAR与左相机之间的旋转矩阵和平移矩阵,
Figure GDA0003541226060000317
表示第i次旋转时,左相机图像的第j个ArUco二维码的第k个顶点在2D-LiDAR坐标系中的对应坐标值,其他符号含义与公式11中相同。
在2D-LiDAR坐标系下,旋转云台绕Y轴的旋转轴为
Figure GDA0003541226060000318
绕轴旋转点为 [XAY,YAY,ZAY]T,与①中公式9类似,可将ArUco特征点(在2D-LiDAR坐标系下) 还原到绕Y轴旋转零角度,其坐标表示为:
Figure GDA0003541226060000321
公式13中,α0表示旋转云台绕Y轴旋转零角度时,传感器记录的旋转角度值,αi则为第i次旋转时传感器记录的旋转角度值;
Figure GDA0003541226060000322
对应着相应特征点还原至绕Y轴旋转零角度时的坐标(2D-LiDAR坐标系下),其余符号含义与公式8中相同。
同样的,在2D-LiDAR坐标系下,旋转云台绕Z轴的旋转轴为
Figure GDA0003541226060000323
绕轴旋转点为[XAZ,YAZ,ZAZ]T,则ArUco二维码各顶点还原到绕Z轴旋转零角度,其坐标表示为:
Figure GDA0003541226060000324
公式14中,θ0表示旋转云台绕Z轴旋转零角度时,传感器记录的旋转角度值,θi则为第i次旋转时传感器记录的旋转角度值;
Figure GDA0003541226060000325
对应着相应特征点还原至绕Z轴旋转零角度时的坐标(2D-LiDAR坐标系下),其余符号与公式13中相同。
最后,利用相机内参数得到ArUco特征点还原至,旋转云台各自由度旋转零角度时的像素坐标,称之为理论像素坐标。具体推导公式如下:
Figure GDA0003541226060000326
公式15中,
Figure GDA0003541226060000331
即为ArUco特征点在旋转云台各自由度旋转零角度时的理论像素坐标,
Figure GDA0003541226060000332
为对应的理论2D-LiDAR坐标值,其余符号含义与公式12和公式14中相同。
仅关注理论像素坐标,与实际旋转云台各自由度旋转零角度时所提取的 ArUco特征点像素坐标(公式10中定义),建立误差函数如下:
Figure GDA0003541226060000333
Figure GDA0003541226060000334
Figure GDA0003541226060000335
公式16中,ΔLu表示左相机图像在图像坐标系u方向的像素差,ΔLv表示左相机图像在图像坐标系v方向的像素差,类似的,ΔRu和ΔRv则表示右相机图像的相应差值,Δu和Δv则反映两个相机在u和v方向的整体像素偏差。其余符号与前述公式中意义相同。
根据公式16,利用最小二乘法迭代求解Y轴与Z轴旋转参数,将①中结果作为初始值,即:
Figure GDA0003541226060000336
Figure GDA0003541226060000337
最终解得旋转云台绕Z轴和Y轴的旋转参数,包括旋转轴
Figure GDA0003541226060000341
及其对应的绕轴旋转点SY,SZ
经过采集末端旋转后,所有在2D-LiDAR坐标系下的点均可按照如下公式还原到采集末端处于零位时的状态:
Figure GDA0003541226060000342
公式17中,LX,LY,LZ表示在2D-LiDAR坐标系下表示的任意一点,α和θ则表示此时传感器记录的旋转云台绕Y轴和Z轴的旋转角度,LZ0,LZ0LZ0则表示将此任意点还原至旋转平台绕Y、Z轴旋转零角度时的坐标值(2D-LiDAR坐标系下)。
按照上述结果与公式17,将(2)中所述激光点云数据还原到Y轴和Z轴旋转零角度时的状态,还原后的激光点云数据能够很好的契合在一起,并且能正确的反应出U型标定板的几何结构,以此证明标定方法和结果的正确性和准确性。
(6)旋转图像空间投影及拼接应用
提出旋转云台一体化坐标系的目的,是将相机的视场、激光测量距离与范围,均在一个坐标系中进行表示,同时也要将,在2D-LiDAR坐标系中表示的采集末端旋转参数,转换到此坐标系中。
具体步骤可以描述为:计算出2D-LiDAR与旋转云台一体化坐标系之间的旋转和平移关系,将旋转参数转换到旋转云台一体化坐标系中;计算左、右相机视场范围关键点,设置为虚拟单目相机视场的关键点,并将这些关键点在2D-LiDAR 坐标系中的空间位置转换到采集末端的旋转云台一体化坐标系中;最后将前述结果用于旋转图像的空间投影及拼接。
①采集末端各自由度旋转参数到旋转云台一体化坐标系的转换
前述采集末端各自由度旋转参数的标定方案均是基于2D-LiDAR坐标系下的,需要将结果转化到旋转云台一体化坐标系下,以便于后续的虚拟单目相机关键点计算,以及旋转图像的拼接处理。
为此,需要将旋转云台的旋转参数结果转换到旋转云台一体化坐标系内(如图2中Omot坐标系)。在认为采集机构的加工精度足够的前提下,按照设计图纸,滑轨坐标系Osld与基座坐标系Obase之间旋转和平移均已知:
Figure GDA0003541226060000351
公式18中,Mbase2sld(d)表示滑轨坐标系与基座坐标系之间的平移矩阵(与基座滑动位移d相关),Xsld,Ysld,Zsld表示滑轨坐标系中点的坐标值,而 Xbase,Ybase,Zbase代表这些点在基座坐标系中的对应位置坐标。
基座坐标系Obase与云台偏转坐标系Oyaw之间旋转和平移均已知:
Figure GDA0003541226060000352
公式19中,Ryaw2base和tyaw2base分别表示基座坐标系与云台偏转坐标系之间的旋转和平移矩阵(通过已知设计尺寸计算得到),Xbase,Ybase,Zbase表示基座坐标系中点的坐标值,而Xyaw,Yyaw,Zyaw则代表这些点在偏转(即绕Z轴旋转)坐标系中的对应位置坐标。
云台偏转坐标系Oyaw与旋转云台一体化坐标系Omot之间旋转和平移也已知:
Figure GDA0003541226060000353
公式20中,Rmot2yaw和tmot2yaw分别表示偏转坐标系与旋转云台一体化坐标系之间的旋转和平移矩阵(通过已知设计尺寸计算得到),Xyaw,Yyaw,Zyaw表示云台偏转坐标系中点的坐标值,
Figure GDA0003541226060000361
则表示这些点在旋转云台Y、Z轴旋转角度均为零角度时,在一体化运动坐标中的对应坐标值。
在旋转云台一体化坐标系中,认为所有旋转均是绕原点进行旋转,即旋转云台绕Y轴旋转的旋转轴为旋转云台一体化坐标系的Y轴方向,旋转点为旋转云台一体化坐标系原点;旋转云台绕Z轴旋转的旋转轴为旋转云台一体化坐标系的Z 轴方向,旋转点也是旋转云台一体化坐标系的原点。于是还原旋转后的旋转云台一体化坐标系中点可表示为:
Figure GDA0003541226060000362
公式21中,Xmot,Ymot,Zmot表示在旋转云台一体化坐标系中点的坐标值,
Figure GDA0003541226060000363
含义见公式20;β对应着此时旋转云台绕旋转云台一体化坐标系Z 轴旋转的传感器记录角度,β0表示Z轴旋转零角度时传感器记录角度,同理α与α0分别表示旋转云台一体化坐标系Y轴旋转的相应角度;
Figure GDA0003541226060000364
Figure GDA0003541226060000365
分别表示旋转云台一体化坐标系Y轴和Z轴的旋转矩阵,矩阵具体形式见公式7。
因为旋转云台一体化坐标系中定义的
Figure GDA0003541226060000366
应与之前所标定旋转云台旋转轴
Figure GDA0003541226060000367
完全重合(
Figure GDA0003541226060000368
Figure GDA0003541226060000369
即为之前标定基于 2D-LiDAR坐标系的归一化旋转方向向量),且
Figure GDA00035412260600003610
的交点应与旋转云台一体化坐标系原点Omot重合,计算
Figure GDA00035412260600003611
最小二乘交点:SAYAZ=[XAYAZ,YAYAZ,ZAYAZ]T,有如下推导:
Figure GDA0003541226060000371
Figure GDA0003541226060000372
公式22中,[XLS,YLS,ZLS]表示2D-LiDAR坐标系内任一点的坐标值,其在旋转云台一体化坐标系中的对应坐标[Xmot,Ymot,Zmot];RLS2mot即为上述两个坐标系之间的旋转矩阵,tLS2mot为平移矩阵,由之前标定的Y、Z轴旋转方向向量及其最小二乘交点确定。
综合公式18-22,即有2D-LiDAR坐标系到滑轨坐标系的转换公式,可以表述如下:
Figure GDA0003541226060000373
结合之前所述相机内参以及2D-LiDAR与相机间关系的标定结果(见公式5),可计算在采集末端进行多自由度旋转后,相机图像坐标中任意点以及2D-LiDAR 坐标系中任意点,至滑轨坐标系中的转换结果坐标。
同样的,对基于旋转云台一体化坐标系的旋转云台Y、Z轴旋转参数标定结果,运用激光点云数据进行校验,其结果原图并无差别。进一步的,随机旋转Y 轴和Z轴,利用2D-LiDAR对办公室环境进行扫描,最后根据基于旋转云台一体化坐标系的旋转参数对激光点云数据进行还原。
②虚拟单目相机视场的关键点计算
根据4.4中所述2D-LiDAR与相机间参数的标定结果,可以将左右相机的视场统一为2D-LiDAR坐标系中的虚拟单目相机模式,以便于后续的左右相机图像拼接,如图11所示。
图中
Figure GDA0003541226060000374
表示左相机图像中的四个顶点在左相机坐标系中的位置,与左相机坐标系原点OLC连线的直线,在2D-LiDAR坐标系ZLS=h平面的角点位置,而
Figure GDA0003541226060000381
则表示右相机图像的四个顶点与此平面的相交位置。 OE则表示左相机图像中心点与左相机坐标系原点OLC连线的直线,与左相机图像中心点与左相机坐标系原点ORC连线的直线,在2D-LiDAR坐标系中的最小二乘交点。取T5
Figure GDA0003541226060000382
Figure GDA0003541226060000383
连线在2D-LiDAR坐标系YLS=0平面的交点,T6
Figure GDA0003541226060000384
Figure GDA0003541226060000385
的连线中点在2D-LiDAR坐标系YLS=0平面的交点,T7
Figure GDA0003541226060000386
Figure GDA0003541226060000387
的连线在 2D-LiDAR坐标系XLS=0平面的交点,T8
Figure GDA0003541226060000388
Figure GDA0003541226060000389
的连线在2D-LiDAR坐标系 XLS=0平面的交点,T0则表示2D-LiDAR坐标系Z轴与ZLS=h交点,即T0(0,0,h)。
上述OE,T0,T5,T6,T7,T8点集,即为在2D-LiDAR坐标系下,用来表示左右相机视场等信息的虚拟单目相机视锥体关键点。并且由于左右相机与2D-LiDAR之间为刚体结构,故这些点在2D-LiDAR坐标系中的相对位置不会随着采集末端的旋转或运动而改变。下面将叙述如何计算这些点在2D-LiDAR坐标系中的位置:
以HL,WL分别表示左相机图像在纵向和横向的像素点个数,HR,WR分别表示右相机图像在横向和纵向的像素点个数,左右相机T0~T4点在相机图像平面坐标系中的坐标分别为:
Figure GDA00035412260600003810
Figure GDA00035412260600003811
图像点转换到相机坐标系中,随后将其转换到2D-LiDAR坐标系中,按照如下公式即可求得左右相机图像内顶点各自在2D-LiDAR坐标系内位置:
Figure GDA00035412260600003812
公式25中,
Figure GDA00035412260600003813
表示前述左相机T0~T4在图像平面中的像素坐标,
Figure GDA00035412260600003814
则表示右相机T0~T4在图像平面中的像素坐标;
Figure GDA00035412260600003815
分别表示左、右相机内参的逆矩阵;
Figure GDA0003541226060000391
分别表示左、右相机与2D-LiDAR之间旋转矩阵的逆矩阵;TLS2LC,TLS2RC分别表示左、右相机与2D-LiDAR之间的平移矩阵。
左、右相机坐标系原点在2D-LiDAR坐标系中转换坐标计算如下:
Figure GDA0003541226060000392
公式26中,
Figure GDA0003541226060000393
为左相机坐标系原点在2D-LiDAR坐标系中的表示,
Figure GDA0003541226060000394
为右相机坐标系原点在2D-LiDAR坐标系中的表示,其余符号含义与公式25中相同。
根据上述各点与相机坐标系原点在2D-LiDAR坐标系中相连所形成射线,与 2D-LiDAR坐标中ZLS=h平面的交点。以左相机为例,其交点坐标为:
Figure GDA0003541226060000395
公式27中符号含义,见公式25、26。
进一步的,T5,T6,T7,T8在2D-LiDAR坐标系位置,可通过上述各点取中间点来计算,这里不再赘述。同样OE在2D-LiDAR坐标系位置可以通过两条直线的最小二乘交点确定,这两条直线分别过左右相机的原点与T0点,即:
Figure GDA0003541226060000396
Figure GDA0003541226060000397
公式28中,X,Y,Z表示需要计算的交点坐标值,其余符号见公式25、26,注意T0点坐标与左、右相机原点坐标的区别。
综上,完成了T5,T6,T7,T8与OE在2D-LiDAR坐标系中位置的计算,结合公式 22,即可将这些点转换到旋转云台一体化坐标系中。最终,实现了基于旋转云台一体化坐标系的虚拟单目相机视场的关键点计算,进一步的可以用T5,T6,T7,T8与 OE来计算虚拟单目相机的垂直于水平方向视场角度等,以达到采集末端的一体化处理效果。
(7)旋转图像的空间投影及拼接
由于相机在拍摄时,其相平面与被拍摄平面并非平行,需要将图像投影到被拍摄物体的平面,随后将左右相机的图像进行拼接,以实现在三维模型表面纹理贴图等目的。
①单相机旋转图像的空间投影变换
本节所述相机图像变换的原理,即将相机图像从其各自的相机坐标系,投影到(6)中②所述的虚拟单目相机平面,下面描述其具体方法:
根据(6)-②中方法,计算左、右相机相平面四个角点在2D-LiDAR坐标系的对应坐标:
Figure GDA0003541226060000401
Figure GDA0003541226060000402
以及左、右相机坐标系原点在2D-LiDAR坐标系中的对应空间坐标:LCOLSRCOLS
此时采集末端传感器记录的旋转云台绕Y轴旋转角度α和绕Z轴旋转角度θ,根据公式17,可将左、右相机相平面角点在2D-LiDAR坐标系的位置还原到旋转云台在旋转零角度状态时的对应坐标,记为
Figure GDA0003541226060000403
Figure GDA0003541226060000404
LCOLS-0RCOLS-0
随后计算这些角点和坐标系原点的连线,在还原至旋转零角度后的2D-LiDAR 坐标系Z=h平面的交点位置。以左相机为例,记这些交点为
Figure GDA0003541226060000405
其计算公式如下:
Figure GDA0003541226060000411
公式29中,
Figure GDA0003541226060000412
表示前述
Figure GDA0003541226060000413
点坐标值的XYZ分量,
Figure GDA0003541226060000414
则表示左相机相平面四个角点还原至旋转零角度后的点
Figure GDA0003541226060000415
的XYZ分量,LCOLS-0(X),LCOLS-0(Y),LCOLS-0(Z)则表示左相机坐标系原点在还原至旋转零角度后的LCOLS-0XYZ分量。
获取这些点在2D-LiDAR坐标系中的最大和最小X、Y坐标值,并以其作为透视变换后图像尺寸的变换依据:
Figure GDA0003541226060000416
Figure GDA0003541226060000417
公式30中,
Figure GDA0003541226060000418
表示在
Figure GDA0003541226060000419
四个点中,X分量最大值,
Figure GDA00035412260600004110
表示X分量的最小值,
Figure GDA00035412260600004111
Figure GDA00035412260600004112
分别为Y分量的最大值和最小值。
右相机的最大最小X、Y坐标计算方法类似。设定透视变换后的图像尺寸与原图像尺寸相同,则原图像四个顶点对应的透视变换后图像坐标,以左相机为例,可以按如下公式来计算:
Figure GDA00035412260600004113
公式31中,
Figure GDA0003541226060000422
Figure GDA0003541226060000423
分别表示原左相机图像四个角点在透视变换后图像中的像素坐标,其余符号含义详见公式29、30。
按照上述步骤,得到了四个控制点在原始图像和透视变换后图像中的坐标位置,计算透视变换矩阵,即可得到左右相机图像的透视变换结果(如图12和图 13所示)。
②旋转图像的拼接
进一步的,采集末端上安装左右两个相机,是为了弥补单相机拍摄视场不足,很多场景下需要将左右相机的图像拼接到一起,形成(6)-②类似单目相机的图像。于是需要将左右相机的图像都变换到同一图像上,再对俩变换后的图像进行权重相加形成一张拼接后的图。于是有:
Figure GDA0003541226060000424
Figure GDA0003541226060000425
公式32中,
Figure GDA0003541226060000426
Figure GDA0003541226060000427
表示在拼接图像四个角点与Z=h平面交点在X 轴方向的最大值和最小值,类似的,
Figure GDA0003541226060000428
Figure GDA0003541226060000429
表示Y轴方向的最大值和最小值,其他符号含义如公式30中所述。
设定拼接后图像纵向的像素个数为H=HL,则横向的像素个数为:
Figure GDA00035412260600004210
公式33中,HL表示左相机原始图像在纵向的像素点数,W为新拼接图像在横向的像素点数,其余符号见公式32。
按照与之前同样的方法(见公式31),可以得到左相机图像四个顶点在新拼接图像中的像素坐标:
Figure GDA00035412260600004211
在公式34中,
Figure GDA0003541226060000431
Figure GDA0003541226060000432
表示左相机图像的四个角点在新拼接图像中的像素坐标,
Figure GDA0003541226060000433
Figure GDA0003541226060000434
含义见公式30-,其他符号公式32中含义相同。经过变换后的左相机图像如图14的左上角区域图所示。
右相机的变换原理与左相机一致,采用透视变换计算即可得到左右相机图像的透视变换到拼接图像的结果,如图14的右上角区域图所示。
随后将两张通过加权相加的方法,实现其拼接,最终结果如图14下部区域图所示,左右相机两张图像拼接在一起,左右相机图像重叠部分重合效果较好。
更进一步的,可以将多组在不同旋转角度位置的图像,按照①与②中方法,根据2D-LiDAR获取的激光点云数据所反映的观测场景真实情况,以提取出不同的相机拍摄平面位置,将相机图像投影到这些平面,并进行拼接,最终实现旋转图像在3D模型中的正确贴图,如图17所示,左右相机拍摄的混泥土桥桥底图像能正确贴合在桥底的多平面模型上。
综上所述,本发明实例仅需多自由度采集末端一次观测特制的U型标定板就可以完成相机内参数、2D-LiDAR与相机关系参数以及采集末端各个旋转轴参数的标定,提高了外参数标定的效率和准确性;通过旋转云台的转动,实现了多传感器几何参数的标定,确定了旋转云台的运动学方程,形成了图像与2D激光点云的有效融合。
尽管本发明的实施方案已公开如上,但并不仅仅限于说明书和实施方案中所列运用,它完全可以被适用于各种适合本发明的领域,对于熟悉本领域的人员而言,可容易地实现另外的修改,因此在不背离权利要求及等同范围所限定的一般概念下,本发明并不限于特定的细节和这里所示出与描述的图例。

Claims (7)

1.一种基于一体化标定参数的图像空间投影拼接方法,其特征在于,包括如下步骤:
步骤一、建立一体化标定***;
所述一体化标定***包括标定板和多自由度采集末端,多自由度采集末端包括旋转云台,旋转云台包括2D-LiDAR和至少两个相机;旋转云台固定连接在Y轴旋转轴,Y轴旋转轴轴接连接有Z轴旋转轴;所述标定板包括ChArUco标定板,ChArUco标定板一侧固定有矩形侧板,另一侧固定有等腰直角三角形侧板,等腰直角三角形侧板的一个直角边为ChArUco标定板的侧边;矩形侧板外侧固定有第一边侧板;等腰直角三角形侧板的斜面上固定有第二边侧板;矩形侧板和等腰直角三角形侧板所在平面均与ChArUco标定板所在平面垂直设置;第一边侧板所在平面与ChArUco标定板坐在平面平行,第二边侧板所在平面与等腰直角三角形侧板所在平面垂直设置;
步骤二、获取2D-LiDAR单次扫描标定板的点云数据,提取标定板的角点信息,根据角点信息确定多自由度采集末端在各个旋转自由度时的零点位置;
步骤三、驱动Y轴旋转轴和Z轴旋转轴预设轨迹旋转,并采集相机和2D-LiDAR数据形成完整的标定数据集;
步骤四、提取完整的标定数据集中每张图像中的ArUco特征点和ChArUco角点坐标序列,所述ArUco特征点包括二维码的ID号和二维码四个顶点的平面坐标;根据采集末端不同旋转角度时相机图像中相应二维码四个顶点的平面坐标信息,标定相机内参;
步骤五、在获取2D-LiDAR扫描标定板的点云数据,提取标定板角点,计算各标定板角点在相应相机图像中的平面坐标信息,结合各个标定板角点的空间和平面坐标信息,确定2D-LiDAR和相机之间的旋转矩阵和平移矩阵;
步骤六:获取2D-LiDAR绕多自由度采集末端的Z轴旋转轴和Y轴旋转轴旋转得到的点云数据和相机图像,进行多自由度采集末端的旋转云台的Y轴和Z轴旋转参数的联合标定,具体步骤如下:
根据完整的标定数据集,进行旋转云台绕Z轴和Y轴的旋转参数标定,标定结果包括代表旋转轴方向的单位向量和绕轴旋转点,旋转云台绕Z轴和Y轴的标定结果最后都统一到同一个坐标系2D-LiDAR坐标系下:
首先,设绕Y轴旋转轴旋转的旋转向量就是2D-LiDAR坐标系中Y轴(0,1,0),通过由粗到细的搜索方法求解绕得到Y轴旋转轴的旋转点,同样通过由粗到细的搜索方法得到Z轴旋转轴的旋转点;利用相机图像中的ArUco特征点作为外部控制点,将旋转时相机图像的ArUco特征点,与旋转到零点角度时相机图像的ArUco特征点之间的偏差作为误差函数,利用最小二乘法,迭代求解误差函数最小时的旋转轴方向与旋转点;求解得到旋转云台的Y轴、Z轴的旋转参数,作为迭代求解时Y、Z轴的旋转轴向量和旋转点的初始值;
其中,由粗到细的搜索方法包括如下步骤:
一)获取Z轴旋转轴旋转不同角度的一系列激光点云数据,读取激光的单次点云数据,记录Z轴旋转轴第i次旋转的角度θi;对第i次转台绕Z轴旋转后采集到的原始点云数据,进行稀疏激光均匀化上采样和滤波处理,得到有效的点云数据,提取得到标定板角点p2、p3;最后得到标定板角点p2、p3与Z轴旋转轴旋转角度的序列{LSPi j,θi}j=2,3,其中,i表示第i次旋转,j表示p2或p3点,LSPi j则表示第i次旋转时p2或p3点在2D-LiDAR坐标系中的坐标值;
二)将2D-LiDAR绕多自由度的旋转简化为空间上一点绕两向量旋转,存在如下关系:
Figure FDA0003541226050000021
Figure FDA0003541226050000022
r0Lnx·Lny·(1-cosθ)+cosθ,r1Lnx·Lny·(1-cosθ)+Lnz·sinθ
r2Lnz·Lnx·(1-cosθ)-Lny·sinθ,r3Lnx·Lny·(1-cosθ)-Lnz·sinθ
r4=(Lny)2·(1-cosθ)+cosθ,r5Lny·Lnz·(1-cosθ)+Lnx·sinθ
r6Lnz·Lnx·(1-cosθ)+Lny·sinθ,r7Lny·Lnz·(1-cosθ)-Lnx·sinθ
r8=(Lnz)2·(1-cosθ)+cosθ
(式7)
式7中,θ为Z轴旋转轴旋转的角度,
Figure FDA0003541226050000023
表示Z轴旋转轴旋转的方向向量,Sz表示Z轴旋转轴的旋转点,
Figure FDA0003541226050000024
则表示旋转云台绕
Figure FDA0003541226050000025
方向旋转θ角度对应的旋转矩阵,P1表示2D-LiDAR坐标系内任意一点,P1
Figure FDA0003541226050000026
轴旋转θ角度对应点为P0;r0、r1、r2...r8分别表示旋转矩阵
Figure FDA0003541226050000031
的各元素XAZ,YAZ,ZAZ分别表示绕Z轴旋转点SZ在2D-LiDAR坐标系的X、Y、Z值,LnxLnyLnz分别表示旋转轴向量
Figure FDA0003541226050000032
在2D-LiDAR坐标系的X、Y、Z轴分量;
确定点P0绕Y轴旋转轴
Figure FDA0003541226050000033
轴旋转α后的几何关系:
Figure FDA0003541226050000034
式8中,α表示旋转云台绕Y轴旋转的角度,
Figure FDA0003541226050000035
表示绕Y轴旋转的方向向量,SY则表示绕Y轴旋转点,
Figure FDA0003541226050000036
表示旋转云台绕
Figure FDA0003541226050000037
方向旋转α角度对应的旋转矩阵;P2则表示P0在绕
Figure FDA0003541226050000038
旋转α角度后的对应点;
将绕Z轴旋转轴旋转后的标定板角点p2、p3还原到旋转零角度的公式即为:
Figure FDA0003541226050000039
公式9中,θ0表示在Z轴旋转轴旋转标定为零角度时,传感器记录的实际角度,而θi则表示在第i次Z轴旋转轴旋转时所记录的角度值;
LSPi j(X),LSPi j(Y),LSPi j(Z)分别表示Z轴旋转轴第i次旋转后,标定板第j个角点在2D-LiDAR坐标系中的坐标值,
Figure FDA00035412260500000310
分别表示标定板角点在还原至绕Z轴旋转零角度后,在2D-LiDAR坐标系中的坐标值;
将所有采集末端绕Z轴旋转采集的激光点云数据均还原到绕Z轴旋转零角度,即形成点云数据;
三)设旋Z轴旋转轴向量为
Figure FDA00035412260500000311
计算Z轴旋转轴的旋转点,根据由粗到细的搜索策略,先利用循环迭代的方法进行粗略的求解:
A.初始时,设旋Z轴旋转轴的旋转点为SZ(0,0,0);
B.Z轴旋转轴旋转采集到的所有激光点云数据,均还原到Z轴旋转轴旋转角度为零时的状态;
C.根据还原后的点云数据,计算得到包含还原后的点云数据的最小外接矩阵,记最小外接矩阵的四个顶点为{Lqk}(k=1,2,3,4);
D.将四个顶点连线,形成四条线段{lk}(k=1,2,3,4),随后给定阈值d和s,得到lk线段各自的外扩矩形;
E.依次判断激光点云落在哪个外扩矩形内,并形成每个外扩矩形的激光点集;根据每个矩形内点集拟合直线,并计算矩形内所有点到对应直线的距离累加值;
F.计算四个外扩矩形的顶点到拟合直线距离累加值之和,如果小于上一次迭代的累加值之和,则根据阈值Δy将旋转点SZ的Y分量加上Δy
Δy=1-3mm;
G.若迭代次数超过设定阈值次数则结束,否则回到步骤B运行循环;
迭代过程停止之后,得到粗略求解的采集末端绕Z轴旋转点,记为
Figure FDA0003541226050000041
其中
Figure FDA0003541226050000042
表示
Figure FDA0003541226050000043
点在2D-LiDAR坐标系中的X分量,
Figure FDA0003541226050000044
表示
Figure FDA0003541226050000045
点在2D-LiDAR坐标系中的Y分量,
Figure FDA0003541226050000046
表示
Figure FDA0003541226050000047
点在2D-LiDAR坐标系中的Z分量;
Figure FDA0003541226050000048
作为细化搜索的迭代初值,设定一个绕Z轴旋转轴旋转点Y分量的迭代阈值dy,dy=0.1*Δy,同时将细化搜索时旋转点Y值搜索范围控制在
Figure FDA0003541226050000049
之内;于是按照步骤A至步骤F的迭代过程,细化搜索得到旋转云台绕Z轴旋转轴的旋转点,最终结果得到的旋转点记为
Figure FDA00035412260500000410
其中
Figure FDA00035412260500000411
表示
Figure FDA00035412260500000412
点在2D-LiDAR坐标系中的X分量,
Figure FDA00035412260500000413
表示
Figure FDA00035412260500000414
点在2D-LiDAR坐标系中的Y分量,
Figure FDA00035412260500000415
表示
Figure FDA00035412260500000416
点在2D-LiDAR坐标系中的Z分量;
根据步骤一)至步骤三),同样得到旋转云台绕Y轴的旋转点,记为
Figure FDA00035412260500000417
其中
Figure FDA00035412260500000418
表示
Figure FDA00035412260500000419
点在2D-LiDAR坐标系中的X分量,
Figure FDA00035412260500000420
表示
Figure FDA00035412260500000421
点在2D-LiDAR坐标系中的Y分量,
Figure FDA00035412260500000422
表示
Figure FDA00035412260500000423
点在2D-LiDAR坐标系中的Z分量;完成对Z轴旋转轴和Y轴旋转轴的标零;
a.在对Z轴旋转轴和Y轴旋转轴标零之后,记录下标零时Z轴旋转轴和Y轴旋转轴的传感器角度值,记为θ0,α0,运用图像处理方法提取此时相机图像中ArUco特征点的像素坐标,记为:
Figure FDA0003541226050000051
其中j表示旋转轴在零角度时,相机图像中第j个ArUco二维码,k则表示第j个ArUco二维码的第k个顶点,LC表示左相机图像,RC则表示右相机图像;
Figure FDA0003541226050000052
表示在旋转零角度时,左相机图像中,第j个ArUco二维码的第k个顶点在图像坐标系中的u轴坐标,
Figure FDA0003541226050000053
表示左相机图像中,第j个ArUco二维码的第k个顶点在图像坐标系中的v轴坐标;
Figure FDA0003541226050000054
分别表示右相机图像中相应点的u轴和v轴坐标;
根据Y轴转动的标定数据和Z轴转动的标定数据集,当第i次Z轴旋转轴和Y轴旋转轴旋转角度分别αii,提取相机图像ArUco特征点:
Figure FDA0003541226050000055
i表示第i次旋转;
根据左相机的内参与外参,将左相机的图像坐标转化为标定板世界坐标:
Figure FDA0003541226050000056
式11中,i表示第i次旋转,j表示图像中的第j个ArUco二维码,k则表示二维码的第k个顶点;RW2C与TW2C分别表示左相机的旋转矩阵和平移矩阵,KLC为左相机内参数矩阵;
Figure FDA0003541226050000057
为左相机图像ArUco特征点的像素坐标,
Figure FDA0003541226050000058
分别为特征点在左相机坐标系中X、Y、Z轴的坐标,
Figure FDA0003541226050000059
分别为特征点在标定板世界坐标系中X、Y、Z轴对应坐标值;
因图像中ArUco二维码的四个顶点世界坐标均已知,故求得
Figure FDA0003541226050000061
根据2D-LiDAR与相机间关系标定结果,得到ArUco二维码各顶点在2D-LiDAR坐标系中的空间坐标表示:
Figure FDA0003541226050000062
式12中,RL2LC与TL2LC分别表示2D-LiDAR与左相机之间的旋转矩阵和平移矩阵,
Figure FDA0003541226050000063
表示第i次旋转时,左相机图像的第j个ArUco二维码的第k个顶点在2D-LiDAR坐标系中的对应坐标值;
同理根据右相机的内参与外参,将右相机的图像坐标转化为2D-LiDAR坐标;在2D-LiDAR坐标系下,旋转云台绕Y轴的旋转轴为
Figure FDA0003541226050000064
绕轴旋转点为[XAY,YAY,ZAY]T,将在2D-LiDAR坐标系下的ArUco特征点还原到绕Y轴旋转轴零角度时的坐标,坐标表示为:
Figure FDA0003541226050000065
式13中,α0表示Y轴旋转轴在零角度时,传感器记录的旋转角度值,αi则为第i次旋转时传感器记录的旋转角度值;
Figure FDA0003541226050000066
分别对应着在2D-LiDAR坐标系下相应特征点还原至Y轴旋转轴在零角度时的坐标;
同样的,在2D-LiDAR坐标系下,旋转云台绕Z轴的旋转轴为
Figure FDA0003541226050000067
绕轴旋转点为[XAZ,YAZ,ZAZ]T,则ArUco二维码各顶点还原到绕Z轴旋转零角度,坐标表示为:
Figure FDA0003541226050000071
式14中,θ0表示旋转云台绕Z轴旋转零角度时,传感器记录的旋转角度值,θi则为第i次旋转时传感器记录的旋转角度值;
Figure FDA0003541226050000072
对应着相应特征点还原至2D-LiDAR坐标系下绕Z轴旋转轴在零角度时的坐标;利用相机内参数得到ArUco特征点还原至Z轴旋转轴和Y轴旋转轴均在零角度时的像素坐标,称之为理论像素坐标;具体推导公式如下:
Figure FDA0003541226050000073
式15中,
Figure FDA0003541226050000074
即为ArUco特征点在Z轴旋转轴和Y轴旋转轴均在零角度时的理论像素坐标,
Figure FDA0003541226050000075
为对应的理论2D-LiDAR坐标值;
关注理论像素坐标,与Z轴旋转轴和Y轴旋转轴均在零角度时所提取的ArUco特征点像素坐标,建立误差函数如下:
Figure FDA0003541226050000076
Figure FDA0003541226050000077
Figure FDA0003541226050000081
式16中,ΔLu表示左相机图像在图像坐标系u方向的像素差,ΔLv表示左相机图像在图像坐标系v方向的像素差,ΔRu和ΔRv则表示右相机图像的相应差值,Δu和Δv则反映两个相机在u和v方向的整体像素偏差;根据式16,利用最小二乘法迭代求解旋转云台绕Z轴和Y轴的旋转参数,将Z轴旋转轴和Y轴旋转轴均在零点时的数据作为初始值,即:
Figure FDA0003541226050000082
Figure FDA0003541226050000083
最终解得旋转云台绕Z轴和Y轴的旋转参数,包括旋转轴
Figure FDA0003541226050000084
及其对应的绕轴旋转点SY,SZ
将所有在2D-LiDAR坐标系下的点按照如下公式还原到Z轴旋转轴和Y轴旋转轴均在零角度时的状态:
Figure FDA0003541226050000085
式17中,LX,LY,LZ表示在2D-LiDAR坐标系下表示的任意一点的X、Y、Z轴坐标,α和θ则表示此时传感器记录的旋转云台绕Y轴和Z轴的旋转角度,LX0LY0LZ0则表示将所述任意一点还原至旋转平台绕Y、Z轴旋转零角度时的坐标值;
b.将旋转云台的旋转参数结果转换到旋转云台一体化坐标系内,其中Y轴旋转轴通过基座滑动连接有直线滑轨;
滑轨坐标系Osld与基座坐标系Obase之间旋转和平移均已知:
Figure FDA0003541226050000091
式18中,Mbase2sld(d)表示滑轨坐标系与基座坐标系之间的平移矩阵,Xsld,Ysld,Zsld表示滑轨坐标系中点的X、Y、Z轴的坐标值,而Xbase,Ybase,Zbase代表滑轨坐标系中的点在基座坐标系中的对应位置坐标;
基座坐标系Obase与云台偏转坐标系Oyaw之间旋转和平移均已知:
Figure FDA0003541226050000092
式19中,Ryaw2base和tyaw2base分别表示基座坐标系与云台偏转坐标系之间的旋转和平移矩阵,Xbase,Ybase,Zbase分别表示基座坐标系中点的X、Y、Z轴坐标值,而Xyaw,Yyaw,Zyaw则代表基座坐标系中的点在绕Z轴旋转的偏转坐标系中的对应位置坐标;
云台偏转坐标系Oyaw与旋转云台一体化坐标系Omot之间旋转和平移也已知:
Figure FDA0003541226050000093
式20中,Rmot2yaw和tmot2yaw分别表示偏转坐标系与旋转云台一体化坐标系之间的旋转和平移矩阵,Xyaw,Yyaw,Zyaw表示云台偏转坐标系中点的X、Y、Z轴坐标值,
Figure FDA0003541226050000094
则表示云台偏转坐标系的点在旋转云台Y、Z轴旋转角度均为零角度时,在一体化运动坐标中的对应坐标值;
还原旋转后的旋转云台一体化坐标系中点表示为:
Figure FDA0003541226050000101
式21中,Xmot,Ymot,Zmot表示在旋转云台一体化坐标系中点的X、Y、Z轴坐标值;β为旋转云台绕旋转云台一体化坐标系Z轴旋转的传感器记录角度,β0表示Z轴旋转零角度时传感器记录角度,α与α0分别表示旋转云台一体化坐标系Y轴旋转的相应角度;
Figure FDA0003541226050000102
Figure FDA0003541226050000103
分别表示旋转云台一体化坐标系Y轴和Z轴的旋转矩阵;
在旋转云台一体化坐标系中
Figure FDA0003541226050000104
与所标定旋转云台三个方向的旋转轴
Figure FDA0003541226050000105
完全重合,其中,
Figure FDA0003541226050000106
Figure FDA0003541226050000107
分别为基于2D-LiDAR坐标系的y轴和z轴旋转方向单位向量,且
Figure FDA0003541226050000108
的交点与旋转云台一体化坐标系原点Omot重合,计算
Figure FDA0003541226050000109
最小二乘交点:SAYAZ=[XAYAZ,YAYAZ,ZAYAZ]T,有如下推导:
Figure FDA00035412260500001010
Figure FDA00035412260500001011
式22中,[XLS,YLS,ZLS]分别表示2D-LiDAR坐标系内任一点的X、Y、Z轴坐标值,在旋转云台一体化坐标系中的对应的X、Y、Z轴坐标为[Xmot,Ymot,Zmot];RLS2mot即为2D-LiDAR坐标系何旋转云台一体化坐标系之间的旋转矩阵,tLS2mot为平移矩阵;XAYAZ、YAYAZ和ZAYAZ分别表示
Figure FDA00035412260500001012
最小二乘交点SAYAZ在2D-LiDAR坐标系中的X、Y、Z轴坐标;
综合式18-22,即有2D-LiDAR坐标系到滑轨坐标系的转换公式,表述如下:
Figure FDA0003541226050000111
Figure FDA0003541226050000112
Xsld、Ysld、Zsld分别表示2D-LiDAR坐标系中任意一点[XLS,YLS,ZLS],在还原至旋转云台旋转零角度时,对应滑轨坐标系中点的X、Y、Z轴坐标;
将相机图像坐标中任意点以及2D-LiDAR坐标系中任意点,至滑轨坐标系中的转换为结果坐标;
步骤七:将Z轴旋转轴和Y轴旋转轴标定结果在2D-LiDAR坐标系内表示;
步骤八:多自由度采集末端采集目标物的图像,并将采集到图像的像素点均转化到2D-LiDAR坐标系内,完成目标物图像的拼接。
2.如权利要求1所述的基于一体化标定参数的图像空间投影拼接方法,其特征在于,所述步骤二包括如下步骤:
相机朝向ChArUco标定板,2D-LiDAR朝向与矩形侧板所在平面平行的方向,采集2D-LiDAR扫描标定板的激光点云;对激光点云进行稀疏激光均匀化采样,激光滤波,得到均匀连续的激光点云;将均匀连续的激光点云转换成图像,对图像进行形态学标定板角点取激光点云清晰连续的骨架,获得获取图像中激光点云的轮廓,当轮廓中包含的点的个数大于或等于预设值时,则为有效有效激光轮廓,否则为无效激光轮廓,重新采集2D-LiDAR扫描标定板的激光点云;
根据有效激光轮廓,进行多边形拟合后,提取有效激光轮廓中关于标定板的多边形轮廓线段,通过外扩矩阵确定标定板角点位置:
Figure FDA0003541226050000121
式1中,K1与K2分别表示多边形轮廓上任一线段的起始点和终止点的坐标;
K3和K4分别表示线段K1K2经过外扩后对应的起始点和终止点的坐标;K3和K4作为更新后的标定板角点位置;
然后根据给定的阈值s,计算外扩矩阵四个顶点的坐标,得到获取包含4个顶点的最小外接矩形;s表示此外接矩形高度的二分之一;
通过给定的阈值s,计算中矩形的4个顶点在2D-LiDAR坐标系中的坐标值,并获取包含4个顶点的最小外接矩形,以判断激光点云在哪个矩形内,将所有激光点云分配到对应的矩形内,然后将矩形内的激光点拟合得到直线,根据两两直线的交点得到有效激光轮廓的角点,分别为p0、p1、p2、p3、p4、p5,其中p1、p2、p3、p4即标定板角点;
其中,d1=3δ,s=2δ,δ表示2D-LiDAR的测量误差;
驱动多自由度采集末端的Y轴旋转轴旋转,采集序列激光剖面,并提取处于ChArUco标定板上p2p3连线与2D-LiDAR坐标系中X轴方向的夹角,取夹角最小的位置为Y轴的零点,记录Y轴为零点时对应的Y轴机械角度θ,并驱动Y轴旋转至Y轴的零点;然后驱动Z轴旋转轴旋转,采集序列激光剖面,并提取线段p2p3的长度,取p2p3长度最小时的位置为Z轴的零点,Z轴的零点对应的Z轴机械角度β,将{θ,β}角度作为多自由度采集末端的零点位置,即完成多自由度采集末端的零点位置标定。
3.如权利要求2所述的基于一体化标定参数的图像空间投影拼接方法,其特征在于,所述步骤三包括如下步骤:
将多自由度采集末端的Y轴旋转轴与Z轴旋转轴均置于零点位置,然后分两步进行标定采集相机和2D-LiDAR数据形成完整的标定数据集:
保持Z轴旋转轴处于Z轴的零点位置不动,驱动多自由度采集末端的Y轴旋转轴转动,每次转动1±0.5°;每次Y轴旋转轴转动后,多自由度采集末端上的两个相机拍摄图像并保存,记录2D-LiDAR一帧点云数据,同时记录多自由度采集末端在Y轴自由度和Z轴自由度的旋转角度,作为一组标定数据;在每个旋转位置保存对应的一组标定数据,作为Y轴转动的标定数据集;
然后保持Y轴旋转轴处于Y轴的零点位置不动,驱动多自由度采集末端的Z轴旋转轴转动,每次转动1±0.5°,将转动过程中每组相机图像、激光点云和多自由度采集末端的姿态数据保存下来,作为Z轴转动的标定数据集;
Y轴转动的标定数据和Z轴转动的标定数据集形成完整的标定数据集。
4.如权利要求3所述的基于一体化标定参数的图像空间投影拼接方法,其特征在于,所述步骤四包括如下步骤:
提取Y轴转动的标定数据中每张图像中的ArUco特征点,所述ArUco特征点包括二维码的ID号和二维码四个顶点的平面坐标;同时提取出每张图像的ChArUco角点坐标序列;基于Y轴旋转轴旋转的多个位置,图像坐标系上的角点序列和相对应的世界坐标序列,根据棋盘格标定方法标定出相机的内参数矩阵K和畸变参数矩阵D,完成相机内参的标定:
Figure FDA0003541226050000131
式2中,fx和fy分别表示x与y方向上像素与实际距离的比例;u0,v0分别为图像中心在像素平面即图像坐标系的u、v轴坐标,u为图像宽度方向为方向,v为图像高度方向;k1,k2,k3表示镜头的径向畸变参数;m1,m2表示镜头的切向畸变参数。
5.如权利要求4所述的基于一体化标定参数的图像空间投影拼接方法,其特征在于,所述步骤五包括如下步骤:
根据Y轴转动的标定数据集,获取激光点云数据,提取有效激光轮廓中的标定板角点,即得到标定板角点p1、p2、p3和p4在2D-LiDAR坐标系中的空间坐标{Pi L|i=1,2,3,4};
根据得到的激光轮廓的角点坐标,通过标定板结构的几何约束,设ChArUco标定板处于标定板世界坐标系Z=0平面,于是通过{Pi L|i=1,2,3,4}以及标定板的几何约束得到所述标定板角点在标定板世界坐标系的坐标{Pi W|i=1,2,3,4},即有:
Figure FDA0003541226050000132
Figure FDA0003541226050000141
Figure FDA0003541226050000142
Figure FDA0003541226050000143
式3中,d为二维码边缘与ChArUco标定板边缘的距离,w为ChArUco标定板的边长;
以相机图像中提取的二维码四个顶点作为控制点,结合标定的相机内参数及畸变参数,通过相机标定方法标定出相机外参,即相机坐标系与标定板世界坐标系之间的旋转矩阵RW2C与平移矩阵TW2C,然后确定标定板角点在图像坐标系中的实际坐标:
Figure FDA0003541226050000147
式4中,WXiWYiWZi分别为标定板角点在标定板世界坐标系中的X、Y、Z分量;
Figure FDA0003541226050000144
分别表示标定板角点在相机坐标系中的X、Y、Z分量;
Figure FDA0003541226050000145
分别表示标定板角点在图像中的u轴和v轴坐标值;
得到标定板角点在图像坐标系中的理论坐标信息:
Figure FDA0003541226050000146
式5中,RL2C和TL2C分别表示2D-LiDAR与相机间的旋转矩阵和平移矩阵;LXi,LYi,LZi分别为标定板角点的X、Y、Z分量;CXi,CYi,CZi分别表示标定板角点在相机坐标系中对应的理论坐标值;Iui,Ivi分别表示标定板角点在图像中的理论像素坐标值;
至此,定义损失函数为:
Figure FDA0003541226050000151
式6中,Δu和Δv分别表示像素点在u和v方向的偏差;
采用最小二乘法,迭代求解得到2D-LiDAR与相机间旋转矩阵和平移矩阵。
6.如权利要求1所述的基于一体化标定参数的图像空间投影拼接方法,其特征在于,所述步骤八包括如下步骤:
Ⅰ.将左右相机的视场统一为2D-LiDAR坐标系中的虚拟单目相机模式;
以HL,WL分别表示左相机图像在纵向和横向的像素点个数,HR,WR分别表示右相机图像在横向和纵向的像素点个数,左右相机t0~t4点在相机图像平面坐标系中的坐标分别为:
Figure FDA0003541226050000152
Figure FDA0003541226050000153
其中,
Figure FDA0003541226050000154
分别表示左相机图像中的四个角点,以及这些点在图像坐标系中u、v轴坐标位置,
Figure FDA0003541226050000155
则表示左相机图像中心点在图像坐标系中u、v轴坐标位置;同样的
Figure FDA0003541226050000156
Figure FDA0003541226050000157
则表示右相机图像中相应点在图像坐标系中u、v轴坐标位置;
将图像点转换到相机坐标系中,随后转换到2D-LiDAR坐标系中,按照如下公式求得左右相机图像四个顶点与中心点在2D-LiDAR坐标系内坐标位置:
Figure FDA0003541226050000158
式25中,
Figure FDA0003541226050000161
分别表示前述左相机
Figure FDA0003541226050000162
在图像平面中的u、v轴像素坐标,
Figure FDA0003541226050000163
则分别表示右相机
Figure FDA0003541226050000164
在图像平面中的u、x轴像素坐标;
Figure FDA0003541226050000165
分别表示左、右相机内参的逆矩阵;
Figure FDA0003541226050000166
分别表示左、右相机与2D-LiDAR之间旋转矩阵的逆矩阵;TLS2LC,TLS2RC分别表示左、右相机与2D-LiDAR之间的平移矩阵;
Figure FDA0003541226050000167
表示左相机
Figure FDA0003541226050000168
点在2D-LiDAR坐标系中对应点的X、Y、Z轴坐标,而
Figure FDA0003541226050000169
则表示右相机
Figure FDA00035412260500001610
点在2D-LiDAR坐标系中对应点的X、Y、Z轴坐标;
左相机坐标系原点OLC和右相机坐标系原点ORC,在2D-LiDAR坐标系中对应点LCOLSRCOLS,其转换坐标计算如下:
Figure FDA00035412260500001611
Figure FDA00035412260500001612
式26中,
Figure FDA00035412260500001613
分别为左相机坐标系原点在2D-LiDAR坐标系中的X、Y、Z轴坐标,
Figure FDA00035412260500001614
为右相机坐标系原点在2D-LiDAR坐标系中X、Y、Z轴坐标;根据
Figure FDA00035412260500001615
Figure FDA00035412260500001616
与左、右相机坐标系原点在2D-LiDAR坐标系中相连所形成射线,射线与2D-LiDAR坐标中ZLS=h平面的交点;ZLS表示2D-LiDAR坐标系的Z轴,h表示相机拍摄平面与2D-LiDAR坐标系原点的垂直距离;ZLS=h平面即虚拟单目相机平面;
射线与ZLS=h平面的交点记为:
Figure FDA00035412260500001617
左相机的交点
Figure FDA00035412260500001618
的X、Y与Z轴坐标计算如下:
Figure FDA0003541226050000171
同理得到右相机的相关交点
Figure FDA0003541226050000172
计算得到射线与ZLS=h平面的交点在2D-LiDAR坐标系中具***置后,得到OE,T5,T6,T7,T8在2D-LiDAR坐标系位置;
OE表示左相机图像中心点与左相机坐标系原点OLC连线的直线与右相机图像中心点与右相机坐标系原点ORC连线的直线,在2D-LiDAR坐标系中的最小二乘交点;T5
Figure FDA0003541226050000173
Figure FDA0003541226050000174
连线在2D-LiDAR坐标系YLS=0平面的交点;T6
Figure FDA0003541226050000175
Figure FDA0003541226050000176
的连线中点在2D-LiDAR坐标系YLS=0平面的交点;T7
Figure FDA0003541226050000177
Figure FDA0003541226050000178
的连线在2D-LiDAR坐标系XLS=0平面的交点;T8
Figure FDA0003541226050000179
Figure FDA00035412260500001710
的连线在2D-LiDAR坐标系XLS=0平面的交点;
OE在2D-LiDAR坐标系中位置,通过两条直线的最小二乘交点确定;所述两条直线分别过左右相机的原点与T0点,具体计算方法:
Figure FDA00035412260500001711
Figure FDA00035412260500001712
式28中,X,Y,Z分别表示需要计算的T0在2D-LiDAR坐标X、Y、Z轴的坐标值;完成基于旋转云台一体化坐标系的虚拟单目相机视场的关键点计算;
Ⅱ.旋转图像的空间投影及拼接
2.1.单相机旋转图像的空间投影变换
相机图像从各自的相机坐标系,投影虚拟单目相机平面:
计算左、右相机相平面四个角点在2D-LiDAR坐标系的对应坐标:
Figure FDA0003541226050000181
Figure FDA0003541226050000182
以及左、右相机坐标系原点在2D-LiDAR坐标系中的对应空间坐标:LCOLSRCOLS
此时传感器记录的旋转云台绕Y轴旋转角度α和绕Z轴旋转角度θ,根据式17,将左、右相机相平面角点和左、右相机坐标系原点在2D-LiDAR坐标系的位置还原到旋转云台在旋转零角度状态时的对应坐标,分别记为
Figure FDA0003541226050000183
Figure FDA0003541226050000184
LCOLS-0RCOLS-0
得到左、右相机相平面角点和坐标系原点的连线,在还原至旋转零角度后的2D-LiDAR坐标系Z=h平面的交点位置;对于左相机交点位置分别为:
Figure FDA0003541226050000185
计算公式如下:
Figure FDA0003541226050000186
式29中,
Figure FDA0003541226050000187
表示
Figure FDA0003541226050000188
点坐标值的XYZ分量,
Figure FDA0003541226050000189
则分别表示左相机相平面四个角点还原至旋转零角度后的点
Figure FDA00035412260500001810
的XYZ分量,LCOLS-0(X),LCOLS-0(Y),LCOLS-0(Z)则分别表示左相机坐标系原点在还原至旋转零角度后的LCOLS-0XYZ分量;
左相机相平面四个角点和左相机坐标系原点在还原至旋转零角度后的点在2D-LiDAR坐标系中的最大和最小X、Y坐标值,并作为透视变换后图像尺寸的变换依据:
Figure FDA0003541226050000191
Figure FDA0003541226050000192
式30中,
Figure FDA0003541226050000193
表示在
Figure FDA0003541226050000194
四个点中,X分量最大值,
Figure FDA0003541226050000195
表示X分量的最小值,
Figure FDA0003541226050000196
Figure FDA0003541226050000197
分别为Y分量的最大值和最小值;
同样得到右相机相平面四个角点和右相机坐标系原点在还原至旋转零角度后的点在2D-LiDAR坐标系中的最大和最小X、Y坐标值;
2.2.设定透视变换后的图像尺寸与原图像尺寸相同,则原图像四个顶点对应的透视变换后图像坐标,对于左相机按如下公式进行计算:
Figure FDA0003541226050000198
式31中,
Figure FDA0003541226050000199
Figure FDA00035412260500001910
分别表示原左相机图像四个角点在透视变换后图像中的像素坐标;
根据四个控制点在原始图像和透视变换后图像中的坐标位置,计算透视变换矩阵,即得到左相机图像的透视变换结果;同理可得右相机图像的透视变换结果;
2.3.旋转图像的拼接
将左右相机的图像都变换到同一图像上,再对两个变换后的图像进行权重相加形成一张拼接后的图,于是有:
Figure FDA00035412260500001911
Figure FDA00035412260500001912
式32中,
Figure FDA00035412260500001913
Figure FDA00035412260500001914
分别表示在拼接图像四个角点与Z=h平面交点在X轴方向的最大值和最小值,
Figure FDA00035412260500001915
Figure FDA00035412260500001916
表示拼接图像四个角点与Z=h平面交点在Y轴方向的最大值和最小值;
设定拼接后图像纵向的像素个数为H=HL,则横向的像素个数为:
Figure FDA00035412260500001917
式33中,HL表示左相机原始图像在纵向的像素点数,W为新拼接图像在横向的像素点数;
得到左相机图像四个顶点在新拼接图像中的像素坐标:
Figure FDA0003541226050000202
在式34中,
Figure FDA0003541226050000203
Figure FDA0003541226050000204
分别表示左相机图像的四个角点在新拼接图像中的像素坐标;
同样得到右相机图像的四个角点在新拼接图像中的像素坐标;
将左相机图像和右相机图像通过加权相加的方法拼接得到新拼接图像。
7.如权利要求6所述的基于一体化标定参数的图像空间投影拼接方法,其特征在于,将多组旋转云台在不同旋转角度位置拍摄的图像投影到ZLS=h平面,得到对应的新拼接图像;在世界坐标系内旋转新拼接图像至旋转云台处于零点时的位置,实现旋转图像在3D模型中的还原贴图。
CN202110169365.7A 2021-02-07 2021-02-07 一种基于一体化标定参数的图像空间投影拼接方法 Active CN112927133B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110169365.7A CN112927133B (zh) 2021-02-07 2021-02-07 一种基于一体化标定参数的图像空间投影拼接方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110169365.7A CN112927133B (zh) 2021-02-07 2021-02-07 一种基于一体化标定参数的图像空间投影拼接方法

Publications (2)

Publication Number Publication Date
CN112927133A CN112927133A (zh) 2021-06-08
CN112927133B true CN112927133B (zh) 2022-04-26

Family

ID=76171105

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110169365.7A Active CN112927133B (zh) 2021-02-07 2021-02-07 一种基于一体化标定参数的图像空间投影拼接方法

Country Status (1)

Country Link
CN (1) CN112927133B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114565689B (zh) * 2022-02-28 2024-02-02 燕山大学 一种面向轴对称类三维模型数据压缩重建方法
CN114758005B (zh) * 2022-03-23 2023-03-28 中国科学院自动化研究所 激光雷达与相机外参标定方法及装置
CN118015446A (zh) * 2022-12-07 2024-05-10 南京农业大学 作物姿态采集检测设备及关键点检测方法
CN115984195B (zh) * 2022-12-16 2023-09-29 秦皇岛燕大滨沅科技发展有限公司 一种基于三维点云的车厢轮廓检测方法及***
CN116718113B (zh) * 2023-06-14 2024-07-23 中国科学院、水利部成都山地灾害与环境研究所 泥石流模拟实验堆积体三维形态实时测量***

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2018049818A1 (zh) * 2016-08-16 2018-03-22 上海汇像信息技术有限公司 一种基于三维测量技术测量物体表面积的***及方法
CN109856642A (zh) * 2018-12-20 2019-06-07 上海海事大学 一种旋转三维激光测量***及其平面标定方法
CN111917984A (zh) * 2020-08-13 2020-11-10 上海航天测控通信研究所 一种虚拟云台及控制方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9124873B2 (en) * 2010-12-08 2015-09-01 Cognex Corporation System and method for finding correspondence between cameras in a three-dimensional vision system
TWI502985B (zh) * 2013-04-08 2015-10-01 Omnivision Tech Inc 一種校準360度照相機系統之系統及方法
CN111536902B (zh) * 2020-04-22 2021-03-09 西安交通大学 一种基于双棋盘格的振镜扫描***标定方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2018049818A1 (zh) * 2016-08-16 2018-03-22 上海汇像信息技术有限公司 一种基于三维测量技术测量物体表面积的***及方法
CN109856642A (zh) * 2018-12-20 2019-06-07 上海海事大学 一种旋转三维激光测量***及其平面标定方法
CN111917984A (zh) * 2020-08-13 2020-11-10 上海航天测控通信研究所 一种虚拟云台及控制方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Mono-Camera based Calibration Method for Two-Axes LRF Measurement System;Zhou Langming, et al.;《11th International Conference on Digital Image Processing (ICDIP)》;20191231;1117948-1至1117948-7 *
旋转平台点云数据的配准方法;周朗明等;《测绘学报》;20130215;第42卷(第01期);73-79 *

Also Published As

Publication number Publication date
CN112927133A (zh) 2021-06-08

Similar Documents

Publication Publication Date Title
CN112927133B (zh) 一种基于一体化标定参数的图像空间投影拼接方法
CN110296691B (zh) 融合imu标定的双目立体视觉测量方法与***
CN111473739B (zh) 一种基于视频监控的隧道塌方区围岩变形实时监测方法
CN112396664B (zh) 一种单目摄像机与三维激光雷达联合标定及在线优化方法
Lavest et al. Three-dimensional reconstruction by zooming
CN110842940A (zh) 一种建筑测量机器人多传感器融合三维建模方法及***
JP4245963B2 (ja) 較正物体を用いて複数のカメラを較正するための方法およびシステム
CN105627948B (zh) 一种大型复杂曲面测量***进行复杂曲面采样的方法
CN107492069B (zh) 基于多镜头传感器的图像融合方法
CN104732577B (zh) 一种基于uav低空航测***的建筑物纹理提取方法
CN106990776B (zh) 机器人归航定位方法与***
CN113850126A (zh) 一种基于无人机的目标检测和三维定位方法和***
CN108629829B (zh) 一种球幕相机与深度相机结合的三维建模方法和***
CN102693543B (zh) 室外环境下变焦云台摄像机全自动标定方法
CN106096207B (zh) 一种基于多目视觉的旋翼无人机抗风评估方法及***
CN101354796B (zh) 基于泰勒级数模型的全向立体视觉三维重建方法
CN113205603A (zh) 一种基于旋转台的三维点云拼接重建方法
CN111079565A (zh) 视图二维姿态模板的构建方法及识别方法、定位抓取***
CN113947638B (zh) 鱼眼相机影像正射纠正方法
Nagy et al. Online targetless end-to-end camera-LiDAR self-calibration
CN111060006A (zh) 一种基于三维模型的视点规划方法
CN111189415A (zh) 一种基于线结构光的多功能三维测量重建***及方法
Kang et al. An automatic mosaicking method for building facade texture mapping using a monocular close-range image sequence
CN116563377A (zh) 一种基于半球投影模型的火星岩石测量方法
Dupont et al. An improved calibration technique for coupled single-row telemeter and ccd camera

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