3.3 基于单应矩阵的摄像机畸变参数估计方法
本节仅讨论基于平面标定模板的畸变参数估计方法。前提条件为：①已知平面标定模板上位置精确的标定点（如棋盘格点）；②已知平面标定模板标定点投影的图像点；③已知摄像机的畸变模型；④已知主点坐标。本节研究如何由上述已知条件确定摄像机畸变参数的过程和方法。从原理上讲，本节的方法也适用于标定物为非平面标定模板的情况，此时需要利用图像主点附近的图像点与对应空间点之间的关系来估计3×4的投影矩阵，其余步骤与平面单应矩阵估计方法相似。
本节基于单应矩阵的畸变参数估计方法的基本思路为：由摄像机畸变公式[见式（3-2）]可知，对任何摄像机畸变模型而言，图像主点附近的图像点畸变效应都不大，因此，利用图像主点附近的图像点与空间模板的对应点估计的平面单应矩阵，应该与针孔模型下摄像机对应的理想的单应矩阵很接近。因此，由主点附近的点估计的单应矩阵，可以比较可靠地估计图像中每个图像点由畸变参数导致的图像偏移量。利用已知的图像偏移量，在给定的摄像机畸变模型下，可以进一步估计摄像机畸变参数的具体数值。
基于上面的基本思路，本节提出的基于单应矩阵的摄像机畸变参数估计的基本算法如下。
算法 基于单应矩阵的摄像机畸变参数估计方法(A Method of Estimating the Camera Distortion Parameters from A Single Image，MECDPSI)
输入：N个空间平面模板点坐标X及其对应的图像平面上的实际图像点坐标Ud，即{Xi↔Udi，i=1,2,3,…,N}，摄像机的畸变模型类型；摄像机的主点坐标为(u0,v0)。输出：畸变参数估计值和理想无畸变下的单应矩阵。
（1）在图像主点(u0,v0)附近选取M≥4个图像点，利用这M个图像点与空间点的对应关系估计初始单应矩阵H0。
（2）将所有空间平面点{Xi，i=1,2,3,…,N}利用第一步计算的H0生成N个虚拟图像点{Uri，i=1,2,3,…,N}。将这N个虚拟图像点{Uri，i=1,2,3,…,N}看作空间平面点在无畸变下的理想图像点。
（3）计算由畸变形成的实际图像点与虚拟图像点之间的图像偏移量{∆Ui=Udi−Uri，i=1,2,3,…,N}。
（4）利用上步估计的图像偏移量{∆Ui=Udi−Uri，i=1,2,3,…,N}估计畸变模型的初始畸变参数P。
（5）将上述初始估计的单应矩阵H0和初始畸变参数P作为初值，通过非线性迭代优化，得到最终的畸变参数和单应矩阵。
下面对关键步骤（1）、（4）与（5）进行详细介绍。
3.3.1 估计初始单应矩阵H0
由3.2.2节的摄像机畸变模型可知，图像点离主点越近，产生的畸变越小，即主点附近的图像点，在针孔模型下的理想图像点与畸变模型下的实际图像点非常接近。因此，依据图像中心点（主点一般在图像中心）附近的点计算的单应矩阵应该接近在无畸变模型下的单应矩阵。MECDPSI算法首先依据选择中心点附近图像点的比例因子t，从所有图像点中选出n（n=t×N，n≥4）个中心点附近的图像点，计算初始的单应矩阵H0。其具体计算过程如下。
在针孔模型下，由图像点与对应空间点之间的成像关系可知，空间平面上的n个点X(xi,yi,1)T（i=1,2,…,n）与图像平面上的n个对应点Udi(ui,vi,1)T（i=1,2,…,n）之间存在以下的单应变换（Homography）。
由于有未知因子λi，因此需要消掉λi，由式（3-3）可以得到下面两个关于hij（i=1,2,3；j=1,2,3）的约束方程，即
展开式（3-4）和式（3-5），可以得到下列方程。
将n组对应得到的约束按上面的形式排在一起，可以得到一个2n×9的齐次矩阵方程，即
式（3-7）中的齐次方程有9个未知数，需要有2n（n≥4）个约束才能求解。因此，由空间平面上的4个非共线点就可以估计该单应矩阵，记为H0。
3.3.2 估计初始摄像机畸变参数P
依据3.3.1节计算出的初始单应矩阵H0与空间平面上的N个点X(xi,yi,1)T（i=1,2,…,N）的信息，可以得出空间平面点任意点X对应的图像中的虚拟点信息{Uri，i=1,2,3,…,N}。记Uri=(uri,vri,1)T（i=1,2,…,N）为针孔模型下的虚拟图像坐标，{Udi=(udi,vdi,1)T，i=1,2,3,…,N}为畸变图像坐标。假设图像存在二阶径向畸变，根据3.2.2节中的畸变模型，每个Ud及其对应点Ur均存在以下等式。
由N个点可得到下列2N个等式。
上面的方程有两个未知数、2n（n≥1）个约束，因此，至少需要一组空间与图像的对应点，便可计算出一组(k1,k2)T，即估计出的畸变参数P。
由已知的畸变图像点{Udi，i=1,2,3,…,N}与畸变参数P，可求解下一步优化初始的单应矩阵H0与畸变参数P。
3.3.3 优化参数
由初始计算的摄像机畸变参数P与已知的畸变图像点{Udi，i=1,2,3,…,N}估计新的单应矩阵与畸变参数，具体算法流程如下。
输入：初始H0、初始畸变参数P、空间平面π上N个点X(xi,yi,1)（i=1,2,…,N）和对应的图像平面I上N个点Udi(ui,vi,1)T（i=1,2,…,n）。
输出：优化后的P与优化后的H0。
优化的步骤如下。
（1）由输入的初始H0与X(xi,yi,1)（i=1,2,…,N）获得图像点{Uri，i=1,2,3,…,N}。
（2）以H0和P为初值，通过LM算法最小化公式
（3）输出新的P与新的H0。
3.3.4 算法说明
上述摄像机畸变参数估计方法的主要特点是：不需要像文献中主流方法那样，首先估计摄像机的线性参数，然后估计摄像机的畸变参数，而是直接通过单应矩阵估计摄像机的畸变参数；并且MECDPSI算法仅需要单幅图像就可以完成摄像机畸变参数估计过程，极大地提高了摄像机标定方法的灵活性。当然，MECDPSI算法要求已知摄像机主点信息，这是基于单幅图像标定需要付出的代价。目前大多数摄像机的主点离图像中心不远，所以“已知主点”信息事实上不是一个很苛刻的条件。
针对MECDPSI算法，有以下几点需要特别说明。
（1）在算法的步骤（1），即在估计初始单应矩阵时，可以利用DLT算法。为了提高估计的可靠性，应尽量利用Hartley[73]提出的数据归一化策略。
（2）在步骤（5）的非线性迭代优化中，目标函数为实际图像点Udi与所估计模型下对应的图像点Uri之间的几何距离，即
（3）最终优化得到的单应矩阵，可以进一步提供关于摄像机内参数、世界坐标系和摄像机坐标系之间的旋转矩阵与平移向量信息，但本节不再进一步讨论这个问题。
（4）假定摄像机仅含两个径向畸变参数和两个切向畸变参数，将对应的理想点与畸变点代入，可得到(k1,k2,p1,p2)T的具体估计方程为
上面的方程中有4个未知数、2×N（N≥2）个约束。当N>2时，可在最小二乘意义下得到一组(k1,k2,p1,p2)的估计值。