由于毕业论文需要，我翻译了张正友教授于1999年发表的论文 Flexible Camera Calibration by Viewing a Plane from Unknown Orientations，顺手发布出来，供大家参考，版权归原作者所有。
译文使用 LaTeX 编辑，全部源代码见下。
也可以点击此处下载Zhang_1999_translation.zip，压缩包内包含文献原文、.tex后缀的译文源代码、编译完成的.pdf后缀的译文等。
\documentclass[UTF8]{ctexart}
\usepackage[left=2.50cm, right=2.50cm, top=2.50cm, bottom=2.50cm]{geometry} %页边距
\usepackage{helvet}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amsfonts, amssymb} % 数学公式、符号
\usepackage{amsthm} % 根据 amsthm 的手册, amsthm 的加载要在 amsmath 之后
% 中文定理环境 % \indent 为了段前空两格
\newtheorem{theorem}{\indent 定理}[section]
\newtheorem{lemma}[theorem]{\indent 引理}
\newtheorem{proposition}[theorem]{\indent 命题}
\newtheorem{corollary}[theorem]{\indent 推论}
\newtheorem{definition}{\indent 定义}[section]
\newtheorem{example}{\indent 例}[section]
\newtheorem{remark}{\indent 注}[section]
\newenvironment{solution}{\begin{proof}[\indent\textbf 解]}{\end{proof}}
\renewcommand{\proofname}{\indent\textbf{证明}}
\renewcommand{\qedsymbol}{$\blacksquare$} % 证毕符号改成黑色正方形
% % English theorem environment
% \newtheorem{theorem}{Theorem}[section]
% \newtheorem{lemma}[theorem]{Lemma}
% \newtheorem{proposition}[theorem]{Proposition}
% \newtheorem{corollary}[theorem]{Corollary}
% \newtheorem{definition}{Definition}[section]
% \newtheorem{remark}{Remark}[section]
% \newtheorem{example}{Example}[section]
% \newenvironment{solution}{\begin{proof}[Solution]}{\end{proof}}
% \renewcommand{\qedsymbol}{$\blacksquare$} % 证毕符号改成黑色正方形
\usepackage{graphicx} % 图片
\usepackage{url} % 超链接
\usepackage{bm} % 加粗方程字体
\usepackage{multirow}
\usepackage{booktabs}
\usepackage{algorithm}
\usepackage{algorithmic}
\bibliographystyle{plain}
\renewcommand{\algorithmicrequire}{ \textbf{Input:}}
\renewcommand{\algorithmicensure}{ \textbf{Initialize:}}
\renewcommand{\algorithmicreturn}{ \textbf{Output:}}
% 算法格式
\usepackage{fancyhdr} % 设置页眉、页脚
\pagestyle{fancy}
\lhead{}
\chead{}
\lfoot{}
\cfoot{\thepage}
\rfoot{}
\usepackage[colorlinks,linkcolor=red,anchorcolor=blue,citecolor=green]{hyperref}
% \usepackage{multicol} % 多栏
\title{\textbf{一种灵活的从未知视角观测平面的相机标定方法}}
\author{\sffamily 张正友}
\date{}
\setlength{\headheight}{15pt}
\begin{document}
\maketitle
\noindent{\bf 摘要：}我们提出了一种用来标定相机的灵活的新方法。
这种方法只需要相机从很少的位置（至少两个）观测平面，而相机和平面都可以自由移动。
本方法已经考虑到了透镜畸变的影响。
本方法先给出了封闭解，然后再基于最大似然估计准则给出了理想的非线性解。
我们已经用计算机模拟数据和真实数据测试了这种方法，并取得了非常理想的结果。
常规方法需要昂贵的设备，比如用两个或者三个正交平面标定的方法，而本方法灵活易用。
它推动了三维计算机视觉从实验室向现实世界的普遍应用迈进了一步。
作者的主页上有本方法相应的软件。
\noindent{\bf 关键词：}相机标定；内参；镜头畸变；灵活的平面标定；运动分析；模型获取
% \begin{multicols}{2} % 双栏 % 由于本文公式较长，使用双栏排版会出现公式越界的问题，无奈只能使用单栏排版
\section{目的}
为了从二维图像中提取几何信息，在三维计算机视觉中，相机标定是必要步骤。
从摄影测量学\cite{duane1971close,faig1975calibration}到最近的计算机视觉\cite{gennery1979stereo,ganapathy1984decomposition,tsai1987versatile,faugeras1986calibration,weng1992camera,wei1993complete,maybank1992theory,faugeras1992camera,triggs1998autocalibration}，前人在这方面已经做过大量工作。
这些技术大致可以分为如下几类：
\begin{description}
\item[摄影测量法标定] 通过观测一个三维尺寸精确已知的物体来标定，这种方法\cite{faugeras1993three}非常高效。
标定的物体通常包含两个或者三个互相正交的平面。
有时，也用一个精确已知变化的平面来标定\cite{tsai1987versatile}。
这些方法需要昂贵的标定设备和复杂的步骤。
\item[自标定] 这种方法不需要任何辅助标定的对象。
只需要在静态场景里移动相机，通过图像信息来完成标定，这种场景对相机的内参有两个约束\cite{maybank1992theory}。
因此，如果用相同的相机和固定的内部参数来获取图像，三张图片足够得到内参和外参，以此我们可以重建出相似的三维结构\cite{luong1997self,hartley1994algorithm}。
虽然这种方法很灵活，但是还不成熟\cite{bougnoux1998projective}，因为有许多参数需要估计，我们无法得到可靠的结果。
\item[其它方法] 现有：用正交方向的灭点进行标定的方法\cite{caprile1990using,liebowitz1998metric}和纯粹的旋转标定方法\cite{hartley1994self,stein1995accurate}。
\end{description}
我们当前的研究集中在桌面视觉系统(DVS)，因为随着相机变得越来越便宜和常见，桌面视觉系统具有极好的应用前景。
桌面视觉系统主要针对那些不擅长计算机视觉的普通大众，典型的计算机用户只是偶尔用一下视觉功能，所以他们不会愿意为昂贵的设备投资。
因此灵活性、鲁棒性、低成本很重要，而本文提出的相机标定方法考虑到了这些事项。
本方法只需要相机从不同方位（至少两个）来观测标定板。
标定板的图案可以是用激光打印机打印出来并贴到一个合适的平面上（比如一个硬的图书封面），相机或标定板可以自由移动。
本文提出的方法介于摄影测量法标定和自标定方法之间，因为我们利用了二维几何信息而不是三维或者纯粹的特征信息。
我们已经用计算机模拟数据和真实数据测试了这种方法，并取得了非常理想的结果。
与经典的摄影测量法标定相比，该方法非常灵活。
与自标定相比，本方法具有较强的鲁棒性。
我们相信本方法使三维计算机视觉从实验室向现实世界的普遍应用迈进了一步。
Bill Trigger\cite{triggs1998autocalibration}最近提出了一种新的自标定方法，这种方法至少需要一个平面图案场景的五个视图。
他的技术比我们的更灵活，但是初始化很困难。
Liebowitz和Zisserman\cite{liebowitz1998metric}描述了一种对标定平面的透视图进行度量矫正的技术。
这种方法利用度量信息，比如一个已知的角度或者两个相等但是未知的角度，和已知的长度比率。
虽然没有算法和实验结果，但是他们仍认为基于这些矫正过的平面，相机内参的标定是可能的。
本文结构如下。
第二部分描述了观测一个标定平面的基本约束条件。
第三部分描述了标定步骤。
我们以一个封闭解开始，然后得到非线性的最优化，同时也考虑到了透镜畸变的影响。
第四部分研究了这个方法的不足之处，而在实践中很容易避免这些缺点。
第五部分提供了实验结果，计算机模拟出的数据和现实数据都用来验证了该方法。
附录中，我们给出了大量详细的说明，包括了估计位图模型与其图像单应性关系的技术。
\section{基本方程}
我们通过观测一个待标定的平面来得到相机内参的约束条件。本文从以下符号开始。
\subsection{符号}
二维点用$\mathbf{m}=\left[ u,v \right] ^T$，三维点用$M=\left[ X,Y,Z \right] ^T$表示。
$\widetilde{x}$表示增广向量，方法是在矩阵中加入1作为最后一个元素：$\widetilde{\mathbf{\mathbf{m}}}=\left[ u,v,1 \right] ^T$，$\widetilde{M}=\left[ X,Y,Z,1 \right] ^T$。
建立一个常见的小孔成像模型，三维点$M$和它图像投影点\textbf{m}关系如下：
\begin{equation}
\label{equation:1}
s\widetilde{\mathbf{m}}=\mathbf{A}\left[ \begin{matrix}
\mathbf{R} & \mathbf{t} \\
\end{matrix} \right] \widetilde{M}\text{，其中}\mathbf{A}=\left[ \begin{matrix}
\alpha & c & u_0 \\
0 & \beta & v_0 \\
0 & 0 & 1 \\
\end{matrix} \right]
\end{equation}
式中$s$是任意的比例因子。
$\left(\mathbf{R},\mathbf{t}\right)$称为外参，$\mathbf{R}$是旋转矩阵，$\mathbf{t}$是平移矩阵。
$\mathbf{A}$是相机内参矩阵，$\left( u_0,v_0 \right)$是坐标的主点，$\alpha$和$\beta$是图像在$u$和$v$轴的比例因子，$c$是描述两个坐标轴倾斜角的参数。
我们简称$\mathbf{A}^{−T}$为$\left(\mathbf{A}^{−1}\right) ^T$或者$\left(\mathbf{A}^{T}\right)^{-1}$。
\subsection{模型平面与其图像间的单应性关系}
一般情况下， 我们假设模型平面在世界坐标系中的$Z$坐标为0。
$\mathbf{R}$的第$i$列旋转矩阵为$\mathbf{r}_i$。
由公式(\ref{equation:1})可知:
\begin{equation*}
s\left[ \begin{array}{l}
u \\
v \\
1 \\
\end{array} \right] =\mathbf{A}\left[ \begin{matrix}
\mathbf{r}_1 & \mathbf{r}_2 & \mathbf{r}_3 & \mathbf{t} \\
\end{matrix} \right] \left[ \begin{array}{l}
X \\
Y \\
0 \\
1 \\
\end{array} \right] =\mathbf{A}\left[ \begin{matrix}
\mathbf{r}_1 & \mathbf{r}_2 & \mathbf{t} \\
\end{matrix} \right] \left[ \begin{array}{l}
X \\
Y \\
1 \\
\end{array} \right]
\end{equation*}
我们仍然用$M$表示位于位图中的一点，当$Z=0$，$M$可以表示为$M=\left[ X,Y \right] ^T$，同样地，$\widetilde{M}=\left[ X,Y,1 \right] ^T$。
因此点$M$和它在图像上映射点$\mathbf{m}$的关系可以用单应矩阵$\mathbf{H}$表示：
\begin{equation}
\label{equation:2}
s\widetilde{\mathbf{m}}=\mathbf{H}\widetilde{M}\text{，其中}\mathbf{H}=\mathbf{A}\left[ \mathbf{r}_1\ \mathbf{r}_2\ \mathbf{t} \right]
\end{equation}
很显然，$\mathbf{H}$是一个$3\times 3$的系数矩阵。
\subsection{内参的约束条件}
对标定平面的一张图片，它的单应矩阵可以估计出来（参见附录\ref{section:A}）。
令$\mathbf{H}=\left[ \mathbf{h}_1\ \mathbf{h}_2\ \mathbf{h}_3 \right]$，由公式(\ref{equation:2})可知：
\begin{equation*}
\left[ \mathbf{h}_1\ \mathbf{h}_2\ \mathbf{h}_3 \right] =\lambda \mathbf{A}\left[ \mathbf{r}_1\ \mathbf{r}_2\ \mathbf{t} \right]
\end{equation*}
式中$\lambda$是任意的标量。
由所学知识可知$\mathbf{r_1}$和$\mathbf{r_2}$是正交的，于是我们有：
\begin{equation}
\label{equation:3}
\mathbf{h}_{1}^{T}\mathbf{A}^{−T}\mathbf{A}^{−1}\mathbf{h}_2=0
\end{equation}
\begin{equation}
\label{equation:4}
\mathbf{h}_{1}^{T}\mathbf{A}^{−T}\mathbf{A}^{−1}\mathbf{h}_1=\mathbf{h}_{2}^{T}\mathbf{A}^{−T}\mathbf{A}^{−1}\mathbf{h}_2
\end{equation}
对一个给定的单应性矩阵，对内参有两个基本的约束条件。
因为这个矩阵有8个自由度和6个外参（其中有3个外参是旋转矩阵的，另外3个外参是平移向量的），这个条件下，我们只能得到两个内参的约束条件。
\section{相机标定的解决方法}
该节详细描述了怎样有效地解决相机的标定问题。
我们先给出一个封闭解，然后根据最大似然估计给出非线性的最优化解，最后再考虑透镜的径向畸变，给出解析解和非线性解。
\subsection{封闭解}\label{section:3.1}
令：
\begin{equation}
\label{equation:5}
\begin{aligned}
\mathbf{B} & =\mathbf{A}^{-T}\mathbf{A}^{-1}\equiv \left[ \begin{matrix}
B_{11} & B_{12} & B_{13} \\
B_{12} & B_{22} & B_{23} \\
B_{13} & B_{23} & B_{33} \\
\end{matrix} \right] \\
& =\left[ \begin{matrix}
\frac{1}{\alpha ^2} & -\frac{c}{\alpha ^2\beta} & \frac{cv_0-u_0\beta}{\alpha ^2\beta} \\
-\frac{c}{\alpha ^2\beta} & \frac{c^2}{\alpha ^2\beta ^2}+\frac{1}{\beta ^2} & -\frac{c\left( cv_0-u_0\beta \right)}{\alpha ^2\beta ^2}-\frac{v_0}{\beta ^2} \\
\frac{cv_0-u_0\beta}{\alpha ^2\beta} & -\frac{c\left( cv_0-u_0\beta \right)}{\alpha ^2\beta ^2}-\frac{v_0}{\beta ^2} & \frac{\left( cv_0-u_0\beta \right) ^2}{\alpha ^2\beta ^2}+\frac{v_{0}^{2}}{\beta ^2}+1 \\
\end{matrix} \right] \\
\end{aligned}
\end{equation}
注意$\mathbf{B}$是对称的，定义一个六维向量：
\begin{equation}
\label{equation:6}
\mathbf{b}=\left[ B_{11},B_{12},B_{22},B_{13},B_{23},B_{33} \right] ^T
\end{equation}
（它实际上描述了绝对二次曲线的图像。）
假设$\mathbf{H}$的第$i$列向量为$\mathbf{h}_i=\left[ h_{i1},h_{i2},h_{i3} \right] ^T$，然后可以得到：
\begin{equation}
\label{equation:7}
\mathbf{h}_{i}^{T}\mathbf{Bh}_j=\mathbf{v}_{ij}^{T}\mathbf{b}
\end{equation}
其中，$\mathbf{v}_{ij}=\left[ h_{i1}h_{j1},h_{i1}h_{j2}+h_{i2}h_{j1},h_{i2}h_{j2}, \right. \left. h_{i3}h_{j1}+h_{i1}h_{j3},h_{i3}h_{j2}+h_{i2}h_{j3},h_{i3}h_{j3} \right] ^T$。
因此，两个基本约束公式(\ref{equation:3})和(\ref{equation:4})可以重写为$\mathbf{b}$中的齐次式：
\begin{equation}
\label{equation:8}
\left[\begin{array}{c}
\mathbf{v}_{12}^{T} \\
\left(\mathbf{v}_{11}-\mathbf{v}_{22}\right)^{T}
\end{array}\right] \mathbf{b}=\mathbf{0}
\end{equation}
如果观测了n张图片，迭代可以得到：
\begin{equation}
\label{equation:9}
\mathbf{Vb}=0
\end{equation}
其中$\mathbf{V}$是一个$2n\times6$的矩阵。
若$n\geqslant 3$，通常会得到一个唯一解$\mathbf{b}$。
若$n = 2$，令倾斜约束参数$c = 0$， 即$\left[ 0,1,0,0,0,0 \right] \mathbf{b}=0$，这作为额外的方程代入公式(\ref{equation:9})。
公式(\ref{equation:9})的解是我们熟知的与最小特征值相关的特征向量$\mathbf{V}^T\mathbf{V}$（等价地，$\mathbf{V}$的右奇异向量与最小奇异值相关）。
一旦$\mathbf{b}$被估计出来，我们可以计算内参矩阵$\mathbf{A}$。
详情请参见附录\ref{section:B}。
已知$\mathbf{A}$后，相应的外部参数也可以被计算出来，由公式(\ref{equation:2})可以得到：
\begin{equation*}
\mathbf{r}_{1}=\lambda \mathbf{A}^{-1} \mathbf{h}_{1}, \mathbf{r}_{2}=\lambda \mathbf{A}^{-1} \mathbf{h}_{2}, \mathbf{r}_{3}=\mathbf{r}_{1} \times \mathbf{r}_{2}, \mathbf{t}=\lambda \mathbf{A}^{-1} \mathbf{h}_{3}
\end{equation*}
其中，$\lambda=1 /\left\|\mathbf{A}^{-1} \mathbf{h}_{1}\right\|=1 /\left\|\mathbf{A}^{-1} \mathbf{h}_{2}\right\|$。
当然，由于噪声数据，因此计算的矩阵$\mathbf{R}=\left[\mathbf{r}_{1}, \mathbf{r}_{2}, \mathbf{r}_{3}\right]$一般不能满足一个旋转矩阵的属性。
附录\ref{section:C}以一般的$3\times3$矩阵为例介绍了一种方法来估计最佳的旋转矩阵。
\subsection{最大似然估计}\label{section:3.2}
上面的解决方案通过最小化代数距离获得的，一般情况下这是没有物理意义的。
我们可以通过最大似然估计理论来完善。
我们得到模型平面的$n$幅图片，模型平面上有$m$个点。
假设图像上像素点的噪声服从独立的同一分布。
最大似然估计可以从通过求以下函数的最小值得到：
\begin{equation}
\label{equation:10}
\sum_{i=1}^{n} \sum_{j=1}^{m}\left\|\mathbf{m}_{i j}-\hat{\mathbf{m}}\left(\mathbf{A}, \mathbf{R}_{i}, \mathbf{t}_{i}, \mathrm{M}_{j}\right)\right\|^{2}
\end{equation}
根据公式(\ref{equation:2})，上式中$\hat{\mathbf{m}}\left(\mathbf{A}, \mathbf{R}_{i}, \mathbf{t}_{i}, \mathrm{M}_{j}\right)$是$\mathrm{M}_{j}$点在第$i$幅图像上的投影。
旋转矩阵$\mathbf{R}$用三个参数的向量$\mathbf{r}$表示，$\mathbf{r}$平行于旋转轴，而且大小等于旋转角度。
$\mathbf{R}$和$\mathbf{r}$关系符合罗德里格斯公式\cite{faugeras1993three}。
求公式(\ref{equation:10})的最小值是一个非线性优化问题，可以通过Levenberg-Marquardt算法\cite{mor1977levenberg}解决。
它需要矩阵$\mathbf{A}$的一个初始的猜测值，$\left\{\mathbf{R}_{i}, \mathbf{t}_{i} \mid i=1 . . n\right\}$可使用在上一小节中描述的方法得到。
\subsection{径向畸变的处理}
到现在为止，我们没有考虑过相机的镜头畸变。
然而，桌面相机的镜头通常会有显著的畸变，尤其是径向畸变。
在本节中，我们只考虑前两个方面的径向畸变。
读者可以参考文献\cite{slama1980manual,duane1971close,faig1975calibration,weng1992camera}获取详细信息。
根据报告\cite{duane1971close,tsai1987versatile,wei1994implicit}，镜头的畸变以径向分量为主，尤其是第一项占主导地位。
同时文献\cite{tsai1987versatile,wei1994implicit}发现任何更复杂的模型不仅没有帮助，而且会导致数值不稳定。
设$(u,v)$是理想（畸变可忽略）像素的图像坐标，$(\breve{u}, \breve{v})$对应的实际观测到的图像坐标。
同样，$(x,y)$和$(\breve{x}, \breve{y})$是理想（无畸变）和实时（含有畸变）情况下的归一化图像坐标。
我们有\cite{duane1971close,wei1994implicit}：
\begin{equation*}
\breve{x}=x+x\left[k_{1}\left(x^{2}+y^{2}\right)+k_{2}\left(x^{2}+y^{2}\right)^{2}\right]
\end{equation*}
\begin{equation*}
\breve{y}=y+y\left[k_{1}\left(x^{2}+y^{2}\right)+k_{2}\left(x^{2}+y^{2}\right)^{2}\right]
\end{equation*}
其中，$k_1$和$k_2$是径向畸变系数。
径向畸变中心的主点是相同的。
从$\breve{u}=u_{0}+\alpha \breve{x}+c \breve{y}$和$\breve{v}=v_{0}+\beta \breve{y}$中可以得到：
\begin{equation}
\label{equation:11}
\breve{u}=u+\left( u-u_0 \right) \left[ k_1\left( x^2+y^2 \right) +k_2\left( x^2+y^2 \right) ^2 \right]
\end{equation}
\begin{equation}
\label{equation:12}
\breve{v}=v+\left( v-v_0 \right) \left[ k_1\left( x^2+y^2 \right) +k_2\left( x^2+y^2 \right) ^2 \right]
\end{equation}
\begin{description}
\item[交替径向畸变的估计]径向畸变一般很小，简单的忽略径向畸变后，可以用\ref{section:3.2}节中的方法估计另外的五个参数。
一种策略是估计其它参数后再估计$k_1$和$k_2$。
然后由公式(\ref{equation:11})和(\ref{equation:12})，对于每幅图的每个点，我们可以得到两个方程：
\begin{equation*}
\left[\begin{array}{cc}
\left(u-u_{0}\right)\left(x^{2}+y^{2}\right) & \left(u-u_{0}\right)\left(x^{2}+y^{2}\right)^{2} \\
\left(v-v_{0}\right)\left(x^{2}+y^{2}\right) & \left(v-v_{0}\right)\left(x^{2}+y^{2}\right)^{2}
\end{array}\right]\left[\begin{array}{l}
k_{1} \\
k_{2}
\end{array}\right]=\left[\begin{array}{l}
\breve{u}-u \\
\breve{v}-v
\end{array}\right]
\end{equation*}
$n$幅图像中共有$m$个点，叠加所有方程可以得到$2mn$个方程，或者用矩阵表示：$\mathbf{Dk}=\mathbf{d}$，其中$\mathbf{k}=\left[k_{1}, k_{2}\right]^{T}$，由最小二乘法得：
\begin{equation}
\label{equation:13}
\mathbf{k}=\left(\mathbf{D}^{T} \mathbf{D}\right)^{-1} \mathbf{D}^{T} \mathbf{d}
\end{equation}
$k_1$和$k_2$估计出来后，可以通过求解公式(\ref{equation:10})来重新估计其它参数，其中可以用公式(\ref{equation:11})和(\ref{equation:12})来代替$\hat{\mathbf{m}}\left(\mathbf{A}, \mathbf{R}_{i}, \mathbf{t}_{i}, \mathrm{M}_{j}\right)$。
交替使用这两个步骤，直到收敛为止。
\item[完整的最大似然估计]实验中，我们发现上述收敛过程比较慢。
根据公式(\ref{equation:10})，可以通过最小化下列函数来估计参数的完整集：
\begin{equation}
\label{equation:14}
\sum_{i=1}^{n} \sum_{j=1}^{m}\left\|\mathbf{m}_{i j}-\breve{\mathbf{m}}\left(\mathbf{A}, k_{1}, k_{2}, \mathbf{R}_{i}, \mathbf{t}_{i}, \mathrm{M}_{j}\right)\right\|^{2}
\end{equation}
根据公式(\ref{equation:2})、(\ref{equation:11})和(\ref{equation:12})，其中,$\breve{\mathbf{m}}\left(\mathbf{A}, k_{1}, k_{2}, \mathbf{R}_{i}, \mathbf{t}_{i}, \mathrm{M}_{j}\right)$是点$M_j$在图像$i$中的投影。
这是一个非线性最小化问题，可以通过Levenberg-Marquardt算法\cite{mor1977levenberg}解决。
旋转矩阵仍然还可用一个\ref{section:3.2}节中所示的三维向量$\mathbf{r}$来表示。
$\mathbf{A}$初始值的选取和$\left\{\mathbf{R}_{i}, \mathbf{t}_{i} \mid i=1 . . n\right\}$可以采用\ref{section:3.1}节或者\ref{section:3.2}节中描述的技术。
$k_1$和$k_2$初始值的选取可采用上一段所述的方法，或着直接将它们设置为0。
\end{description}
\subsection{小结}
建议校准程序如下：
\begin{enumerate}
\item 打印图案，并将它贴到一个平面；
\item 移动相机或者平面，从不同角度拍下几张图片；
\item 在图像中检测特征点；
\item 利用\ref{section:3.1}节中所述方法，估计五个内参和全部的外参；
\item 通过公式(\ref{equation:13})估计径向畸变系数;
\item 通过最小化公式(\ref{equation:14})改进全部参数。
\end{enumerate}
\section{退化情形}
本节我们研究不满足相机内参约束的退化情形。
因为公式(\ref{equation:3})和(\ref{equation:4})是从旋转矩阵的性质推导出来的，如果$\mathbf{R}_2$和$\mathbf{R}_1$不相互独立，那么第二幅图对于相机标定来说不起作用。
下面我们将考虑更复杂的情况。
\begin{proposition}
如果第二个位置的模型平面平行于第一个位置的模型平面，那么第二幅图片不提供额外的约束条件。
\end{proposition}
由于篇幅原因证明过程省略，具体细节可以参见技术文档\cite{zhang1998flexible}。
实际过程中，非常容易避免退化：只需要改变平面与相机拍照点的相对位置即可。
\section{实验结果}
本文提出的算法测试了计算机模拟数据和真实数据。
封闭解是来分解一个$2n\times 6$的矩阵从而求得奇异解，其中$n$是图片数。
用Levenberg-Marquardt算法求非线性解需要3到5次迭代才收敛。
\subsection{计算机模拟}
模拟的相机具有如下性质：$\alpha=1250$，$\beta=900$，$c= 1.09083$（相当于$89.95°$），$u_0=255$，$v_0=255$。
图像分辨率为$512\times 512$像素。
标定板包含$10\times 14 = 140$个角点，尺寸是18cm$\times $25cm。平面的旋转矩阵用三维向量$\mathbf{r}$来表示，平移矩阵用一个三维向量$\mathbf{t}$（米制单位）表示。
\begin{description}
\item[关于噪声水平的性能] 实验中，我们用三个平面$\mathbf{r}_{1}=\left[20^{\circ}, 0,0\right]^{T}$，$\mathbf{t}_{1}=\left[-9, -12.5,500\right]^{T}$，$\mathbf{r}_{2}=\left[0, 20^{\circ},0\right]^{T}$，$\mathbf{t}_{2}=\left[-9, -12.5,510\right]^{T}$，$\mathbf{r}_{3}=\frac{1}{\sqrt{5}}\left[-30^{\circ},-30^{\circ},-15^{\circ}\right]^{T}$，$\mathbf{t}_{3}=\left[-10.5, -12.5,525\right]^{T}$。
将期望为0、标准差为$\sigma$的高斯噪声加入投影图像，对比估计出的相机参数与真值。
我们测试了$\alpha$和$\beta$的相对误差、$u_0$和$v_0$的绝对误差。
加入的噪声水平分布在0.1像素到1.5像素之间，对每个噪声等级，做了100次独立试验，并得到平均结果。
从图\ref{figure:1}可以看到，误差随噪声水平呈线性增长。
\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figure/figure1.jpg}
\caption{不同噪声水平标定的误差对比} \label{figure:1}
\end{figure}
当$\sigma=0.5$（大于实际噪声的标准差）时，$\alpha$和$\beta$的误差小于0.3\%，$u_0$和$v_0$误差在1像素附近。
$u_0$的误差大于$v_0$的误差，主要原因是$u$轴方向上的数据比$v$方向上数据少。
\item[不同平面数的性能测试] 该实验探索了不同模型平面数与性能的关系。
前三幅图像的方向和位置和上一段的相同。
第四幅图是在均匀的球体中任选一个旋转轴，然后再旋转$30^{\circ}$得到。
图片数从2到16不等，对每个图片数，做了100次独立实验，实验中噪声相互独立，并服从期望为0，标准差为0.5像素的高斯分布。
图\ref{figure:2}显示了平均结果，图示可以看出，平面数从2到3，误差水平显著减小。
\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figure/figure2.jpg}
\caption{不同平面数标定误差对比} \label{figure:2}
\end{figure}
\item[模型平面方向的性能测试] 实验研究了位面方向对于标定结果的影响。
实验用了三幅图片。
平面方位选取如下：开始模型平面与图像平面平行，旋转轴从均匀球体中随机选取，然后绕轴旋转$\theta$角， 图像中加入期望为0、标准差为0.5 像素的高斯噪声。
重复测试100次，计算平均误差。
$\theta$变化范围为$5^{\circ}$到$75^{\circ}$，结果如图\ref{figure:3}所示。
当$\theta$为$0.5^{\circ}$时，40\%的实验因为平面之间几乎是相互平行而失败，实验结果中已经包含了这种情况。
角度为$45^{\circ}$左右时，效果似乎最佳。
值得注意的是，实践中当角度增加时，透影缩减现象让角点检测变得困难，但是实验中并未考虑这点。
\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figure/figure3.jpg}
\caption{不同角度标定的误差对比} \label{figure:3}
\end{figure}
\end{description}
\subsection{真实数据测试}
本文提出的方法已在我们的视觉组和图像组得到普遍应用。
在这里，我们以一实际例子来说明标定过程。
标定PULNiX牌镜头焦距为6mm的CCD相机，图像分辨率为$640\times 480$。
标定板含$8\times 8$的正方形，因此有256个角点。
标定板尺寸为17cm$\times$17cm，得到五幅不同方位的图像，如图\ref{figure:4}所示，图中可以看到明显的畸变。
\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figure/figure4.jpg}
\caption{标定平面的五幅图片以及提取的角点（角点用十字表示，但太小而无法观察）} \label{figure:4}
\end{figure}
实验中，分别采用不同数量的图像来标定，即2、3、4、5幅图像。
结果如表\ref{table:1}所示，对每种情形，第一列是对封闭解的估计，第二列是最大似然估计的结果，第三列是标准差的估计，以此表示结果的不确定性。
可以清楚的看到，封闭解比较合理，最终的最大似然估计值也具有很好的一致性。
同时可以看出随着图片数的增多，估计结果的不一致性降低了。
表\ref{table:1}的最后一行是检测出的角点和投影点的均方根距离，单位是像素，可知最大似然估计使结果提高了很多。
\begin{table}[htbp]
\centering
\caption{用2幅图像到5幅的标定数据}
\label{table:1}
\resizebox{\textwidth}{!}{
\begin{tabular}{c|ccc|ccc|ccc|ccc}
\hline
\multirow{2}{*}{nb} & \multicolumn{3}{c|}{2 images} & \multicolumn{3}{c|}{3 images} & \multicolumn{3}{c|}{4 images} & \multicolumn{3}{c}{5 images} \\ \cline{2-13}
& initial & final & $\sigma$ & initial & final & $\sigma$ & initial & final & $\sigma$ & initial & final & $\sigma$ \\ \hline
$\alpha$ & 825.59 & 830.47 & 4.74 & 917.65 & 830.80 & 2.06 & 876.62 & 831.81 & 1.56 & 877.16 & 832.50 & 1.41 \\
$\beta$ & 825.26 & 830.24 & 4.85 & 920.53 & 830.69 & 2.10 & 876.22 & 831.82 & 1.55 & 876.80 & 832.53 & 1.38 \\
$c$ & 0 & 0 & 0 & 2.2956 & 0.1676 & 0.109 & 0.0658 & 0.2867 & 0.095 & 0.1752 & 0.2045 & 0.078 \\
$u_0$ & 295.79 & 307.03 & 1.37 & 277.09 & 305.77 & 1.45 & 301.31 & 304.53 & 0.86 & 301.04 & 303.96 & 0.71 \\
$v_0$ & 217.69 & 206.55 & 0.93 & 223.36 & 206.42 & 1.00 & 220.06 & 206.79 & 0.78 & 220.41 & 206.56 & 0.66 \\
$k_1$ & 0.161 & -0.227 & 0.006 & 0.128 & -0.229 & 0.006 & 0.145 & -0.229 & 0.005 & 0.136 & -0.228 & 0.003 \\
$k_2$ & -1.955 & 0.194 & 0.032 & -1.986 & 0.196 & 0.034 & -2.089 & 0.195 & 0.028 & -2.042 & 0.190 & 0.025 \\ \hline
RMS & 0.761 & \multicolumn{2}{c|}{0.295} & 0.987 & \multicolumn{2}{c|}{0.393} & 0.927 & \multicolumn{2}{c|}{0.361} & 0.881 & \multicolumn{2}{c}{0.335} \\ \hline
\end{tabular}}
\end{table}
细心的读者可能发现了封闭解和最大似然估计解中的$k_1$和$k_2$的一致性不好。
原因在于解封闭解时，相机内参估计的前提是假设没有畸变，预测点更靠近图像中心而不是检测值，而随后对畸变的估计试着延伸外部点，并且增大系数以此减小距离，使畸变形状和真实畸变不一致。
最大似然估计最终覆盖了畸变形态。
估计出的畸变参数让我们校正原图像，图\ref{figure:5}展示了校正后的第1、2幅图像，可以看出方框之间的线条变直了。
\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figure/figure5.jpg}
\caption{标定畸变后的第一幅和第二幅图片} \label{figure:5}
\end{figure}
\begin{description}
\item[标定结果的变化] 表\ref{table:1}中，我们展示了用2到5张图片标定的结果，而且我们发现结果的一致性非常好。
为了进一步探究该算法的稳定性，我们从中任取四张图片，再运用该方法标定。
结果如表\ref{table:2}所示。
例如第三列，实现了用第1、2、3、5幅图片标定的结果。
最后两列代表均值和标准差。
代表倾斜度的参数$c$在0附近，偏差系数$0.086/0.1401 = 0.6$，是一个相对较大的数。
$c = 0.1401$，$\alpha = 832.85$对应$89.99°$，非常接近$90°$，我们同样计算了比率$\alpha/\beta $，平均值为0.99995，偏差为0.00012，非常接近1，即标定后图案是正方形的。
\item[基于图像的建模] 用上文中提到的相机拍摄一个茶叶罐的两幅图片，如图\ref{figure:6}所示。
\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figure/figure6.jpg}
\caption{茶叶罐的两幅图片} \label{figure:6}
\end{figure}
茶叶罐的两边是可见的，我们在每边上手动挑选8个点，然后用我们之前开发的运动建模软件得到该茶叶罐的部分模型。
该重构模型是用虚拟现实建模语言(VRML)实现的，三维渲染细节如图\ref{figure:7}所示，每边重构出的点确实共面，计算重构的两个平面的夹角为$94.7°$。
\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{figure/figure7.jpg}
\caption{茶叶罐重构后模拟的三面} \label{figure:7}
\end{figure}
虽然我们没有确切值，但是茶叶罐的两面确实几乎正交。
\end{description}
真实的数据、结果和软件可以从下列网页获得：
\url{http://research.microsoft.com/p~zhang/Calib/}
\begin{table}[htbp]
\centering
\caption{用四幅不同图片标定结果对比}
\label{table:2}
\resizebox{0.8\textwidth}{!}{
\begin{tabular}{c|ccccc|cc}
\hline
quadruple & (1234） & (1235） & (1245） & (1345） & (2345） & mean & deviation \\ \hline
$\alpha$ & 831.81 & 832.09 & 837.53 & 829.69 & 833.14 & 832.85 & 2.90 \\
$\beta$ & 831.82 & 832.10 & 837.53 & 829.91 & 833.11 & 832.90 & 2.84 \\
$c$ & 0.2867 & 0.1069 & 0.0611 & 0.1363 & 0.1096 & 0.1401 & 0.086 \\
$u_0$ & 304.53 & 304.32 & 304.57 & 303.95 & 303.53 & 304.18 & 0.44 \\
$v_0$ & 206.79 & 206.23 & 207.30 & 207.16 & 206.33 & 206.76 & 0.48 \\
$k_1$ & -0.229 & -0.228 & -0.230 & -0.227 & 0.229 & 0.229 & 0.001 \\
$k_2$ & 0.195 & 0.191 & 0.193 & 0.179 & 0.190 & 0.190 & 0.006 \\ \hline
RMS & 0.361 & 0.357 & 0.262 & 0.358 & 0.334 & 0.334 & 0.04 \\ \hline
\end{tabular}}
\end{table}
\section{结论}
本文提出了一种新的灵活标定相机的方法。
本方法只需要相机从不同方向（至少两个）观测标定板，可以任意移动相机或者标定板，同时，本方法已经考虑到了透镜畸变的影响。
本方法先给出封闭解，然后根据最大似然估计给出非线性优化结果。
我们已经用计算机模拟数据和真实数据测试了这种方法，并取得了非常理想的结果。
与经典方法相比，该方法相当灵活。
\begin{description}
\item[致谢] 感谢Brian Guenter的关于角点提取的软件以及Bill Triggs具有洞察力的评论。
感谢Andrew Zisserman的CVPR98成果\cite{liebowitz1998metric}，它使用相同的约束但不同的形式。
感谢Bill Triggs和Gideon Stein给实验提出的宝贵意见\cite{zhang1998flexible}。
感谢Anandan和Charles Loop校对了本文英语。
\end{description}
\appendix
\section{标定平面与其图像的单应性估计}
\label{section:A}
有多种方法来估计标定平面与其图像的单应性关系，本文我们用最大似然估计来处理。
假设$M_i$和$\mathbf{m}_i$分别是标定平面和图像上的点，理想情形下它们符合公式(\ref{equation:2})，但实际上，因为图像上存在噪声，所以不符合上述假设。
假设$\mathbf{m}_i$上加均值为0、协方差矩阵为$\Lambda _{\mathbf{m}_i}$的高斯噪声，通过求下式的最小值，得到$\mathbf{H}$。
\begin{equation*}
\sum_i{\left( \mathbf{m}_i-\mathbf{\hat{m}}_i \right)}^T\Lambda _{\mathbf{m}_i}^{-1}\left( \mathbf{m}_i-\mathbf{\hat{m}}_i \right) \text{，}
\end{equation*}
\begin{equation*}
\text{其中}\mathbf{\hat{m}}_i=\frac{1}{\overline{\mathbf{h}}_{3}^{T}\text{M}_i}\left[ \begin{array}{l}
\overline{\mathbf{h}}_{1}^{T}\text{M}_i \\
\overline{\mathbf{h}}_{2}^{T}\text{M}_i \\
\end{array} \right] \text{，}\overline{\mathbf{h}}_i\text{是}\mathbf{H}\text{的第}i\text{行}
\end{equation*}
在工程实际中，我们简单的假设对于所有的$i$有$\Lambda_{\mathbf{m}_{i}}=\sigma^{2} \mathbf{I}$。
如果图像的点在同一过程中被独立的提取出来，这种假设就不合理了。
在这种情况下，以上的问题便转化成一个非线性最小二乘法的问题，例如$\min {\mathbf{H}} \sum_{i}\left\|\mathbf{m}_{i}-\hat{\mathbf{m}}_{i}\right\|^{2}$，这种非线性最小化可以用Levenberg-Marquardt算法\cite{mor1977levenberg}来实现，这需要一个初始猜测值，它可用如下方法得到：
令$\mathbf{x}=\left[\overline{\mathbf{h}}_{1}^{T}, \overline{\mathbf{h}}_{2}^{T}, \overline{\mathbf{h}}_{3}^{T}\right]^{T}$，公式(\ref{equation:2})可以写为：
\begin{equation*}
\left[\begin{array}{ccc}
\widetilde{\mathrm{M}}^{T} & \mathbf{0}^{T} & -u \widetilde{\mathrm{M}}^{T} \\
0^{T} & \widetilde{\mathrm{M}}^{T} & -v \widetilde{\mathrm{M}}^{T}
\end{array}\right] \mathbf{x}=0
\end{equation*}
当已知$n$个点，就可得到$n$个上述方程，它们可以写成矩阵方程，如：$\mathbf{L x}=\mathbf{0}$，这里$\mathbf{L}$是一个$2n\times 9$的矩阵。
定义$\mathbf{x}$为一个比例因子，这种方法被称作伴有最小奇异值的$\mathbf{L}$矩阵的右奇异向量（等价于矩阵$\mathbf{L}^T \mathbf{L}$最小特征值对应的特征向量）。
$\mathbf{L}$中的元素一部分是常数1，一部分在图像坐标系，一部分在世界坐标系里，一部分是两者的相乘，这就使得$\mathbf{L}$矩阵没有充分的数字化。
可以在运行以上程序之前将数据归一化，可以得到很好的结果，例如在\cite{hartley1995defence}中提出的方法。
\section{从\texorpdfstring{$\mathbf{B}$}{}矩阵提取内部参数}
\label{section:B}
在\ref{section:3.1}节中描述的矩阵$\mathbf{B}$被估计为一个比例因子，即$\mathbf{B}=\lambda \mathbf{A}^{-T} \mathbf{A}$，$\lambda$为一个任意的因子。
我们可以很容易的从矩阵$\mathbf{B}$中提取出内部参数。
\begin{equation*}
\begin{aligned}
v_{0} & =\left(B_{12} B_{13}-B_{11} B_{23}\right) /\left(B_{11} B_{22}-B_{12}^{2}\right) \\
\lambda & =B_{33}-\left[B_{13}^{2}+v_{0}\left(B_{12} B_{13}-B_{11} B_{23}\right)\right] / B_{11} \\
\alpha & =\sqrt{\lambda / B_{11}} \\
\beta & =\sqrt{\lambda B_{11} /\left(B_{11} B_{22}-B_{12}^{2}\right)} \\
c & =-B_{12} \alpha^{2} \beta / \lambda \\
u_{0} & =c v_{0} / \alpha-B_{13} \alpha^{2} / \lambda
\end{aligned}
\end{equation*}
\section{用旋转矩阵近似\texorpdfstring{$3\times3$}{}矩阵}
\label{section:C}
本节是为了构造一个最佳的旋转矩阵$\mathbf{R}$来近似一个$3\times3$的矩阵$\mathbf{Q}$。
这里的最佳是基于$\mathbf{R}-\mathbf{Q}$差值的最小F范数。
解决方法参见我们给出的技术报告\cite{zhang1998flexible}。
\bibliography{Zhang_1999_translation.bib}
% \end{multicols}
\end{document}
2 条评论
优秀！
我家另一只狗子来了