1.相机标定：坐标系转换学习
2.相机的畸变介绍
3.相机标定原理(不考虑径向畸变)
A Flexible New Technique for Camera Calibration 张氏标定论文原文
在第一篇博客中已经对坐标系转换做了详细介绍。在张氏标定原理中便不再介绍，直接使用其结论。请配合论文原文进行食用。
[ P x i ∗ Z c P y i ∗ Z c Z c ] = [ M ∗ f − M ∗ f ∗ c o t θ M ∗ P x 0 N ∗ f / s i n θ N ∗ P y 0 0 1 ] ∗ ( R ∣ − R C ) ∗ [ x w y w z w 1 ] left[begin{array}{ccccc} P_{x_i}*Z_c\ P_{y_i}*Z_c\ Z_c end{array}right]= left[ begin{array}{ccccc} M*f&-M*f*cotθ&M*P_x\ 0&N*f/sinθ&N*P_y\ 0&0&1 end{array} right]*(R|-RC)* left[begin{array}{ccccc} x_w\ y_w\ z_w\ 1 end{array}right] ⎣⎡Pxi∗ZcPyi∗ZcZc⎦⎤=⎣⎡M∗f00−M∗f∗cotθN∗f/sinθ0M∗PxN∗Py1⎦⎤∗(R∣−RC)∗⎣⎢⎢⎡xwywzw1⎦⎥⎥⎤
如上式是第一篇博客中得出来的结论。为了和张氏标定论文中的符号能对照起来，这里做如下变换。
令 m ~ = [ P x i , P y i , 1 ] T = [ u , v , 1 ] T widetilde{m}=[P_{x_i},P_{y_i},1]^T=[u,v,1]^T m =[Pxi,Pyi,1]T=[u,v,1]T
令 s = Z c s=Z_c s=Zc
令 M ~ = [ x w , y w , z w , 1 ] T = [ X , Y , Z , 1 ] T widetilde{M}=[x_w,y_w,z_w,1]^T=[X,Y,Z,1]^T M =[xw,yw,zw,1]T=[X,Y,Z,1]T。
令内参矩阵 K K K转换成如下矩阵：
A = K = [ M ∗ f − M ∗ f ∗ c o t θ M ∗ P x 0 N ∗ f / s i n θ N ∗ P y 0 0 1 ] = [ α γ u 0 0 β v 0 0 0 1 ] A=K= left[begin{array}{ccccc} M*f&-M*f*cotθ&M*P_x\ 0&N*f/sinθ&N*P_y\ 0&0&1 end{array}right]= left[begin{array}{ccccc} α&γ&u_0\ 0&β&v_0\ 0&0&1 end{array}right] A=K=⎣⎡M∗f00−M∗f∗cotθN∗f/sinθ0M∗PxN∗Py1⎦⎤=⎣⎡α00γβ0u0v01⎦⎤
令平移矩阵 − R C = t -RC=t −RC=t。
则整个坐标转换公式转化如下：
s ∗ m ~ = A ∗ [ R ∣ t ] ∗ M ~ s ∗ [ u v 1 ] = [ α γ u 0 0 β v 0 0 0 1 ] [ R ∣ t ] [ X Y Z 1 ] s*widetilde{m}=A*[R|t]*widetilde{M}\ s* left[ begin{array}{ccccc} u\v\1 end{array} right]= left[ begin{array}{ccccc} α&γ&u_0\ 0&β&v_0\ 0&0&1 end{array} right][R|t] left[ begin{array}{ccccc} X\Y\Z\1 end{array} right] s∗m =A∗[R∣t]∗M s∗⎣⎡uv1⎦⎤=⎣⎡α00γβ0u0v01⎦⎤[R∣t]⎣⎢⎢⎡XYZ1⎦⎥⎥⎤
这里称呼 s s s为一个尺度因子， A A A为内参矩阵， [ R ∣ t ] [R|t] [R∣t]为外参矩阵； m ~ widetilde{m} m 为2D点的齐次坐标， M ~ widetilde{M} M 为3D点的齐次坐标。
这里需要指出的是：
s s s 作为尺度因子相对于一个定点是一个定值，而并非任意值。
张氏标定法是目前求解相机内、外参数和畸变系数的常用方法。
我们在3.相机标定原理(不考虑径向畸变)中讲到的标定方法，需要一个比较复杂的标定装置，要求三个面互相垂直。而张氏标定法的优点就在于，使用一个打印的图片就可以完成相机的内参数和畸变系数的标定。
张氏标定是基于一个假设而进行的，即假设3D点 [ X , Y , Z ] T [X,Y,Z]^T [X,Y,Z]T的 Z Z Z坐标为0。
至于这个假设是否会影响到标定得出来的真实结果(显然不会，因为在大量实验的基础上已经证实了张氏标定的准确性)，抑或是为什么这个假设不会影响到标定的真实结果？暂且按下不表，这是后话。
这个冗长乏味的推导，将分为如下几个简短的部分，以方便大家观看。
于是在这个假设之下，坐标转化公式就变成了如下形式， r 1 , r 2 , r 3 , t r_1,r_2,r_3,t r1,r2,r3,t均为三维列向量：
s ∗ [ u v 1 ] = [ α γ u 0 0 β v 0 0 0 1 ] [ r 1 r 2 r 3 t ] [ X Y 0 1 ] s*left[ begin{array}{ccccc} u\v\1 end{array} right]= left[ begin{array}{ccccc} α&γ&u_0\ 0&β&v_0\ 0&0&1 end{array} right] left[ begin{array}{ccccc} r_1&r_2&r_3&t end{array} right] left[ begin{array}{ccccc} X\Y\0\1 end{array} right] s∗⎣⎡uv1⎦⎤=⎣⎡α00γβ0u0v01⎦⎤[r1r2r3t]⎣⎢⎢⎡XY01⎦⎥⎥⎤
很容易知道， r 3 r_3 r3由于 Z = 0 Z=0 Z=0的影响，在计算中起不到任何作用，于是对上述公式进行简化：
s ∗ [ u v 1 ] = [ α γ u 0 0 β v 0 0 0 1 ] [ r 1 r 2 t ] [ X Y 1 ] s*left[ begin{array}{ccccc} u\v\1 end{array} right]= left[ begin{array}{ccccc} α&γ&u_0\ 0&β&v_0\ 0&0&1 end{array} right] left[ begin{array}{ccccc} r_1&r_2&t end{array} right] left[ begin{array}{ccccc} X\Y\1 end{array} right] s∗⎣⎡uv1⎦⎤=⎣⎡α00γβ0u0v01⎦⎤[r1r2t]⎣⎡XY1⎦⎤
更新 M ~ = [ X , Y , 1 ] T widetilde{M}=[X,Y,1]^T M =[X,Y,1]T。
不妨令：
H = [ α γ u 0 0 β v 0 0 0 1 ] [ r 1 r 2 t ] = A [ r 1 r 2 t ] = [ h 1 h 2 h 3 ] H= left[ begin{array}{ccccc} α&γ&u_0\ 0&β&v_0\ 0&0&1 end{array} right] left[ begin{array}{ccccc} r_1&r_2&t end{array} right]=A left[ begin{array}{ccccc} r_1&r_2&t end{array} right] =left[ begin{array}{ccccc} h_1&h_2&h_3 end{array} right] H=⎣⎡α00γβ0u0v01⎦⎤[r1r2t]=A[r1r2t]=[h1h2h3]
易知 H H H为3×3的矩阵，令 H = [ h 1 h 2 h 3 ] H=left[ begin{array}{ccccc} h_1&h_2&h_3 end{array} right] H=[h1h2h3]，称为单应性矩阵，我更愿意称之为映射变换矩阵。其反映了一个点到另外一个点之间的映射关系。
对照原文的时候我们会发现，在这里多出了一个任意标量 λ λ λ。而且在论文的后文中 λ λ λ出现的位置又不对，总之挺费解的，目前尚未搞清楚这个 λ λ λ的作用，因此在本文中，默认 λ = 1 λ=1 λ=1，这并不影响我们对张氏标定方法的理解。
已知我们的旋转矩阵 R R R是一个单位正交矩阵，即其行与行(列与列)之间互相正交即点乘为0，且每行(列)的模为1。用数学公式表达如下：
r 1 T ⋅ r 2 = 0 r 1 T . r 1 = ∣ ∣ r 1 ∣ ∣ = 1 r_1^T·r_2=0\ r_1^T.r_1=||r_1||=1 r1T⋅r2=0r1T.r1=∣∣r1∣∣=1
在做 H H H的代换的时候我们得到如下的一个方程：
A [ r 1 r 2 t ] = [ h 1 h 2 h 3 ] A left[ begin{array}{ccccc} r_1&r_2&t end{array} right] =left[ begin{array}{ccccc} h_1&h_2&h_3 end{array} right] A[r1r2t]=[h1h2h3]
应用旋转矩阵的两个性质，我们做如下计算：
[ r 1 r 2 t ] = A − 1 [ h 1 h 2 h 3 ] = [ A − 1 h 1 A − 1 h 2 A − 1 h 3 ] left[ begin{array}{ccccc} r_1&r_2&t end{array} right] =A^{-1} left[ begin{array}{ccccc} h_1&h_2&h_3 end{array} right]= left[ begin{array}{ccccc} A^{-1}h_1&A^{-1}h_2&A^{-1}h_3 end{array} right] [r1r2t]=A−1[h1h2h3]=[A−1h1A−1h2A−1h3]
则：
r 1 T ⋅ r 2 = ( A − 1 h 1 ) T A − 1 h 2 = h 1 T A − T A − 1 h 2 = 0 r 1 T . r 1 = ∣ ∣ r 1 ∣ ∣ = h 1 T A − T A − 1 h 1 = 1 r 2 T . r 2 = ∣ ∣ r 2 ∣ ∣ = h 2 T A − T A − 1 h 2 = 1 r_1^T·r_2=(A^{-1}h_1)^TA^{-1}h_2=h_1^TA^{-T}A^{-1}h_2=0\ r_1^T.r_1=||r_1||=h_1^TA^{-T}A^{-1}h_1=1\ r_2^T.r_2=||r_2||=h_2^TA^{-T}A^{-1}h_2=1 r1T⋅r2=(A−1h1)TA−1h2=h1TA−TA−1h2=0r1T.r1=∣∣r1∣∣=h1TA−TA−1h1=1r2T.r2=∣∣r2∣∣=h2TA−TA−1h2=1
令 B = A − T A − 1 B=A^{-T}A^{-1} B=A−TA−1，则 h 1 T B h 2 = 0 , h 1 T B h 1 = h 2 T B h 2 h_1^TBh_2=0,h_1^TBh_1=h_2^TBh_2 h1TBh2=0,h1TBh1=h2TBh2。
接下来我们要计算 A − 1 A^{-1} A−1。
求取矩阵的逆通常有三种方法，第一种待定系数法，第二种伴随矩阵法，第三种初等变换法。由于我们的 A A A中含有多个未知量，使用待定系数法和初等变换法都不是很方便，无奈之下选择伴随矩阵法来进行计算。
这里就不再介绍伴随矩阵的求法了，直接给出结果！
A ∗ = [ β − γ γ v 0 − β u 0 0 α − α v 0 0 0 α β ] A^*= left[begin{array}{ccccc} β&-γ&γv_0-βu_0\ 0&α&-αv_0\ 0&0&αβ end{array}right] A∗=⎣⎡β00−γα0γv0−βu0−αv0αβ⎦⎤
伴随矩阵法求逆还需要求 ∣ ∣ A ∣ ∣ ||A|| ∣∣A∣∣。容易求得 ∣ ∣ A ∣ ∣ = α β ||A||=αβ ∣∣A∣∣=αβ。
于是乎：
A − 1 = [ 1 α − γ α β γ v 0 − β u 0 α β 0 1 β − v 0 β 0 0 1 ] ， A − T = [ 1 α 0 0 − γ α β 1 β 0 γ v 0 − β u 0 α β − v 0 β 1 ] A^{-1}= left[begin{array}{ccccc} frac{1}α&-frac{γ}{αβ}&frac{γv_0-βu_0}{αβ}\ 0&frac{1}β&-frac{v_0}β\ 0&0&1 end{array}right]， A^{-T}= left[begin{array}{ccccc} frac{1}α&0&0\ -frac{γ}{αβ}&frac{1}β&0\ frac{γv_0-βu_0}{αβ}&-frac{v_0}β&1 end{array}right] A−1=⎣⎡α100−αβγβ10αβγv0−βu0−βv01⎦⎤，A−T=⎣⎡α1−αβγαβγv0−βu00β1−βv0001⎦⎤
于是乎：
B = A − T A − 1 = [ 1 α 2 − γ α 2 β γ v 0 − β u 0 α 2 β − γ α 2 β γ 2 α 2 β 2 + 1 β 2 − γ ( γ v 0 − β u 0 ) α 2 β 2 − v 0 β 2 γ v 0 − β u 0 α 2 β − γ ( γ v 0 − β u 0 ) α 2 β 2 − v 0 β 2 ( γ v 0 − β u 0 ) 2 ( α β ) 2 + v 0 2 β 2 + 1 ] = [ B 11 B 12 B 13 B 21 B 22 B 23 B 31 B 32 B 33 ] B=A^{-T}A^{-1}= left[begin{array}{ccccc} frac{1}{α^2}&-frac{γ}{α^2β}&frac{γv_0-βu_0}{α^2β}\ -frac{γ}{α^2β}&frac{γ^2}{α^2β^2}+frac{1}{β^2}&-frac{γ(γv_0-βu_0)}{α^2β^2}-frac{v_0}{β^2}\ frac{γv_0-βu_0}{α^2β}&-frac{γ(γv_0-βu_0)}{α^2β^2}-frac{v_0}{β^2}&frac{(γv_0-βu_0)^2}{(αβ)^2}+frac{v_0^2}{β^2}+1 end{array}right]= left[begin{array}{ccccc} B_{11}&B_{12}&B_{13}\ B_{21}&B_{22}&B_{23}\ B_{31}&B_{32}&B_{33} end{array}right] B=A−TA−1=⎣⎢⎡α21−α2βγα2βγv0−βu0−α2βγα2β2γ2+β21−α2β2γ(γv0−βu0)−β2v0α2βγv0−βu0−α2β2γ(γv0−βu0)−β2v0(αβ)2(γv0−βu0)2+β2v02+1⎦⎥⎤=⎣⎡B11B21B31B12B22B32B13B23B33⎦⎤
对比上述式子，我们容易发现， B B B是一个对称矩阵，即 B 12 = B 21 , B 13 = B 31 , B 23 = B 32 B_{12}=B_{21},B_{13}=B_{31},B_{23}=B_{32} B12=B21,B13=B31,B23=B32
已知 B B B是一个对称矩阵，其未知量共有6个，为 B 11 , B 21 , B 22 , B 31 , B 32 , B 33 B_{11},B_{21},B_{22},B_{31},B_{32},B_{33} B11,B21,B22,B31,B32,B33。而这6个未知量又是使用5个参数表示的，因此准确说来，未知量的个数为5个。
我们做如下的代换，用向量 b b b来表示矩阵 B B B的未知量，同样的向量 b b b的自由度为5。：
b = [ B 11 B 21 B 22 B 31 B 32 B 33 ] b= left[begin{array}{ccccc} B_{11}\B_{21}\B_{22}\B_{31}\B_{32}\B_{33} end{array}right] b=⎣⎢⎢⎢⎢⎢⎢⎡B11B21B22B31B32B33⎦⎥⎥⎥⎥⎥⎥⎤
我们在4.2.3中已知 h 1 T B h 2 = 0 , h 1 T B h 1 = h 2 T B h 2 h_1^TBh_2=0,h_1^TBh_1=h_2^TBh_2 h1TBh2=0,h1TBh1=h2TBh2。
令 h 1 = [ h 11 h 12 h 13 ] h_1=left[begin{array}{ccccc}h_{11}\h_{12}\h_{13}end{array}right] h1=⎣⎡h11h12h13⎦⎤， h 2 = [ h 21 h 22 h 23 ] h_2=left[begin{array}{ccccc}h_{21}\h_{22}\h_{23}end{array}right] h2=⎣⎡h21h22h23⎦⎤
带入计算，结果用向量 b b b来表示：
h 1 T B h 2 = [ h 11 h 12 h 13 ] [ B 11 B 12 B 13 B 21 B 22 B 23 B 31 B 32 B 33 ] [ h 21 h 22 h 23 ] h 1 T B h 2 = [ h 11 h 21 h 12 h 21 + h 11 h 22 h 12 h 22 h 13 h 21 + h 11 h 23 h 13 h 23 + h 12 h 23 h 13 h 23 ] [ B 11 B 21 B 22 B 31 B 32 B 33 ] h_1^TBh_2=left[begin{array}{ccccc}h_{11}&h_{12}&h_{13}end{array}right] left[begin{array}{ccccc} B_{11}&B_{12}&B_{13}\ B_{21}&B_{22}&B_{23}\ B_{31}&B_{32}&B_{33} end{array}right] left[begin{array}{ccccc}h_{21}\h_{22}\h_{23}end{array}right]\ h_1^TBh_2= left[begin{array}{ccccc} h_{11}h_{21}&h_{12}h_{21}+h_{11}h_{22}&h_{12}h_{22}&h_{13}h_{21}+h_{11}h_{23}&h_{13}h_{23}+h_{12}h_{23}&h_{13}h_{23} end{array}right] left[begin{array}{ccccc} B_{11}\B_{21}\B_{22}\B_{31}\B_{32}\B_{33} end{array}right] h1TBh2=[h11h12h13]⎣⎡B11B21B31B12B22B32B13B23B33⎦⎤⎣⎡h21h22h23⎦⎤h1TBh2=[h11h21h12h21+h11h22h12h22h13h21+h11h23h13h23+h12h23h13h23]⎣⎢⎢⎢⎢⎢⎢⎡B11B21B22B31B32B33⎦⎥⎥⎥⎥⎥⎥⎤
将上述公式简记为：
h 1 T B h 2 = v 12 T b v 12 T = [ h 11 h 21 h 12 h 21 + h 11 h 22 h 12 h 22 h 13 h 21 + h 11 h 23 h 13 h 23 + h 12 h 23 h 13 h 23 ] h_1^TBh_2=v_{12}^Tb\ v_{12}^T= left[begin{array}{ccccc} h_{11}h_{21}&h_{12}h_{21}+h_{11}h_{22}&h_{12}h_{22}&h_{13}h_{21}+h_{11}h_{23}&h_{13}h_{23}+h_{12}h_{23}&h_{13}h_{23} end{array}right] h1TBh2=v12Tbv12T=[h11h21h12h21+h11h22h12h22h13h21+h11h23h13h23+h12h23h13h23]
同理，可计算得到：
h 1 T B h 1 − h 1 T B h 1 = ( v 11 T − v 22 T ) b = 0 h_1^TBh_1-h_1^TBh_1=(v_{11}^T-v_{22}^T)b=0 h1TBh1−h1TBh1=(v11T−v22T)b=0
我们可以将上述二式列为矩阵的形式，我们姑且称该矩阵中的前项为二矩阵，该二矩阵是2×6维矩阵：
[ v 12 T ( v 11 − v 22 ) T ] b = 0 left[begin{array}{ccccc} v_{12}^T\ (v_{11}-v_{22})^T end{array}right]b=0 [v12T(v11−v22)T]b=0
我们从推导中可以知道，上述的齐次矩阵方程所包含的未知量，就是相机的内参数和外参数。相对于一张图片上的不同点，相机的内参数和外参数是不发生变化的，因此，一张图片上我们只能够得到一组上述的齐次矩阵方程。但是这远远不足以让我们求出所有的未知量，于是我们将该公式推广到n张图中，以期获得n组上述的齐次矩阵方程。
现在我们将其推广到n张图中，于是我们得到一个如下的矩阵：
V b = 0 Vb=0 Vb=0式中的 V V V是由每张图中的二矩阵组成的，因此 V V V的维数是2n×6维的矩阵。
在上述的公式中，我们假设 V V V是已知的，我们的未知量是包含有5个自由度的向量 b b b，而一张图可提供两个方程。因此我们需要针对n的个数进行一个讨论。
当n≥3时，我们有5个自由度，等于或多于6个方程，可能出现超定方程组问题。
当n=2时，我们有5个自由度，4个方程，求不出解，因此我们额外施加一个约束： γ = 0 γ=0 γ=0。
在第一篇的坐标转换中我们讲到 ”由于制作工艺等的误差，可能会使像素坐标系的x轴和y轴之间的夹角θ并非一定是90°“，这里我们施加的约束为夹角 θ = 90 ° θ=90° θ=90°，则 c o s θ = 0 ， γ = 0 cosθ=0，γ=0 cosθ=0，γ=0。因此该约束被称为：不倾斜约束。于是我们就得到 B 21 = 0 B_{21}=0 B21=0，未知量个数减1，可求解。
当n=1时，我们无法求解所有的未知量，当知道相机内参中的3个参数时，可求解剩余的两个。
关于超定方程求解的问题，我们在第三篇文章中已经提起过，我们将任务转化为求解一个 b b b使得 ∣ ∣ V b ∣ ∣ ||Vb|| ∣∣Vb∣∣的值最小且 ∣ ∣ b ∣ ∣ = 1 ||b||=1 ∣∣b∣∣=1。求解方法就是对 V V V进行奇异值分解， b b b的解就是 V V V的最小奇异值所对应的右奇异向量。
于是乎，我们求解得到了 b b b，我们可以根据 b b b的值来获得矩阵 B B B。在得知矩阵 B B B后，我们可以根据其和内参数之间的对应关系获得内参数的值。
B = A − T A − 1 = [ 1 α 2 − γ α 2 β γ v 0 − β u 0 α 2 β − γ α 2 β γ 2 α 2 β 2 + 1 β 2 − γ ( γ v 0 − β u 0 ) α 2 β 2 − v 0 β 2 γ v 0 − β u 0 α 2 β − γ ( γ v 0 − β u 0 ) α 2 β 2 − v 0 β 2 ( γ v 0 − β u 0 ) 2 ( α β ) 2 + v 0 2 β 2 + 1 ] = [ B 11 B 12 B 13 B 21 B 22 B 23 B 31 B 32 B 33 ] B=A^{-T}A^{-1}= left[begin{array}{ccccc} frac{1}{α^2}&-frac{γ}{α^2β}&frac{γv_0-βu_0}{α^2β}\ -frac{γ}{α^2β}&frac{γ^2}{α^2β^2}+frac{1}{β^2}&-frac{γ(γv_0-βu_0)}{α^2β^2}-frac{v_0}{β^2}\ frac{γv_0-βu_0}{α^2β}&-frac{γ(γv_0-βu_0)}{α^2β^2}-frac{v_0}{β^2}&frac{(γv_0-βu_0)^2}{(αβ)^2}+frac{v_0^2}{β^2}+1 end{array}right]= left[begin{array}{ccccc} B_{11}&B_{12}&B_{13}\ B_{21}&B_{22}&B_{23}\ B_{31}&B_{32}&B_{33} end{array}right] B=A−TA−1=⎣⎢⎡α21−α2βγα2βγv0−βu0−α2βγα2β2γ2+β21−α2β2γ(γv0−βu0)−β2v0α2βγv0−βu0−α2β2γ(γv0−βu0)−β2v0(αβ)2(γv0−βu0)2+β2v02+1⎦⎥⎤=⎣⎡B11B21B31B12B22B32B13B23B33⎦⎤
由上述的对应关系可以求得：
α = 1 B 11 β = B 11 B 11 B 22 − B 12 2 γ = − B 12 β B 11 v 0 = B 12 B 13 − B 11 B 23 B 11 B 22 − B 12 2 u 0 = B 12 2 B 13 − B 11 B 12 B 23 B 11 2 B 22 − B 11 B 12 2 − B 13 B 11 α=sqrt{frac{1}{B_{11}}}\ β=sqrt{frac{B_{11}}{B_{11}B_{22}-B_{12}^2}}\ γ=-frac{B_{12}β}{B_{11}}\ v_0=frac{B_{12}B_{13}-B_{11}B_{23}}{B_{11}B_{22}-B_{12}^2}\ u_0=frac{B_{12}^2B_{13}-B_{11}B_{12}B_{23}}{B_{11}^2B_{22}-B_{11}B_{12}^2}-frac{B_{13}}{B_{11}} α=B111 β=B11B22−B122B11 γ=−B11B12βv0=B11B22−B122B12B13−B11B23u0=B112B22−B11B122B122B13−B11B12B23−B11B13
如下是论文原文中的结果：
但需要注意的是，论文原文中的 u 0 = γ v 0 / α − B 13 α 2 / λ u_0=γv_0/α-B_{13}α^2/λ u0=γv0/α−B13α2/λ有误，正确的结果应当是 u 0 = γ v 0 / β − B 13 α 2 / λ u_0=γv_0/β-B_{13}α^2/λ u0=γv0/β−B13α2/λ
本文在这里并没有考虑 λ λ λ的取值问题，默认其值为1，但这里也着实没办法从矩阵 B B B中推断出 λ λ λ的取值。
如果我们得到了向量 b b b，我们就可以得到相机的内参数，进而得到矩阵 A A A，我们就可以求得 A − 1 A^{-1} A−1
通过在求导过程中得到的如下公式，我们可以求得旋转和平移矩阵：
r 1 T ⋅ r 2 = ( A − 1 h 1 ) T A − 1 h 2 = h 1 T A − T A − 1 h 2 = 0 r 1 T . r 1 = ∣ ∣ r 1 ∣ ∣ = h 1 T A − T A − 1 h 1 = 1 r 2 T . r 2 = ∣ ∣ r 2 ∣ ∣ = h 2 T A − T A − 1 h 2 = 1 r_1^T·r_2=(A^{-1}h_1)^TA^{-1}h_2=h_1^TA^{-T}A^{-1}h_2=0\ r_1^T.r_1=||r_1||=h_1^TA^{-T}A^{-1}h_1=1\ r_2^T.r_2=||r_2||=h_2^TA^{-T}A^{-1}h_2=1 r1T⋅r2=(A−1h1)TA−1h2=h1TA−TA−1h2=0r1T.r1=∣∣r1∣∣=h1TA−TA−1h1=1r2T.r2=∣∣r2∣∣=h2TA−TA−1h2=1
我们可以得到：
r 1 = A − 1 h 1 ， r 2 = A − 1 h 2 ， r 3 = r 1 × r 2 ， t = A − 1 h 3 r_1=A^{-1}h_1，r_2=A^{-1}h_2，r_3=r_1×r_2，t=A^{-1}h_3 r1=A−1h1，r2=A−1h2，r3=r1×r2，t=A−1h3
到目前为止，我们已经求得了相机的内参矩阵和旋转矩阵与平移矩阵。但是，我们目前所得到的结果准确来说其实是基于两个假设而得到：
第一个假设其实很好理解，我们来讨论为什么假设 Z = 0 Z=0 Z=0对标定结果没有影响。我们在3.相机标定原理(不考虑径向畸变)中已经对相机标定原理做了一个简单的介绍。我们在通用相机标定原理里特别指出了，选取的6个点不应在同一个平面上，不然会出现退化解的情况。而张氏标定法假设 Z = 0 Z=0 Z=0其实就是做了如下的一件事，假设我选取的点就是在一个平面上。那么张氏标定法为什么没有出现退化解的情况？（在写上一篇文章的时候，我以为我可以避开这个问题o(╥﹏╥)o天道好轮回啊，如下的介绍请配合3.相机标定原理(不考虑径向畸变)食用，当然其实你也可以跳过去）
我们来介绍一下退化解的概念：
百度百科：退化解是在进行线性规划问题时，其基本可行解出现一个或多个基变量等于零时，或者按最小比值来确定换出基的变量时，存在两个以上相同最小比值的线性规划问题，出现的原因是模型中存在多余的约束，使多个基本可行解对应同一顶点。
(⊙o⊙)…看懂了吗？我是不太懂的。
但是我在知乎上找到了一个比较通俗的解释，我们大概可以这么理解：
退化类似“降维”，在 P m = 0 Pm=0 Pm=0的求解过程中，由于存在过多的约束，我们看似要求解 m m m中的11个自由度，其实我们真正需要求解的自由度可能只有10或者8个甚至更少。
那么我们来看，如果将所有的点都选在一个平面上，其带来的一个约束是什么，要么这些点的坐标 Z = 0 Z=0 Z=0，要么这些点的坐标 X = 0 X=0 X=0，要么这些点的坐标 Y = 0 Y=0 Y=0。如若如此，那么我们就到了张氏标定法，因此我们的 m m m中的自由度其中的一个或多个其实已经确定，例如我们求解的单应性矩阵 H H H的自由度只有8（下文会讲到）。
而张氏标定法，按照我的个人理解，其实就是对退化解的一个应对策略，其假设 Z = 0 Z=0 Z=0，也可响应的假设 X = 0 X=0 X=0或 Y = 0 Y=0 Y=0都是一样的。
我们重点来讨论第二个假设"假设 V V V是已知的"。
我们在4.2.5节中讨论，为了求解方程组 V b = 0 Vb=0 Vb=0,假设 V V V是已知的。且我们在求解旋转矩阵和平移矩阵的时候认为 h 1 , h 2 , h 3 h_1,h_2,h_3 h1,h2,h3是已知的。然而这些值真的是已知的吗？我们来看 V V V是由什么组成的。
v 12 T = [ h 11 h 21 h 12 h 21 + h 11 h 22 h 12 h 22 h 13 h 21 + h 11 h 23 h 13 h 23 + h 12 h 23 h 13 h 23 ] v_{12}^T= left[begin{array}{ccccc} h_{11}h_{21}&h_{12}h_{21}+h_{11}h_{22}&h_{12}h_{22}&h_{13}h_{21}+h_{11}h_{23}&h_{13}h_{23}+h_{12}h_{23}&h_{13}h_{23} end{array}right] v12T=[h11h21h12h21+h11h22h12h22h13h21+h11h23h13h23+h12h23h13h23]
V V V是由众多的形如 v 12 T v_{12}^T v12T的式子组成的，更深一层可以认为 V V V是由众多的 h 1 , h 2 , h 3 h_1,h_2,h_3 h1,h2,h3组成的，再深一层可以认为 V V V是由相机内参数和多种外参数(每一张图片都会对应一个旋转和平移矩阵)组成的。于是就有了一个疑惑，这不就形成了“我证明我自己”的问题了？很显然，这样是不能够解决问题的。那么在这个过程中，肯定有某种办法可以是可以求解得到单应性矩阵 H H H的。
关于求解单应性矩阵 H H H原理的讲解，个人觉得相机标定之张正友标定法数学原理详解（含python源码）这篇文章讲的还是比较清楚的。单应性矩阵Homography计算和优化以及SLAM入门之视觉里程计(5)：单应矩阵两篇文章中关于单应性矩阵的计算也给了详尽的解释，不过这两篇文章是通过两张图像之间的对应点来求解单应性矩阵的，原理是相似的。
我们最初获得的单应性矩阵如下，(为了与上文中的 h 1 , h 2 h_1,h_2 h1,h2对应，请读者留意矩阵 H H H各行各列的编码变化)：
H = [ α γ u 0 0 β v 0 0 0 1 ] [ r 1 r 2 t ] = A [ r 1 r 2 t ] = [ h 1 h 2 h 3 ] = [ h 11 h 21 h 31 h 12 h 22 h 32 h 13 h 23 h 33 ] s m ~ = H M ~ H= left[ begin{array}{ccccc} α&γ&u_0\ 0&β&v_0\ 0&0&1 end{array} right] left[ begin{array}{ccccc} r_1&r_2&t end{array} right]=A left[ begin{array}{ccccc} r_1&r_2&t end{array} right] =left[ begin{array}{ccccc} h_1&h_2&h_3 end{array} right]= left[ begin{array}{ccccc} h_{11}&h_{21}&h_{31}\ h_{12}&h_{22}&h_{32}\ h_{13}&h_{23}&h_{33} end{array} right]\ swidetilde{m}=Hwidetilde{M} H=⎣⎡α00γβ0u0v01⎦⎤[r1r2t]=A[r1r2t]=[h1h2h3]=⎣⎡h11h12h13h21h22h23h31h32h33⎦⎤sm =HM
通过上式，我们已知图像上的已知2D点 m ~ widetilde{m} m 和世界坐标系下的3D点 M ~ widetilde{M} M 之间存在一个 H H H的变换，接下来，我们需要思考如下的问题：
能否通过图像上的多个点和世界坐标系下的多个点之间的相对关系求解单应性矩阵？
根据这些对应关系，我们如果已知2D点 m ~ widetilde{m} m 和3D点 M ~ widetilde{M} M ，我们可以列出如下的足够多的方程来求解：
s 1 m 1 ~ = H M 1 ~ s 2 m 2 ~ = H M 2 ~ … … s n m n ~ = H M n ~ s_1widetilde{m_1}=Hwidetilde{M_1}\ s_2widetilde{m_2}=Hwidetilde{M_2}\ ……\ s_nwidetilde{m_n}=Hwidetilde{M_n} s1m1 =HM1 s2m2 =HM2 ……snmn =HMn
2D点 m ~ widetilde{m} m 我们可以在图像上找到，那么3D点 M ~ widetilde{M} M 呢？我们在进行标定时，所打印的图片是一个2维的平面，由于不考虑 Z Z Z轴的影响(我们已经假设其为0了)，那么该2维平面上我们可以设置一个坐标原点，作为世界坐标系的原点，该平面则为x-y平面，这样我们就可以得到该平面上的所有点的坐标了 ，于是我们就知道了 M ~ widetilde{M} M 的坐标。利用上述方程，我们就可以求解得到单应性矩阵 H H H。
我们需要多少个点才能够求解单应性矩阵 H H H呢？
在考虑这个问题之前，我们还需要考虑一个问题， H H H的自由度是多少？确定了自由度，我们就可以确定方程的个数，进而确定所需要点的个数。
H H H是一个3×3的矩阵，理论上来说其自由度应该是9。
我们来做一个简单的推导。
s [ u 1 v 1 1 ] = [ h 11 h 21 h 31 h 12 h 22 h 32 h 13 h 23 h 33 ] [ X 1 Y 1 1 ] sleft[begin{array}{ccccc} u_1\ v_1\ 1 end{array}right]= left[ begin{array}{ccccc} h_{11}&h_{21}&h_{31}\ h_{12}&h_{22}&h_{32}\ h_{13}&h_{23}&h_{33} end{array} right] left[ begin{array}{ccccc} X_1\ Y_1\ 1 end{array} right] s⎣⎡u1v11⎦⎤=⎣⎡h11h12h13h21h22h23h31h32h33⎦⎤⎣⎡X1Y11⎦⎤
将上述式子展开：
s u 1 = h 11 X 1 + h 21 Y 1 + h 31 s v 1 = h 12 X 1 + h 22 Y 1 + h 32 s = h 13 X 1 + h 23 Y 1 + h 33 su_1=h_{11}X_1+h_{21}Y_1+h_{31}\ sv_1=h_{12}X_1+h_{22}Y_1+h_{32}\ s=h_{13}X_1+h_{23}Y_1+h_{33} su1=h11X1+h21Y1+h31sv1=h12X1+h22Y1+h32s=h13X1+h23Y1+h33
我们可以得到：
u 1 = h 11 X 1 + h 21 Y 1 + h 31 h 13 X 1 + h 23 Y 1 + h 33 v 1 = h 12 X 1 + h 22 Y 1 + h 32 h 13 X 1 + h 23 Y 1 + h 33 u_1=frac{h_{11}X_1+h_{21}Y_1+h_{31}}{h_{13}X_1+h_{23}Y_1+h_{33}}\ v_1=frac{h_{12}X_1+h_{22}Y_1+h_{32}}{h_{13}X_1+h_{23}Y_1+h_{33}} u1=h13X1+h23Y1+h33h11X1+h21Y1+h31v1=h13X1+h23Y1+h33h12X1+h22Y1+h32
展开：
( h 13 X 1 + h 23 Y 1 + h 33 ) u 1 − ( h 11 X 1 + h 21 Y 1 + h 31 ) = 0 ( h 13 X 1 + h 23 Y 1 + h 33 ) v 1 − ( h 12 X 1 + h 22 Y 1 + h 32 ) = 0 (h_{13}X_1+h_{23}Y_1+h_{33})u_1-(h_{11}X_1+h_{21}Y_1+h_{31})=0\ (h_{13}X_1+h_{23}Y_1+h_{33})v_1-(h_{12}X_1+h_{22}Y_1+h_{32})=0 (h13X1+h23Y1+h33)u1−(h11X1+h21Y1+h31)=0(h13X1+h23Y1+h33)v1−(h12X1+h22Y1+h32)=0
已知 h i j h_{ij} hij是我们待求的未知量，将上式转换成矩阵形式：
[ − X 1 − Y 1 − 1 0 0 0 X 1 u 1 Y 1 u 1 u 1 0 0 0 − X 1 − Y 1 − 1 X 1 v 1 Y 1 v 1 v 1 ] [ h 11 h 21 h 31 h 12 h 22 h 32 h 13 h 23 h 33 ] = 0 left[begin{array}{ccccc} -X_1&-Y_1&-1&0&0&0&X_1u_1&Y_1u_1&u_1\ 0&0&0&-X_1&-Y_1&-1&X_1v_1&Y_1v_1&v_1 end{array}right] left[begin{array}{ccccc} h_{11}\h_{21}\h_{31}\h_{12}\h_{22}\h_{32}\h_{13}\h_{23}\h_{33} end{array}right]=0 [−X10−Y10−100−X10−Y10−1X1u1X1v1Y1u1Y1v1u1v1]⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡h11h21h31h12h22h32h13h23h33⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=0
我们不难发现这是一个齐次方程组。举个例子，已知一个齐次方程 x + y = 0 x+y=0 x+y=0 我们需要知道几个约束条件来求解得到 x , y x,y x,y 呢？很显然我们只需要再知道一个约束方程就可以求解到 x , y x,y x,y 。这是因为该齐次方程本身就已经是一个约束条件了。
同理，我们也可以发现，我们只需要知道 H H H矩阵中的8个未知量，带入到上述的齐次方程组中，就可以求得第9个未知量。
因此，我们可以得出结论： H H H的自由度是8。也就是说我们至少需要8个方程才可以求解得到 H H H。
而且，从上述推导中，我们可以发现，每找到图像上的一个点，我们就可以得到两个方程。
于是，我们可以得出结果，我们至少需要在图像上找到四个点，这样我们就可以求解 H H H。得到了 H H H之后，我们就可以 V V V，之后就可以求解向量 b b b，继而得到矩阵 B B B，再求解得到内参和外参。
求解单应性矩阵 H H H
在实际的标定过程中，使用的点数通常多于四个点了，借此增加标定结果的鲁棒性。
于是我们得到这样一个方程组 N h = 0 Nh=0 Nh=0：
[ − X 1 − Y 1 − 1 0 0 0 X 1 u 1 Y 1 u 1 u 1 0 0 0 − X 1 − Y 1 − 1 X 1 v 1 Y 1 v 1 v 1 − X 2 − Y 2 − 1 0 0 0 X 2 u 2 Y 2 u 2 u 2 0 0 0 − X 2 − Y 2 − 1 X 2 v 2 Y 2 v 2 v 2 … … … … − X n − Y n − 1 0 0 0 X n u n Y n u n u n 0 0 0 − X n − Y n − 1 X n v n Y n v n v n ] [ h 11 h 21 h 31 h 12 h 22 h 32 h 13 h 23 h 33 ] = 0 left[begin{array}{ccccc} -X_1&-Y_1&-1&0&0&0&X_1u_1&Y_1u_1&u_1\ 0&0&0&-X_1&-Y_1&-1&X_1v_1&Y_1v_1&v_1\ -X_2&-Y_2&-1&0&0&0&X_2u_2&Y_2u_2&u_2\ 0&0&0&-X_2&-Y_2&-1&X_2v_2&Y_2v_2&v_2\ ……&……\ -X_n&-Y_n&-1&0&0&0&X_nu_n&Y_nu_n&u_n\ 0&0&0&-X_n&-Y_n&-1&X_nv_n&Y_nv_n&v_n\ end{array}right] left[begin{array}{ccccc} h_{11}\h_{21}\h_{31}\h_{12}\h_{22}\h_{32}\h_{13}\h_{23}\h_{33} end{array}right]=0 ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡−X10−X20……−Xn0−Y10−Y20……−Yn0−10−10−100−X10−X20−Xn0−Y10−Y20−Yn0−10−10−1X1u1X1v1X2u2X2v2XnunXnvnY1u1Y1v1Y2u2Y2v2YnunYnvnu1v1u2v2unvn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡h11h21h31h12h22h32h13h23h33⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=0
这是就出现了上文4.2.6节曾经出现过的一个问题：超定方程组的问题。将任务转化为求解一个 h h h使得 ∣ ∣ N h ∣ ∣ ||Nh|| ∣∣Nh∣∣的值最小且 ∣ ∣ h ∣ ∣ = 1 ||h||=1 ∣∣h∣∣=1。求解方法就是对 N N N进行奇异值分解， h h h的解就是 N N N的最小奇异值所对应的右奇异向量。
自此我们已经解出了单应性矩阵 H H H。
由于求解精度的问题，我们得到旋转矩阵和内参矩阵是存在一定误差的，我们认为，误差自然是越小越好。
于是我们定义一个新的约束条件，称为最大似然估计(说人话就是，我们通过多次计算，得到多个内参和外参，然后我们用这些内参和外参计算世界坐标下的点在图像上的投影位置，挑选出来计算结果和真实结果最相近的那一组内参和外参作为标定的结果)：
E = ∑ i = 1 n ∑ j = 1 m ∣ ∣ m i j − m ^ ( A , R i , t i , M j ) ∣ ∣ 2 E=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}||m_{ij}-widehat{m}(A,R_i,t_i,M_j)||^2 E=i=1∑nj=1∑m∣∣mij−m (A,Ri,ti,Mj)∣∣2
该方程中 m i j m_{ij} mij表示第 i i i 张图片中的第 j j j 个点，是在图像中真实得到的。 m ^ widehat{m} m 则表示我们在已知世界坐标、内参矩阵、外参矩阵的情况下，计算得到的第 i i i 张图像中的第 j j j 个点。 A A A表示相机的内参数， R i , t i R_i,t_i Ri,ti表示第 i i i 个世界坐标系(标定板的图片)相对于相机的外参矩阵， M j M_j Mj则表示第 i i i 个世界坐标系中的第 j j j 个点。
该方程的实际意义就是，计算每张图片中每个点的真实值和计算值之间的误差的平方和。
最大似然估计就是说，要让方程中 E E E 的值最小，即误差平方和最小，那么对应的内参矩阵和外参矩阵就是误差最小，精度最高的。
而这种最小化的实现通常采用的是列文伯格—马夸尔特优化算法或者是牛顿优化算法。(这里不做介绍)
我们来进行一个对照，将通用的标定方法和张氏标定法进行一个对比：
通用标定方法：
张氏标定法的标定过程：
整体来看，通用标定方法和张氏标定方法是大体上是一致的，二者之间的区别就在于，前者从三个平面上取点，后者从一个平面上取点。这就更验证了我的猜测：张氏标定是对通用标定法将所有点取在一个平面上时出现的退化问题的一个应对策略，这反而简化了标定的过程。
本文发布于:2024-01-28 10:37:39，感谢您对本站的认可！
本文链接：https://www.4u4v.net/it/17064094646825.html
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系，我们将在24小时内删除。
| 留言与评论（共有 0 条评论） |
