计算机视觉大型攻略 —— 立体视觉（2）张正友标定

最新推荐文章于 2023-10-16 21:08:07 发布

linusyue

最新推荐文章于 2023-10-16 21:08:07 发布

阅读量2.1k

点赞数 6

分类专栏：立体视觉文章标签：立体视觉张正友标定校准 undistort

本文链接：https://blog.csdn.net/plateros/article/details/102564202

版权

立体视觉专栏收录该内容

5 篇文章 20 订阅

订阅专栏

这篇文章是立体视觉系列的第二部分，讲解相机标定和消除畸变。

双目相机校正一般分为两步，

消除图像失真 (标定矫正，Calibration and undistortion)
立体校正 (Stereo Rectification)

在求深度或者做3D重建算法的时候，我们都会假设相机符合一定的数学模型，但是由于种种原因导致的失真，使我们最终呈现的图像偏离了这些模型。这就需要我们首先去标定和矫正。之后，立体匹配算法为了性能和准确性，需要对图片做极线平行，就需要做第二步。

这一篇说说标定和矫正算法，这类算法最经典的就是张正友标定法。本文大量图片和公式推导出自他的原始论文及Learning Opencv3。

原始论文及出处：A Flexible New Technique for Camera Calibration, IEEE Transactions on Pattern Analysis and Machine Intelligence ( Volume: 22 , Issue: 11 , Nov 2000 )

针孔相机模型

上图是我们熟知的针孔模型，光线从远处物体发射过来，通过小孔，物体被投影到成像平面(Image plane)。物体高度X, 投影高度x, 相机到物体的距离Z，焦距f，通过相似三角形，可知，

$-x=f\frac{X}{Z}$

负号是因为目标图像是倒立的。该模型可以进一步转换成另一种模式，将成像平面(image plane)放到针孔平面(Pinhole)和物体中间。

针孔被理解为投影中心或光心，每一条光线，从远处某个点出发，到达投影中心。光轴与成像平面(Image plane)的交点称为主点(principal point)。这样形成了新的三角形相似关系，

$\frac{x}{f} = \frac{X}{Z}$

由于图像不再倒立，因此去掉了负号。令(u, v)为图像坐标系上点q的坐标，(X, Y, Z)为相机坐标系下的点Q的坐标，fx, fy分别为x轴，y轴的焦距，引入了新的公式，

$u=f_{x}\frac{X}{Z}+u_{0}, v=f_{y}\frac{Y}{Z}+v_{0}$

u0. v0是图像坐标系下Principle Point的坐标，由于装配精度的问题，数值上一般不等于图像的中心。

上述公式也可以写做，

$\vec{q} = M\cdot \vec{Q}$ ，

其中，

$\vec{q}=w\begin{bmatrix} u\\ v\\ 1 \end{bmatrix}, M = \begin{bmatrix} f_{x} & 0 &u_{0} \\ 0& f_{y} &v_{0}\\ 0& 0 & 1 \end{bmatrix}, \vec{Q} = \begin{bmatrix} X\\ Y\\ Z \end{bmatrix}$

以上摘自Learning Opencv3第18章,很好理解，与原文不同的是，这里把图像坐标系的坐标全改成了u, v，以示区分。下文也均采用此命名规则。

下面开始解释张正友标定法，直接借用原作者的图。

与LearningOpenCV3的图类似，假设点M在世界坐标系下的坐标为(X,Y,Z), m是从光心C发出的射线到点M在图像上的交点，图像坐标系下的坐标为(u,v)，PrinciplePoint坐标在图像坐标系下为(u0,v0)。注意，与上面的模型稍有不同，这里的(X,Y,Z)是世界坐标系的坐标而不是相机坐标系。

按照向量的表示方法，令

$\widetilde{m}= [u,v,1]^{T},\widetilde{M}=[X,Y,Z,1]^{T}$ ，M与m的关系可以定义为，

其中，s为任意缩放系数。(R, t)是世界坐标系到相机坐标系的旋转平移矩阵，~~通过RT将世界坐标系下的坐标转到相机坐标系，~~称为外参(extrinsic parameters)，虽然R是一个3*3的矩阵，t是3*1的矩阵，但其实只有6个参数(3个角度，3个平移）。

A是相机内参，一共有五个参数(u0,v0,α,β,γ)

其中u0, v0是principal point的坐标，一般情况下不等于width/2, height/2。
α，β分别是u方向和v方向上的缩放系数，即我们常说的fx, fy。如果α，β不一致，就会导致，

γ~~(gamma)~~是倾斜系数。定义为γ=cotθ, θ为u,v两个轴的夹角。如果像素是方的且θ=90度，那么γ=0。否则，

校正工作的第一步就是要获得这6个外参和5个内参。抛开γ(gamma)不谈，如果令RT等于单位阵，~~也就是说让世界坐标系等于相机坐标系~~，这个公式和上述引自LearningOpencv3的公式完全相同。

畸变模型

畸变模型引用张正友的相关文献，常见畸变成因有两种，径向畸变(radial distortion)，偏心畸变。

径向畸变指的是光线在远离透镜中心的地方比靠近中心的地方更加弯曲。比如说本来是直线，成像后弯曲了。一般是由于镜头形状的误差导致。

偏心畸变一般是由于镜片装配误差导致，会产生径向畸变和切向畸变。

关于这方面的详细内容请参考下面这个论文。这里只放结论了。

Camera Calibration with Distortion Models and Accuracy Evaluation. J Weng，P Cohen，M Herniou - 《Pattern Analysis & Machine Intelligence IEEE Transactions on》

首先对误差建模，

$\bar{x}=x+\delta _{x}, \bar{y}=y+\delta _{y},$

$(\bar{x}, \bar{y})$ ， (x, y) 分别表示相机坐标系下点P未发生畸变的和发生畸变的坐标。 $(\delta _{x}, \delta _{y})$ 为畸变系数。

令 $r=\sqrt{x^{2}+y^{y}}$ , 畸变系数可定义为，

其中，ki为径向畸变的系数，pi为偏心畸变的系数。

张正友标定法

从上面的公式可以看到，只要我们已知足够的(u,v)和(X,Y,Z)的对应关系，联立方程组就可以把P求出来，再通过一些数学手段就可以求出A, R, t。所以，传统的校准方式采用的是3D治具精确的获得治具的3D坐标。不过他的这样做缺点也很多，首先，3D治具本身价格不菲，其次，对治具的精度要求也很高。张正友在98年提出了一种基于2D平面的标定方法（与Z无关），治具和操作简单，结果准确。基本过程就是准备一张棋盘格，在相机前面以不同的姿态，位置拍摄几十张照片。通过这些照片估算出相机的内外参。

先引入一个概念，单应矩阵。单应矩阵有很多含义这里只是简单的表达平面上的点到图像上点的映射关系。存在一个结论，给定一个平面的图像，总是可以估计出他的单应矩阵(Homography)。原论文的附录A给出了极大似然法估计的证明。

定义一个平面到图像的单应矩阵(Homography)

既然和Z无关，那么可以令Z=0，根据上面的公式，可以得到

令， H即为单应矩阵。由于给定一个平面(棋盘格)的图像，他的单应矩阵一定可以求出来，那么H就可以当做已知变量。问题就变成了已知H，求r1, r2, t, A。

求解内外参

将H写成列向量的形式，

，

可得， $r_{1} = \lambda A^{-1}h_{1}$ , $r_{1} = \lambda A^{-1}h_{1}$

由于旋转矩阵r1和r2正交，即

$r_{1}^{T}r_{2} = 0$ ， $\left \| r_{1} \right \|=\left \| r_{2} \right \| = 1$

可得约束条件，

原始论文中同时给出了这两个约束几何解释，这里略过。接下来

已知，

，

令

注意到B是一个对称阵，只有6个元素生效。因此可以写成，

，

令hi表示H的第i个列向量， $h_{i} = \begin{bmatrix} h_{i1} & h_{i2} &h_{i3} \end{bmatrix}^{T}$

可得，

且，

上面两个约束公式可以进一步写成，

因此，给定一个H（或者说一幅棋盘格），我们可以给出上面两个公式，而b总共6个参数，因此当H的个数>=3时，我们可以把六个未知数求出来，从而求出B。通过B进一步求出所有内参，

求出A后，根据以下公式求外参。

极大似然估计优化

以上的算法并没有考虑实际中的存在的误差，因此，作者采用极大似然法做了进一步优化。

假设有n个棋盘格，每个棋盘格上有m个点。假设这些点符合同样分布的噪声，且彼此独立。那么极大似然估计等效于求下面公式的最小值，

，这是一个非线性最小化的问题，作者采用了Levenberg-Marquardt算法求解。这种方法可以看成是牛顿法的一种改良。具体细节不在这里说明了。

由于成像畸变主要由径向畸变导致，张氏标定法使用极大似然的思想对k1, k2进行了优化。优化方式同样采用可以采用Levenverg-Marquardt算法。

矫正

将校准得到的内参和畸变系数代入畸变模型，就可以得到未发生畸变的新的像素坐标。当畸变图像M转向非畸变图像M'时，M'上并不是所有像素都有值，这里还需要差值处理。

张正友标定法只对径向畸变的k1, k2做了矫正。令 $(\breve{x} \; \breve{y})$ 为未发生畸变的点在相机坐标系下的点，(x y)为发生了畸变的点的坐标。 $(\bar{u}, \bar{v})$ , (u, v)为相应的图像坐标系上的坐标。

转换成图像坐标系

代码实例

Matlab与Opencv均实现了相机校准和矫正的功能。Opencv4中的代码实例是sample/cpp/calibration。我跑了一下，与数据集一起放到了github上: https://github.com/linusyue/calibMono.git

该程序运行参数，

 ./calibMono -w=9 -h=6 -s=0.02423 -o=camera.yml -op -oe -su imagelist.xml

w是宽度方向角点个数，h为高度方向角点个数。注意这里角点指的是内角点，数据集棋盘格是10*7，内角点就是9*6，s是单个棋盘格的边长，o是生成的相机内外参数，imagelist.xml指明了棋盘图像的路径。

程序大体流程，

从视频流或者文件读入图像。
对每一个图像调用findChessboardCorners，这个函数通过Harris角点算法获取内角点坐标，

503             case CHESSBOARD:
504                 found = findChessboardCorners( view, boardSize, pointbuf,
505                     CALIB_CB_ADAPTIVE_THRESH | CALIB_CB_FAST_CHECK | CALIB_CB_NORMALIZE_IMAGE);
506                 break;

再调用cornerSubPix获取亚像素坐标值，

517        // improve the found corners' coordinate accuracy
518         if( pattern == CHESSBOARD && found) cornerSubPix( viewGray, pointbuf, Size(winSize,winSize),
519             Size(-1,-1), TermCriteria( TermCriteria::EPS+TermCriteria::COUNT, 30, 0.0001 ));

为每一幅棋盘图计算世界坐标系坐标。这里Z=0，且每幅图的世界坐标系都可以从各自最边上的格子开始，也就是说，每幅图的每个角点的世界坐标是相同的。

118 static void calcChessboardCorners(Size boardSize, float squareSize, vector<Point3f>& corners, Pattern patternType = CHESSBOARD)
119 {
120     corners.resize(0);
121 
122     switch(patternType)
123     {
124       case CHESSBOARD:
125       case CIRCLES_GRID:
126         for( int i = 0; i < boardSize.height; i++ )
127             for( int j = 0; j < boardSize.width; j++ )
128                 corners.push_back(Point3f(float(j*squareSize),
129                                           float(i*squareSize), 0));
130         break;

使用calibrateCameraRO校准返回内参矩阵cameraMatrix, 畸变系数distCoeffs, 外参rvecs, tvecs。

171     rms = calibrateCameraRO(objectPoints, imagePoints, imageSize, iFixedPoint,
172                             cameraMatrix, distCoeffs, rvecs, tvecs, newObjPoints,
173                             flags | CALIB_FIX_K3 | CALIB_USE_LU);
174     printf("RMS error reported by calibrateCamera: %g\n", rms);
175

最后通过initUndistortRectifyMap和remap对图片做矫正。initUndistortRectifyMap根据畸变模型的公式，计算出两个map，大小等同于图片大小，存储了畸变像素和未畸变像素的关系。remap通过这个关系，将畸变图像转换成无畸变图像。INTER_LINEAR为差值方法。

587     if( !capture.isOpened() && showUndistorted )
588     {
589         Mat view, rview, map1, map2;
590         initUndistortRectifyMap(cameraMatrix, distCoeffs, Mat(),
591                                 getOptimalNewCameraMatrix(cameraMatrix, distCoeffs, imageSize, 1, imageSize, 0),
592                                 imageSize, CV_16SC2, map1, map2);
593 
594         for( i = 0; i < (int)imageList.size(); i++ )
595         {
596             view = imread(imageList[i], 1);
597             if(view.empty())
598                 continue;
599             //undistort( view, rview, cameraMatrix, distCoeffs, cameraMatrix );
600             remap(view, rview, map1, map2, INTER_LINEAR);
601             imshow("Image View", rview);
602             char c = (char)waitKey();
603             if( c == 27 || c == 'q' || c == 'Q' )
604                 break;
605         }
606     }

Opencv的实现其实是在Bouguet的MATLAB tool box的基础上直接翻译的。 http://www.vision.caltech.edu/bouguetj/calib_doc/

Bouguet在张正友标定法的基础上做了大量的改进，比如，先通过vanish ponit法粗略求内外参，再通过非线性优化算法LM优化求解，不仅对径向畸变求解，也求解了切向畸变。具体的代码在modules/calib3d/src/calibration.cpp中。有一个PPT讲解了Opencv的实现，上传到CSDN https://download.csdn.net/download/plateros/11890853

获得了相机内参和畸变系数后，就可以进行双目标定。下一篇会详细说明。