前言

前面的几篇博客中介绍了有关相机标定的基础知识（三维重建学习(1)：基础知识：旋转矩阵与旋转向量、三维重建学习(2)：相机标定基础）。这次介绍一个十分经典的单目相机标定方法——张正友标定，并给出数学理论推导。

基本方程

模型

我们首先约定如下表示：
二维点坐标：m=[uv]m = \begin{bmatrix} u \\ v \end{bmatrix}，三维点坐标：M=⎡⎣⎢XYZ⎤⎦⎥M = \begin{bmatrix} X \\ Y \\ Z \end{bmatrix}
将上面的二维点和三维点的坐标表示成齐次形式：
二维点坐标：m~=⎡⎣⎢uv1⎤⎦⎥\tilde{m} = \begin{bmatrix} u \\ v \\ 1 \end{bmatrix}，，三维点坐标：M~=⎡⎣⎢⎢⎢XYZ1⎤⎦⎥⎥⎥\tilde{M} = \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix}

参考前面的博客：三维重建学习(2)：相机标定基础，我们可以建立一个常见的相机小孔成像模型：

s⋅m~=A[Rt]M~

s \cdot \tilde{m} = A \begin{bmatrix} R & t \end{bmatrix} \tilde{M}
其中：

A=⎡⎣⎢α00cβ0u0v01⎤⎦⎥

A = \begin{bmatrix} \alpha & c & u_0 \\ 0 & \beta & v_0 \\ 0 & 0 & 1 \end{bmatrix}
ss为任意比例因子，RR、 tt都是相机外参，RR为旋转矩阵， tt对应三个轴的平移。AA是相机内参， (u0,v0)(u_0, v_0)是坐标的主点， α\alpha和 β\beta是图像在 uu轴和vv轴的比例因子， cc是描述两个坐标轴亲倾斜角的参数（注：如果两个坐标轴相互垂直，则c=0c=0，即默认情况下这个 cc都是为0的）。

模型平面与图像之间的单应性关系

设RR的第 ii列旋转矩阵为rir_i，那么 RR可以表示为：

R=[r1r2r3]

R = \begin{bmatrix} r_1 & r_2 & r_3 \end{bmatrix}
代入原方程中：

s⋅⎡⎣⎢uv1⎤⎦⎥=A[r1r2r3t]⎡⎣⎢⎢⎢XYZ1⎤⎦⎥⎥⎥

s \cdot \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = A \begin{bmatrix} r_1 & r_2 & r_3 & t \end{bmatrix} \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix}
假设模型平面在世界坐标系中 ZZ坐标为0，那么ZZ坐标的值都为0：

s⋅⎡⎣⎢uv1⎤⎦⎥=A[r1r2r3t]⎡⎣⎢⎢⎢XY01⎤⎦⎥⎥⎥

s \cdot \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = A \begin{bmatrix} r_1 & r_2 & r_3 & t \end{bmatrix} \begin{bmatrix} X \\ Y \\ 0 \\ 1 \end{bmatrix}
乘出来都还是0，省略掉这一部分：

s⋅[uv]=A[r1r2t]⎡⎣⎢XY1⎤⎦⎥

s \cdot \begin{bmatrix} u \\ v \end{bmatrix} = A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix} \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix}
此时坐标表示也要变换一下：
二维点坐标： m~=[uv]\tilde{m} = \begin{bmatrix} u \\ v \end{bmatrix}，，三维点坐标： M~=⎡⎣⎢XY1⎤⎦⎥\tilde{M} = \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix}

因此点MM和它在图像上的映射点mm之间的关系可以使用单应矩阵HH来表示：

s⋅m~=H⋅M~

s \cdot \tilde{m} = H \cdot \tilde{M}
其中：

H=A[r1r2t]

H = A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix}
很显然， HH是一个3×33 \times 3的矩阵。 AA对应着内参，[r1r2t]\begin{bmatrix} r_1 & r_2 & t \end{bmatrix}对应着外参。

内参的约束条件

令H=[h1h2h3]H = \begin{bmatrix} h_1 & h_2 & h_3 \end{bmatrix}，h1h_1、h2h_2、h3h_3都是3×13 \times 1的矩阵，各自对应一列。
则有[h1h2h3]=λA[r1r2t]\begin{bmatrix} h_1 & h_2 & h_3 \end{bmatrix} = \lambda A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix}，式中λ\lambda是任意的标量。

我们知道旋转矩阵的每一列两两正交，即r1r_1与r2r_2正交。这里不对基础知识赘述，如果有疑问，请查看：旋转矩阵(Rotate Matrix)的性质分析。（吐槽一下：旋转矩阵真的是一个完美的矩阵）
根据r1r_1与r2r_2正交，我们能得到条件：

rT1r2=0,∥r1∥=∥r2∥=1

r_1^Tr_2 = 0,\|r_1\|=\|r_2\|=1
上面这个条件先放在这里，后面再用。
由前面公式： [h1h2h3]=λA[r1r2t]\begin{bmatrix} h_1 & h_2 & h_3 \end{bmatrix} = \lambda A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix}一一对应得到方程组：

⎧⎩⎨h1=λAr1h2=λAr2h3=λAt

\begin{cases} h_1 = \lambda A r_1 \\ h_2 = \lambda A r_2 \\ h_3 = \lambda A t \end{cases}
推出：

{r1=λ−1A−1h1r2=λ−1A−1h2

\begin{cases} r_1 = \lambda^{-1} A^{-1} h_1 \\ r_2 = \lambda^{-1} A^{-1} h_2 \end{cases}
代入前面的条件：

rT1r2=0,∥r1∥=∥r2∥=1

r_1^Tr_2 = 0,\|r_1\|=\|r_2\|=1
得到：（ λ\lambda是常数）

rT1r2=λ−ThT1A−TA−1h2λ−1=λ−2hT1A−TA−1h2=0

r_1^Tr_2 = \lambda^{-T} h_1^{T} A^{-T} A^{-1} h_2 \lambda^{-1} = \lambda^{-2} h_1^{T} A^{-T} A^{-1} h_2=0

{∥r1∥2=λ−2hT1A−TA−1h1=1∥r2∥2=λ−2hT2A−TA−1h2=1⇒λ−2hT1A−TA−1h1=λ−2hT2A−TA−1h2

\begin{cases} \|r_1\|^2 = \lambda^{-2} h_1^{T} A^{-T} A^{-1} h_1 = 1 \\\|r_2\|^2 = \lambda^{-2} h_2^{T} A^{-T} A^{-1} h_2 = 1 \end{cases} \\ \Rightarrow \lambda^{-2} h_1^{T} A^{-T} A^{-1} h_1 = \lambda^{-2} h_2^{T} A^{-T} A^{-1} h_2

由上面的式子我们可以发现，对于一个给定的单应性矩阵HH，对于内参有2个基本的约束条件：λ−2hT1A−TA−1h2=0\lambda^{-2} h_1^{T} A^{-T} A^{-1} h_2=0和λ−2hT1A−TA−1h1=λ−2hT2A−TA−1h2\lambda^{-2} h_1^{T} A^{-T} A^{-1} h_1 = \lambda^{-2} h_2^{T} A^{-T} A^{-1} h_2。
对于HH矩阵来说，它是一个3×33\times3的矩阵，有9个参数，那么就有8个自由度。对应的外参有6个（注：旋转矩阵RR有3个，平移向量tt有3个）。到这一步，我们只能得到内参的约束条件，却没法解出来。

解决相机标定

利用约束条件求出内参矩阵AA

好的，还是看一下前面给出的约束条件，注意到中间都有这么一个东西：A−TA−1A^{-T} A^{-1}。
我们试着将它表示出来：令B=A−TA−1B=A^{-T} A^{-1}。
我们在最前面给出了AA的定义：

A=⎡⎣⎢α00cβ0u0v01⎤⎦⎥

A = \begin{bmatrix} \alpha & c & u_0 \\ 0 & \beta & v_0 \\ 0 & 0 & 1 \end{bmatrix}
计算A的逆矩阵，步骤就省略了，很简单。得到：

A=⎡⎣⎢⎢⎢⎢⎢1α00−cαβ1β0v0c−u0βαβ−v0β1⎤⎦⎥⎥⎥⎥⎥

A = \begin{bmatrix} \frac{1}{\alpha} &- \frac{c}{\alpha \beta} & \frac{v_0 c - u_0 \beta}{\alpha \beta} \\ 0 & \frac{1}{\beta} & - \frac{v_0}{\beta} \\ 0 & 0 & 1 \end{bmatrix}
计算出BB来：

B=⎡⎣⎢⎢⎢⎢⎢⎢⎢1α−cαβv0c−u0βαβ01β−v0β001⎤⎦⎥⎥⎥⎥⎥⎥⎥⋅⎡⎣⎢⎢⎢⎢⎢1α00−cαβ1β0v0c−u0βαβ−v0β1⎤⎦⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢1α2−cα2βv0c−u0βα2β−cα2βc2α2β+1β2−c(v0c−u0β)α2β2−v0β2v0c−u0βα2β−c(v0c−u0β)α2β2−v0β2(v0c−u0β)2α2β2+v20β2+1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥

\begin{align} B &= \begin{bmatrix} \frac{1}{\alpha} & 0 & 0 \\ - \frac{c}{\alpha \beta} & \frac{1}{\beta} & 0 \\ \frac{v_0 c - u_0 \beta}{\alpha \beta} & - \frac{v_0}{\beta} & 1 \end{bmatrix} \cdot \begin{bmatrix} \frac{1}{\alpha} &- \frac{c}{\alpha \beta} & \frac{v_0 c - u_0 \beta}{\alpha \beta} \\ 0 & \frac{1}{\beta} & - \frac{v_0}{\beta} \\ 0 & 0 & 1 \end{bmatrix} \\&= \begin{bmatrix} \frac{1}{\alpha^2} & -\frac{c}{\alpha^2 \beta} & \frac{v_0 c - u_0 \beta}{\alpha^2 \beta} \\ -\frac{c}{\alpha^2 \beta} & \frac{c^2}{\alpha^2 \beta} + \frac{1}{\beta^2} & -\frac{c(v_0 c - u_0 \beta)}{\alpha^2 \beta^2} - \frac{v_0}{\beta^2} \\ \frac{v_0 c - u_0 \beta}{\alpha^2 \beta} & -\frac{c(v_0 c - u_0 \beta)}{\alpha^2 \beta^2} - \frac{v_0}{\beta^2} & \frac{(v_0 c - u_0 \beta)^2}{\alpha^2 \beta^2} + \frac{v_0^2}{\beta^2} + 1 \end{bmatrix} \end{align}
不难发现BB是对称的，我们可以使用6个变量来表示出BB。定义一个6维向量：

b=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢B11B12B22B13B23B33⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

b = \begin{bmatrix} B_{11} \\ B_{12} \\ B_{22} \\ B_{13} \\ B_{23} \\ B_{33} \end{bmatrix}
BB可以表示为：

B=⎡⎣⎢B11B12B13B12B22B23B13B23B33⎤⎦⎥

B = \begin{bmatrix} B_{11} & B_{12} & B_{13} \\ B_{12} & B_{22} & B_{23} \\ B_{13} & B_{23} & B_{33} \end{bmatrix}
与前面表示旋转矩阵时类似，我们假设HH的第ii列为：hi=⎡⎣⎢hi1hi2hi3⎤⎦⎥h_i = \begin{bmatrix} h_{i1} \\ h_{i2} \\ h_{i3} \end{bmatrix}，那么H=[h1h2h3]H = \begin{bmatrix} h_1 & h_2 & h_3 \end{bmatrix}。
对于hTiBhjh_i^T B h_j，我们把前面的公式代入看看：

hTiBhj=[hi1hi2hi3]⋅⎡⎣⎢B11B12B13B12B22B23B13B23B33⎤⎦⎥⋅⎡⎣⎢hj1hj2hj3⎤⎦⎥=(hi1B11+hi2B12+hi3B13)hj1+(hi1B12+hi2B22+hi3B23)hj2+(hi3B13+hi2B23+hi3B33)hj3=hi1hj1B11+(hi2hj1+hi1hj2)B12+hi2hj2B22+(hi3hj1+hi1hj3)B13+(hi3hj2+hi2hj3)B23+hi3hj3B33

\begin{align} h_i^T B h_j &= \begin{bmatrix} h_{i1} & h_{i_2} & h_{i3} \end{bmatrix} \cdot \begin{bmatrix} B_{11} & B_{12} & B_{13} \\ B_{12} & B_{22} & B_{23} \\ B_{13} & B_{23} & B_{33} \end{bmatrix} \cdot \begin{bmatrix} h_{j1} \\ h_{j2} \\ h_{j3} \end{bmatrix} \\ &= (h_{i1} B_{11} + h_{i2}B_{12} + h_{i3}B_{13}) h_{j1} + (h_{i1} B_{12} + h_{i2} B_{22} + h_{i3} B_{23})h_{j2} + (h_{i3} B_{13} + h_{i2} B_{23} + h_{i3} B_{33}) h_{j3} \\ &= h_{i1} h_{j1} B_{11} + (h_{i2} h_{j1} + h_{i1} h_{j2})B_{12} + h_{i2} h_{j2} B_{22} + (h_{i3} h_{j1} + h_{i1} h_{j3})B_{13} + (h_{i3} h_{j2} + h_{i2} h_{j3})B_{23} + h_{i3} h_{j3} B_{33} \end{align}
我们前面已经给出了一个6维向量bb：

b=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢B11B12B22B13B23B33⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

b = \begin{bmatrix} B_{11} \\ B_{12} \\ B_{22} \\ B_{13} \\ B_{23} \\ B_{33} \end{bmatrix}
那么上面那一堆东西可以表示为：

hTiBhj=vTijb

h_i^T B h_j = v_{ij}^T b
其中：

vij=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hi1hj1(hi2hj1+hi1hj2)hi2hj2(hi3hj1+hi1hj3)(hi3hj2+hi2hj3)hi3hj3⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥

v_{ij} = \begin{bmatrix} h_{i1} h_{j1} \\ (h_{i2} h_{j1} + h_{i1} h_{j2}) \\ h_{i2} h_{j2} \\ (h_{i3} h_{j1} + h_{i1} h_{j3}) \\ (h_{i3} h_{j2} + h_{i2} h_{j3}) \\ h_{i3} h_{j3} \end{bmatrix}
回到我们前面推导出的约束条件：

{λ−2hT1A−TA−1h2=0λ−2hT1A−TA−1h1=λ−2hT2A−TA−1h2

\begin{cases} \lambda^{-2} h_1^{T} A^{-T} A^{-1} h_2=0 \\ \lambda^{-2} h_1^{T} A^{-T} A^{-1} h_1 = \lambda^{-2} h_2^{T} A^{-T} A^{-1} h_2 \end{cases}
这两个约束条件可以改写为齐次形式：

{vT12b=0(v11−v22)Tb=0

\begin{cases} v_{12}^T b = 0 \\ (v_{11} - v_{22})^T b = 0 \end{cases}
我们用一个新的矩阵VV来表示这两个式子：

V=[vT12(v11−v22)T]

V = \begin{bmatrix} v_{12}^T \\ (v_{11} - v_{22})^T \end{bmatrix}
这里的VV是一个2×62 \times 6的矩阵。约束条件变为：V⋅b=0V \cdot b = 0
如果观察了n张图片，那么可以得到nn个方程V⋅b=0V \cdot b = 0。我们想要解出bb，bb是一个6维向量，要求出唯一解，则至少需要6个方程。一个V⋅b=0V \cdot b = 0有2个约束条件，那么要求出唯一解，至少需要3个V⋅b=0V \cdot b = 0，即至少需要3张图片（n≥3n \geq 3）。
如果我们求出了唯一解bb，那么就可以得到BB，那就也可以求出相机内参AA（注：使用cholesky分解）。
下面直接给出结果：（已知BB，求解出AA中的各个参数）

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪v0=(B12B13−B11B23)/(B11B22−B212)λ=B33−[B213+v0(B12B13−B11B23)]/B11α=λ/B11−−−−−√β=λB11/(B11B22−B212)−−−−−−−−−−−−−−−−−√c=−B12α2β/λu0=cv0/α−B13α2/λ

\begin{cases} v_0 = (B_{12}B_{13} - B_{11} B_{23}) / (B_{11} B_{22} - B_{12}^2) \\ \lambda = B_{33} - [B_{13}^2 + v_0(B_{12} B_{13} - B_{11} B_{23})] / B_{11} \\ \alpha = \sqrt{ \lambda / B_{11} } \\ \beta = \sqrt{ \lambda B_{11} / (B_{11} B_{22} - B_{12}^2)} \\ c = -B_{12} \alpha^2 \beta / \lambda \\ u_0 = c v_0 / \alpha - B_{13} \alpha^2 / \lambda \end{cases}

利用内参矩阵AA求解外参矩阵

通过前面的方法，假设我们已经求到了内参矩阵AA。
我们使用下面的式子表示点MM和它在图像上的映射点mm之间的关系：

s⋅m~=H⋅M~

s \cdot \tilde{m} = H \cdot \tilde{M}
如果我们已知图片中的点mm的坐标，以及点MM在三维空间中的坐标，s又只是个常数，那么我们可以求出单应性矩阵HH。
对于已知的单应性矩阵H=[h1h2h3]H = \begin{bmatrix} h_1 & h_2 & h_3 \end{bmatrix}，我们有：

[h1h2h3]=λA[r1r2t]

\begin{bmatrix} h_1 & h_2 & h_3 \end{bmatrix} = \lambda A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix}
解出外参：

⎧⎩⎨⎪⎪r1=λ−1A−1h1r2=λ−1A−1h2t=λ−1A−1h3

\begin{cases} r_1 = \lambda^{-1} A^{-1} h_1 \\ r_2 = \lambda^{-1} A^{-1} h_2 \\ t = \lambda^{-1} A^{-1} h_3 \end{cases}
由旋转矩阵的性质得到：

{r3=r1×r2∥r1∥=∥r2∥=1

\begin{cases} r_3 = r_1 \times r_2 \\ \|r_1\| = \|r_2\| = 1 \end{cases}
这些东西整合一下：

⎧⎩⎨⎪⎪⎪⎪⎪⎪r1=λ−1A−1h1r2=λ−1A−1h2r3=r1×r2t=λ−1A−1h3∥λ−1A−1h1∥=∥λ−1A−1h2∥=1⇒λ=∥A−1h1∥=∥A−1h2∥

\begin{cases} r_1 = \lambda^{-1} A^{-1} h_1 \\ r_2 = \lambda^{-1} A^{-1} h_2 \\ r_3 = r_1 \times r_2 \\ t = \lambda^{-1} A^{-1} h_3 \end{cases} \\ \| \lambda^{-1} A^{-1} h_1 \| = \| \lambda^{-1} A^{-1} h_2 \| = 1 \\ \Rightarrow \lambda = \| A^{-1} h_1 \| = \| A^{-1} h_2 \|
这样就求出了外参的各项参数。

极大似然估计

普通形式

张正友在论文中提到，前面的这些数学原理和推导并没有太多的物理意义，仅仅是为后面的极大似然优化提供了一个初值。
假设我们得到了模型平面的n幅图片，模型平面上有m个点，假设图像上像素点的噪声服从独立的同一分布，下面给出极大似然优化问题：

∑i=1n∑j=1m∥mij−m^(A,Ri,ti,Mj)∥

\sum_{i=1}^{n} \sum_{j=1}^{m} \| m_{ij} - \hat{m}(A,R_i,t_i,M_j) \|
其中： m^(A,Ri,ti,Mj)\hat{m}(A,R_i,t_i,M_j)表示的是点 MjM_j在第 ii幅图像上的投影，旋转矩阵RR使用有三个参数的向量 rr表示，rr平行于旋转轴，模长为旋转角度。这里的 RR和rr关系符合Rodrigue公式（参考前面的博客：三维重建学习(1)：基础知识：旋转矩阵与旋转向量）。
好了，这是一个非线性优化问题，他衡量的是点的实际坐标与估计值的误差，我们期望它尽可能地小。解决这个优化问题有很多种工具，但不是这里的重点，所以不做赘述。只是有一点需要注意，它需要一个 AA的初值，我们使用前面那一大堆公式导出一个猜测值赋给它。为什么这么大费周章？很简单，因为这比随机初始化初值收敛的快得多，有助于加快算法。

径向畸变的处理

如前面的博客：三维重建学习(2)：相机标定基础中所说，通常畸变有两种：径向畸变、切向畸变。通常切向畸变可以忽略不计，我们只考虑径向畸变。
令(u,v)(u,v)为点在理想的（无畸变的）像素坐标系中的坐标，令 (u~,v~)(\tilde{u}, \tilde{v})为实际观察到的坐标（有畸变）。同样地，令 (x,y)(x,y)表示理想的（无畸变）点在图像坐标系中的坐标， x~,y~\tilde{x}, \tilde{y}为实际观察到的坐标（有畸变）。
下面给出张正友在论文中给出的公式：

x~=x[1+k1(x2+y2)+k2(x2+y2)2]y~=y[1+k1(x2+y2)+k2(x2+y2)2]

\tilde{x} = x[1 + k_1(x^2 + y^2) + k_2(x^2 + y^2)^2] \\ \tilde{y} = y[1 + k_1(x^2 + y^2) + k_2(x^2 + y^2)^2]
k1k_1、 k2k_2是径向畸变参数。
从公式： u~=u0+αx~+cy~\tilde{u} = u_0 + \alpha \tilde{x} + c \tilde{y}、 v~=v0+βy~\tilde{v} = v_0 + \beta \tilde{y}（注：由内参矩阵的参数导出），得到：

u~=u+(u−u0)[k1(x2+y2)+k2(x2+y2)2]v~=v+(v−v0)[k1(x2+y2)+k2(x2+y2)2]

\tilde{u} = u + (u - u_0)[k_1(x^2 + y^2) + k_2(x^2 + y^2)^2] \\ \tilde{v} = v + (v - v_0)[k_1(x^2 + y^2) + k_2(x^2 + y^2)^2]
从上面的式子我们可以得到两个方程：

u~−u=(u−u0)[k1(x2+y2)+k2(x2+y2)2]v~−v=(v−v0)[k1(x2+y2)+k2(x2+y2)2]

\tilde{u} - u = (u - u_0)[k_1(x^2 + y^2) + k_2(x^2 + y^2)^2] \\ \tilde{v} - v = (v - v_0)[k_1(x^2 + y^2) + k_2(x^2 + y^2)^2]
表示成矩阵形式：

[(u−u0)(x2+y2)(v−v0)(x2+y2)(u−u0)(x2+y2)2(v−v0)(x2+y2)2]⋅[k1k2]=[u~−uv~−v]

\begin{bmatrix} (u - u_0)(x^2 + y^2) & (u - u_0) (x^2 + y^2)^2 \\ (v - v_0)(x^2 + y^2) & (v - v_0) (x^2 + y^2)^2 \end{bmatrix} \cdot \begin{bmatrix} k_1 \\ k_2 \end{bmatrix} = \begin{bmatrix} \tilde{u} - u \\ \tilde{v} - v \end{bmatrix}
nn幅图像中各有mm个点，迭代所有的方程会得到一个有 2mn2mn个方程的方程组。
用矩阵表示来表示这个方程组： D⋅k=dD \cdot k = d，其中 k=[k1k2]k = \begin{bmatrix} k_1 \\ k_2 \end{bmatrix}。
使用最小二乘法求解出 kk，直接套公式：

k=(DTD)−1DTd

k = (D^T D)^{-1} D^T d
如此求出了畸变参数后，我们可以回到前面的优化问题中再求解其他参数。

完整形式

下面给出完整的极大似然估计优化问题：

∑i=1n∑j=1m∥mij−m^(A,k1,k2,Ri,ti,Mj)∥

\sum_{i=1}^{n} \sum_{j=1}^{m} \| m_{ij} - \hat{m}(A,k_1,k_2,R_i,t_i,M_j) \|
其中： m^(A,ki,kj,Ri,ti,Mj) \hat{m}(A,k_i,k_j,R_i,t_i,M_j)表示的是点 MjM_j在第 ii幅图像上的投影。
与前面类似，这里依然回一个非线性优化问题。其中内参矩阵AA和 [Rt]\begin{bmatrix} R & t \end{bmatrix}的初始值选取采用前面的方法求出， k1k_1和 k2k_2的选取可以采用上一段的方法。
到这里，整个算法的数学原理已经推导完毕了。
在论文中，张正友还给出了相机标定的程序流程：

个人感觉用中文翻译意思上总是会有点偏差，姑且还是给出翻译后的流程：

建议采用如下校准步骤：
1. 打印一个图案，并把它贴到一个平面上；
2. 通过移动平面或者相机，从不同的方向拍摄一些图片；
3. 检测图片中的特征点；
4. 采用3.1节所述的方法，估计5个内参和全部的外参；
5. 使用最小二乘法(13)估计径向畸变系数；
6. 最小化优化问题(14)，改进所有的参数

这里用到的数学公式全部都在前面介绍了，理应很清楚了，不做赘述。下一篇博客再进行程序实现。

参考资料：

[图像]摄像机标定(2) 张正友标定推导详解
张正友标定法的真实理解
Flexible Camera Calibration By Viewing a Plane From Unknown Orientations（Zhengyou Zhang）

三维重建学习(3)：张正友相机标定推导相关推荐

张正友相机标定程序实现
版权声明:本文为博主原创文章,未经博主允许不得转载. https://blog.csdn.net/hongbin_xu/article/details/78988450 前言在前面的博客中( 三维重 ...
[毕设系列] 一、张正友相机标定
张正友相机标定预备知识 0.1 刚体.仿射.线性.旋转变换 0.2 什么是标定? 0.3 为什么要进行标定? 0.4 什么是畸变? 一.张正友标定法 1.1 简介 1.2 流程 1.3 畸变公式二 ...
OpenCV实现张正友相机标定源代码
本源代码基于VC++和opencv Opencv2.4.13.6版本开发,实现张正友相机标定源代码,资源包括完整源代码和12张棋盘图片,完美运行. Opencv2.4.13.6安装包下载地址:http ...
张正友相机标定Opencv实现以及标定流程标定结果评价图像矫正流程解析（附标定程序和棋盘图）
from:https://blog.csdn.net/dcrmg/article/details/52939318 使用Opencv实现张正友法相机标定之前,有几个问题事先要确认一下,那就是相机为什么 ...
张正友相机标定算法详解
张正友相机标定算法详解 1.齐次表示法与一些基本结论 1.1 点与直线的齐次表示在射影几何中,通常用齐次方式来表达点与直线.比如p=(u,v)p=(u,v)p=(u,v)被表示成p^=(x1,x ...
opencv C++ 张正友相机标定
//张正友相机标定 //https://blog.csdn.net/u010925447/article/details/77997735 #include <opencv2/core/core ...
张正友相机标定（全流程，含畸变，matlab源代码解析）
张正友标定的具体原理很多文章已经介绍,这里主要结合源代码对其中的基本原理及本人遇到的问题进行介绍.(仅介绍基本原理供本人复习,同时方便他人,如有问题,请及时指正勿喷) 1. 标定基本思路介绍相机标定 ...
张正友相机标定算法原理与源代码（OpenCV+C++）
摄像机的标定问题是机器视觉领域的入门问题,可以分为传统的摄像机定标方法和摄像机自定标方法.定标的方法有很多中常见的有:Tsai(传统)和张正友(介于传统和自定标)等, 摄像机成像模型和四个坐标系(通用 ...
【计算机视觉】张正友相机标定Calibration原理过程结果
一.数学原理 1.将世界坐标转换为图像坐标: 2.图像坐标转换为像素坐标系: 3.针孔成像下的透视投影矩阵: 4.世界坐标系转换为像素坐标系 5.张氏标定法的畸变模型: 以上参考:https://bl ...

三维重建学习(3)：张正友相机标定推导

前言