CN113284196B - 一种摄像机畸变逐像素标定方法 - Google Patents

一种摄像机畸变逐像素标定方法 Download PDF

Info

Publication number
CN113284196B
CN113284196B CN202110819872.0A CN202110819872A CN113284196B CN 113284196 B CN113284196 B CN 113284196B CN 202110819872 A CN202110819872 A CN 202110819872A CN 113284196 B CN113284196 B CN 113284196B
Authority
CN
China
Prior art keywords
pixel
picture
calibration plate
camera
mask
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
CN202110819872.0A
Other languages
English (en)
Other versions
CN113284196A (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.)
Hangzhou Xianao Technology Co ltd
Original Assignee
Hangzhou Xianao 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 Hangzhou Xianao Technology Co ltd filed Critical Hangzhou Xianao Technology Co ltd
Priority to CN202110819872.0A priority Critical patent/CN113284196B/zh
Publication of CN113284196A publication Critical patent/CN113284196A/zh
Application granted granted Critical
Publication of CN113284196B publication Critical patent/CN113284196B/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
    • 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

Landscapes

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

Abstract

本发明公开了一种摄像机畸变逐像素标定方法，该方法如下：对拍摄标定板图片、标准标定板图片及其掩膜进行数字图像相关计算，并根据计算结果提取控制点，对摄像机内参和外参进行初始估计；将虚拟标定板和虚拟掩膜板投影到图像上，对投影结果和标定板图片进行数字图像相关计算，以计算得到的映射矩阵中非零值的均方差作为目标函数，以摄像机内参和外参作为优化变量进行优化计算；利用摄像机内参和外参最优值将虚拟标定板和虚拟掩膜板投影到图像上，对投影结果和标定板图片进行数字图像相关计算，取计算得到的映射矩阵中每个元素在各组映射矩阵中非零值的平均值作为每一点的畸变量，得到畸变矩阵。本方法可显著提高摄像机畸变和内参数标定的准确性。

Description

一种摄像机畸变逐像素标定方法
技术领域
本发明属于计算机视觉技术领域，涉及一种摄像机畸变逐像素标定方法。
背景技术
在双目视觉测量及三维重建技术中，摄像机标定的准确性对测量、匹配及重建的准确性有决定性的影响。目前广泛使用小孔成像模型来标定相机内参，对于畸变的标定使用径向和切向畸变模型或其他畸变模型，其特点是用高次方程组来近似描述成像平面不同位置的畸变大小。
摄像机畸变产生的原因是在相机的前方加了透镜，由透镜形状引起的畸变称为径向畸变，光线在远离透镜中心的地方比靠近透镜中心的地方更加弯曲，主要分为两大类：桶形畸变和枕形畸变。如图1所示，从左至右分别为正常图像、桶形畸变和枕形畸变。
径向畸变由公式（2）描述：
Figure 147915DEST_PATH_IMAGE001
(1)
Figure 533897DEST_PATH_IMAGE002
(2)
其中(x,y)为归一化图像平面中某点坐标，(xr,yr)为进行畸变校正后归一化图像平面中对应点坐标，k1、k2、k3为畸变系数，r为点(x,y)与归一化图像平面中心点的距离，如公式（1）所示。
由机械组装过程中透镜和成像平面没有完全平行引起的畸变称为切向畸变，由公式（3）所示：
Figure 782476DEST_PATH_IMAGE003
(3)
其中p1、p2为畸变系数。
径向和切向畸变系数的值，通常依照张正友标定法计算得到。
因为目前常用的畸变模型是摄像机真实畸变关系的近似，使用这种模型进行标定所能达到的精度是十分有限的。如何更为准确的对摄像机进行标定，是亟待解决的问题。
发明内容
本发明的目的在于针对现有技术的不足，提供一种误差小于通常方法的摄像机畸变标定方法，可以对摄像机畸变进行逐像素标定，标定结果以两张映射表的形式给出。
本发明采用的技术方案如下：
准备标定所使用的图像，并根据此图像制作特定的标定板。由公式（4）生成散斑图 案，作为标定所使用的图像。其中I(x,y)表示散斑图案中像素坐标（x,y）处的值，n表示散斑 的数量，k表示散斑序号，D表示散斑的半径（以像素为单位），(xk,yk)为随机数，表示第k个散 斑的位置，
Figure 646527DEST_PATH_IMAGE004
也为随机数，表示第k个散斑的峰值大小。设置散斑图的尺寸为4000×4000像 素，n=15000，D=60。此参数设置不是唯一方案，可有所增减，其原则是散斑图黑、白区域面积 大致相等，印刷在标定板表面并由相机拍摄后（拍摄方法见后文拍摄相机标定所需图片部 分），散斑的细节清晰可见。散斑区域周围留有黑色边缘，如图2所示，黑色边缘的宽度大于 后文所述DIC检测步骤中设置的子区半径。
Figure 929740DEST_PATH_IMAGE005
(4)
用上一步所述方法生成的散斑图片制作特定的标定板。特定的标定板的大小由所标定摄像机的焦距和视场角决定。特定的标定板大小的选择原则是，在用需要标定的摄像机拍摄标定板时，画面中的标定板约占整个画面的2/3大小。制作特定的标定板时，可以使用高精度打印机打印散斑图片，并粘贴在同样大小的正方形玻璃板上。
准备一张分辨率大于待标定相机垂直方向分辨率2倍的标定板标准图片（以下称为“高像素标准标定板图片”），以及一张像素数约等于待标定相机垂直方向分辨率3/4的标定板标准图片（以下称为“标准标定板图片”）。高像素标准标定板图片和标准标定板图片由上一步所述方法生成的散斑图片进行降采样得到。
另外准备两个掩膜图片，即高像素标准标定板图片的掩膜图片和标准标定板图片的掩膜图片。掩膜图片用来在DIC检测中标记感兴趣区域（散斑区域）和其他区域（边框区域），其外观如图3所示。掩膜图片由逻辑0和逻辑1的区域组成。高像素标准标定板图片的掩膜图片大小与高像素标准标定板图片相同，其逻辑1的区域范围与高像素标准标定板图片中的散斑区域范围相同，逻辑0的区域范围与高像素标准标定板图片中的黑色边框区域范围相同。标准标定板图片的掩膜图片大小与标准标定板图片相同，其逻辑1的区域范围与标准标定板图片中的散斑区域范围相同，逻辑0的区域范围与标准标定板图片中的黑色边框区域范围相同。
首先用待标定的摄像机拍摄标定板，得到n张图片（以下称为“拍摄标定板图片”）。拍摄标定板图片的数量n≥3且不是唯一确定，为达到良好的标定效果，n通常可以取15-30。拍摄的时候注意画面中的标定板要占整个画面大约2/3面积。拍摄时对标定板的姿态有所要求。以标定板平面正对摄像机镜头为基准，如图4所示，在笛卡尔坐标系内，标定板在平面x-y内的旋转角α要尽量小，在y-z平面内的旋转角β和z-x平面内的旋转角γ可以任意设置，通常β和γ取0°~30°，如果角度太大可能导致DIC检测失败。拍摄时对标定板姿态的另一要求是，确保相机画面的所有位置都被至少一张图片中的标定板散斑图案覆盖。
以标准标定板图片的掩膜作为感兴趣区域，以标准标定板图片作为参考图片，以拍摄标定板图片作为当前图片，采用数字图像相关法（digital image correlation DIC），为标准标定板图片中每个像素，找到拍摄标定板图片中对应像素相对于此像素的位移（此过程以下简称“DIC检测”）。“对应像素”的定义是，如果对于图片A中位于某个特征的中心位置的像素(u1,v1)，如能在图片B中找到相同特征，则称其中心位置的像素(u2,v2)是图片A中像素(u1,v1)在图片B中的对应像素。将DIC方法中得到的所有像素x方向和y方向的位移分别保存成两个映射矩阵。对所有n张拍摄标定板图片应用此方法，得到n组共2n个映射矩阵（以下称为“标准图片映射矩阵”）。
每组映射矩阵中存储的是每幅拍摄标定板图片中每个像素相对于其在标准标定 板图片中对应像素(xr,yr)的位移
Figure 221044DEST_PATH_IMAGE006
，由此可以计算出标准标定板图片上每个像素在 第k幅拍摄标定板图片上的对应像素
Figure 455454DEST_PATH_IMAGE007
，如下所示。
Figure 123196DEST_PATH_IMAGE008
Figure 526496DEST_PATH_IMAGE009
在映射矩阵中提取s行s列的点阵列，确保s≥2。点阵列的提取方式是在标准标定板图片的散斑区域中等间隔的取s行s列的像素点，确保s≥2。并按照标准图片映射矩阵，取这些像素点在每幅拍摄标定板图片上对应的像素坐标，保存这n×s×s个像素坐标。相应的取标定板所在平面上同样间隔的s×s点阵列，并保存它们在世界坐标系中的坐标。如图5所示，图中黑点为标定板上所取的7×7控制点阵列，实际阵列的行列数量并不固定。在标定板上取的点阵列在标定板上的散斑图案中的位置与在标准标定板图片上取的点阵列在标准标定板图片上的散斑图案中的位置是相同的。保存的拍摄标定板图片上的像素坐标和标定板上的世界坐标以下称为“控制点”。
利用张正友标定法，对摄像机的内参K和n组外参Rk、tk（每幅图片对应标定板相对于相机的旋转和平移）进行初始估计。
根据摄像机内外参数的初始估计，将高像素标准标定板图片和它的掩膜从世界坐标系投影到像素坐标系，投影过程参见公式（9）-公式（12）及其描述。因为有n组Rk和tk，对于每组Rk和tk得到1张投影后的高像素标准标定板图片（以下称为“高像素标准标定板投影图片”）和1张投影后的掩膜（以下称为“掩膜投影图片”），共有n张高像素标准标定板投影图片和n张掩膜投影图片。第k张高像素标准标定板投影图片中的标定板姿态与第k张拍摄标定板图片中的标定板姿态几乎相同，其不同之处在于内外参数初始估计误差以及待标定的畸变量。
对所有n张拍摄标定板图片，n张高像素标准标定板投影图片和n张掩膜投影图片 进行n次DIC计算。将每次检测结果中所有像素x方向和y方向的位移分别保存成两个映射矩 阵
Figure 723122DEST_PATH_IMAGE010
和
Figure 211872DEST_PATH_IMAGE011
，因此得到n组共2n个映射矩阵（以下称为“投影图片映射矩阵”）。
将这2n个映射矩阵中非零值的均方差作为目标函数, 摄像机内外参数K和Rk、tk作为优化变量，摄像机内外参数的初始估计作为初值，用优化方法进行计算，如公式（5）-公式（12）所示。找到使公式（5）达到最小值的摄像机内外参数值，即优化方法计算得到的最优值。
Figure 417725DEST_PATH_IMAGE012
(5)
Figure 409952DEST_PATH_IMAGE013
(6)
Figure 308638DEST_PATH_IMAGE014
(7)
Figure 254991DEST_PATH_IMAGE015
(8)
Figure 998956DEST_PATH_IMAGE016
(9)
Figure 376847DEST_PATH_IMAGE017
(10)
Figure 180855DEST_PATH_IMAGE018
(11)
Figure 378619DEST_PATH_IMAGE019
(12)
Figure 926275DEST_PATH_IMAGE020
(13)
公式（5）-公式（13）中，i、j表示映射矩阵的行列序号，同时i、j也是高像素标准标 定板投影图片中的像素坐标。k表示第k组映射矩阵的序号，同时也是高像素标准标定板投 影图片、掩膜投影图片、拍摄标定板图片的图片序号。ik、jk表示第k组高像素标准标定板投 影图片中的像素坐标。公式（5）-公式（7）中，公式（5）表示2n个映射矩阵中非零值的均方差，
Figure 893094DEST_PATH_IMAGE021
为第k组映射矩阵中的x方向映射矩阵中第i行第j列的值，表示第k幅拍摄标定板图片 中像素(ic,jc)相对于其在第k幅高像素标准标定板投影图片中对应像素(i,j)的x方向位 移，
Figure 133582DEST_PATH_IMAGE022
为第k组映射矩阵中的y方向映射矩阵中第i行第j列的值，表示第k幅拍摄标定板图 片中像素(ic,jc)相对于其在第k幅高像素标准标定板投影图片中对应像素(i,j)的y方向位 移。公式（6），即
Figure 553062DEST_PATH_IMAGE023
表示所有在第i行第j列有非零值的mi,j,x个x方向映射矩阵中
Figure 402944DEST_PATH_IMAGE024
的平 均值。mi,j,x表示所有x方向映射矩阵中在第i行第j列有非零值的矩阵数量。公式（7），即
Figure 489849DEST_PATH_IMAGE025
表示所有在第i行第j列有非零值的mi,j,y个y方向映射矩阵中
Figure 901239DEST_PATH_IMAGE026
的平均值。mi,j,y表示所有 y方向映射矩阵中在第i行第j列有非零值的矩阵数量。
公式（8）中的fdic表示DIC计算，其计算结果为公式（5）-公式（7）中用到的n组共2n 个映射矩阵的元素
Figure 808015DEST_PATH_IMAGE021
、
Figure 697473DEST_PATH_IMAGE027
。公式（8）中i、j表示高像素标准标定板投影图片中的像素坐 标。k表示高像素标准标定板投影图片、掩膜投影图片、拍摄标定板图片的图片序号，
Figure 638884DEST_PATH_IMAGE028
表 示第k幅高像素标准标定板投影图片，
Figure 221176DEST_PATH_IMAGE029
表示第k幅拍摄标定板图片，
Figure 615248DEST_PATH_IMAGE030
表示第k幅掩膜投 影图片。高像素标准标定板投影图片
Figure 573977DEST_PATH_IMAGE028
和掩膜投影图片
Figure 605780DEST_PATH_IMAGE030
是由投影过程，即公式（10）-公 式（13）得到。在进行投影过程之前，需要先根据高像素标准标定板图片
Figure 358972DEST_PATH_IMAGE031
，在世界坐标系中 定义一虚拟标定板BS，并根据高像素标准标定板图片的掩膜
Figure 240341DEST_PATH_IMAGE032
，在世界坐标系中定义一虚 拟掩膜板Bm。设标定板的尺寸为D×D，高像素标准标定板图片及其掩膜图片的像素数为W× W。虚拟标定板BS和虚拟掩膜板Bm都由W×W个坐标点组成的点云描述，坐标点间的距离为D/ W。公式（9）表示虚拟标定板BS中坐标点(u,v,0)的值BS(u,v,0)与高像素标准标定板图片中 像素
Figure 737181DEST_PATH_IMAGE033
的值
Figure 387605DEST_PATH_IMAGE034
相同，同样的，对于虚拟掩膜板Bm，其坐标点(u,v,0) 的值Bm (u,v,0)与高像素标准标定板图片的掩膜中像素
Figure 46120DEST_PATH_IMAGE033
的值
Figure 414784DEST_PATH_IMAGE035
相同。
公式（10）中的
Figure 449736DEST_PATH_IMAGE036
表示第k幅高像素标准标定板投影图片中像素(i,j)的像素 值，BS(u,v,0)表示虚拟标定板中点(u,v,0)的值。
Figure 718781DEST_PATH_IMAGE037
表示第k幅掩膜投影图片中像素 (i,j)的像素值，Bm (u,v,0)表示虚拟掩膜板中点(u,v,0)的值。公式（10）表示高像素标准 标定板投影图片和虚拟标定板以及掩膜投影图片和虚拟掩膜板之间的映射，其中的映射关 系，即第k幅高像素标准标定板投影图片中像素(i,j)和虚拟标定板中点(u,v,0)的对应关 系，以及第k幅掩膜投影图片中像素(i,j)和虚拟掩膜板中点(u,v,0)的对应关系，由公式 （11）-公式（13）得到。
已知第k幅高像素标准标定板投影图片中像素(i,j)和虚拟标定板中点(u,v,0)的 对应关系，根据公式（9）-公式（10），可由高像素标准标定板图片中像素
Figure 813776DEST_PATH_IMAGE033
的值
Figure 669737DEST_PATH_IMAGE034
得到第k幅高像素标准标定板投影图片中像素(i,j)的像素值
Figure 508380DEST_PATH_IMAGE036
。已知 第k幅掩膜投影图片中像素(i,j)和虚拟掩膜板中点(u,v,0)的对应关系，根据公式（9）-公 式（10），可由高像素标准标定板图片的掩膜中像素
Figure 867817DEST_PATH_IMAGE033
的值
Figure 133713DEST_PATH_IMAGE035
得到 第k幅掩膜投影图片中像素(i,j)的像素值
Figure 211391DEST_PATH_IMAGE037
。通常(u,v,0)对应的
Figure 119304DEST_PATH_IMAGE033
为亚像 素坐标，即坐标值为小数的坐标，因此，第k幅高像素标准标定板投影图片中的像素(i,j)的 像素值
Figure 598827DEST_PATH_IMAGE036
由高像素标准标定板图片中像素
Figure 537089DEST_PATH_IMAGE033
近邻的四个整数坐标像素的像 素值进行双线性插值得到，第k幅掩膜投影图片中的像素(i,j)的像素值
Figure 102063DEST_PATH_IMAGE037
由高像素标 准标定板图片的掩膜中像素
Figure 282508DEST_PATH_IMAGE033
近邻的四个整数坐标像素的像素值进行双线性插 值得到。
公式（11）-公式（13）表示第k幅高像素标准标定板投影图片中像素(i,j)和虚拟标 定板中点(u,v,0)间的坐标变换关系。公式（11）表示由图像坐标系变换到像素坐标系，（ik, jk）表示第k幅虚拟标定板投影图片中的像素坐标，K表示相机内参，
Figure 616538DEST_PATH_IMAGE038
表示点（ik,jk） 在图像坐标系中的坐标。公式（12）表示由摄像机坐标系变换到图像坐标系，
Figure 224236DEST_PATH_IMAGE038
表示 摄像机坐标系的坐标（xk,yk,zk）进行归一化得到的像平面坐标，即图像坐标系坐标。公式 （13）表示由世界坐标系变换到摄像机坐标系。Rk、tk表示第k组外参。(u,v,0)表示虚拟标定 板（或虚拟掩膜板）上某点的坐标。因为虚拟标定板（或虚拟掩膜板）平面与世界坐标系的x- y平面重合，所以虚拟标定板（或虚拟掩膜板）平面上的点的第三个维度坐标值为0。（xk,yk, zk）表示 (u,v,0)在摄像机坐标系下的坐标。
找K和Rk、tk的最优值后，根据这些最优值，利用公式（10）-公式（13）将虚拟标定板 和虚拟掩膜板投影到图像上，并进行如公式（8）所示的DIC计算，得到n组投影图片映射矩 阵。根据公式（6）和公式（7），取每一个像素点(i,j)在n组投影图片映射矩阵中的mi,j,x个x方 向非0位移的平均值
Figure 276506DEST_PATH_IMAGE023
和mi,j,y个y方向非0位移的平均值
Figure 729484DEST_PATH_IMAGE025
，即代表计算得到的逐像素的 畸变量，把它们按像素顺序排列即得到畸变矩阵。
畸变矩阵中的元素表示含畸变的图片中每个像素相对于去畸变后的图片中对应像素的位移。在对新拍摄图片进行畸变校正时，按照畸变矩阵，依次为校正后图像的每个像素点查找其在含畸变图像中的对应位置，因为此位置坐标一般是小数，取含畸变图像中此位置近邻的四个整数坐标像素的像素值的双线性插值来作为校正后图片中此像素点的值。
本发明中将数字图像相关计算纳入到优化计算过程中，所用数字图像相关算法可以是FA-GN（Forward Additive Gauss-Newton）方法、IC-GN（Inverse CompositionalGauss-Newton）方法或它们的变体，所述的优化方法包括但不限于数学优化方法（梯度法、Gauss-Newton法、Levenberg-Marquardt方法、线性规划方法等）或进化计算方法（粒子群算法、蚁群算法、遗传算法、模拟退火算法等）或它们的结合；当采用Levenberg-Marquardt方法时，优化计算的效果更好。
本发明的有益效果是：
本发明方法将数字图像相关法应用到优化计算过程中，并将多张拍摄图片中每个像素点和其根据当前内外参估计计算得到的理想像素坐标之间的差异的均方差作为目标函数，脱离了常规畸变模型的限制，对摄像机畸变进行逐点标定，标定结果能够完整描述相机畸变，而不再是一个近似描述，可以有效提高标定精度，从而可有效提高实际应用中视觉测量和三维重建的准确性。经检验，本发明方法的重投影误差达到0.0786，小于张正友标定法的重投影误差0.1085。
附图说明
图1是径向畸变的示意图；
图2是标准标定板图片的示意图；
图3是标准标定板图片的掩膜的示意图；
图4是拍摄标定板图片时的标定板旋转方向示意图；
图5标定板散斑区域的控制点阵列示意图；
图6是本发明方法的流程示意图。
具体实施方式
实施例
如图6是本发明方法的流程示意图，具体方法如下：
准备标定所使用的图像，并根据此图像制作标定板。由公式（4）生成散斑图案，作为标定所使用的图像。设置散斑图的尺寸为4000 ×4000像素，n=15000，D=60。此参数设置不是唯一方案，可有所增减，其原则是散斑图黑、白区域面积大致相等，印刷在标定板表面并由相机拍摄后，散斑的细节清晰可见。散斑区域周围留有500像素宽度的黑色边缘，如图2所示。
标定板实物的大小由所标定摄像机的焦距和视场角决定。标定板大小的选择原则是，在用需要标定的摄像机拍摄标定板时，画面中的标定板约占整个画面的2/3大小。制作尺寸为60mm×60mm的标定板。制作标定板时，使用高精度打印机打印上一步中生成的带边框的散斑图案，并粘贴在同样大小的正方形玻璃板上。
准备一张像素数为2000×2000，纹理与标定板实物完全相同的高精度标准标定板图片，其边框宽度200像素，散斑区域尺寸1600×1600像素。并准备一张像素数为1000×1000，纹理与标定板实物完全相同的标准标定板图片，其边框宽度100像素，散斑区域尺寸800×800像素。
准备高像素标准标定板图片和标准标定板图片的掩膜，如图3所示。其中高像素标准标定板图片的掩膜像素数为2000×2000，边缘200像素宽度的区域为逻辑0，其他区域为逻辑1，标准标定板图片的掩膜像素数为1000×1000，边缘100像素宽度的区域为逻辑0，其他区域为逻辑1。
使用像素数为1920×1080的相机拍摄标定板，得到20张拍摄标定板图片。拍摄标定板图片的数量n不是唯一确定，n≥3，n通常可以取15-30。拍摄的时候注意画面中的标定板要占整个画面大约2/3面积。拍摄时对标定板的姿态有所要求。以标定板平面正对摄像机镜头为基准，如图4所示，在笛卡尔坐标系内，标定板在平面x-y内的旋转角α要尽量小，在y-z平面内的旋转角β和z-x平面内的旋转角γ可以任意设置，通常β和γ取0°~30°，如果角度太大可能导致DIC检测失败。拍摄时对标定板姿态的另一要求是，确保相机画面的所有位置都被至少一张图片中的标定板散斑图案覆盖。
对20张拍摄标定板图片、标准标定板图片及其掩膜图片进行20次DIC检测。设置DIC检测的参数，子区半径设为40，子区间隔设为0，收敛精度设为1e-6，最大迭代次数设为50，此参数设置不是唯一方案，可根据实际情况进行调整。DIC检测为标准标定板图片中每个像素，找到拍摄标定板图片中对应像素相对于此像素的位移。将这些像素x方向和y方向的位移分别保存成两个映射矩阵。对所有20张拍摄标定板图片应用此方法，得到20组共40个映射矩阵。
每组映射矩阵中存储的是每幅拍摄标定板图片中每个像素相对于其在标准标定板图片中对应像素的位移，由此可以计算出标准标定板图片上每个像素在第k幅拍摄标定板图片上的对应像素的坐标，如公式（3）、公式（4）所示。
在映射矩阵中提取7行7列的点阵列。点阵列的提取方式是在标准标定板图片的散斑区域中等间隔的取像素点，在标准标定板图片的第[200:100:800]行和第[200:100:800]列提取7×7个控制点，并根据映射矩阵计算出标准标定板图片上的这7×7个控制点在第k幅拍摄标定板图片上的对应像素的坐标，保存这20×7×7个像素坐标。相应的在标定板实物上取间距6mm的7×7点阵列，并保存它们在世界坐标系中的坐标，具体取点位置如图5所示。保存的拍摄标定板图片上的像素坐标和标定板上的世界坐标以下称为“控制点”。
根据像素坐标系和世界坐标系中的控制点位置，利用张正友标定法，对摄像机的内参K和20组外参Rk、tk（每幅图片对应标定板相对于相机的旋转和平移）进行初始估计，得到一个内参矩阵K和20组外参矩阵Rk、tk。
根据摄像机内外参矩阵的初始估计，将高像素标准标定板图片和它的掩膜从世界坐标系投影到像素坐标系，投影过程参见公式（9）-公式（13）及其描述。因为有20组Rk和tk，对于每组Rk和tk得到1张投影后的高像素标准标定板图片和1张投影后的掩膜，共有20张高像素标准标定板投影图片和20张掩膜投影图片。第k张高像素标准标定板投影图片中的标定板姿态与第k张拍摄标定板图片中的标定板姿态几乎相同，其不同之处在于内外参数初始估计误差以及待标定的畸变量。
对所有20张拍摄标定板图片，20张高像素标准标定板投影图片和20张掩膜投影图片进行20次DIC检测。将每次检测结果中所有像素x方向和y方向的位移分别保存成两个映射矩阵，因此得到n组共2n个映射矩阵。
将这2n个映射矩阵中非零值的均方差作为目标函数，如公式（5）-公式（13）所示,摄像机内外参数K和Rk、tk作为优化变量，摄像机内外参数的初始估计作为初值，用列文伯格-马夸尔特优化方法进行计算。找到使（5）达到最小值的摄像机内外参数值，即优化计算得到的最优值。
找到K和Rk、tk的最优值后，根据这些最优值，利用公式（9）-公式（13）得到高像素标准标定板图片及其掩膜图片的投影，并进行如公式（8）所示的DIC计算，得到n组投影图片映射矩阵。根据公式（6）和公式（7）即代表计算得到的逐像素的畸变量，把它们按像素顺序排列即得到畸变矩阵。畸变矩阵中的元素表示含畸变的图片中每个像素相对于去畸变后的图片中对应像素的位移。在对新拍摄图片进行畸变校正时，按照畸变矩阵，依次为校正后图像的每个像素点查找其在含畸变图像中的对应位置，因为此位置坐标一般是小数，取含畸变图像中此位置近邻的四个整数坐标像素的像素值的双线性插值来作为校正后图片中此像素点的值。
表1列出了使用本发明的摄像机畸变逐像素标定方法和张正友标定法利用不同组图片重复十次标定实验得到的数据，每组图片数量为20张。在张正友标定法中，用于特征点提取的标定板包括棋盘格、三角、圆点和散斑标定板。在表1中，第一列为不同标定方法的重投影误差平均值，和其他对照实验的结果相比，本发明的摄像机畸变逐像素标定方法重投影误差平均值最小。后面四列为标定结果中对焦距和主点位置估计值的均方根误差，表示参数估计的稳定性。重复多次标定中每次估计的参数值差异更小说明更稳定，可以看出和其他对照实验的结果相比，本发明的方法参数估计是最稳定的。
表1 标定误差评价
Figure 682134DEST_PATH_IMAGE039

Claims (8)

1.一种摄像机畸变逐像素标定方法，其特征在于，包括如下：
准备特定的标定板、高像素标准标定板图片及其掩膜、标准标定板图片及其掩膜；利用待标定摄像机拍摄所述标定板，得到n张拍摄标定板图片，其中n≥3；
以标准标定板图片的掩膜作为感兴趣区域，以标准标定板图片作为参考图片，以拍摄标定板图片作为当前图片，采用数字图像相关法，为标准标定板图片中每个像素找到拍摄标定板图片中对应点相对于此像素的位移，将所有像素x方向和y方向的位移分别保存成两个映射矩阵；即对所有n张拍摄标定板图片得到n组共2n个标准图片映射矩阵；
提取控制点，利用张正友标定法对摄像机的内参和外参进行初始估计；得到一组内参及n组外参；
根据得到的摄像机内外参数的初始估计，将高像素标准标定板图片及其掩膜从世界坐标系投影到像素坐标系，得到n张标准标定板投影图片及其掩膜；
以标准标定板投影图片的掩膜作为感兴趣区域，以标准标定板投影图片作为参考图片，以拍摄标定板图片作为当前图片，用数字图像相关法为标准标定板投影图片中每个像素找到拍摄标定板图片中对应点相对于此像素的位移，将所有像素x方向和y方向的位移分别保存成两个映射矩阵，即对所有n张拍摄标定板图片得到n组共2n个投影图片映射矩阵；
将这2n个投影图片映射矩阵的均方差作为目标函数,摄像机内外参数作为优化变量，摄像机内外参数的初始估计作为初值，用优化方法进行计算，找到使所述目标函数达到最小值的摄像机内外参数值，即最优值；
找到最优值后，在对应的n组投影图片映射矩阵中，取x和y方向位移各自的平均值，得到畸变矩阵，其中的元素表示含畸变的图片中每个像素相对于去畸变后的图片中对应像素的位移；
所述的提取控制点，具体方法如下：在映射矩阵中提取s行s列的点阵列，确保s≥2；点阵列的提取方式是在标准标定板图片的散斑区域中等间隔的取s行s列的像素点，确保s≥2；并按照标准图片映射矩阵，取这些像素点在每幅拍摄标定板图片上对应的像素坐标，保存这n×s×s个像素坐标；相应的取标定板所在平面上同样间隔的s×s点阵列，并保存它们在世界坐标系中的坐标；保存的拍摄标定板图片上的像素坐标和标定板上的世界坐标称为控制点。
2.根据权利要求1所述的摄像机畸变逐像素标定方法，其特征在于，所述的特定的标定板的图案为散斑图案，高像素标准标定板图片、标准标定板图片为散斑图片。
3.根据权利要求1所述的摄像机畸变逐像素标定方法，其特征在于，所述的n为15-30。
4.根据权利要求1所述的摄像机畸变逐像素标定方法，其特征在于，所述的目标函数具体为：
Figure FDA0003255365170000021
Figure FDA0003255365170000022
Figure FDA0003255365170000023
Figure FDA0003255365170000024
Figure FDA0003255365170000025
Figure FDA0003255365170000026
Figure FDA0003255365170000027
Figure FDA0003255365170000031
Figure FDA0003255365170000032
公式(1)-公式(9)中，i、j表示映射矩阵的行、列序号；k表示第k组映射矩阵的序号；ik、jk表示第k组高像素标准标定板投影图片中的像素坐标；公式(1)-公式(3)中，公式(1)表示2n个映射矩阵中非零值的均方差，
Figure FDA0003255365170000033
为第k组映射矩阵中的x方向映射矩阵中第i行第j列的值，表示第k幅拍摄标定板图片中像素(ic,jc)相对于其在第k幅高像素标准标定板投影图片中对应像素(i,j)的x方向位移，
Figure FDA0003255365170000034
为第k组映射矩阵中的y方向映射矩阵中第i行第j列的值，表示第k幅拍摄标定板图片中像素(ic,jc)相对于其在第k幅高像素标准标定板投影图片中对应像素(i,j)的y方向位移；公式(2)，即
Figure FDA0003255365170000035
表示所有在第i行第j列有非零值的mi,j,x个x方向映射矩阵中
Figure FDA0003255365170000036
的平均值；mi,j,x表示所有x方向映射矩阵中在第i行第j列有非零值的矩阵数量；公式(3)，即
Figure FDA0003255365170000037
表示所有在第i行第j列有非零值的mi,j,y个y方向映射矩阵中
Figure FDA0003255365170000038
的平均值；mi,j,y表示所有y方向映射矩阵中在第i行第j列有非零值的矩阵数量；
公式(4)中的fdic表示DIC计算，其计算结果为公式(1)-公式(3)中用到的n组共2n个映射矩阵的元素
Figure FDA0003255365170000039
公式(4)中
Figure FDA00032553651700000310
表示第k幅高像素标准标定板投影图片，
Figure FDA00032553651700000311
表示第k幅拍摄标定板图片，
Figure FDA00032553651700000312
表示第k幅掩膜投影图片；高像素标准标定板投影图片
Figure FDA0003255365170000041
和掩膜投影图片
Figure FDA0003255365170000042
是由投影过程，即公式(6)-公式(9)得到；在进行投影过程之前，需要先根据高像素标准标定板图片
Figure FDA0003255365170000043
在世界坐标系中定义一虚拟标定板BS，并根据高像素标准标定板图片的掩膜
Figure FDA0003255365170000044
在世界坐标系中定义一虚拟掩膜板Bm；设标定板的尺寸为D×D，高像素标准标定板图片及其掩膜图片的像素数为W×W；虚拟标定板BS和虚拟掩膜板Bm都由W×W个坐标点组成的点云描述，坐标点间的距离为D/W；公式(5)表示虚拟标定板BS中坐标点(u,v,0)的值BS(u,v,0)与高像素标准标定板图片中像素
Figure FDA0003255365170000045
的值
Figure FDA0003255365170000046
相同，同样的，对于虚拟掩膜板Bm，其坐标点(u,v,0)的值Bm(u,v,0)与高像素标准标定板图片的掩膜中像素
Figure FDA0003255365170000047
的值
Figure FDA0003255365170000048
相同；
公式(6)中的
Figure FDA0003255365170000049
表示第k幅高像素标准标定板投影图片中像素(i,j)的像素值，BS(u,v,0)表示虚拟标定板中点(u,v,0)的值；
Figure FDA00032553651700000410
表示第k幅掩膜投影图片中像素(i,j)的像素值，Bm(u,v,0)表示虚拟掩膜板中点(u,v,0)的值；公式(6)表示高像素标准标定板投影图片和虚拟标定板以及掩膜投影图片和虚拟掩膜板之间的映射；
已知第k幅高像素标准标定板投影图片中像素(i,j)和虚拟标定板中点(u,v,0)的对应关系，根据公式(5)-公式(6)，可由高像素标准标定板图片中像素
Figure FDA00032553651700000411
的值
Figure FDA00032553651700000412
得到第k幅高像素标准标定板投影图片中像素(i,j)的像素值
Figure FDA00032553651700000413
已知第k幅掩膜投影图片中像素(i,j)和虚拟掩膜板中点(u,v,0)的对应关系，根据公式(5)-公式(6)，可由高像素标准标定板图片的掩膜中像素
Figure FDA0003255365170000051
的值
Figure FDA0003255365170000052
得到第k幅掩膜投影图片中像素(i,j)的像素值
Figure FDA0003255365170000053
第k幅高像素标准标定板投影图片中的像素(i,j)的像素值
Figure FDA0003255365170000054
由高像素标准标定板图片中像素
Figure FDA0003255365170000055
近邻的四个整数坐标像素的像素值进行双线性插值得到，第k幅掩膜投影图片中的像素(i,j)的像素值
Figure FDA0003255365170000056
由高像素标准标定板图片的掩膜中像素
Figure FDA0003255365170000057
近邻的四个整数坐标像素的像素值进行双线性插值得到；
公式(7)-公式(9)表示第k幅高像素标准标定板投影图片中像素(i,j)和虚拟标定板中点(u,v,0)间的坐标变换关系；公式(7)表示由图像坐标系变换到像素坐标系(ik,jk)，表示第k幅虚拟标定板投影图片中的像素坐标，K表示相机内参，
Figure FDA0003255365170000058
表示点(ik,jk)在图像坐标系中的坐标；公式(8)表示由摄像机坐标系变换到图像坐标系，
Figure FDA0003255365170000059
表示摄像机坐标系的坐标(xk,yk,zk)进行归一化得到的像平面坐标，即图像坐标系坐标；公式(9)表示由世界坐标系变换到摄像机坐标系；Rk、tk表示第k组外参；(u,v,0)表示虚拟标定板或虚拟掩膜板上某点的坐标；(xk,yk,zk)表示(u,v,0)在摄像机坐标系下的坐标；
在此优化过程中，摄像机内外参数K和Rk、tk作为优化变量。
5.根据权利要求4所述的摄像机畸变逐像素标定方法，其特征在于，在生成畸变矩阵时，根据K和Rk、tk的最优值，利用公式(6)-公式(9)将虚拟标定板和虚拟掩膜板投影到图像上，并进行如公式(4)所示的DIC计算，得到n组投影图片映射矩阵；根据公式(2)和公式(3)，取每一个像素点(i,j)在n组投影图片映射矩阵中的mi,j,x个x方向非0位移的平均值
Figure FDA00032553651700000510
和mi,j,y个y方向非0位移的平均值
Figure FDA00032553651700000511
即代表计算得到的逐像素的畸变量，把它们按像素顺序排列即得到畸变矩阵。
6.根据权利要求1所述的摄像机畸变逐像素标定方法，其特征在于，所述的数字图像相关法是Forward Additive Gauss-Newton方法、Inverse CompositionalGauss-Newton方法或它们的变体。
7.根据权利要求1所述的摄像机畸变逐像素标定方法，其特征在于，所述的优化方法采用数学优化方法或进化计算方法或它们的结合，所述的数学优化方法为梯度法、Gauss-Newton法、Levenberg-Marquardt方法或线性规划方法，所述进化计算方法为粒子群算法、蚁群算法、遗传算法或模拟退火算法。
8.根据权利要求1所述的摄像机畸变逐像素标定方法，其特征在于，在对新拍摄图片进行畸变校正时，按照畸变矩阵，依次为校正后图像的每个像素点查找其在含畸变图像中的对应像素点的位置，取含畸变图像中此位置近邻的四个整数坐标像素的像素值的双线性插值来作为校正后图片中此像素点的值。
CN202110819872.0A 2021-07-20 2021-07-20 一种摄像机畸变逐像素标定方法 Active CN113284196B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110819872.0A CN113284196B (zh) 2021-07-20 2021-07-20 一种摄像机畸变逐像素标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110819872.0A CN113284196B (zh) 2021-07-20 2021-07-20 一种摄像机畸变逐像素标定方法

Publications (2)

Publication Number Publication Date
CN113284196A CN113284196A (zh) 2021-08-20
CN113284196B true CN113284196B (zh) 2021-10-22

Family

ID=77286774

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110819872.0A Active CN113284196B (zh) 2021-07-20 2021-07-20 一种摄像机畸变逐像素标定方法

Country Status (1)

Country Link
CN (1) CN113284196B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117576228B (zh) * 2024-01-16 2024-04-16 成都合能创越软件有限公司 基于实时场景的相机坐标标定方法及***

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106803273A (zh) * 2017-01-17 2017-06-06 湖南优象科技有限公司 一种全景摄像机标定方法
CN109272570A (zh) * 2018-08-16 2019-01-25 合肥工业大学 一种基于立体视觉数学模型的空间点三维坐标求解方法
CN109598762A (zh) * 2018-11-26 2019-04-09 江苏科技大学 一种高精度双目相机标定方法
CN110517202A (zh) * 2019-08-30 2019-11-29 的卢技术有限公司 一种车身摄像头标定方法及其标定装置
CN112258588A (zh) * 2020-11-13 2021-01-22 江苏科技大学 一种双目相机的标定方法、***及存储介质

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102750697B (zh) * 2012-06-08 2014-08-20 华为技术有限公司 一种参数标定方法及装置
EP3200148B1 (en) * 2014-10-31 2019-08-28 Huawei Technologies Co., Ltd. Image processing method and device
EP3742399A1 (en) * 2019-05-21 2020-11-25 Sportlogiq Inc. Systems and methods for image registration and camera calibration using learned error functions

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106803273A (zh) * 2017-01-17 2017-06-06 湖南优象科技有限公司 一种全景摄像机标定方法
CN109272570A (zh) * 2018-08-16 2019-01-25 合肥工业大学 一种基于立体视觉数学模型的空间点三维坐标求解方法
CN109598762A (zh) * 2018-11-26 2019-04-09 江苏科技大学 一种高精度双目相机标定方法
CN110517202A (zh) * 2019-08-30 2019-11-29 的卢技术有限公司 一种车身摄像头标定方法及其标定装置
CN112258588A (zh) * 2020-11-13 2021-01-22 江苏科技大学 一种双目相机的标定方法、***及存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A New Solution for Camera Calibration and;Rui Melo等;《IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING》;20120217;第634-643页 *
面成型3D打印工艺中投影图像畸变校正技术研究;贾红帅;《中国优秀硕士学位论文全文数据库 信息科技辑》;20170115;第1-42页 *

Also Published As

Publication number Publication date
CN113284196A (zh) 2021-08-20

Similar Documents

Publication Publication Date Title
CN109598762B (zh) 一种高精度双目相机标定方法
CN106780388B (zh) 一种线阵相机光学畸变矫正方法
CN109272570B (zh) 一种基于立体视觉数学模型的空间点三维坐标求解方法
CN112013792B (zh) 一种复杂大构件机器人面扫描三维重建方法
CN111351446B (zh) 一种用于三维形貌测量的光场相机校准方法
CN109272574B (zh) 基于投影变换的线阵旋转扫描相机成像模型构建方法和标定方法
CN115861445B (zh) 一种基于标定板三维点云的手眼标定方法
CN109827521B (zh) 一种快速多线结构光视觉测量***标定方法
CN112258587B (zh) 一种基于灰狼粒子群混合算法的相机标定方法
CN112614075B (zh) 一种面结构光3d***的畸变校正方法及设备
CN111047649A (zh) 一种基于最优偏振角的相机高精度标定方法
CN113920205B (zh) 一种非同轴相机的标定方法
CN112258588A (zh) 一种双目相机的标定方法、***及存储介质
CN110349257B (zh) 一种基于相位伪映射的双目测量缺失点云插补方法
CN112947885B (zh) 一种曲面屏展平图像的生成方法及装置
CN113298886B (zh) 一种投影仪的标定方法
CN114283203A (zh) 一种多相机***的标定方法及***
CN115457147A (zh) 相机标定方法、电子设备及存储介质
CN113284196B (zh) 一种摄像机畸变逐像素标定方法
CN113920206A (zh) 透视移轴相机的标定方法
CN113610929B (zh) 一种相机与多线激光的联合标定方法
CN110751692B (zh) 一种摄像机成像误差标定方法与矫正方法
CN115100078B (zh) 一种曲面屏图像中点阵坐标修正与填补的方法及相关装置
CN113865514B (zh) 一种线结构光三维测量***标定方法
CN116277979B (zh) 一种dlp打印机光机畸变矫正方法

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