双目立体视觉之单相机标定（一、相机模型）

近段时间，由于项目需要做了部分拼接的工作，虽然是使用opencv的sticher类进行实现的，但是出于好奇依然对这类多视角图像的融合技术进行了比较浅显的调研和了解，深入查阅资料后其实发现其本质依然是图像与图像间的变换关系，这部分知识在双目立体视觉中有比较深入的理论研究和数学推导。在笔者看来，技术可以是千变万化的，但是其核心思想与数理方法是具有很强的迁移性的，因此接下来几篇博客就深入讨论一下双目立体视觉的数学原理。

提到双目立体视觉，我们就不得不花很大的功夫去学习一下张正友教授在1998年提出的张氏标定法【a flexible new technique for camera calibration】，引用量达到1.1w的大牛文章，是做cv的必读经典文章，因此第一篇博客就花一些功夫去总结一下这篇文章，同时这也是做双目的一些基础。本博客中很多图直接引用文章中内容，在此由衷地向前辈大牛表示深深的敬意和感谢，哈哈，扯远了。另外，对于需要深入研究计算机视觉的原理方法的人，我在这里强烈推荐两本英文教材：《computer vision: algrithms and applications》以及《multiple-view geometry》。

1.针孔相机模型

图1.1 针孔相机模型

图1.1是笔者自己画的一个针孔相机模型的示意图，其实针孔模型非常好理解，就是我们小时候玩的小孔成像，而这个小孔的位置就在X_c和Y_c坐标系的原点。根据相似三角形原理，真实像面（左边）可以等效对称到右边的对称面上去，因此我们一般看到的模型是省略了CCD的真实位置的，如图1.2所示。

图1.2 针孔相机模型及其参数意义

抛开像差等问题，针孔模型的局限性之一就在于它忽略了实际镜头成像的景深，不能真实反应长工作距、大相对孔径镜头的浅景深系统的成像特性【都是笔者的猜想，欢迎大家留言批评指正】，虽然不影响系统的标定，但是对于我们在建立这类问题的成像模型用于复原等情况时，这个估计就没法直接使用了。我曾经纠结了很久这个focal length到底是不是焦距，按照实际推导应为像距。我们生活中常用的广角相机，由于其短焦特性，一般拍摄稍微远一点的物体时，根据高斯成像公式，可知此时像面离镜头的焦平面及其近，此时其远景一般都是到无穷远的，这种情况非常满足针孔模型......

另外，细心的读者应该已经发现图1.1我把相机坐标系中的三维空间点表示为（X,Y,Z,1），这就是齐次坐标，至于什么是齐次坐标，大家感兴趣的可以好好去翻翻射影变换的书。使用齐次坐标最重要的目的是合并矩阵运算中的乘法和加法，就细节而言有以下好处：

（1）点、线、面的统一化表达，线和面都可以直接使用齐次坐标进行表达。

（2）在（1）的基础上，点线面的运算就极其简洁了，比如点x在线l上就可以直接表达为x‘l=0，这些性质是数学描述对极几何的基础。

（3）无穷远点（X,Y,Z,0）的引入，这是一个很自然的数学联想

2.成像建模

标定的过程说简单一点就是探索某些已知的真实世界和图像世界之间的联系，这些联系的数理表征就是我们常用的建模了。

2.1 内参（intrinsic matrix）

内参，顾名思义就是相机内部参数，即由相机的配置决定，它包含了图像的缩放（仿射变换），像素坐标与像面坐标的偏移（）

相机坐标系中的点X到其像点x的转换关系可以由图1.2获得：

用矩阵表示为：

我们在计算时，图像都是对图像平面采样量化后的离散矩阵，其坐标原点在左上角，而模型中的图像平面的坐标原点在中心，这就意味着图像坐标与像素坐标存在一个偏移，同时像素坐标存在采样量化间隔（像素有大小），考虑这些因素，即内参矩阵表示为：

其中：

（px，py）表示主点即相机中心在图像坐标系中的位置

mx和my分别表示单位长度（一般采用m）内x和y方向的像素的个数。

一般我们通过标定后获得的内参就是这四个参数，其具体的物理意义其上已经稍作解释。

2.2 外参（extrinsic matrix）

以上的X点的坐标位置是相对于相机坐标系的，而X点是世界坐标中的点，因此需要将其坐标值变换到相机坐标内。以上操作可以通过一个旋转和平移矩阵做到，即先平移世界坐标系让其原点和相机坐标系原点重合，再绕着原点进行旋转变换。其对应点之间的变换关系表示为：

    

2.3 单映矩阵（P）

单映矩阵表示世界坐标中的点在图像像素坐标的成像对应关系，因此其是内参和外参的共同作用体。

3.标定过程

有了上面的模型，我们就很自然而然知道了标定其实就是求取P矩阵中的各个参数，并将其分解为K，R，t。由2.3我们建立了一个映射关系，我们可以设计已知场景（本文采用棋盘格）进行拍摄，通过上述映射方程求取P中未知参数。这里引用文章中的经验结论，一般需要拍摄10-15幅以上不同姿态棋盘格的图片就行拍摄。这里贴上matlab的标定代码：

% Define images to process
clear;
clc;
close all;
imageFileNames = {};

% Detect checkerboards in images
[imagePoints, boardSize, imagesUsed] = detectCheckerboardPoints(imageFileNames);
imageFileNames = imageFileNames(imagesUsed);

% Read the first image to obtain image size
originalImage = imread(imageFileNames{1});
[mrows, ncols, ~] = size(originalImage);

% Generate world coordinates of the corners of the squares
squareSize = ;  % in units of 'centimeters'
worldPoints = generateCheckerboardPoints(boardSize, squareSize);

% Calibrate the camera
[cameraParams, imagesUsed, estimationErrors] = estimateCameraParameters(imagePoints, worldPoints, ...
    'EstimateSkew', false, 'EstimateTangentialDistortion', false, ...
    'NumRadialDistortionCoefficients', 2, 'WorldUnits', 'centimeters', ...
    'InitialIntrinsicMatrix', [], 'InitialRadialDistortion', [], ...
    'ImageSize', [mrows, ncols]);

% View reprojection errors
h1=figure; showReprojectionErrors(cameraParams);

% Visualize pattern locations
h2=figure; showExtrinsics(cameraParams, 'CameraCentric');

% Display parameter estimation errors
displayErrors(estimationErrors, cameraParams);

% For example, you can use the calibration data to remove effects of lens distortion.
undistortedImage = undistortImage(originalImage, cameraParams);

由于方程数大于未知数个数，即可以采用最小二乘进行拟合，文中使用的LM算法。我们在下一章博客结合详细源代码深入讲解其求解算法和矩阵特性约束。