A Flexible New Technique for Camera Calibration----张正友相

?????? 做相机校正时用到了张正友相机标定的方法。了解到Opencv和Matlab中的相机标定都是用的张正友相机标定方法，很经典的方法，但是我在网上竟然没有找到它的中文翻译文章。我回头仔细看了张正友相机标定的英文原文，大致是看明白了过程是什么。于是就想着翻译一下，逐字逐句的来真的翻译时，才发现自己的英文水平真的很有限

~~，这里本着不放弃的精神，坚持翻译完算法的讲解部分。有些地方真的不明白怎么翻译合适，希望各位同仁们赶快来纠正我错误的地方，我想能有个像样的翻译版本供大家查看，我将感激不尽。

????? 张正友相机标定--A Flexible New Technique for Camera Calibration。原文免费下载地址。

? ? 从原文的第3面开始：即从“2 Basic Equations“开始 到“4 Degenerate Configurations”那部分结束。

??? 现在开始我的翻译：

2基本方程

? ? 我们来测试相机内参矩阵的约束条件，内参矩阵是通过观察一个平面获得的。我们现在从本文中使用的符号开始。

2.1符号

? ? 一个二维点表示为：m=[u,v]^T。一个三维点表示为：M=[X,Y,Z]^T。我们使用

来表示增广向量，方法是在向量中加一个1作为最后一个元素：

和

。相机模型使用常用的针孔模型建模：即一个3D点M和它的投影图像之间的关系，如下：

??????

????? ?????????? （1）

????其中s是一个任意的缩放因子，（R,T）叫做外参，是连接世界坐标系系统与相机坐标系系统的旋转与平移矩阵，A叫做内参矩阵，如下所示：

? ? ? ?(u0,v0)是主点位置，

是图像在u轴和v轴的尺寸因子，

参数描述的是图像两个坐标系的扭曲因子。

我们使用缩写方式A^-T代表(A^-1)^T或者(A^T)^-1。

2.2模板平面和它对应图像之间的单应性矩阵

???

? ? 通常，我们假定模板平面在世界坐标系系统的Z=0处。我们将旋转矩阵R的第i^th列表示为r_i。由（1）式知道：

??? 我们仍用M表示模板上的一点，但是因为Z总是等于0。而

。因此，一个模板上的点M和它在图像上点m之间的关系用单应性矩阵表示：

????????????????????

???????????????? (2)

? ? 很明显，3x3的矩阵H等价于一个缩放因子。

2.3内参矩阵的约束条件

? ? 给定一个平面模板的图像，那么单应性矩阵可以被估计出来（见附录A）。我们将它表示为H=[h₁h₂ h₃]。从（2）式，我们知道：

? ? 其中

是一个任意的尺寸因子。基于r₁和r₂是正交的原理，我们得到：

???????

?????????????????????? (3)?
??

???????????? (4)

? ? 这就是内参矩阵的两个基本约束条件，已给定一个单应性矩阵。因为一个单应性矩阵有8个自由度和6个外参数（3个旋转的，3个平移的），我们只能找到2个内参约束条件。注意到A^-TA^-1实际上描述的是绝对二次曲线的图像。在下一节中，我们将给出几何解释。

2.4几何解释

? ? 我们现在将（3）和（4）联系到绝对二次曲线。

? ? 不难发现，在我们前面约定好的前提下，模板在相机坐标系系统中描述为如下方程：

? ? 其中w=0表示点在无穷远处，否则w=1。这个平面和另一个平面相交于一条直线的无穷远处，我们可以很容易的发现，

和

是那条直线上的两个特殊点。这条直线上的任何一个点都可以表示为这两个点的线性组合。例如：

? ??? 现在，我们用绝对二次曲线来计算上式的交点。根据定义，

点处，我们通常所说的圆形点（循环点），满足：

，例如：

??? 方程的解是

，其中

。也就是，方程的两个解是：

? ? 它们在图像平面的投影就可以知道了，加上一个缩放因子，为：

???

点是在二次曲线图像上的，描述为

。这可以得到：

? ? 要求实部和虚部都是0在（3）和（4）范围内。

3? 实现相机标定

? ? ? ?这一部分提供如何有效解决相机标定问题的细节。我们根据一种基于最大似然估计的非线性优化技术，开始分析求解。最后将镜头畸变考虑进去，给出分析和非线性解。

3.1封闭型解

令：

??

????????????????? (5)?

? ? ?注意B是对称的，可以定义为一个6D向量：

???????????????????????????? ?

?????????????? ????? （6）

? ? ?令H的第i^th列向量为：

。然后我们就得到：

???????????????????????????????????????????????????????????????????????????????????????????????????????

? ?? ? 因此，由给定的单应性矩阵得到的，这两个基础的约束条件(3)和(4)，可被重写为用b表示的2个齐次方程：

????????????????????????????????????????????????????????????? ??????????

??????????? （8）

? ?? ?如果有n各模板的观测图片，通过迭代n个（8）式的方程组，我们可得到：

???????????????????????

???????????? （9）

? ? ?其中V是一个2nx6的矩阵。如果n>=3，我们通常将有一个独特的解B，被定义为缩放因子。如果n=2，我们将利用扭曲因子

，例如：[0,1,0]b=0，被作为（9）式的一个附加方程。（如果n=1，我们只能解出两个相机的内参矩阵，例如：

，假定

已经知道（例如：在图像的中心）并且

，且这也确实是我们在[19]中为检测头部姿势所做的方法，头部检测是基于人眼和嘴是共面的事实）。（9）式的解就是我们熟知的

最小特征值对应的特征向量（等价的，V的右奇异向量对应最小奇异值）。

? ? ? ?一旦b被估算出来了，我们可以出计算相机内参矩阵A的所有元素，了解更多细节参考附录B。

? ? ? ?一旦A知道了，每个图像的外参系数很容易就可计算出来。从（2）式可知：

? ? ? 其中

。当然，因为数据中的噪声，所谓的计算矩阵

通常不能满足旋转矩阵的特性。附录C描述了从一个3x3的矩阵估计最优旋转矩阵的方法。

3.2? 最大似然估计

? ? ? ?上述解是通过最小化一个没有实际物理意义的代数距离获得的。我们能通过最大似然方法推导细化。

?????? 假设给定n个模板的图片且每个模板上有m各点。假定图像上点集被独立同分布的噪声损坏了。最大似然估计值可以通过最小化以下方程获得：

??????????????

????????（10）

? ? ? 根据方程（2），其中

是点M_j在图像i上的投影。一个旋转矩阵R是由一个3个参数的向量参数化的，用r表示，r是平行于旋转轴的，r的大小是和旋转角等价的。R和r是由Rodrigues公式联系起来的。最小化（10）是一个非线性最小化的问题，它可通过Levenberg-Marquardt算法解决，在MInpack[18]中有实现方法。它需要一个初始猜测的A，{R_i，t_i|i=1..n}，可以使用前面部分描述的技术获得。

3.3? 处理径向畸变

? ? ? ? 现在为止，我们还没有考虑相机的透镜的畸变。然而，一个桌面相机通常呈现明显的透镜畸变，特别是径向畸变。在这一部分中，我们仅考虑前两阶径向畸变。读者可以参考[20,2,4,26]获取更多的阐述。基于[2,23,25]文章中的报告，很可能失真函数完全是由径向分量主导的，尤其是由一阶系数主导。同时也发现，更细的建模不仅不会有帮助（与传感器数量级相比可忽略），还会引起数量不稳定[23,25]。

? ? ? ? 令（u，v）为理想（无失真）图像像素坐标，

是相应的实际观测的图像坐标。理想的点是模板上的点经过针孔模型的投影。相似的，

是理想（无失真）和实际（畸变）归一化的图像坐标系。我们有：

，

? ? ? ? 其中k₁和k₂是径向畸变系数。径向畸变的光学中心和主点是相同的。由

和

以及

我们得到：

?????????

????????????????????(11)?

???? ?

? ???????? ?????? ? ?(12)

? ? ? ? 迭代估计径向失真。径向失真一般认为较小，人们估计其他五个参数时一般认为是使用3.2节描述的技术，所以忽略失真也是认为合理的。估计k₁和k₂的一种策略是估计出其他参数，这样就给出我们理想的像素坐标（u，v）。然后，从（11）和（12）式，对每个图片的每个点我们有两个方程：

。

? ? ? ? 已知在n幅图片中有m个点，我们可以将所有的方程堆在一起以获得总共2mn个方程，或者以矩阵形式表示为Dk=d,其中k=[k₁,k₂]^T.线性最小平方解（最小二乘法）如下：

?????????????????????????????????????????? ??

??????? (13)

? ? ? ? 一旦k₁和k₂被估算出来，你就可以细化其他参数的估计值，通过求解（10），用(11)和(12)代替

。我们迭代这两个步骤直到收敛。

完成最大似然估计。通过实验，我们发现上面方法的收敛较慢。（10）的一个自然延伸接下来通过最小化下面功能来估计整套参数。

???????????????????????????? ????????

?????（14）

? ? ? ? 其中

是点M_j通过方程（2）在图像i中的投影，并接着经过（11）和（12）式的畸变。这是一个非线性的最小化问题，它可通过Levenberg-Marquardt算法解决，在MInpack[18]中有实现方法。一个旋转矩阵仍是被一个3维向量r参数化，见3.2节。一个A和{R_i，t_i| i=1..n}初始化值可通过3.1或3.2节中描述的方法获得。k₁和k₂的初始化猜测值可以通过上一段描述的方法获得，或者简单的设置为0.

3.4 总结

???????推荐的标定步骤如下：

1.??????打印一个棋盘格图案，将它贴到一个平坦的表面；

2.??????在不同方向下拍摄一些模板平面的图像。可以通过移动模板或者移动相机改变方向；

3.??????检测图像的特征点；

4.??????使用3.1节描述的封闭式方程估算五个内参矩阵的参数以及外参。

5.??????通过求解线性最小二乘式（13），估算径向畸变系数。

6.??????通过最小化（14）式，细化所有参数。

4?? 不好的配置

? ? ? ?我们在这一部分的配置中学习附加的图片不会提供更多相机内参的约束条件的情况。因为（3）和（4）是从旋转矩阵的性质中提取的，如果R₂不是与R₁相互独立的，那么图片2不提供额外的约束。尤其，如果一个平面经过一个纯粹的平移，然后R₂=R₁这样图2对相机标定时是无帮助的。接下来，我们考虑一个更复杂的配置。

命题1.如果模板平面在第二个位置与第一个位置是平行的，这时第二个图像的单应性矩阵不提供额外的约束条件。

? ? ? ? 证明。在我们的约定下，R₂和R₁是通过绕z轴旋转关联起来的。也就是，

，

? ? ? ? ?其中

是相对旋转角度。我们将分别使用角标^（1）和^（2）来表示图1和图2相关的向量。很清楚我们有：

?????那么，从第2幅图得到的第一个约束条件（3）式变成，

? ? ? ? 上式是H₁提供的两个约束条件的线性组合。相似的，我们可以表明从图2获得的第二个约束条件，也是H₁提供的两个约束条件的线性组合。因此，我们不能从H₂获得任何约束条件。

? ? ? ? 这个结果明显是因为平行的模板交平面于无穷远处在相同的圆环点处，并且这样根据2.4节它们提供了相同的约束条件。

? ? ? ? 实践中，非常容易就可以避免不好的配置：我们只需要从一次拍照到另一次拍照之间改变模板的方向。

? ? ? ? 尽管如果模板只是经过纯平移时提出的方法将不能有效，但如果知道平移量相机标定仍然可能计算出来。请参考附录D。