目录

• 1 相机为什么要标定？
• 2 反解相机矩阵
• 3 反解畸变系数
• 4 去畸变原理与实战

1 相机为什么要标定？

直接给结论：相机标定的目的是为了获得相机矩阵和畸变系数。这部分内容可以参考计算机视觉教程0-3：为何拍照会有死亡视角？详解相机矩阵与畸变。获得了相机矩阵和畸变系数后就能进行更高级的应用，主要应用于立体视觉领域，例如视觉SLAM、三维导航、双目测距等，总之，相机标定相当重要。

相机标定的经典方法是张正友标定法，需要下面的张正友标定板。本节将详细地从底层数学原理推导张正友标定法。

2 反解相机矩阵

不考虑畸变，对图像中标定板角点进行检测，得到角点像素坐标值[ u v 1 ] T \left[ \begin{matrix} u& v& 1\\\end{matrix} \right] ^T [u​v​1​]T，根据已知棋盘格大小，计算世界坐标系下标定板角点的世界坐标值[ W  ⁣  ⁣ X W  ⁣  ⁣ ⁣  ⁣  Y W  ⁣  ⁣ Z 1 ] T \left[ \begin{matrix} ^{\boldsymbol{W}\!}\!X& ^{\boldsymbol{W}\!}\!\!\:\!\:Y& ^{\boldsymbol{W}\!}\!Z& 1\\\end{matrix} \right] ^T [WX​WY​WZ​1​]T，此时从世界坐标到像素坐标的映射为

[ u ^ v ^ 1 ] = K[ W C    ⁣ ⁣ ⁣ R C  ⁣ pw0] [ W  ⁣  ⁣ X W  ⁣  ⁣ ⁣  ⁣  Y W  ⁣  ⁣ Z 1 ] \left[ \begin{array}{c} \hat{u}\\ \hat{v}\\ 1\\\end{array} \right] =\boldsymbol{K}\left[ \begin{matrix} _{\boldsymbol{W}}^{\boldsymbol{C}}\;\!\!\!\boldsymbol{R}& ^{\boldsymbol{C}}\!\boldsymbol{p}_{w_0}\\\end{matrix} \right] \left[ \begin{array}{c} ^{\boldsymbol{W}\!}\!X\\ ^{\boldsymbol{W}\!}\!\!\:\!\:Y\\ ^{\boldsymbol{W}\!}\!Z\\ 1\\\end{array} \right] ⎣⎡​u^v^1​⎦⎤​=K[WC​R​Cpw0​​​]⎣⎢⎢⎡​WXWYWZ1​⎦⎥⎥⎤​

通常将世界坐标系与标定板固连使 W  ⁣  ⁣ Z = 0 ^{\boldsymbol{W}\!}\!Z=0 WZ=0，如图所示

此时映射关系简化为

[ u ^ v ^ 1 ] = K[ r 1 r 2 p ] [ W  ⁣  ⁣ X W  ⁣  ⁣ ⁣  ⁣  Y 1 ] \left[ \begin{array}{c} \hat{u}\\ \hat{v}\\ 1\\\end{array} \right] =\boldsymbol{K}\left[ \begin{matrix} \boldsymbol{r}_1& \boldsymbol{r}_2& \boldsymbol{p}\\\end{matrix} \right] \left[ \begin{array}{c} ^{\boldsymbol{W}\!}\!X\\ ^{\boldsymbol{W}\!}\!\!\:\!\:Y\\ 1\\\end{array} \right] ⎣⎡​u^v^1​⎦⎤​=K[r1​​r2​​p​]⎣⎡​WXWY1​⎦⎤​

其中 r 1 \boldsymbol{r}_1 r1​、 r 2 \boldsymbol{r}_2 r2​是旋转矩阵 W C    ⁣ ⁣ ⁣ R _{\boldsymbol{W}}^{\boldsymbol{C}}\;\!\!\!\boldsymbol{R} WC​R的前两列。

设从标定板到成像平面的单应性矩阵为

H =[ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] =[ h 1 h 2 h 3 ] \boldsymbol{H}=\left[ \begin{matrix} h_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\\\end{matrix} \right] =\left[ \begin{matrix} \boldsymbol{h}_1& \boldsymbol{h}_2& \boldsymbol{h}_3\\\end{matrix} \right] H=⎣⎡​h11​h21​h31​​h12​h22​h32​​h13​h23​h33​​⎦⎤​=[h1​​h2​​h3​​]

这部分内容可以参考计算机视觉教程1-2：单应性矩阵估计

下面基于已知的$\mathbf{H}$反解内参矩阵$\mathbf{K}$，代入 h i = Kr i ( i = 1 , 2 ) \boldsymbol{h}_i=\boldsymbol{Kr}_i\left( i=1,2 \right) hi​=Kri​(i=1,2)可得

{ h 1 T Q h 2 = 0 h 1 T Q h 1 = h 2 T Q h 2 = 1 \begin{cases} \boldsymbol{h}_{1}^{T}\boldsymbol{Qh}_2=0\\ \boldsymbol{h}_{1}^{T}\boldsymbol{Qh}_1=\boldsymbol{h}_{2}^{T}\boldsymbol{Qh}_2=1\\\end{cases} {h1T​Qh2​=0h1T​Qh1​=h2T​Qh2​=1​

其中对称矩阵 Q = K−T K−1 \boldsymbol{Q}=\boldsymbol{K}^{-T}\boldsymbol{K}^{-1} Q=K−TK−1

考虑到

hiT Q h j =[ h i 1 h i 2 h i 3 ] [ q 1 q 2 q 3 q 2 q 4 q 5 q 3 q 5 q 6 ] [ h j 1 h j 2 h j 3 ]=v i j T q \boldsymbol{h}_{i}^{T}\boldsymbol{Qh}_j=\left[ \begin{matrix} h_{i1}& h_{i2}& h_{i3}\\\end{matrix} \right] \left[ \begin{matrix} q_1& q_2& q_3\\ q_2& q_4& q_5\\ q_3& q_5& q_6\\\end{matrix} \right] \left[ \begin{array}{c} h_{j1}\\ h_{j2}\\ h_{j3}\\\end{array} \right] \\=\boldsymbol{v}_{ij}^{T}\boldsymbol{q} hiT​Qhj​=[hi1​​hi2​​hi3​​]⎣⎡​q1​q2​q3​​q2​q4​q5​​q3​q5​q6​​⎦⎤​⎣⎡​hj1​hj2​hj3​​⎦⎤​=vijT​q

可以改写为

[ v 12 T v 11 T − v 22 T ] q = 0 \left[ \begin{array}{c} \boldsymbol{v}_{12}^{T}\\ \boldsymbol{v}_{11}^{T}-\boldsymbol{v}_{22}^{T}\\\end{array} \right] \boldsymbol{q}=0 [v12T​v11T​−v22T​​]q=0

其中矩阵$\mathbf{V}$完全由已知的$\mathbf{H}$导出，$\mathbf{q}$为$\mathbf{Q}$的向量形式。由于一个$\mathbf{H}$最多只能贡献两个线性方程，因此确定内参矩阵至少要3张标定板(工程上一般取15~20张为宜)，这些方程共同构成$\mathbf{V}\mathbf{q}=0$，可用计算机视觉教程1-2：单应性矩阵估计的基本DLT算法解算，解出来就是相机内参矩阵。

3 反解畸变系数

张氏标定法仅考虑畸变中影响较大的径向畸变而忽略切向畸变，用求得的内参矩阵反解畸变参数。

设无畸变像素坐标、成像平面坐标分别为[ u ^ v ^ ] T \left[ \begin{matrix} \hat{u}& \hat{v}\\\end{matrix} \right] ^T [u^​v^​]T、[ x ^ y ^ ] T \left[ \begin{matrix} \hat{x}& \hat{y}\\\end{matrix} \right] ^T [x^​y^​​]T；有畸变像素坐标、成像平面坐标分别为[ u v ] T \left[ \begin{matrix} u& v\\\end{matrix} \right] ^T [u​v​]T、[ x y ] T \left[ \begin{matrix} x& y\\\end{matrix} \right] ^T [x​y​]T；设内参$s=0$，则

{ u = f u x + c u v = f v y + c v    { u ^ = f u x ^ + c u v ^ = f v y ^ + c v \begin{cases} u=f_ux+c_u\\ v=f_vy+c_v\\\end{cases}\,\, \begin{cases} \hat{u}=f_u\hat{x}+c_u\\ \hat{v}=f_v\hat{y}+c_v\\\end{cases} {u=fu​x+cu​v=fv​y+cv​​{u^=fu​x^+cu​v^=fv​y^​+cv​​

两式相减，并代入 f u x = u − c u f_ux=u-c_u fu​x=u−cu​、 f v y = v − c v f_vy=v-c_v fv​y=v−cv​以及径向畸变公式

{ x ^ = x ( 1 + κ 1 r 2 + κ 2 r 4 ) y ^ = y ( 1 + κ 1 r 2 + κ 2 r 4 ) \begin{cases} \hat{x}=x\left( 1+\kappa _1r^2+\kappa _2r^4 \right)\\ \hat{y}=y\left( 1+\kappa _1r^2+\kappa _2r^4 \right)\\\end{cases} {x^=x(1+κ1​r2+κ2​r4)y^​=y(1+κ1​r2+κ2​r4)​

即得

[ ( u −cu ) r 2 ( u −cu ) r 4 ( v −cv ) r 2 ( v −cv ) r 4 ] [ κ 1 κ 2 ] =[ u ^ − u v ^ − v ] \left[ \begin{matrix} \left( u-c_u \right) r^2& \left( u-c_u \right) r^4\\ \left( v-c_v \right) r^2& \left( v-c_v \right) r^4\\\end{matrix} \right] \left[ \begin{array}{c} \kappa _1\\ \kappa _2\\\end{array} \right] =\left[ \begin{array}{c} \hat{u}-u\\ \hat{v}-v\\\end{array} \right] [(u−cu​)r2(v−cv​)r2​(u−cu​)r4(v−cv​)r4​][κ1​κ2​​]=[u^−uv^−v​]

考虑到

[ u ^ v ^ 1 ] = K[ W C    ⁣ ⁣ ⁣ R C  ⁣ pw0] [ W  ⁣  ⁣ X W  ⁣  ⁣ ⁣  ⁣  Y W  ⁣  ⁣ Z 1 ] \left[ \begin{array}{c} \hat{u}\\ \hat{v}\\ 1\\\end{array} \right] =\boldsymbol{K}\left[ \begin{matrix} _{\boldsymbol{W}}^{\boldsymbol{C}}\;\!\!\!\boldsymbol{R}& ^{\boldsymbol{C}}\!\boldsymbol{p}_{w_0}\\\end{matrix} \right] \left[ \begin{array}{c} ^{\boldsymbol{W}\!}\!X\\ ^{\boldsymbol{W}\!}\!\!\:\!\:Y\\ ^{\boldsymbol{W}\!}\!Z\\ 1\\\end{array} \right] ⎣⎡​u^v^1​⎦⎤​=K[WC​R​Cpw0​​​]⎣⎢⎢⎡​WXWYWZ1​⎦⎥⎥⎤​

对于$M$幅标定板图片，每张图片取$N$个角点，则可构造非齐次线性方程

D 2 M N × 2 κ =d 2 M N × 1\boldsymbol{D}_{2MN\times 2}\boldsymbol{\kappa }=\boldsymbol{d}_{2MN\times 1} D2MN×2​κ=d2MN×1​

其中 κ = [ κ 1 κ 2 ] T \boldsymbol{\kappa }=\left[ \begin{matrix} \kappa _1& \kappa _2\\\end{matrix} \right] ^T κ=[κ1​​κ2​​]T即为畸变系数。通过最小二乘法进行回归

κ = ( D TD )− 1 DT d {\boldsymbol{\kappa }=\left( \boldsymbol{D}^T\boldsymbol{D} \right) ^{-1}\boldsymbol{D}^T\boldsymbol{d}} κ=(DTD)−1DTd

即得畸变系数

4 去畸变原理与实战

综合上面两节，给出去畸变的原理

• 准备张氏标定法的标准棋盘格，用相机对其进行不同角度拍摄，得到一组图像
• 忽略畸变，每张标定板可通过世界坐标与像素坐标的关系，求解二者间的单应性矩阵，用于反解相机内参矩阵；忽略切向畸变与内参矩阵误差，反解径向畸变参数
• 利用L-M(Levenberg Marquardt)等非线性迭代优化算法对上述参数进行优化，得到最终的内参矩阵与畸变参数。

按照上面的原理，我们将相机标定与去畸变的代码写出来

/* * @breif:采用张氏标定法标定相机并生成标定文件.txt * @param[in]:folderPath->标定板图像集的文件夹路径; * @param[in]:boardSize->标定板上每行、列的内角点数; * @param[in]:realSizePerBox->实际测量的标定板棋盘格几何尺寸，与相机矩阵单位pix/mm对应，Unit:mm; * @param[in]:savePath->标定文件.txt的保存路径; * @retval:None */void cameraCalibration(String folderPath, Size boardSize, Size realSizePerBox, String savePath){std::vector<cv::String> fileNames;// 存放标定图片文件名的序列cv::glob(folderPath, fileNames);// 将标定文件夹中的文件名读取到序列Mat cameraMatrix = Mat(3, 3, CV_32FC1, Scalar::all(0)); // 像机内参矩阵Mat distCoeffs = Mat(1, 5, CV_32FC1, Scalar::all(0));// 像机畸变系数：k1,k2,p1,p2,k3vector<Mat> tvecsMat;// 相机外参->旋转向量vector<Mat> rvecsMat;// 相机外参->平移向量int imgCnt = 0;// 图像数量vector<int> ptCnt;// 图像中角点的数量Size imgSize;// 图像尺寸vector<Point2f> imgPtBuffer;// 缓存每幅图像上检测到的角点vector<vector<Point2f>> imgPtSequence;// 检测到的所有角点序列，Unit:pixvector<vector<Point3f>> realImgPtSequence;// 检测到的所有角点序列，Unit:mmdouble totalErr = 0.0;// 所有图像的重投影误差和/*========================================= 像素平面角点坐标提取 =========================================================*/for (size_t i = 0; i < fileNames.size(); ++i){imgCnt++;printf("read the number %d image\n", i+1);Mat imgInput = imread(fileNames[i]);ptCnt.push_back(boardSize.width * boardSize.height);// 初始化每幅图像中的角点数量，假定每幅图像中都可看到完整标定板if (imgCnt == 1)//读入第一张图片时获取图像宽高信息{imgSize.width = imgInput.cols;imgSize.height = imgInput.rows;}if (findChessboardCorners(imgInput, boardSize, imgPtBuffer) == 0)// 提取ChessBoard类型角点{printf("can't find the corner point of image number %d, please check!\n", i);return;}else{Mat imgGray;cvtColor(imgInput, imgGray, CV_RGB2GRAY);//图像灰度化find4QuadCornerSubpix(imgGray, imgPtBuffer, Size(11, 11));//对粗提取的角点进行亚像素精确化imgPtSequence.push_back(imgPtBuffer);}}/*========================================= 几何平面角点坐标提取 =========================================================*/for (int i = 0; i < imgCnt; i++){vector<Point3f> tempPtSetPerImg;for (int j = 0; j < boardSize.height; j++){for (int k = 0; k < boardSize.width; k++){Point3f realPoint;realPoint.x = j * realSizePerBox.width;realPoint.y = k * realSizePerBox.height;realPoint.z = 0;// 张氏标定法假设标定板放在世界坐标系中z=0的平面上tempPtSetPerImg.push_back(realPoint);}}realImgPtSequence.push_back(tempPtSetPerImg);}/*======================================== 相机标定与重投影误差计算 ======================================================*/calibrateCamera(realImgPtSequence, imgPtSequence, imgSize, cameraMatrix, distCoeffs, rvecsMat, tvecsMat);totalErr = calReprojectErr(realImgPtSequence, imgPtSequence, rvecsMat, tvecsMat, cameraMatrix, distCoeffs, ptCnt);}

看看效果，只看标定格：

🚀 计算机视觉基础教程说明

章号                                    内容
  0                              色彩空间与数字成像
  1                              计算机几何基础
  2                              图像增强、滤波、金字塔
  3                              图像特征提取
  4                              图像特征描述
  5                              图像特征匹配
  6                              立体视觉
  7                              项目实战

🔥 更多精彩专栏：

• 《机器人原理与技术》
• 《计算机视觉教程》
• 《机器学习》
• 《嵌入式系统》
• 《数值优化方法》
• …

🏠 欢迎加入社区和更多志同道合的朋友交流：AI 技术社

👇配套代码 · 优质体验 · 系统知识 请关注👇