CN103886595B - 一种基于广义统一模型的折反射相机自标定方法 - Google Patents

Abstract

本发明公开了一种基于广义统一模型的折反射相机自标定方法。该方法包括以下步骤：在同一场景不同位置使用折反射相机获取两幅左右成像清晰的折反射全景图；提取左右两图的SIFT特征点，并进行特征点匹配；初始化广义统一模型的内参；使用极几何计算出本质矩阵E，在计算E的过程中使用RANSAC提高精度，同时去除外点；使用非线性优化Levenberg Marquardt方法，最小化特征点的重投影误差；不断的优化内参，直到重投影误差小于所设定的阈值，优化出鲁棒和精确的结果。本发明仅需两幅不同位置的图即可标定广义统一模型的内参和两个位置之间的关系，标定过程完全自动，具有方便快捷，操作简单，精度较高的特点。

Description

一种基于广义统一模型的折反射相机自标定方法

技术领域

本发明涉及图像数据、处理分析的方法，尤其是涉及一种基于广义统一模型的折反射相机自标定方法。

背景技术

传统的透视相机视场较小，在某些场景中使用较为受限制。随着图像传感器技术和机器视觉技术的进步，大视场图的应用越来越广泛。折反射相机由传统的透视相机和一个反射镜构成，能够采集水平方向360度视场的全景图，并且系统构成简单，功耗和成本都较小。通过精心设计不同类型的反射镜面可以构成不同特性的折反射相机系统，在此基础上可以精确得到空间的度量信息。由于其大视场和图像信息结构化等特点，折反射相机系统被广泛的使用于机器人导航、监控、视频会议、场景重建等等应用中

折反射全向相机根据入射光线是否都通过一个虚拟点可以分为：单一视点折反射相机和非单视点折反射相机。单一视点折反射相机成像原理简单，标定和使用方便，但是很难得到，实际当中不存在完美的单视点折反射相机。当镜面和透视相机之间的位姿关系不满足特定条件的时候，如镜面发生旋转或者平移，使得镜面对称轴和相机光轴不再对准，折反射相机发生失准（misalignment），单一视点折反射相机转换为非单视点折反射相机。现实当中的折反射相机都或多或少的存在失准。为了描述非单视点折反射相机，不同模型和标定方法被提出。文献1（C.Mei and P.Rives,"Single View Point OmnidirectionalCamera Calibration from Planar Grids,"in Robotics and Automation,2007IEEEInternational Conference on,2007）在统一模型的基础上加入了切向畸变来弥补失准，但是切向畸变只能补偿较小的失准，失准较大时切向畸变无法补偿，使得在失准较大时，标定结果并不理想。文献2（Z.Xiang,B.Sun,and X.Dai,"The Camera Itself as aCalibration Pattern:A Novel Self-Calibration Method for Non-CentralCatadioptric Cameras,"Sensors12,7299-7317(2012).）没有将非单视点折反射相机看成一个整体，而是将镜面和透视相机看成两个部分，标定折反射相机就是标定镜面和透视相机直接的位姿关系，该方法利用镜面边缘和透视相机镜头成像来标定位姿关系，但是该方法的成像模型复杂，透视相机内参需要离线标定。文献3（Xiang,Z.,X.Dai,and X.Gong,Noncentral catadioptric camera calibration using a generalized unifiedmodel.Optics letters,2013.38(9):p.1367-1369）提出了一种针对非单视点折反射相机的广义统一模型，该模型成像结构简单，并且能够非常好的弥补失准，即使在失准较大的情况下，也能很好的对非单视点折反射相机进行建模。

使用标定板标定广义统一模型虽然精度较高，但是需要在场景中不同位置放置标定板，并且需要手工提取标定板在图像中的网格点，在折反射相机配置经常发生配置的场景下，使用较为不便。

发明内容

针对背景技术中标定方法的不足，本发明的目的在于提出一种基于广义统一模型的折反射相机自标定方法。使用折反射相机广义统一模型对非单视点折反射相机进行建模，使用折反射相机在空间采集两幅图，提取两幅图上匹配的特征点，利用配对的特征点自动标定广义统一模型的内参和两幅图之间的外参的方法。该方法不需要在场景当中加入任何标记物，只需要采集场景中不同位置的两幅图，即可对折反射相机进行建模，并标定出相应的广义统一模型的参数。

本发明采用的技术方案的步骤如下：

步骤1)在同一场景不同位置使用折反射相机获取两幅左右成像清晰的折反射全景图,要求两图中有场景重复的部分，不需要对环境做特殊的设置，也不需要任何标定模板；

步骤2)提取左右两图的SIFT特征点，并进行特征点匹配，得到图像上n对匹配的特征点对和i=1..n，n对特征点对中存在误配点；

步骤3)初始化广义统一模型的内参；

步骤4)使用极几何计算出本质矩阵E，在计算E的过程中使用RANSAC提高精度，同时去除外点；

步骤5)使用非线性优化Levenberg Marquardt方法，最小化特征点的重投影误差，从而优化出内参；

步骤6)重新回到步骤4)，循环步骤4)和步骤5)，不断的优化内参，同时得到更多的内点，直到重投影误差小于所设定的阈值，优化出鲁棒和精确的结果。

步骤3)中所述要初始化广义统一模型的内参，首先使用广义统一模型对折反射相机进行建模；

折反射相机广义统一模型由虚拟的单位球，虚拟视点，归一化平面和成像平面П_p组成；以单位球中心C_m为原点建立坐标系F_m，以虚拟光心C_p为原点建立坐标系F_p，两坐标系相互平行，归一化平面为F_p中的z＝1平面，归一化平面上的点通过径向失真和透视变换投影到成像平面П_p；

折反射相机广义统一模型，有两个坐标系，以虚拟光心C_p为原点建立坐标系F_p，和以单位球中心C_m为原点建立坐标系F_m，原点C_p在F_m坐标系中的坐标为ξ＝(ξ1,ξ2,ξ3)^T；

空间点X_m在世界坐标系F_w中的表示为(X_m)_W＝(x_w,y_w,z_w)，通过以下前向投影过程，即成像过程如下：

(1)空间点从世界坐标系F_w转换到F_m坐标系，其中R_w和T_w为世界坐标系和F_m坐标系的转换关系，此时空间点在F_m坐标系中的坐标为 ${\left(X_m\right)}_{F_m}={(x_m,y_m,z_m)}^T;$

(2)将投影到单位球面上坐标为(x_s,y_s,z_s)^T，

(3)转换到F_p坐标系下其中ξ＝(ξ1,ξ2,ξ3)^T， ${\left(X_s\right)}_{F_p}={\left(X_s\right)}_{F_m}+\mathrm{\ξ}={(x_s+\mathrm{\ξ}1,y_s+\mathrm{\ξ}2,z_s+\mathrm{\ξ}3)}^T;$

(4)将投影到归一化平面的m_u坐标为(x,y,1)^T， $m_u=\frac{1}{\mathrm{\λ}}{\left(X_s\right)}_{F_p}={(\frac{x_s+\mathrm{\ξ}1}{z_s+\mathrm{\ξ}3},\frac{y_s+\mathrm{\ξ}2}{z_s+\mathrm{\ξ}3},1)}^T=h\left({\left(X_s\right)}_{F_m}\right),$ 其中λ＝z_s+ξ3，h()为折反射投影函数；

(5)在归一化平面上加入径向失真以弥补真实系统中透视相机镜头的失真，径向失真函数为径向畸变函数L(ρ)＝1+k₁ρ²+k₂ρ⁴+k₃ρ⁶，其中，k₁,k₂,k₃为失真系数，失真后的点m'_u＝(x',y',1)^T，其中x'＝x*L(ρ),y'＝y*L(ρ)；

(6)通过透视变换将归一化平面上的点m'_u投影到成像平面П_p上点p坐标为(u，v，1)^T， $p={(u,v,1)}^T=K{m}_u^{\′}=\left[\begin{array}{ccc}{f}_1\mathrm{\η}& f_1\mathrm{\η\α}& u_0\\ 0& f_2\mathrm{\η}& v_0\\ 0& 0& 1\end{array}\right]m_u^{\′}=k\left(m_u^{\′}\right);$ 其中k()为透视投影函数，f₁，f₂是相机的焦距，(u₀,v₀)是相机的主点，α是偏斜，η是镜面参数；此处的K为归一化平面到成像平面的透视变换矩阵，不再将折反射系统看成是一个独立的镜面和一个独立的相机，而是将之视为一个整体，内参矩阵K包含了相机和镜面的参数，f₁，f₂和η无法单独得到，只能得到γ₁＝f₁η，γ₂＝f₂η；

从图像上的点p反向推导出相应的单位球面上的点的方向的过程称为后向投影过程：

(1)从图像点p＝(u,v,1)^T映射到归一化平面上的点m'_u，m′_u＝K^-1p＝(x,y,1)^T；

(2)对m′_u去除径向畸变，径向畸变函数L(ρ′)＝1+k₁ρ^′2+k₂ρ^′4+k₃ρ^′6，k₁,k₂,k₃为失真系数，去失真后的点m_u＝(x,y,1)^T，其中x=x'/L(ρ′),y＝y′/L(ρ′)

(3)从归一化平面上的点m_u映射到单位球面上的点 其中h()是投影函数，h^-1()是投影函数的逆函数；

由前向投影过程中的步骤(4)(3)(2)解出从图像点到空间射线的后向投影方程 ${\left(X_s\right)}_{F_m}=\frac{{\mathrm{\ξ}}^T{K}^{-1}p+\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}p\right)}^2-\left(p^T{K}^{-T}{K}^{-1}p\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p^T{K}^{-T}{K}^{-1}p}{K}^{-1}p-\mathrm{\ξ},$ 此处未加入去除径向畸变，因为径向畸变非线性，加入径向畸变以后后向投影公式过于复杂，实际使用当中在初始化阶段忽略径向畸变，在优化过程中分步进行后向投影的计算而不是利用后向投影方程一次性得出。

通过移动折反射相机，在不同的位置采集空间点的成像，空间点X_m在世界坐标系F_w中的表示为(X_m)_W＝(x_w,y_w,z_w)，在不同位置折反射相机中成像，成像过程如前面前向投影所述，在两幅图上形成的成像点分别为p₁＝(u₁,v₁,1)^T和p₂＝(u₂,v₂,1)^T，根据后向投影方程计算出p₁和p₂对应的和射线和都经过空间点X_m，得到：其中本质矩阵E＝T₁₂×R₁₂，其中R₁₂和T₁₂为第一个相机和第二个相机之间的坐标转换关系；接着得到广义统一模型的极几何公式：

$\begin{array}{c}(\frac{{\mathrm{\ξ}}^T{K}^{-1}{p}_2+\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}{p}_2\right)}^2-\left(p_2^T{K}^{-T}{K}^{-1}{p}_2\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p_2^T{K}^{-T}{K}^{-1}{p}_2}{K}^{-1}{p}_2-\mathrm{\ξ})E\\ (\frac{{\mathrm{\ξ}}^T{K}^{-1}{p}_1+\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}{p}_1\right)}^2-\left(p_1^T{K}^{-T}{K}^{-1}{p}_1\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p_1^T{K}^{-T}{K}^{-1}{p}_1}{K}^{-1}{p}_1-\mathrm{\ξ})=0\end{array};$

广义统一模型含有需要标定的参数为，k₁,k₂,k₃，ξ＝(ξ1,ξ2,ξ3)^T，γ₁，γ₂，[u₀,v₀]和α；

其次，为了获取模型的初值，先假设折反射全向相机为单一视点模型，设k₁≈k₂≈k₃≈α≈ξ1≈ξ2≈0，γ₁≈γ₂＝γ，所需要估计的内参为ξ3、γ、u₀、v₀，实验表明ξ3对系统参数的初值计算影响不大，设ξ3＝1，主点[u₀,v₀]的初始值估计，先提取图像中镜面边缘，再取图像中镜面的中心为主点的初值；

将 $K=\left[\begin{array}{ccc}\mathrm{\γ}& 0& u_0\\ 0& \mathrm{\γ}& v_0\\ 0& 0& 1\end{array}\right]$ 和ξ＝[0,0,1]^T带入极几何公式得到 $(\frac{2}{p_2^T{K}^{-T}{K}^{-1}{p}_2}{K}^{-1}{p}_2-\mathrm{\ξ})E(\frac{2}{p_1^T{K}^{-T}{K}^{-1}{p}_1}{K}^{-1}{p}_1-\mathrm{\ξ})=0,$ 将 $K^{-1}{p}_j={[\frac{u_j-u_0}{\mathrm{\γ}},\frac{v_j-v_0}{\mathrm{\γ}},1]}^T={[\frac{u_{\mathrm{cj}}}{\mathrm{\γ}},\frac{v_{\mathrm{cj}}}{\mathrm{\γ}},1]}^T(j=\mathrm{1,2})$ 带入得到 $[\mathrm{\γ}{u}_{c2},\mathrm{\γ}{v}_{c2},\frac{{\mathrm{\γ}}^2}{2}-\frac{u_{c2}^2+v_{c2}^2}{2}]E{[\mathrm{\γ}{u}_{c1},\mathrm{\γ}{v}_{c1},\frac{{\mathrm{\γ}}^2}{2}-\frac{u_{c1}^2+v_{c2}^2}{2}]}^T=0,$ 设 $E=\left[\begin{array}{ccc}{e}_1& e_2& e_3\\ e_4& e_5& e_6\\ e_7& e_8& e_9\end{array}\right],$ 重新整理后得到多项式特征值问题[D₀+γD₁+γ²D₂+γ³D₃+γ⁴D₄]e＝0，其中

$D_0=[\mathrm{0,0,0,0,0,0,0,0},u_{c2}^2{u}_{c1}^2+u_{c2}^2{v}_{c1}^2+u_{c1}^2{v}_{c2}^2+v_{c2}^2{v}_{c1}^2]$

$D_1=[\mathrm{0,0},-2u_{c2}{v}_{c1}^2-2u_{c2}{u}_{c1}^2,\mathrm{0,0},-{2u}_{c1}^2{v}_{c2}-{2v}_{c2}{v}_{c1}^2,-{2u}_{c2}^2{u}_{c1}-{2u}_{c1}{v}_{c2}^2,-{2u}_{c2}^2{v}_{c1}-{2v}_{c2}^2{v}_{c1},0]$

$D_2=[{4u}_{c2}{u}_{c1},{4u}_{c2}{v}_{c1},0,{4u}_{c1}{v}_{c2},{4v}_{c2}{v}_{c1},\mathrm{0,0,0},-u_{c2}^2-u_{c1}^2-v_{c2}^2-v_{c1}^2]$

D₃=[0,0,2u_c2,0,0,2v_c2,2u_c1,2v_c1,0]

D₄=[0,0,0,0,0,0,0,0,1]

e＝[e₁,e₂,e₃,e₄,e₅,e₆,e₇,e₈,e₉]

该多项式特征值问题通过9个特征值点对解出e和γ，使用RANSAC来提高精度,求得γ的初值。

所述步骤（4）具体为：

在内参初值K和ξ已知的情况下,每一对特征点带入极几何公式，得到一个关于E的线性方程其中从后向投影方程得到， ${(x_j^i,y_j^i,z_j^i)}^T=(\frac{2}{p_j^{\mathrm{iT}}{K}^{-T}{K}^{-1}{p}_j^i}{K}^{-1}{p}_j^i-\mathrm{\ξ}),$ j＝1,2，i＝1..n；若设 $E=\left[\begin{array}{ccc}{e}_1& e_2& e_3\\ e_4& e_5& e_6\\ e_7& e_8& e_9\end{array}\right],$ 线性方程转换为aⁱe＝0，其中 $a^i=[x_2^i{x}_1^i,x_2^i{y}_1^i,x_2^i{z}_1^i,y_2^i{x}_1^i,y_2^i{y}_1^i,y_2^i{z}_1^i,z_2^i{x}_1^i,z_2^i{y}_1^i,z_2^i{z}_1^i],$ e＝[e₁,e₂,e₃,e₄,e₅,e₆,e₇,e₈,e₉]；当特征点对数目大于等于8时，求出e的最小二乘解，包含一个尺度因子；在计算E的过程中使用RANSAC提高精度，排除误差过大的外点；

奇异值分解E，E＝U diag(σ₁,σ₂,σ₃)VT，其中σ₁≥σ₂≥σ₃，U、V是单位正交矩阵；得到R₁₂和T₁₂， $\left\{\begin{array}{c}{R}_{12}={\mathrm{UR}}_z^T(\±\frac{\mathrm{\π}}{2})V^T\\ {\hat{T}}_{12}={\mathrm{UR}}_z(\±\frac{\mathrm{\π}}{2})\mathrm{\Σ}{U}^T\end{array}\right.,$ 其中Σ＝diag(1,1,0)，T₁₂包含一尺度因子；选择具有物理意义的那一组R₁₂和T₁₂，即两幅图对应相机朝向相近的那一组。

所述步骤（5）具体为：

使用非线性优化方法Levenberg Marquardt方法,通过最小化特征点的重投影误差来优化出内参；

要计算特征点的重投影误差，要对特征点进行三维重构；在内参K和ξ已知的情况下，并且在得到外参R₁₂和T₁₂的情况下，对空间点进行重构；对任意特征点对p₁＝(u₁,v₁,1)^T和p₂＝(u₂,v₂,1)^T，根据反向投影方程，得到对应的空间入射光线的方向，即和通过外参R₁₂和T₁₂，可将两光线变换到同一坐标系下和理想情况下两光线相交，实际当中由于存在误差，两光线为异面直线，取异面直线中垂线的中点X'_m为重构点；通过前向投影得到X'_m在左右两图中的重投影p'₁＝(u'₁,v'₁,1)^T和p'₂＝(u'₂,v'₂,1)^T；这对特征点的重投影误差就为

$e_{\mathrm{reprojection}}=\frac{\left|(p_1^{\′}-p_1)\right|+\left|(p_2^{\′}-p_2)\right|}{2};$

使用非线性优化Levenberg Marquardt方法，最小化所有特征点的重投影误差的和，从而优化出内参。

所述的步骤（6）具体为：

步骤(4)中所用的内参初值从步骤(3)中得到，初始时初值误差较大，导致步骤（4）去除外点过多，进而影响了步骤（5）中的优化结果的精度；在经过步骤（5）的优化后，内参的值比初始时要精确，使用该内参去重新进行步骤（4）和步骤（5）；如此重复，不断优化得到更精确的内参，直到重投影误差小于说要求的阈值。

本发明具有的有益效果是：

本发明提出一种基于广义统一模型的折反射相机自标定方法。使用广义统一模型对非单视点折反射相机进行建模。使用折反射相机在空间不同位置采集的两幅图，完全自动的标定广义统一模型的内参和两幅图之间的位置关系。该方法不需要在场景中做特殊的标记，完全自动，不需要任何人工干预，使用方便且精度较高。

附图说明

图1是基于广义统一模型的折反射相机自标定方法步骤图。

图2是折反射相机的广义统一模型。

图3是折反射相机广义统一模型的极几何图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步的描述。

如图1所示，是基于广义统一模型的折反射相机自标定方法的步骤图。该自标定方法的步骤为：

步骤1)在同一场景不同位置使用折反射相机获取两幅左右成像清晰的折反射全景图,要求两图中有场景重复的部分，不需要对环境做特殊的设置，也不需要任何标定模板；

步骤2)提取左右两图的SIFT特征点，并进行特征点匹配，得到图像上n对匹配的特征点对和i=1..n，n对特征点对中存在误配点；文献4(Lowe,D.G.,Distinctive imagefeatures from scale-invariant keypoints.International journal of computervision,2004.60(2):p.91-110.)详细介绍了SIFT特征点的提取与匹配，得到图像上n对匹配的特征点对和i=1..n。n对特征点对中存在误配点。

步骤3)初始化广义统一模型的内参；

步骤4)使用极几何计算出本质矩阵E，在计算E的过程中使用RANSAC提高精度，同时去除外点；

步骤5)使用非线性优化Levenberg Marquardt方法，最小化特征点的重投影误差，从而优化出内参；

步骤6)重新回到步骤4)，循环步骤4)和步骤5)，不断的优化内参，同时得到更多的内点，直到重投影误差小于所设定的阈值，优化出鲁棒和精确的结果。

要使用基于广义统一模型的自标定方法来标定折反射相机，首先要使用广义统一模型来对折反射进行建模，图2所示为折反射相机的广义统一模型。同时还要提出广义统一模型的极几何，图3显示了折反射相机广义统一模型的极几何图。

广义统一模型的前向投影是空间点成像的过程，图2中有两个坐标系F_p（以虚拟光心C_p为原点）和F_m（以单位球中心C_m为原点），原点C_p在F_m坐标系中的坐标为ξ＝(ξ1,ξ2,ξ3)^T，在单位球模型中ξ1＝0，ξ2＝0，ξ3就是单位球模型的ξ。

空间点X_m(x_w,y_w,z_w)（世界坐标系F_w中的表示）的成像过程：

(1)从世界坐标转换到F_m坐标系，其中R_W和T_w为世界坐标系和F_m坐标系的转换关系，此时 ${\left(X_m\right)}_{F_m}={(x_m,y_m,z_m)}^T$

${\left(X_m\right)}_{F_m}=R_w{\left(X_m\right)}_w+T_w---\left(1\right)$

(2)将投影到单位球面上

${\left(X_s\right)}_{F_m}=\frac{{\left(X_m\right)}_{F_m}}{\left|{\left(X_m\right)}_{F_m}\right|}={(x_s,y_s,z_s)}^T---\left(2\right)$

(3)转换到F_p坐标系下其中ξ＝(ξ1,ξ2,ξ3)^T:

${\left(X_s\right)}_{F_p}={\left(X_s\right)}_{F_m}+\mathrm{\ξ}={(x_s+\mathrm{\ξ}1,y_s+\mathrm{\ξ}2,z_s+\mathrm{\ξ}3)}^T---\left(3\right)$

(4)将投影到归一化平面其中λ＝z_s+ξ3：

$m_u=\frac{1}{\mathrm{\λ}}{\left(X_s\right)}_{F_p}={(\frac{x_s+\mathrm{\ξ}1}{z_s+\mathrm{\ξ}3},\frac{y_s+\mathrm{\ξ}2}{z_s+\mathrm{\ξ}3},1)}^T=h\left({\left(X_s\right)}_{F_m}\right)---\left(4\right)$

(5)可以在m_u上加入径向畸变

(6)通过一个普通的透视变换将归一化平面上的点m_u投影到成像平面上点p＝(u,v,1)^T

$p=K{m}_u=\left[\begin{array}{ccc}{f}_1\mathrm{\η}& f_1\mathrm{\η\α}& u_0\\ 0& f_2\mathrm{\η}& v_0\\ 0& 0& 1\end{array}\right]m_u=k\left(m_u\right)---\left(5\right)$

其中f₁，f₂是相机的焦距，(u₀,v₀)是相机的主点，α是偏斜（skew），η是镜面参数。此处的K为归一化平面到成像平面的透视变换矩阵，不再将折反射系统看成是一个独立的镜面和一个独立的相机，而是将之视为一个整体，内参矩阵K包含了相机和镜面的参数，f₁，f₂和η无法单独得到，只能得到γ₁＝f₁η，γ₂＝f₂η。

后向投影是图像点到空间射线的投影过程，从图像点p＝(u,v,1)^T到空间射线的反向映射过程：

(1)从图像点p＝(u,v,1)^T映射到归一化平面上的点m_u：

m_u＝K^-1p＝(x,y,1)^T （6）

(2)如果正向投影加入径向畸变,需要对m_u去除径向畸变。

(3)从归一化平面上的点m_u映射到单位球面上的点：

${\left(X_s\right)}_{F_m}={(x_s,y_s,z_s)}^T=h^{-1}\left(m_u\right)---\left(7\right)$

由公式（4）可以得到其中λ＝z_s+ξ3，然而并未知道z_s的值由公式（3）可以得到：

${\left(X_s\right)}_{F_m}={\left(X_s\right)}_{F_p}-\mathrm{\ξ}=\mathrm{\λ}{K}^{-1}p-\mathrm{\ξ}---\left(8\right)$

再由公式（2）得到：

$\left|{\left(X_s\right)}_{F_m}\right|=1\⇒|\mathrm{\λ}{K}^{-1}p-\mathrm{\ξ}|=1$

展开得到：

(λK^-1p-ξ)^T(λK^-1p-ξ)＝1

λ²p^TK^-TK^-1p-2λξ^TK^-1p+ξ^Tξ＝1

可以解出:

$\mathrm{\λ}=\frac{{\mathrm{\ξ}}^T{K}^{-1}p\±\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}p\right)}^2-\left(p^T{K}^{-T}{K}^{-1}p\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p^T{K}^{-T}{K}^{-1}p}$

根据λ的定义及模型，明显λ取较大值

$\mathrm{\λ}=\frac{{\mathrm{\ξ}}^T{K}^{-1}p+\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}p\right)}^2-\left(p^T{K}^{-T}{K}^{-1}p\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p^T{K}^{-T}{K}^{-1}p}---\left(9\right)$

带入到公式（8），可以得到反向投影方程：

${\left(X_s\right)}_{F_m}=\frac{{\mathrm{\ξ}}^T{K}^{-1}p+\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}p\right)}^2-\left(p^T{K}^{-T}{K}^{-1}p\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p^T{K}^{-T}{K}^{-1}p}{K}^{-1}p-\mathrm{\ξ}---\left(10\right)$

传统透视相机极几何是空间同一点在两个透视相机中或者同一透视相机不同位置的成像的关系。通过移动折反射相机系统，也可以得到同一场景的不同位置的图，找出这两图之间的关系，也就是推算出折反射相机的极几何。

图3显示了空间点X_m(x_w,y_w,z_w)（世界坐标系F_w中的表示）在不同位置折反射相机（折反射相机用广义统一模型进行建模）中成像的过程，在两幅图上形成的成像点分别为p₁＝(u₁,v₁,1)^T和p₂＝(u₂,v₂,1)^T，在内参已知的情况下，根据反向投影方程（10）可以计算出对应的和射线和都经过空间点X_m，那么明显有：

${\mathrm{\λ}}_1{\left(X_{s1}\right)}_{F_{m1}}=R_{w1}{X}_w+T_{w1}---\left(11\right)$

${\mathrm{\λ}}_2{\left(X_{s2}\right)}_{F_{m2}}=R_{w2}{X}_w+T_{w2}---\left(12\right)$

其中，R_w1和T_w1为世界坐标系F_w和第一个相机F_m1之间的坐标转换关系，R_w2和T_w2为世界坐标系F_w和第二个相机F_m2之间的坐标转换关系，公式（11）和公式（12）联立消去X_m得到：

$R_{w1}^{-1}({\mathrm{\λ}}_1{\left(X_{s1}\right)}_{F_{m1}}-T_{w1})=R_{w2}^{-1}({\mathrm{\λ}}_2{\left(X_{s2}\right)}_{F_{m2}}-T_{w2})$

继续推导可以得到：

${\mathrm{\λ}}_1{R}_{12}{\left(X_{s1}\right)}_{F_{m1}}+T_{12}={\mathrm{\λ}}_2{\left(X_{s2}\right)}_{F_{m2}}---\left(13\right)$

其中R₁₂和T₁₂为第一个相机和第二个相机之间的坐标转换关系。公式（13）左右同乘

${\mathrm{\λ}}_1\widehat{T_{12}}{R}_{12}{\left(X_{s1}\right)}_{F_{m1}}={\mathrm{\λ}}_2\widehat{T_{12}}{\left(X_{s2}\right)}_{F_{m2}}$

再同乘

${\left(X_{s2}\right)}_{F_{m2}}^T{\hat{T}}_{12}{R}_{12}{\left(X_{s1}\right)}_{F_{m1}}=0---\left(14\right)$

定义 $E={\hat{T}}_{12}{R}_{12}:$

${\left(X_{s2}\right)}_{F_{m2}}^{T}E{\left(X_{s1}\right)}_{F_{m1}}=0---\left(15\right)$

公式（10）带入到公式（15）得到广义统一模型的极几何：

$\begin{array}{c}(\frac{{\mathrm{\ξ}}^T{K}^{-1}{p}_2+\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}{p}_2\right)}^2-\left(p_2^T{K}^{-T}{K}^{-1}{p}_2\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p_2^T{K}^{-T}{K}^{-1}{p}_2}{K}^{-1}{p}_2-\mathrm{\ξ})E\\ (\frac{{\mathrm{\ξ}}^T{K}^{-1}{p}_1+\sqrt{{\left({\mathrm{\ξ}}^T{K}^{-1}{p}_1\right)}^2-\left(p_1^T{K}^{-T}{K}^{-1}{p}_1\right)({\mathrm{\ξ}}^T\mathrm{\ξ}-1)}}{p_1^T{K}^{-T}{K}^{-1}{p}_1}{K}^{-1}{p}_1-\mathrm{\ξ})=0\end{array}---\left(16\right)$

在内参K和ξ已知的情况下，每一个点对根据公式（15）可以构造一个关于E的线性方程：

$(x_2^i,y_2^i,z_2^i)E{(x_1^i,y_1^i,z_1^i)}^T=0---\left(17\right)$

若设 $E=\left[\begin{array}{ccc}{e}_1& e_2& e_3\\ e_4& e_5& e_6\\ e_7& e_8& e_9\end{array}\right],$ 公式（17）可以转化为：

aⁱe＝0 （18）

其中 $a^i=[x_2^i{x}_1^i,x_2^i{y}_1^i,x_2^i{z}_1^i,y_2^i{x}_1^i,y_2^i{y}_1^i,y_2^i{z}_1^i,z_2^i{x}_1^i,z_2^i{y}_1^i,z_2^i{z}_1^i],$ e＝[e₁,e₂,e₃,e₄,e₅,e₆,e₇,e₈,e₉]，当对应点对的数目大于8时，可以求出e的最小二乘解，当然是包含一个尺度因子的。在计算E的过程中使用RANSAC提高精度，文献5（Fischler,M.A.and R.C.Bolles,Randomsample consensus:a paradigm for model fitting with applications to imageanalysis and automated cartography.Communications of the ACM,1981.24(6):p.381-395.）介绍了RANSAC方法，排除误差过大的外点。

因此知道E后可以分解出R₁₂和T₁₂(归一化T₁₂，T₁₂包含一尺度因子，这个尺度因子会使最终重构结果也包含一尺度因子)。首先将E进行奇异值分解：

E＝U diag(σ₁,σ₂,σ₃)V^T （19）

其中σ₁≥σ₂≥σ₃，U、V是单位正交矩阵。通过

$\left\{\begin{array}{c}{R}_{12}={\mathrm{UR}}_z^T(\±\frac{\mathrm{\π}}{2})V^T\\ {\hat{T}}_{12}={\mathrm{UR}}_z(\±\frac{\mathrm{\π}}{2})\mathrm{\Σ}{U}^T\end{array}\right.---\left(20\right)$

其中Σ＝diag(1,1,0)。

在内参K和ξ已知的情况下，并且在得到外参R₁₂和T₁₂的情况下，可以对空间点进行重构。在两幅图上两个像点分别为p₁＝(u₁,v₁,1)^T和p₂＝(u₂,v₂,1)^T，根据反向投影公式（10）可以得到对应的和先通过外参R₁₂和T₁₂将和C_m1转换到F_m2坐标系下：

${\left(C_{m1}\right)}_{F_{m2}}=R_{12}{C}_{m1}+T_{12}$

${\left(X_{s1}\right)}_{F_{m2}}=R_{12}{\left(X_{s1}\right)}_{F_{m1}}$

理想情况下射线和相交于一点，这点是X_m在F_m2下的位置，然而实际上由于误差的存在，射线和是异面直线，以异面直线和垂线的中点作为的估计值。在重构出空间点以后就可以得到空间点的重投影误差，通过前向投影可以得到X_m在左右图中的重投影p'₁＝(u'₁,v'₁,1)^T和p'₂＝(u'₂,v'₂,1)^T，可以定义重投影误差为：

$e_{\mathrm{reprojection}}=\frac{\left|(p_1^{\′}-p_1)\right|+\left|(p_2^{\′}-p_2)\right|}{2}---\left(21\right)$

内参初始化:

自标定过程是基于非线性优化的，精确的初值能够使得优化变快同时得到更好的优化结果，在估计内参初值时，为了较为简单的计算出初值，将折反射相机退化为单视点。在这种情况下，设k₁≈k₂≈k₃≈α≈ξ1≈ξ2≈0，γ₁≈γ₂≈γ。剩余所需要估计的内参为ξ3，γ，以及(u₀,v₀)。实验表明ξ3的初值对单位球模型的最终优化结果影响不大,设ξ3＝1。剩余所需要估计的内参为γ，以及(u₀,v₀)。(u₀,v₀)的初值使用提取镜面边缘的成像椭圆，以椭圆的中心为(u₀,v₀)的初值。首先使用canny算子提取出全向图的边缘，再根据边缘线的长短去除大部分非镜面边缘点，最后使用基于RANSAC的5点算法提取出镜面边缘椭圆，文献6（Xiang,Z.,B.Sun,and X.Dai,The Camera Itself as a Calibration Pattern:A NovelSelf-Calibration Method for Non-Central Catadioptric Cameras.Sensors,2012.12(6):p.7299-7317.）中详细阐述了镜面边缘的提取方法。

沿用前面的假设 $K=\left[\begin{array}{ccc}\mathrm{\γ}& 0& u_0\\ 0& \mathrm{\γ}& v_0\\ 0& 0& 1\end{array}\right],$ ξ＝[0,0,1]^T带入到极几何限制公式（16），可以简化为：

$(\frac{2}{p_2^T{K}^{-T}{K}^{-1}{p}_2}{K}^{-1}{p}_2-\mathrm{\ξ})E(\frac{2}{p_1^T{K}^{-T}{K}^{-1}{p}_1}{K}^{-1}{p}_1-\mathrm{\ξ})=0---\left(22\right)$

将 $K^{-1}{p}_j={[\frac{u_j-u_0}{\mathrm{\γ}},\frac{v_j-v_0}{\mathrm{\γ}},1]}^T={[\frac{u_{\mathrm{cj}}}{\mathrm{\γ}},\frac{v_{\mathrm{cj}}}{\mathrm{\γ}},1]}^T(j=\mathrm{1,2})$ 带入公式（22）：

$\frac{1}{(\frac{u_{c2}^2+v_{c2}^2}{{\mathrm{\γ}}^2}+1)}[\frac{u_{c2}}{\mathrm{\γ}},\frac{v_{c2}}{\mathrm{\γ}},\frac{1}{2}-\frac{u_{c2}^2+v_{c2}^2}{2{\mathrm{\γ}}^2}]E{[\frac{u_{c1}}{\mathrm{\γ}},\frac{v_{c1}}{\mathrm{\γ}},\frac{1}{2}-\frac{u_{c1}^2+v_{c1}^2}{2{\mathrm{\γ}}^2}]}^T\frac{1}{(\frac{u_{c1}^2+v_{c1}^2}{{\mathrm{\γ}}^2}+1)}=0---\left(23\right)$

简化后可以得到：

$[\mathrm{\γ}{u}_{c2},\mathrm{\γ}{v}_{c2},\frac{{\mathrm{\γ}}^2}{2}-\frac{u_{c2}^2+v_{c2}^2}{2}]E{[\mathrm{\γ}{u}_{c1},\mathrm{\γ}{v}_{c1},\frac{{\mathrm{\γ}}^2}{2}-\frac{u_{c1}^2+v_{c2}^2}{2}]}^T=0---\left(24\right)$

设 $E=\left[\begin{array}{ccc}{e}_1& e_2& e_3\\ e_4& e_5& e_6\\ e_7& e_8& e_9\end{array}\right],$ 重新整理24后可以得到：

[D₀+γD₁+γ²D₂+γ³D₃+γ⁴D₄]e＝0 （25）

其中

$D_0=[\mathrm{0,0,0,0,0,0,0,0},u_{c2}^2{u}_{c1}^2+u_{c2}^2{v}_{c1}^2+u_{c1}^2{v}_{c2}^2+v_{c2}^2{v}_{c1}^2]$

$D_1=[\mathrm{0,0},-2u_{c2}{v}_{c1}^2-2u_{c2}{u}_{c1}^2,\mathrm{0,0},-{2u}_{c1}^2{v}_{c2}-{2v}_{c2}{v}_{c1}^2,-{2u}_{c2}^2{u}_{c1}-{2u}_{c1}{v}_{c2}^2,-{2u}_{c2}^2{v}_{c1}-{2v}_{c2}^2{v}_{c1},0]$

$D_2=[{4u}_{c2}{u}_{c1},{4u}_{c2}{v}_{c1},0,{4u}_{c1}{v}_{c2},{4v}_{c2}{v}_{c1},\mathrm{0,0,0},-u_{c2}^2-u_{c1}^2-v_{c2}^2-v_{c1}^2]$

D₃=[0,0,2u_c2,0,0,2v_c2,2u_c1,2v_c1,0]

D₄=[0,0,0,0,0,0,0,0,1]

e＝[e₁,e₂,e₃,e₄,e₅,e₆,e₇,e₈,e₉]

方程（25）为多项式特征值问题(Polynomial Eigenvalue Problem)，文献7（Bai,Z.,Templates for the solution of algebraic eigenvalue problems:a practicalguide.Vol.11.2000:Siam.），任意9个点对可以解出e和γ,使用RANSAC来提高精度,求得γ的初值。

内参初值已知，并且知道E，可以分解出两幅图之间的外参R₁₂和T₁₂。在已知内参初值和外参的情况下，SIFT特征点对可以进行重构。由于内参并非真实值和噪声的存在，重构的点存在误差，将重构点通过前向投影再投影回两个成像平面，并与原始的特征点进行比较，就能得到重投影误差。公式（26）表示了重投影误差，j=1表示左图，j=2表示右图，函数X表示通过两个点重构的空间点，函数R表示重构点重投影回图像平面得到的图像点，函数H表示重投影误差。

$\begin{array}{c}H(\mathrm{\ξ}1,\mathrm{\ξ}2,\mathrm{\ξ}3,k_1,k_2,\mathrm{\α},{\mathrm{\γ}}_1,{\mathrm{\γ}}_2)=\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{2}\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{[R_j(\mathrm{\ξ}1,\mathrm{\ξ}2,\mathrm{\ξ}3,k_1,k_2,\mathrm{\α},{\mathrm{\γ}}_1,{\mathrm{\γ}}_2,X(p_{i1},p_{i2}))-p_{\mathrm{ij}}]}^2\end{array}---\left(26\right)$

使用非线性优化方法，最小化公式（26）可以优化出内参。

实施例：

在真实图像的测试中，使用了由NEOVISION制造的折反射相机系统。它原本是由H3S双曲镜面和Sony XCD-SX910CR相机组合成的单视点折反射相机。透视相机Sony XCD-SX910CR的内参为：f_x＝1455.07，f_y＝1459.51，k_s＝0，u₀＝639.2，v₀＝482.2。图像的分辨率为1280×960，但是该分辨率过大,导致处理的速度很慢,因此将图像分辨率缩减为800×600，对应的相机参数变为：f_x＝909.42，f_y＝912.2，k_s＝0，u₀＝399.5，v₀＝301.38。镜面的参数为：a＝0.0281m，b＝0.0234m，c＝0.0366m，d＝2c＝0.0731m，在镜面与相机构成单一视点模型的情况下其统一球面模型对应参数为： $\mathrm{\ξ}=\frac{d}{\sqrt{d^2+4p^2}}=0.9665,\mathrm{\η}=\frac{-2p}{\sqrt{d^2+4p^2}}=0.2565,$

γ₁＝f_xη＝233.27，γ₂＝f_yη＝233.98。这里的参数是镜面和相机构成单一视点模型的情况下，而实际当中相机焦点未知，镜面焦点未知，安装存在误差，所装配折反射相机必然存在失准。

使用上述方法进行标定，得到自标定结果为[ξ1,ξ2,ξ3]＝[0.0068,0.10132,1.3378]，[u₀,v₀]＝[384.3743,293.8175]，[γ₁,γ₂]＝[339.2043,336.1176]，[k₁,k₂]＝[-0.1965,0.2895]，重投影误差为0.1518像素，重投影误差的方差为0.1411，分解得到的外参 $R=\left[\begin{array}{ccc}0.9943& 0.1063& 0.0027\\ -0.1063& 0.9943& -0.0035\\ -0.0030& 0.0032& 1.0000\end{array}\right],$ 归一化以后的T＝[-0.7294 -0.6839 0.0192]^T。重投影误差小于0.5个像素，表明标定结果精确。

Claims (4)

1.一种基于广义统一模型的折反射相机自标定方法，其特征在于，该方法包括以下步骤：

步骤1)在同一场景不同位置使用折反射相机获取两幅左右成像清晰的折反射全景图,要求两图中有场景重复的部分，不需要对环境做特殊的设置，也不需要任何标定模板；

步骤2)提取左右两图的SIFT特征点，并进行特征点匹配，得到图像上n对匹配的特征点对和i＝1..n，n对特征点对中存在误配点；

步骤3)初始化广义统一模型的内参；

步骤4)使用极几何计算出本质矩阵E，在计算E的过程中使用RANSAC提高精度，同时去除外点；

步骤5)使用非线性优化Levenberg Marquardt方法，最小化特征点的重投影误差，从而优化出内参；

步骤6)重新回到步骤4)，循环步骤4)和步骤5)，不断的优化内参，同时得到更多的内点，直到重投影误差小于所设定的阈值，优化出鲁棒和精确的结果；

步骤3)中所述初始化广义统一模型的内参，首先使用广义统一模型对折反射相机进行建模；

折反射相机广义统一模型由虚拟的单位球，虚拟视点，归一化平面和成像平面П_p组成；以单位球中心C_m为原点建立坐标系F_m，以虚拟光心C_p为原点建立坐标系F_p，两坐标系相互平行，归一化平面为F_p中的z＝1平面，归一化平面上的点通过径向失真和透视变换投影到成像平面П_p；

折反射相机广义统一模型，有两个坐标系，以虚拟光心C_p为原点建立坐标系F_p，和以单位球中心C_m为原点建立坐标系F_m，原点C_p在F_m坐标系中的坐标为ξ＝(ξ1,ξ2,ξ3)^T；

空间点X_m在世界坐标系F_w中的表示为(X_m)_W＝(x_w,y_w,z_w)，通过以下前向投影过程，即成像过程如下：

(1)空间点从世界坐标系F_w转换到F_m坐标系，其中R_w和T_w为世界坐标系和F_m坐标系的转换关系，此时空间点在F_m坐标系中的坐标为

(2)将投影到单位球面上坐标为(x_s,y_s,z_s)^T，

(3)转换到F_p坐标系下其中ξ＝(ξ1,ξ2,ξ3)^T，

(4)将投影到归一化平面的m_u坐标为(x,y,1)^T，其中λ＝z_s+ξ3，h( )为折反射投影函数；

(5)在归一化平面上加入径向失真以弥补真实系统中透视相机镜头的失真，径向失真函数为径向畸变函数L(ρ)＝1+k₁ρ²+k₂ρ⁴+k₃ρ⁶，其中，k₁,k₂,k₃为失真系数，失真后的点m'_u＝(x',y',1)^T，其中x'＝x*L(ρ),y'＝y*L(ρ)；

(6)通过透视变换将归一化平面上的点m'_u投影到成像平面П_p上点p坐标为(u,v,1)^T，其中k( )为透视投影函数，f₁，f₂是相机的焦距，(u₀,v₀)是相机的主点，α是偏斜，η是镜面参数；此处的K为归一化平面到成像平面的透视变换矩阵，不再将折反射系统看成是一个独立的镜面和一个独立的相机，而是将之视为一个整体，内参矩阵K包含了相机和镜面的参数，f₁，f₂和η无法单独得到，只能得到γ₁＝f₁η，γ₂＝f₂η；

从图像上的点p反向推导出相应的单位球面上的点的方向的过程称为后向投影过程：

(1)从图像点p＝(u,v,1)^T映射到归一化平面上的点m'_u，m′_u＝K^-1p＝(x,y,1)^T；

(2)对m′_u去除径向畸变，径向畸变函数L(ρ′)＝1+k₁ρ′²+k₂ρ′⁴+k₃ρ′⁶，k₁,k₂,k₃为失真系数，去失真后的点m_u＝(x,y,1)^T，其中x＝x'/L(ρ′),y＝y′/L(ρ′)

(3)从归一化平面上的点m_u映射到单位球面上的点 其中h( )是投影函数，h^-1( )是投影函数的逆函数；

由前向投影过程中的步骤(4)(3)(2)解出从图像点到空间射线的反向投影 方程此处未加入去除径向畸变，因为径向畸变非线性，加入径向畸变以后后向投影公式过于复杂，实际使用当中在初始化阶段忽略径向畸变，在优化过程中分步进行后向投影的计算而不是利用反向投影方程一次性得出；

通过移动折反射相机，在不同的位置采集空间点的成像，空间点X_m在世界坐标系F_w中的表示为(X_m)_W＝(x_w,y_w,z_w)，在不同位置折反射相机中成像，成像过程如前面前向投影所述，在两幅图上形成的成像点分别为p₁＝(u₁,v₁,1)^T和p₂＝(u₂,v₂,1)^T，根据反向投影方程计算出p₁和p₂对应的和射线和都经过空间点X_m，得到：其中本质矩阵E＝T₁₂×R₁₂，其中R₁₂和T₁₂为第一个相机和第二个相机之间的坐标转换关系；接着得到广义统一模型的极几何公式：

广义统一模型含有需要标定的参数为，k₁,k₂,k₃，ξ＝(ξ1,ξ2,ξ3)^T，γ₁，γ₂，[u₀,v₀]和α；

其次，为了获取模型的初值，先假设折反射全向相机为单一视点模型，设k₁≈k₂≈k₃≈α≈ξ1≈ξ2≈0，γ₁≈γ₂＝γ，所需要估计的内参为ξ3、γ、u₀、v₀，实验表明ξ3对系统参数的初值计算影响不大，设ξ3＝1，主点[u₀,v₀]的初始值估计，先提取图像中镜面边缘，再取图像中镜面的中心为主点的初值；

将和ξ＝[0,0,1]^T带入极几何公式得到将带入得到设重新整理后得到多项式特征值问题[D₀+γD₁+γ²D₂+γ³D₃+γ⁴D₄]e＝0，其中

D₃＝[0,0,2u_c2,0,0,2v_c2,2u_c1,2v_c1,0]

D₄＝[0,0,0,0,0,0,0,0,1]

e＝[e₁,e₂,e₃,e₄,e₅,e₆,e₇,e₈,e₉]

该多项式特征值问题通过9个特征值点对解出e和γ，使用RANSAC来提高精度,求得γ的初值。

2.根据权利要求1所述的一种基于广义统一模型的折反射相机自标定方法，其特征在于：所述步骤4)具体为：

在内参初值K和ξ已知的情况下,每一对特征点带入极几何公式，得到一个关于E的线性方程其中从反向投影方程得到，j＝1,2，j＝1表示左图，j＝2表示右图,i＝1..n；若设所述E的线性方程转换为aⁱe＝0，其中e＝[e₁,e₂,e₃,e₄,e₅,e₆,e₇,e₈,e₉]；当特征点对数目大于等于8时，求出e的最小二乘解，包含一个尺度因子；在计算E的过程中使用RANSAC提高精度，排除误差过大的外点；

奇异值分解E，E＝U diag(σ₁,σ₂,σ₃)V^T，其中σ₁≥σ₂≥σ₃，U、V是单位正交矩阵； 得到R₁₂和T₁₂，其中Σ＝diag(1,1,0)，T₁₂包含一尺度因子；选择具有物理意义的那一组R₁₂和T₁₂，即两幅图对应相机朝向相近的那一组。

3.根据权利要求1所述的一种基于广义统一模型的折反射相机自标定方法，其特征在于：所述步骤5)具体为：

使用非线性优化方法Levenberg Marquardt方法,通过最小化特征点的重投影误差来优化出内参；

要计算特征点的重投影误差，要对特征点进行三维重构；在内参K和ξ已知的情况下，并且在得到外参R₁₂和T₁₂的情况下，对空间点进行重构；对任意特征点对p₁＝(u₁,v₁,1)^T和p₂＝(u₂,v₂,1)^T，根据反向投影方程，得到对应的空间入射光线的方向，即和通过外参R₁₂和T₁₂，可将两光线变换到同一坐标系下和理想情况下两光线相交，实际当中由于存在误差，两光线为异面直线，取异面直线中垂线的中点X'_m为重构点；通过前向投影得到X'_m在左右两图中的重投影p'₁＝(u'₁,v'₁,1)^T和p'₂＝(u'₂,v'₂,1)^T；这对特征点的重投影误差就为

使用非线性优化Levenberg Marquardt方法，最小化所有特征点的重投影误差的和，从而优化出内参。

4.根据权利要求1所述的一种基于广义统一模型的折反射相机自标定方法，其特征在于：所述的步骤6)具体为：

步骤4)中所用的内参初值从步骤3)中得到，初始时初值误差较大，导致步骤4)去除外点过多，进而影响了步骤5)中的优化结果的精度；在经过步骤5)的优化后，内参的值比初始时要精确，使用该内参去重新进行步骤4)和步骤5)；如此重复，不断优化得到更精确的内参，直到重投影误差小于说要求的阈值。