相机标定（具体过程详解）张正友、单应矩阵、B、R、T

首先先声明一点，本文介绍的方法并非很多文章中介绍的方法只是应用，直接调用matlab或者opencv的内部原有相机标定函数，本文主要介绍求解相机参数的一个大体流程。

相机标定分：内参标定、外参标定、c矩阵

内参标定：采用棋盘图或者点阵图，将相机固定（焦距、光圈等内部元素）。用相机拍摄多组不同角度的棋盘图，这点和普通的matlab的棋盘图标定类似。matlab利用棋盘图标定，它也给出了标定的外参，但是给定的并非我们所求，而是将世界坐标系定在了棋盘图上，我们需要的世界坐标系是固定不变的或者至少是与相机的相对位置不会发生改变。
内参标定输入为：像素坐标u、v和世界坐标x、y、z，该z没有定义具体的数值，因为在此处的求解不会应用到，并且我们此步也无法求得。那么这些参数怎么来的呢？u、v为图像中棋盘格的角点坐标，与其对应的x、y表示该角点在整个棋盘图下的坐标，该坐标系定在棋盘图上，因此我们拍摄5次图像，这些图像相应的坐标应该相同。为了计算时数据的变换较为方便，这里z也给出，值给0或1都行，不影响计算。
内参标定输出为：内部参数+相机畸变参数。内部参数：fu、fv、u、v；相机畸变参数：k1，k2

内参求解步骤：

1.求单应矩阵

这里的λ只是1/s的另一种表示，并非内参矩阵的λ，谨记

这里的h33设置为1……课本中给出的理由是它可以设置为任意数，二郎也和同事讨论过，没有定论，有一种理解有道理：计算中h33能都提到右边，那么在根据公式进行求解时，它变成几，都是对所有公式进行操作的，因此公式之间的关系不会发生改变，所求出来代表关系的h参数也不会受到影响，只是按h33进行了缩放。
比如：两边先同时除以h33（由于没有找到这里的h33=1的解释，因此二郎做了这种假设）
每个特征点，可以求得两个方程。这里有8个未知数，因此可以用4个特征点完成求解。
这里用到了棋盘格，在一张棋盘格图中，我们得到了很多（u，v）与（x，y）的对应，因此可以解出我们的单应矩阵。
这里每个棋盘格图对应了一个单应矩阵，这些单应矩阵主要的不同在于（R，t），每个棋盘图都有自己的世界坐标系，原点在棋盘图的一个角点上，这个matlab内部程序已经指定，一般情况下为左上角。此处单应矩阵的求解用到的最小二乘法。下面的推导和求解，均是建立在单应矩阵已知的情况下进行的。
本文直接调用的matlab的单应矩阵求解代码。

H = homography(img(:,tnum:tnum+snum(i)-1),world(:,tnum:tnum+snum(i)-1));

其中img代表图像坐标的uv，world代表世界坐标的xy。tnum:tnum+snum(i)-1)则是将我们拍的五组图分散开来。

2.单应矩阵构造v，进行B的求解

这里有人会不理解，r1=和r2=是怎么来的？其实自己可以推导一下，把A设成[3X3]的矩阵，乘进去，应该能够看到，组成了[Ar1 Ar2 Ar3]，它们互相还是独立的。至于为什么A到A逆，这个是最基本的转换，左右同时左乘A逆，看一下，是不是变过来了。

这里的A矩阵为内参矩阵，上面两个式子其实不难证明，
第一个式子：其实是$A^{-1}{h}_1$和$A^{-1}{h}_2$的点乘关系，如果他们垂直那么点乘便为0。
$A^{-1}H=A^{-1}\lambda A\left[\begin{array}{ccc}{\text{r}}_1& {\text{r}}_2& \text{t}\end{array}\right]\text{=}\left[\begin{array}{ccc}{\text{r}}_1& {\text{r}}_2& \text{t}\end{array}\right]$
所以$h_1^T{A}^{-T}{A}^{-1}{h}_2={\text{r}}_1^T{\text{r}}_2=0$
第二个式子是
${\text{r}}_1^T{\text{r}}_1={\text{r}}_2^T{\text{r}}_2$由此便可以看出这两步都是在说明r1和r2是单位正交的关系。

由于B为对称矩阵，因此B中元素可以写为

大家可以从上面看到B是由A矩阵组成的，那们A的组合满足的关系，B也满足



V的由来……大部分文献都没有提到，其实比较简单，直接计算hBh，计算完构造vb，b已经在前面定义。
后面的vb等于0的两项也来自前面约束条件

用奇异值分解最小奇异值对应的特征向量便可以表示为b

[u s v]=svd(V);        %用正交分解解出b
b=v(:,6);

图片来自：https://www.cnblogs.com/litaotao-doctor/p/5324160.html

3.利用B，求出内参矩阵A

$A=\left[\begin{array}{ccc}{f}_u& s& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right]$

4.优化——提高精度关键步骤LM算法

注意，这里的λ和B求出的不是一个概念，这个λ是1/s

上边的λ为了精确，选用了两个的平均值。
求解出来的旋转矩阵需要进行正交化处理。

RL=[r1 r2 r3]; 
[u s v]=svd(RL);
RL=u*v';    

利用上述步骤求解出了外参矩阵。
因此我们可以通过求解出的内参矩阵和外参矩阵，输入外参[X,Y,Z]求得内参[u’,v’]。

x=A*ML*M(:,tnum:tnum+snum(i)-1);

其中M矩阵为外参数据，ML为拥有R和t的外参矩阵，x求得的内参[u’,v’]。
这里需要强调，为了归一化处理保证Z为1，因此X=X/Z;Y=Y/Z，即[λu’,λv’,λ]/λ

x=[x(1,:)./x(3,:) ; x(2,:)./x(3,:); x(3,:)./x(3,:)];    

最后利用LM方法，优化x的[u’,v’,1]与实际的[u,v,1]的差值，使其达到最小，进而优化参数。

5.利用理想的[u’,v’]与实际的[u,v]求解内参

这里仅考虑了主要畸变原因：径向畸变，因此求解的畸变参数仅为k1、k2
本节需要强调一下，下面出现的[u,v]其实是设成了理想的图像坐标，而像[u’,v’]的一项为实际图像坐标，这里要注意。至于理想的图像坐标来源是通过实际坐标通过修正的外参矩阵和内参矩阵的转变而来的。
首先我们研究理想的外部坐标系有畸变时

又因为其与图像坐标是对应关系

因此图像畸变表示为

下面的矩阵形式是我们主要用到的

最终公式

k=inv(D'*D)*D'*d;   

6.把畸变系数k1、k2代入矩阵，如上一次优化一样

比较求得和实际图像的是否一样，需要注意这里的[u,v]虽说理想的，但是他们并不直接来自图像坐标，而是通过外部数据[X,Y,Z]变化的到的。

7.优化后通过优化的参数，做最后优化

优化完后，我们便可以用我们优化过的内参外参以及畸变系数，求取外部坐标对应的“实际内参”，然后

再求取一次k1,k2

此时我们的内部参数便全部求完了

这中间可以加入一步初始化，类似于F初始化求解。