当前位置:   article > 正文

二维平面成像模型(Models for transformations)(一)_transformation model image

transformation model image


  本节将介绍利用针孔相机对场景中2D平面成像的相关知识。在这种情况下,场景中2D平面上的点和像平面上的点满足一一对应的关系。

1 2D转换模型

转换模型:transformation model

1.1 Euclidean transformation model

  Euclidean transformation model描述场景中与像平面正面平行、距光心距离为 D D D的平面上的点与像平面对应点的数学关系。
  首先假设世界坐标系的 u u u- v v v平面与该平面重合,也即 w w w轴与之垂直,则该平面上的点可表示为 w = [ u , v , 0 ] T w=[u,v,0]^T w=[u,v,0]T。回想小孔成像的数学模型:
λ x ~ = Λ [ Ω , τ ] w ~ (1) \lambda \widetilde{x}=\Lambda[\Omega,\tau]\widetilde{w}\tag{1} λx =Λ[Ω,τ]w (1)
  应用此时的场景,便有:
λ [ x y 1 ] = [ ϕ x γ δ x 0 ϕ y δ y 0 0 1 ] [ ω 11 ω 12 0 τ x ω 21 ω 22 0 τ y 0 0 1 D ] [ u v 0 1 ] = [ ϕ x γ δ x 0 ϕ y δ y 0 0 1 ] [ ω 11 ω 12 τ x ω 21 ω 22 τ y 0 0 D ] [ u v 1 ] = [ ϕ x γ δ x 0 ϕ y δ y 0 0 D ] [ ω 11 ω 12 τ x ω 21 ω 22 τ y 0 0 1 ] [ u v 1 ] (2) \lambda \left[

xy1
\right]= \left[
ϕxγδx0ϕyδy001
\right] \left[
ω11ω120τxω21ω220τy001D
\right]\left[
uv01
\right]\\ =\left[
ϕxγδx0ϕyδy001
\right] \left[
ω11ω12τxω21ω22τy00D
\right]\left[
uv1
\right]\\ =\left[
ϕxγδx0ϕyδy00D
\right] \left[
ω11ω12τxω21ω22τy001
\right]\left[
uv1
\right]\tag{2} λxy1=ϕx00γϕy0δxδy1ω11ω210ω12ω220001τxτyDuv01=ϕx00γϕy0δxδy1ω11ω210ω12ω220τxτyDuv1=ϕx00γϕy0δxδyDω11ω210ω12ω220τxτy1uv1(2)
  上式中由倒数第二行到最后一行的变化为将距离 D D D移到了 Λ \Lambda Λ中,注意这样处理有利于后续的处理而且并不改变计算结果。之后我们再等式两边同时左乘 Λ − 1 \Lambda^{-1} Λ1可得到:
λ [ x ′ y ′ 1 ] = [ ω 11 ω 12 τ x ω 21 ω 22 τ y 0 0 1 ] [ u v 1 ] (3) \lambda \left[
xy1
\right]= \left[
ω11ω12τxω21ω22τy001
\right]\left[
uv1
\right]\tag{3}
λxy1=ω11ω210ω12ω220τxτy1uv1(3)

  上式便是Euclidean transformation. 还可表示为:
[ x ′ y ′ ] = [ ω 11 ω 12 ω 21 ω 22 ] [ u v ] + [ τ x τ y ] (4) \left[
xy
\right]= \left[
ω11ω12ω21ω22
\right]\left[
uv
\right]+\left[
τxτy
\right]\tag{4}
[xy]=[ω11ω21ω12ω22][uv]+[τxτy](4)

  或者简记为:
x ′ = e u c [ w , Ω , τ ] (5) x'=euc[w,\Omega,\tau]\tag{5} x=euc[w,Ω,τ](5)
  注意,Euclidean transformation只能刚性旋转和平移,虽然它表面上有6个未知参数,实际上只有3个参数(旋转角度 θ \theta θ、平移参数 τ x \tau_x τx τ y \tau_y τy), Ω \Omega Ω可表示为:
[ ω 11 ω 12 ω 21 ω 22 ] = [ cos ⁡ θ sin ⁡ θ − sin ⁡ θ cos ⁡ θ ] (6) \left[
ω11ω12ω21ω22
\right]=\left[
cosθsinθsinθcosθ
\right]\tag{6}
[ω11ω21ω12ω22]=[cosθsinθsinθcosθ](6)

1.2 Similarity transformation model

  Similarity transformation model描述场景中与像平面正面平行、到光心距离未知的平面上的点与像平面对应点的数学关系。此时平面上坐标点 w = [ u , v ] T w=[u,v]^T w=[u,v]T和对应投影点 x = [ x , y ] T x=[x,y]^T x=[x,y]T的关系为:
λ [ x ′ y ′ 1 ] = [ ω 11 ω 12 τ x ω 21 ω 22 τ y 0 0 D ] [ u v 1 ] (7) \lambda \left[

xy1
\right]= \left[
ω11ω12τxω21ω22τy00D
\right]\left[
uv1
\right]\tag{7} λxy1=ω11ω210ω12ω220τxτyDuv1(7)
  等式两边同时乘以 ρ = 1 / D \rho=1/D ρ=1/D得:
ρ λ [ x ′ y ′ 1 ] = [ ρ ω 11 ρ ω 12 ρ τ x ρ ω 21 ρ ω 22 ρ τ y 0 0 1 ] [ u v 1 ] (8) \rho\lambda \left[
xy1
\right]= \left[
ρω11ρω12ρτxρω21ρω22ρτy001
\right]\left[
uv1
\right]\tag{8}
ρλxy1=ρω11ρω210ρω12ρω220ρτxρτy1uv1(8)

  再把 λ \lambda λ τ x \tau_x τx τ y \tau_y τy前面的 ρ \rho ρ并入可得:
λ [ x ′ y ′ 1 ] = [ ρ ω 11 ρ ω 12 τ x ρ ω 21 ρ ω 22 τ y 0 0 1 ] [ u v 1 ] (9) \lambda \left[
xy1
\right]= \left[
ρω11ρω12τxρω21ρω22τy001
\right]\left[
uv1
\right]\tag{9}
λxy1=ρω11ρω210ρω12ρω220τxτy1uv1(9)

  上式便是Similarity transformation. 还可表示为:
[ x ′ y ′ ] = [ ρ ω 11 ρ ω 12 ρ ω 21 ρ ω 22 ] [ u v ] + [ τ x τ y ] (10) \left[
xy
\right]= \left[
ρω11ρω12ρω21ρω22
\right]\left[
uv
\right]+\left[
τxτy
\right]\tag{10}
[xy]=[ρω11ρω21ρω12ρω22][uv]+[τxτy](10)

  或者简记为:
x ′ = s i m [ w , Ω , τ , ρ ] (11) x'=sim[w,\Omega,\tau,\rho]\tag{11} x=sim[w,Ω,τ,ρ](11)
  注意,Similarity transformation只能刚性旋转、平移和缩放,它比Euclidean transformation多了一个参数 ρ \rho ρ

1.3 Affine transformation model

  Affine transformation model被希望(后面将会看到这不可行)用来描述场景中任意位置的平面上的点与像平面对应点的数学关系:
λ [ x ′ y ′ 1 ] = [ ϕ 11 ϕ 12 τ x ϕ 21 ϕ 22 τ y 0 0 1 ] [ u v 1 ] (12) \lambda \left[

xy1
\right]= \left[
ϕ11ϕ12τxϕ21ϕ22τy001
\right]\left[
uv1
\right]\tag{12} λxy1=ϕ11ϕ210ϕ12ϕ220τxτy1uv1(12)
  上式便是Affine transformation. 还可表示为:
[ x ′ y ′ ] = [ ϕ 11 ϕ 12 ϕ 21 ϕ 22 ] [ u v ] + [ τ x τ y ] (13) \left[
xy
\right]= \left[
ϕ11ϕ12ϕ21ϕ22
\right]\left[
uv
\right]+\left[
τxτy
\right]\tag{13}
[xy]=[ϕ11ϕ21ϕ12ϕ22][uv]+[τxτy](13)

  或者简记为:
x ′ = a f f [ w , Φ , τ ] (14) x'=aff[w,\Phi,\tau]\tag{14} x=aff[w,Φ,τ](14)
  注意,由于相机内部参数矩阵 Λ \Lambda Λ也是一个仿射矩阵(最后一行前两个元素为0),而两个仿射矩阵的乘积仍然是一个仿射矩阵,所以可以用仿射变换来直接描述像平面上点和场景中平面上对应点的数学关系。
  注意,仿射变换适用的条件为:场景中平面的深度变换比其到相机的平均距离小。这种情况一般表现为观察角度不是很倾斜、相机距离平面较远而且视角相对较小。除此以外的情况不适合用仿射变换来描述,这是因为仿射变换属于线性变换,线性变换的特点之一是平行的线变换后仍保持平行;而实际生活中,例如我们观察较远处的火车轨道时,它往往表现为将要交汇于某一点。

1.4 Projective transformation model

  透视变换可以描述从任意角度、任意距离观察任意大小物体时的坐标变化关系。此时,场景中平面上的点 w = [ u , v , 0 ] T w=[u,v,0]^T w=[u,v,0]T和与之对应的像平面上的投影点的关系为:
λ [ x y 1 ] = [ ϕ x γ δ x 0 ϕ y δ y 0 0 1 ] [ ω 11 ω 12 ω 13 τ x ω 21 ω 22 ω 23 τ y ω 31 ω 32 ω 33 τ z ] [ u v 0 1 ] = [ ϕ x γ δ x 0 ϕ y δ y 0 0 1 ] [ ω 11 ω 12 τ x ω 21 ω 22 τ y ω 31 ω 32 τ z ] [ u v 1 ] (15) \lambda \left[

xy1
\right]= \left[
ϕxγδx0ϕyδy001
\right] \left[
ω11ω12ω13τxω21ω22ω23τyω31ω32ω33τz
\right]\left[
uv01
\right]\\ =\left[
ϕxγδx0ϕyδy001
\right] \left[
ω11ω12τxω21ω22τyω31ω32τz
\right]\left[
uv1
\right]\tag{15} λxy1=ϕx00γϕy0δxδy1ω11ω21ω31ω12ω22ω32ω13ω23ω33τxτyτzuv01=ϕx00γϕy0δxδy1ω11ω21ω31ω12ω22ω32τxτyτzuv1(15)
  将两个 3 × 3 3\times 3 3×3的矩阵相乘可得:
λ [ x y 1 ] = [ ϕ 11 ϕ 12 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ] [ u v 1 ] (16) \lambda \left[
xy1
\right]= \left[
ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]\left[
uv1
\right]\tag{16}
λxy1=ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32ϕ13ϕ23ϕ33uv1(16)

  上式便是透视变换(projective transformation),还被称为单应变换(homography). 在笛卡尔坐标系下表示为:
x = ϕ 11 u + ϕ 12 v + ϕ 13 ϕ 31 u + ϕ 32 v + ϕ 33 y = ϕ 21 u + ϕ 22 v + ϕ 23 ϕ 31 u + ϕ 32 v + ϕ 33 (17) x=\frac{\phi_{11}u+\phi_{12}v+\phi_{13}}{\phi_{31}u+\phi_{32}v+\phi_{33}}\\ y=\frac{\phi_{21}u+\phi_{22}v+\phi_{23}}{\phi_{31}u+\phi_{32}v+\phi_{33}}\tag{17} x=ϕ31u+ϕ32v+ϕ33ϕ11u+ϕ12v+ϕ13y=ϕ31u+ϕ32v+ϕ33ϕ21u+ϕ22v+ϕ23(17)
  或者简记为:
x = h o m [ w , Φ ] (18) x=hom[w,\Phi]\tag{18} x=hom[w,Φ](18)
  注意,透视变换可以将平面上一个四边形变化为任意形状的四边形。而且虽然矩阵 Φ \Phi Φ中有9个未知参数,实际上它只有8个自由度。

2 转换模型参数的学习

  实际应用中,场景平面成像过程中容易受到各种噪声的影响,因此成像点坐标带有不确定性,因此我们考虑用概率模型进一步描述。对于透视投影(Homography),概率模型为:
P r ( x ∣ w ) = N o r m x [ h o m [ w , Φ ] , δ 2 I ] (19) Pr(x|w)=Norm_x[hom[w,\Phi],\delta^2I]\tag{19} Pr(xw)=Normx[hom[w,Φ],δ2I](19)

2.1 问题描述

  已知场景中2D平面上 I I I个不同的2D坐标点 w i = [ u i , v i ] T w_i=[u_i,v_i]^T wi=[ui,vi]T以及与之对应的图像平面2D投影点 x i = [ x i , y i ] T x_i=[x_i,y_i]^T xi=[xi,yi]T,对于给定的转换模型 t r a n s [ w i , θ ] trans[w_i,\theta] trans[wi,θ],求解模型的参数 θ \theta θ。应用极大似然法可得:
θ ^ = argmax ⁡ θ   [ ∏ i = 1 I N o r m x [ t r a n s [ w , θ ] , δ 2 I ] ] = argmax ⁡ θ   [ ∑ i = 1 I log ⁡ [ N o r m x [ t r a n s [ w , θ ] , δ 2 I ] ] ] = argmin ⁡ θ   [ ∑ i = 1 I ( x i − t r a n s [ w i , θ ] ) T ( x i − t r a n s [ w i , θ ] ) ] (20) \widehat{\theta}={\underset {\theta}{\operatorname {argmax} }}\,\Bigg[\prod_{i=1}^I Norm_x[trans[w,\theta],\delta^2I]\Bigg]\\ ={\underset {\theta}{\operatorname {argmax} }}\,\Bigg[\sum_{i=1}^I\log[Norm_x[trans[w,\theta],\delta^2I]]\Bigg]\\ ={\underset {\theta}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-trans[w_i,\theta])^T(x_i-trans[w_i,\theta]) \Bigg]\tag{20} θ =θargmax[i=1INormx[trans[w,θ],δ2I]]=θargmax[i=1Ilog[Normx[trans[w,θ],δ2I]]]=θargmin[i=1I(xitrans[wi,θ])T(xitrans[wi,θ])](20)

2.2 学习Euclidean参数

  Enclidean transformation的参数包括一个 2 × 2 2\times 2 2×2的旋转矩阵 Ω \Omega Ω和一个 2 × 1 2\times 1 2×1的位移向量 τ = [ τ x , τ y ] T \tau=[\tau_x,\tau_y]^T τ=[τx,τy]T。所以问题变为:
Ω ^ , τ ^ = argmin ⁡ Ω , τ   [ ∑ i = 1 I ( x i − e u c [ w i , Ω , τ ] ) T ( x i − e u c [ w i , Ω , τ ] ) ] = argmin ⁡ Ω , τ   [ ∑ i = 1 I ( x i − Ω w i − τ ) T ( x i − Ω w i − τ ) ] (21) \widehat{\Omega},\widehat{\tau} ={\underset {\Omega,\tau}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-euc[w_i,\Omega,\tau])^T(x_i-euc[w_i,\Omega,\tau]) \Bigg]\\ ={\underset {\Omega,\tau}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-\Omega w_i-\tau)^T(x_i-\Omega w_i-\tau) \Bigg] \tag{21} Ω ,τ =Ω,τargmin[i=1I(xieuc[wi,Ω,τ])T(xieuc[wi,Ω,τ])]=Ω,τargmin[i=1I(xiΩwiτ)T(xiΩwiτ)](21)
  将目标函数对 τ \tau τ求导并等于零可以求解得到:
τ ^ = ∑ i = 1 I x i − Ω w i I = μ x − Ω μ w (22) \widehat{\tau} =\frac{\sum_{i=1}^I x_i-\Omega w_i}{I}=\mu_x-\Omega\mu_w \tag{22} τ =Ii=1IxiΩwi=μxΩμw(22)
  代入式(21)可得:
Ω ^ = argmin ⁡ Ω   [ ∑ i = 1 I ( ( x i − μ x ) − Ω ( w i − μ w ) ) T ( ( x i − μ x ) − Ω ( w i − μ w ) ) ] (23) \widehat{\Omega} ={\underset {\Omega}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I ((x_i-\mu_x)-\Omega (w_i-\mu_w))^T((x_i-\mu_x)-\Omega (w_i-\mu_w)) \Bigg] \tag{23} Ω =Ωargmin[i=1I((xiμx)Ω(wiμw))T((xiμx)Ω(wiμw))](23)
  定义矩阵 B = [ x 1 − μ x , x 2 − μ x , ⋯   , x I − μ x ] B=[x_1-\mu_x,x_2-\mu_x,\cdots,x_I-\mu_x] B=[x1μx,x2μx,,xIμx]和矩阵 A = [ w 1 − μ w , w 2 − μ w , ⋯   , w I − μ w ] A=[w_1-\mu_w,w_2-\mu_w,\cdots,w_I-\mu_w] A=[w1μw,w2μw,,wIμw],此时问题转化为如下的Orthogonal Procrustes Problem:
Ω ^ = argmin ⁡ Ω   [ ∣ Ω A − B ∣ F ] (24) \hat \Omega={\underset {\Omega}{\operatorname {argmin} }}\,[|\Omega A-B|_F]\tag{24} Ω^=Ωargmin[ΩABF](24)
  问题的解为:
B A T = U L V T Ω ^ = U V T (25) BA^T=ULV^T\\ \widehat{\Omega}=UV^T\tag{25} BAT=ULVTΩ =UVT(25)

2.3 学习Similarity参数

  Similarity transformation的参数包括一个 2 × 2 2\times 2 2×2的旋转矩阵 Ω \Omega Ω、一个 2 × 1 2\times 1 2×1的位移向量 τ = [ τ x , τ y ] T \tau=[\tau_x,\tau_y]^T τ=[τx,τy]T和一个标量 ρ \rho ρ。所以问题变为:
Ω ^ , τ ^ , ρ ^ = argmin ⁡ Ω , τ , ρ   [ ∑ i = 1 I ( x i − s i m [ w i , Ω , τ , ρ ] ) T ( x i − s i m [ w i , Ω , τ , ρ ] ) ] = argmin ⁡ Ω , τ   [ ∑ i = 1 I ( x i − ρ Ω w i − τ ) T ( x i − ρ Ω w i − τ ) ] (26) \widehat{\Omega},\widehat{\tau},\widehat{\rho} ={\underset {\Omega,\tau,\rho}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-sim[w_i,\Omega,\tau,\rho])^T(x_i-sim[w_i,\Omega,\tau,\rho]) \Bigg]\\ ={\underset {\Omega,\tau}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-\rho\Omega w_i-\tau)^T(x_i-\rho\Omega w_i-\tau) \Bigg] \tag{26} Ω ,τ ,ρ =Ω,τ,ρargmin[i=1I(xisim[wi,Ω,τ,ρ])T(xisim[wi,Ω,τ,ρ])]=Ω,τargmin[i=1I(xiρΩwiτ)T(xiρΩwiτ)](26)
  和2.2一样可以求解得到 Ω ^ \widehat{\Omega} Ω ;
  代入可得 ρ ^ \widehat{\rho} ρ τ ^ \widehat{\tau} τ 为:
ρ ^ = ∑ i = 1 I ( x i − μ x ) T Ω ^ ( w i − μ w ) ∑ i = 1 I ( w i − μ w ) T ( w i − μ w ) τ ^ = ∑ i = 1 I ( x i − ρ ^ Ω ^ w i ) I = μ x − Ω μ w (27) \widehat{\rho}=\frac{\sum_{i=1}^I (x_i-\mu_x)^T\widehat{\Omega}(w_i-\mu_w)}{\sum_{i=1}^I (w_i-\mu_w)^T(w_i-\mu_w)}\\ \widehat{\tau}=\frac{\sum_{i=1}^I (x_i-\widehat{\rho}\widehat{\Omega} w_i)}{I}=\mu_x-\Omega\mu_w \tag{27} ρ =i=1I(wiμw)T(wiμw)i=1I(xiμx)TΩ (wiμw)τ =Ii=1I(xiρ Ω wi)=μxΩμw(27)

2.4 学习Affine参数

  Affine transformation的参数包括一个 2 × 2 2\times 2 2×2的旋转矩阵 Ω \Omega Ω、一个 2 × 1 2\times 1 2×1的位移向量 τ = [ τ x , τ y ] T \tau=[\tau_x,\tau_y]^T τ=[τx,τy]T。所以问题变为:
Φ ^ , τ ^ = argmin ⁡ Φ , τ   [ ∑ i = 1 I ( x i − a f f [ w i , Φ , τ ] ) T ( x i − a f f [ w i , Φ , τ ] ) ] = argmin ⁡ Ω , τ   [ ∑ i = 1 I ( x i − Φ w i − τ ) T ( x i − Φ w i − τ ) ] (28) \widehat{\Phi},\widehat{\tau} ={\underset {\Phi,\tau}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-aff[w_i,\Phi,\tau])^T(x_i-aff[w_i,\Phi,\tau]) \Bigg]\\ ={\underset {\Omega,\tau}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-\Phi w_i-\tau)^T(x_i-\Phi w_i-\tau) \Bigg] \tag{28} Φ ,τ =Φ,τargmin[i=1I(xiaff[wi,Φ,τ])T(xiaff[wi,Φ,τ])]=Ω,τargmin[i=1I(xiΦwiτ)T(xiΦwiτ)](28)
  我们对 Φ w i + τ \Phi w_i+\tau Φwi+τ做如下处理:
Φ w i + τ = [ u i v i 1 0 0 0 0 0 0 u i v i 1 ] [ ϕ 11 ϕ 12 τ x ϕ 21 ϕ 22 τ y ] = A i b (29) \Phi w_i+\tau= \left[

uivi1000000uivi1
\right] \left[
ϕ11ϕ12τxϕ21ϕ22τy
\right]=A_ib\tag{29} Φwi+τ=[ui0vi0100ui0vi01]ϕ11ϕ12τxϕ21ϕ22τy=Aib(29)
  这样,式(28)变为:
b ^ = argmin ⁡ b   [ ∑ i = 1 I ( x i − A i b ) T ( x i − A i b ) ] (30) \widehat{b} ={\underset {b}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-A_ib)^T(x_i-A_ib) \Bigg] \tag{30} b =bargmin[i=1I(xiAib)T(xiAib)](30)
  这样便可以按照Least squares problem求解。

2.5 学习Projective参数

  Projective transformation的参数包括一个 3 × 3 3\times 3 3×3的矩阵 Φ \Phi Φ,因此问题变为:
Φ ^ = argmin ⁡ Φ   [ ∑ i = 1 I ( x i − h o m [ w i , Φ ] ) T ( x i − h o m [ w i , Φ ] ) ] = argmin ⁡ Φ   [ ∑ i = 1 I ( x i − ϕ 11 u i + ϕ 12 v i + ϕ 13 ϕ 31 u i + ϕ 32 v i + ϕ 33 ) 2 + ( y i − ϕ 21 u i + ϕ 22 v i + ϕ 23 ϕ 31 u i + ϕ 32 v i + ϕ 33 ) 2 ] (31) \widehat{\Phi} ={\underset {\Phi}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-hom[w_i,\Phi])^T(x_i-hom[w_i,\Phi]) \Bigg]\\ ={\underset {\Phi}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x_i-\frac{\phi_{11}u_i+\phi_{12}v_i+\phi_{13}}{\phi_{31}u_i+\phi_{32}v_i+\phi_{33}})^2+(y_i-\frac{\phi_{21}u_i+\phi_{22}v_i+\phi_{23}}{\phi_{31}u_i+\phi_{32}v_i+\phi_{33}})^2 \Bigg] \tag{31} Φ =Φargmin[i=1I(xihom[wi,Φ])T(xihom[wi,Φ])]=Φargmin[i=1I(xiϕ31ui+ϕ32vi+ϕ33ϕ11ui+ϕ12vi+ϕ13)2+(yiϕ31ui+ϕ32vi+ϕ33ϕ21ui+ϕ22vi+ϕ23)2](31)
  这个问题没有封闭解,只能应用基于梯度的非线性优化求解方法。像之前一样,我们可以求得一个合适的初始解。在齐次坐标表示下:
λ [ x i y i 1 ] = [ ϕ 11 ϕ 12 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ] [ u i v i 1 ] (32) \lambda \left[

xiyi1
\right]= \left[
ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]\left[
uivi1
\right]\tag{32} λxiyi1=ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32ϕ13ϕ23ϕ33uivi1(32)
  每个齐次坐标表示一条通过光心的直线,因此上式表明等号左边的向量与等号右边的向量相等,也即同方向,因此有:
x ~ × Φ w ~ = 0 (33) \widetilde{x}\times \Phi\widetilde{w}=0\tag{33} x ×Φw =0(33)
  展开后可以得到:
[ y ( ϕ 31 u + ϕ 32 v + ϕ 33 ) − ( ϕ 21 u + ϕ 22 v + ϕ 23 ) ( ϕ 11 u + ϕ 12 v + ϕ 13 ) − x ( ϕ 31 u + ϕ 32 v + ϕ 33 ) x ( ϕ 21 u + ϕ 22 v + ϕ 23 ) − y ( ϕ 11 u + ϕ 12 v + ϕ 13 ) ] = 0 \left[
y(ϕ31u+ϕ32v+ϕ33)(ϕ21u+ϕ22v+ϕ23)(ϕ11u+ϕ12v+ϕ13)x(ϕ31u+ϕ32v+ϕ33)x(ϕ21u+ϕ22v+ϕ23)y(ϕ11u+ϕ12v+ϕ13)
\right]=0
y(ϕ31u+ϕ32v+ϕ33)(ϕ21u+ϕ22v+ϕ23)(ϕ11u+ϕ12v+ϕ13)x(ϕ31u+ϕ32v+ϕ33)x(ϕ21u+ϕ22v+ϕ23)y(ϕ11u+ϕ12v+ϕ13)=0

  将 I I I个点对应的等式全部整合为下式:
[ 0 0 0 − u 1 − v 1 − 1 u 1 y 1 v 1 y 1 y 1 u 1 v 1 1 0 0 0 − u 1 x 1 − v 1 x 1 − x 1 0 0 0 − u 2 − v 2 − 1 u 2 y 2 v 2 y 2 y 2 u 2 v 2 1 0 0 0 − u 2 x 2 − v 2 x 2 − x 2 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 0 0 − u I − v I − 1 u I y I v I y I y I u I v I 1 0 0 0 − u I x I − v I x I − x I ] [ ϕ 11 ϕ 12 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ] = 0 (34) \left[
000u1v11u1y1v1y1y1u1v11000u1x1v1x1x1000u2v21u2y2v2y2y2u2v21000u2x2v2x2x2000uIvI1uIyIvIyIyIuIvI1000uIxIvIxIxI
\right]\left[
ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]=0\tag{34}
0u10u20uI0v10v20vI010101u10u20uI0v10v20vI0101010u1y1u1x1u2y2u2x2uIyIuIxIv1y1v1x1v2y2v2x2vIyIvIxIy1x1y2x2yIxIϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33=0(34)

  此时,问题就变为形如 A b = 0 Ab=0 Ab=0的minimum direction problem (请阅读参考资料),即求解奇异值分解 A = U L V T A=ULV^T A=ULVT并且设 b ^ \hat{b} b^ V V V的最后一列。

3 利用转换模型进行推理

  已知一个转换模型 t r a n s [ w , θ ] trans[w,\theta] trans[w,θ]和像平面上的一个点 x = [ x , y ] T x=[x,y]^T x=[x,y]T,求解场景中平面上的对应点:
w ^ = argmax ⁡ w   [ ∑ i = 1 I log ⁡ [ N o r m x [ t r a n s [ w , θ ] , δ 2 I ] ] ] = argmin ⁡ w   [ ∑ i = 1 I ( x − t r a n s [ w , θ ] ) T ( x − t r a n s [ w , θ ] ) ] (35) \widehat{w}={\underset {w}{\operatorname {argmax} }}\,\Bigg[\sum_{i=1}^I\log[Norm_x[trans[w,\theta],\delta^2I]]\Bigg]\\ ={\underset {w}{\operatorname {argmin} }}\,\Bigg[\sum_{i=1}^I (x-trans[w,\theta])^T(x-trans[w,\theta]) \Bigg]\tag{35} w =wargmax[i=1Ilog[Normx[trans[w,θ],δ2I]]]=wargmin[i=1I(xtrans[w,θ])T(xtrans[w,θ])](35)
x = t r a n s [ w , θ ] (36) x=trans[w,\theta]\tag{36} x=trans[w,θ](36)
λ [ x y 1 ] = [ ϕ 11 ϕ 12 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ] [ u v 1 ] (37) \lambda \left[

xy1
\right]= \left[
ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]\left[
uv1
\right]\tag{37} λxy1=ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32ϕ13ϕ23ϕ33uv1(37)
λ ′ [ u v 1 ] = [ ϕ 11 ϕ 12 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ] − 1 [ x y 1 ] (38) \lambda' \left[
uv1
\right]= \left[
ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]^{-1}\left[
xy1
\right]\tag{38}
λuv1=ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32ϕ13ϕ23ϕ331xy1(38)

4 学习外部参数

  已知场景中平面上 I I I个不同的3D坐标点 { w i } i = 1 I \{w_i\}_{i=1}^I {wi}i=1I(假设 w i = 0 w_i=0 wi=0),与之对应的图像平面2D投影点 { x i } i = 1 I \{x_i\}_{i=1}^I {xi}i=1I和内部参数 Λ \Lambda Λ,求解外部参数 { Ω , τ } \{\Omega,\tau\} {Ω,τ}
Ω ^ ,   τ ^ = argmax ⁡ Ω ,   τ   [ ∑ i = 1 I log ⁡ [ P r ( x i ∣ w i , Λ , Ω , τ ) ] ] = argmax ⁡ Ω ,   τ [ ∑ i = 1 I log ⁡ [ N o r m x i [ p i n h o l e [ w i , Λ , Ω , τ ] , σ 2 I ] ] ] (39) \widehat{\Omega},\,\widehat{\tau}={\underset {\Omega,\,\tau}{\operatorname {argmax} }}\,\Bigg[\sum_{i=1}^I\log[Pr(x_i|w_i,\Lambda,\Omega,\tau)]\Bigg]\\ ={\underset{\Omega,\,\tau}{\operatorname{argmax}}}\Bigg[\sum_{i=1}^I \log\big[Norm_{x_i}[pinhole[w_i,\Lambda,\Omega,\tau],\sigma^2I]\big]\Bigg]\tag{39} Ω ,τ =Ω,τargmax[i=1Ilog[Pr(xiwi,Λ,Ω,τ)]]=Ω,τargmax[i=1Ilog[Normxi[pinhole[wi,Λ,Ω,τ],σ2I]]](39)
  和之前的情况一样,这个问题仍然没有解析解,但可以在齐次坐标的表示下求得一个合适的初始解:
λ [ x y 1 ] = λ ′ [ ϕ x γ δ x 0 ϕ y δ y 0 0 1 ] [ ω 11 ω 12 τ x ω 21 ω 22 τ y ω 31 ω 32 τ z ] [ u v 1 ] = [ ϕ 11 ϕ 12 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ] [ u v 1 ] (40) \lambda \left[

xy1
\right]=\lambda' \left[
ϕxγδx0ϕyδy001
\right] \left[
ω11ω12τxω21ω22τyω31ω32τz
\right]\left[
uv1
\right]=\left[
ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]\left[
uv1
\right] \tag{40} λxy1=λϕx00γϕy0δxδy1ω11ω21ω31ω12ω22ω32τxτyτzuv1=ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32ϕ13ϕ23ϕ33uv1(40)
  在2.5节中我们已经学习了如何求解homography Φ \Phi Φ,所以接下来我们只需要根据式(40)拆分矩阵 Φ \Phi Φ来获得旋转矩阵 Ω \Omega Ω和位移矩阵 τ \tau τ。具体为首先式(40)后两项左乘 Λ − 1 \Lambda^{-1} Λ1以消除内部参数的影响,可得:
Φ ′ = Λ − 1 Φ = λ ′ [ ω 11 ω 12 τ x ω 21 ω 22 τ y ω 31 ω 32 τ z ] = [ ϕ 11 ′ ϕ 1 2 ′ ϕ 13 ′ ϕ 21 ′ ϕ 22 ′ ϕ 23 ′ ϕ 31 ′ ϕ 32 ′ ϕ 33 ′ ] (41) \Phi'=\Lambda^{-1}\Phi= \lambda' \left[
ω11ω12τxω21ω22τyω31ω32τz
\right]=\left[
ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]\tag{41}
Φ=Λ1Φ=λω11ω21ω31ω12ω22ω32τxτyτz=ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32ϕ13ϕ23ϕ33(41)

  接下来估计旋转矩阵 Ω \Omega Ω中的前两列(这里的求解方法和Orthogonal Procrustes Problem类似):
[ ϕ 11 ′ ϕ 1 2 ′ ϕ 21 ′ ϕ 22 ′ ϕ 31 ′ ϕ 32 ′ ] = U L V T [ ω 11 ω 12 ω 21 ω 22 ω 31 ω 32 ] = U [ 1 0 0 1 0 0 ] V T (42) \left[
ϕ11ϕ12ϕ21ϕ22ϕ31ϕ32
\right]=ULV^T\\ \left[
ω11ω12ω21ω22ω31ω32
\right]=U\left[
100100
\right]V^T\tag{42}
ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32=ULVTω11ω21ω31ω12ω22ω32=U100010VT(42)

  为了求出旋转矩阵最后一列,我们只需要求前两列的叉积(旋转矩阵是正交矩阵),并保证最终的旋转矩阵行列式等于1.式(41)中的标量 λ ′ \lambda' λ可以按下式求解:
λ ′ = ∑ m = 1 3 ∑ n = 1 2 ϕ m n ′ / w m n 6 (43) \lambda'=\frac{\sum_{m=1}^3\sum_{n=1}^2 \phi_{mn}'/w_{mn}}{6}\tag{43} λ=6m=13n=12ϕmn/wmn(43)
  最后,位移矩阵为:
τ = [ ϕ 13 ′ , ϕ 23 ′ , ϕ 33 ′ ] T / λ ′ (44) \tau=[\phi_{13}',\phi_{23}',\phi_{33}']^T/\lambda'\tag{44} τ=[ϕ13,ϕ23,ϕ33]T/λ(44)

5 学习内部参数

  如果希望利用场景中一个平面上的点去标定相机,需要相机在多个位姿下对平面成像。基于此,此处的相机标定问题的描述为:已知场景中某个平面上 I I I个不同的3D坐标点 { w i } i = 1 I \{w_i\}_{i=1}^I {wi}i=1I和与之对应的 J J J个不同位姿下的图像平面2D投影点 { x i } i = 1 , j = 1 I , J \{x_i\}_{i=1,j=1}^{I,J} {xi}i=1,j=1I,J,利用最大似然方法求解内部参数 { Ω , τ } \{\Omega,\tau\} {Ω,τ}
Λ ^ = argmax ⁡ Λ [ max ⁡ Ω 1... J , τ 1... J   [ ∑ i = 1 I ∑ j = 1 J log ⁡ [ P r ( x i j ∣ w i , Λ , Ω j , τ j ) ] ] ] (45) \widehat{\Lambda}={\underset{\Lambda}{\operatorname{argmax}}}\Bigg[{\underset {\Omega_{1...J},\tau_{1...J}}{\operatorname {max} }}\,\Bigg[\sum_{i=1}^I \sum_{j=1}^J \log[Pr(x_{ij}|w_i,\Lambda,\Omega_j,\tau_j)]\Bigg]\Bigg]\tag{45} Λ =Λargmax[Ω1...J,τ1...Jmax[i=1Ij=1Jlog[Pr(xijwi,Λ,Ωj,τj)]]](45)
  和之前一样,这里介绍的方法采用了Coordinate ascent method,即迭代执行下面两步:

  • 求解相机 J J J个位姿下的外部参数:
    Ω ^ j ,   τ ^ j = argmax ⁡ Ω j ,   τ j   [ ∑ i = 1 I log ⁡ [ P r ( x i j ∣ w i , Λ , Ω j , τ j ) ] ] (46) \widehat{\Omega}_j,\,\widehat{\tau}_j={\underset {\Omega_j,\,\tau_j}{\operatorname {argmax} }}\,\Bigg[\sum_{i=1}^I\log[Pr(x_{ij}|w_i,\Lambda,\Omega_j,\tau_j)]\Bigg]\tag{46} Ω j,τ j=Ωj,τjargmax[i=1Ilog[Pr(xijwi,Λ,Ωj,τj)]](46)
  • 求解当前迭代步中的内部参数:
    Λ ^ = argmax ⁡ Λ   [ ∑ i = 1 I ∑ j = 1 J log ⁡ [ P r ( x i j ∣ w i , Λ , Ω j , τ j ) ] ] (47) \widehat{\Lambda}={\underset{\Lambda}{\operatorname{argmax}}}\,\Bigg[\sum_{i=1}^I \sum_{j=1}^J \log[Pr(x_{ij}|w_i,\Lambda,\Omega_j,\tau_j)]\Bigg]\tag{47} Λ =Λargmax[i=1Ij=1Jlog[Pr(xijwi,Λ,Ωj,τj)]](47)
      类似的,该方法的效率也很低,实际应用中还有更好的方法。

6 推理3D坐标点

  已知一个(注意:由于此处一一对应的关系,只需要一个相机便能推理)标定好并且已知位姿的相机(也即已知 Λ , Ω , τ \Lambda, \Omega, \tau Λ,Ω,τ)和场景中某个平面上一个未知3D坐标点 w = [ u , v , 0 ] T w=[u,v,0]^T w=[u,v,0]T在相机上的2D投影点坐标 x x x,求该3D坐标点的位置:
T = [ ϕ 11 ϕ 12 ϕ 13 ϕ 21 ϕ 22 ϕ 23 ϕ 31 ϕ 32 ϕ 33 ] = [ ϕ x γ δ x 0 ϕ y δ y 0 0 1 ] [ ω 11 ω 12 τ x ω 21 ω 22 τ y ω 31 ω 32 τ z ] (48) T=\left[

ϕ11ϕ12ϕ13ϕ21ϕ22ϕ23ϕ31ϕ32ϕ33
\right]=\left[
ϕxγδx0ϕyδy001
\right] \left[
ω11ω12τxω21ω22τyω31ω32τz
\right]\tag{48} T=ϕ11ϕ21ϕ31ϕ12ϕ22ϕ32ϕ13ϕ23ϕ33=ϕx00γϕy0δxδy1ω11ω21ω31ω12ω22ω32τxτyτz(48)
  在世界坐标系下为:
w ~ = T − 1 x ~ (49) \widetilde{w}=T^{-1}\widetilde{x}\tag{49} w =T1x (49)
  转换到相机坐标系下为:
w ′ = Ω w + τ (50) w'=\Omega w+\tau\tag{50} w=Ωw+τ(50)

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/weixin_40725706/article/detail/793964
推荐阅读
相关标签
  

闽ICP备14008679号