[发明专利]一种基于畸变校正的像机焦距求解方法

摘要

本发明提出了一种基于畸变校正的像机焦距求解方法，包括图像畸变线性模型的建立和像机焦距的求解两个基本步骤。步骤一，首先建立像机成像模型，进一步通过公式推导，求解畸变前后像点间的关系式，进而建立图像畸变线性模型；步骤二，首先建立像机标定层次化模型，进而求取像机最优焦距。本发明通过建立图像畸变线性模型，进而层次化的求解像机焦距，可灵活有效的解决像机自标定过程中同时求解像机内部参数和畸变系数的问题。

主权项

一种基于畸变校正的像机焦距求解方法，其特征在于，包括图像畸变线性模型的建立和像机焦距的求解两个步骤，所述图像畸变线性模型的建立包括：(1)建立像机成像模型：将世界坐标系中的点X＝(X,Y,Z)^T投影到图像平面上，表示为x＝(x,y)^T，将其写成齐次坐标的形式X＝(X,Y,Z,1)^T、x＝(x,y,1)^T，则满足投影方程：x～PX   (1)其中，P为3×4的投影矩阵，通常表示为：P＝K[R|t]   (2)其中，R、t分别为像机坐标系相对于世界坐标系的旋转矩阵和平移向量，K为像机的内部参数矩阵，K可以写成如下形式：$K=\left[\begin{array}{ccc}\mathrm{\τ}f& c& u_0\\ 0& f& v_0\\ 0& 0& 1\end{array}\right]---\left(3\right)$其中，f为像机的焦距，(u₀,v₀)为主点坐标，c为歪斜量(通常c＝0)；实际上，镜头的成像都带有不同程度的畸变，因此，理论成像点x＝(x,y)^T在受到镜头畸变影响后的实际像点为x′＝(x′,y′)^T，二者间的关系见式(4)：x＝x′+δx              (4)y＝y′+δy其中，δx和δy为非线性畸变值，理论上镜头会同时存在径向畸变和切向畸变，但切向畸变比较小，可忽略，径向畸变的表达式见式(5)：δx＝(x′‑u₀)(k₁r²+k₂r⁴+…)                           (5)δy＝(y′‑v₀)(k₁r²+k₂r⁴+…)其中，r²＝(x′‑u₀)²+(y′‑v₀)²；一般而言，一阶径向畸变已足够描述非线性畸变模型，因此，式(5)可写成：δx＝(x′‑u₀)kr²                 (6)δy＝(y′‑v₀)kr²像机非线性模型的内部参数由线性模型参数f、τ、(u₀,v₀)和非线性畸变参数k共同构成；(2)畸变参数的求解：假定(u,v)是真实的像点坐标，和分别为主点坐标校正到图像中心、且纵横轴坐标归一化到相同比例的理想像点和畸变后的像点，(x,y)和是归一化到理想像机下的理想像点和畸变后的像点；通常只考虑一阶径向畸变，存在如下畸变公式：其中，k为径向畸变系数；对于主点坐标校正到图像中心、且纵横轴坐标归一化到相同比例的畸变后的像点坐标存在以下关系：则有：其中，$\tilde{k}=k/f^2;$考虑极坐标系：$\tilde{u}=\tilde{r}cos\widetilde{\mathrm{\θ}}$$\tilde{v}=\tilde{r}sin\widetilde{\mathrm{\θ}}$                    (10)对于径向畸变，存在式(11)：代入上述畸变公式可得：整理可得：已知和通常情况下，可用如下求根公式计算当时，则归一化后的坐标可由畸变后的坐标求得：将式(15)写成如下形式：其中：则分别对式(18)和式(19)求偏导可得：因此：同理，将式(16)写成如下形式：则有：其中，并且，已知：因此，利用泰勒公式展开可得：上式省略二阶无穷小可得：将式(33)、式(34)写成矩阵形式：所述像机焦距的求解包括：(1)建立层次化模型：由式(3)可知，像机的内部参数矩阵可写为：$K=\left[\begin{array}{ccc}\mathrm{\τ}f& 0& u_0\\ 0& f& v_0\\ 0& 0& 1\end{array}\right]---\left(36\right)$基本矩阵为：$\begin{array}{c}F=K^{-T}{\[t\]}_{\×}{\mathrm{RK}}^{-1}=K={\left[\begin{array}{ccc}\mathrm{\τ}f& 0& u_0\\ 0& f& v_0\\ 0& 0& 1\end{array}\right]}^{-T}{\[t\]}_{\×}R{\left[\begin{array}{ccc}\mathrm{\τ}f& 0& u_0\\ 0& f& v_0\\ 0& 0& 1\end{array}\right]}^{-T}\\ ={\left[\begin{array}{ccc}\mathrm{\τ}& 0& u_0\\ 0& 1& v_0\\ 0& 0& 1\end{array}\right]}^{-T}{\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]}^{-T}{\[t\]}_{\×}R{\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]}^{-T}{\left[\begin{array}{ccc}\mathrm{\τ}& 0& u_0\\ 0& 1& v_0\\ 0& 0& 1\end{array}\right]}^{-T}\end{array}---\left(37\right)$其中：t＝[t_x t_y t_z]^T   (38)${\[t\]}_{\×}=\left[\begin{array}{ccc}0& -t_z& t_y\\ t_z& 0& -t_x\\ -t_y& t_x& 0\end{array}\right]---\left(39\right)$由基本矩阵的基本性质可知：x^TFx'＝0   (40)其中，是两幅图像的任意一对理想匹配点；将式(37)代入式(40)可得：$x^T{\left[\begin{array}{ccc}\mathrm{\τ}& 0& u_0\\ 0& 1& v_0\\ 0& 0& 1\end{array}\right]}^{-T}{\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]}^{-T}{\[t\]}_{\×}R{\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]}^{-1}{\left[\begin{array}{ccc}\mathrm{\τ}& 0& u_0\\ 0& 1& v_0\\ 0& 0& 1\end{array}\right]}^{-1}{x}^{\′}=0---\left(41\right)$上式可简写为：${\tilde{x}}^{T}G{\tilde{x}}^{\′}=0---\left(42\right)$其中，$\tilde{x}={\left[\begin{array}{ccc}\mathrm{\τ}& 0& u_0\\ 0& 1& v_0\\ 0& 0& 1\end{array}\right]}^{-1}x,$${\tilde{x}}^{\′}={\left[\begin{array}{ccc}\mathrm{\τ}& 0& u_0\\ 0& 1& v_0\\ 0& 0& 1\end{array}\right]}^{-1}{x}^{\′},$$G={\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]}^{-T}{\[t\]}_{\×}R{\left[\begin{array}{ccc}f& 0& 0\\ 0& f& 0\\ 0& 0& 1\end{array}\right]}^{-1};$G中只包含了像机的外部参数和焦距；根据上述对畸变的分析，式可重写为：${\left(\begin{array}{c}{u}_0\\ v_0\\ 1\end{array}\right)}^T{M}^T(u,v,\mathrm{\τ},\tilde{k}){\mathrm{GM}}^{\′}(u,v,\mathrm{\τ},\tilde{k})\left(\begin{array}{c}{u}_0\\ v_0\\ 1\end{array}\right)=0---\left(43\right)$式中：由归一化后的图像及其匹配点，可通过规范化的八点法求解G矩阵；(2)求取像机最优焦距：首先，对G进行奇异值分解可得：G＝PΛW^T   (46)其中，Λ＝diag(a,b,0)，(a,b>0)P和W为两个正交矩阵；$P=\left[\begin{array}{ccc}{p}_{11}& p_{12}& p_{13}\\ p_{21}& p_{22}& p_{23}\\ p_{31}& p_{32}& p_{33}\end{array}\right],$$W=\left[\begin{array}{ccc}{W}_{11}& W_{12}& W_{13}\\ W_{21}& W_{\mathrm{22}}& W_{\mathrm{23}}\\ W_{31}& W_{\mathrm{32}}& W_{\mathrm{33}}\end{array}\right]$将上式代入Kruppa方程可得：Gdiag(f²,f²,1)G^T∝[e']_×diag(f²,f²,1)[e']_×^T   (47)其中，[e']_×为e'的反对称阵；进一步上式可写为：PΛW^Tdiag(f²,f²,1)WΛP^T∝[p₃]_×diag(f²,f²,1)[p₃]_×^T   (48)其中，p₃＝[p₁₃ p₂₃ p₃₃]，[p₃]_×为p₃的反对称阵；式(48)可分解为如下三个方程：f²(ap₁₃p₂₃(1‑w₁₃²)+bw₁₃w₂₃(1‑p₂₃²))+p₂₃w₁₃(ap₁₃w₁₃+bp₂₃w₂₃)＝0   (49)f²(aw₁₃w₂₃(1‑p₁₃²)+bp₁₃p₂₃(1‑w₂₃²))+p₁₃w₂₃(ap₁₃w₁₃+bp₂₃w₂₃)＝0   (50)f⁴(a²(1‑p₁₃²)(1‑w₁₃²)‑b²(1‑p₂₃²)(1‑w₂₃²))+f²(a²(p₁₃²+w₁₃²‑2p₁₃²w₁₃²)   (51)‑b²(p₂₃²+w₂₃²‑2p₂₃²w₂₃²))+(a²p₁₃²w₁₃²‑b²p₂₃²w₂₃²)＝0上三式为基本的标定方程；通常，由式(51)即可求得焦距，式(49)和式(50)作为验证标定结果的方程；从左右两个视角可得到n幅图像，进而可得到个基本矩阵，个G矩阵；并且，式(51)可写为：Af⁴+Bf²+C＝0   (52)其中，A＝a²(1‑p₁₃²)(1‑w₁₃²)‑b²(1‑p₂₃²)(1‑w₂₃²)，C＝(a²p₁₃²w₁₃²‑b²p₂₃²w₂₃²)，B＝a²(p₁₃²+w₁₃²‑2p₁₃²w₁₃²)‑b²(p₂₃²+w₂₃²‑2p₂₃²w₂₃²)；对个G矩阵进行奇异值分解，可得到下列方程组：$\left\{\begin{array}{c}{A}_1{f}^4+B_1{f}^2+C_1=0\\ A_2{f}^4+B_2{f}^2+C_2=0\\ \·\\ \·\\ \·\\ A_n{f}^4+B_n{f}^2+C_n=0\end{array}\right.n=1,2...C_n^2---\left(53\right)$进而，通过最小二乘法求得最优的焦距。