Direct Linear Transform アルゴリズムを使って、投影変換の行列Hを計算する。

Multiple View Geometry in Computer Vision Second Edition4章から。
対応する4つの点を使って、投影変換を表す行列Hの計算します。Direct Linear Transform (DLT) とたいそうな名前がついていますが、点同士の制約から最小二乗法を用いて行列を求めるアルゴリズムです。

入力データ

Gauche-CVで投影変換Hを求めるのデータをそのまま使います。

xi	x'i
(5,5)	(20,15)
(15,5)	(25,20)
(15,15)	(25,25)
(5,15)	(15,20)

点の数は計算が一番楽な4点を使います。

DLTアルゴリズム

STEP1 それぞれの対応する点について、Ai行列を計算する。
STEP2 n個の2x9行列を一つの2nx9の行列Aにする。
STEP3 AのSVD (singular value decomposition, 特異値分解) を行う。A=UDV^Tの形になった時、一番小さな固有値に対応する単位固有ベクトルが求めるHの要素である。
STEP4 SVDの結果からHを決定する。

STEP1 それぞれの対応する点について、Ai行列を計算する

投影変換の式(i=1,2,3,4)
$x'_{i}=Hx_{i}$

この式は外積を使った表現と同じである。
$x'_{i}\times Hx_{i}=0$

xi, Hxiをこのように表現する。（h^1T の1は添字でTは転置の意味）
$Hx_{i}=\begin{pmatrix} h^{1T}x_{i} \\ h^{2T}x_{i} \\ h^{3T}x_{i} \end{pmatrix}$
$x'_{i}=\begin{pmatrix} x'_{i} \\ y'_{i} \\ w'_{i} \end{pmatrix}$

Hxiに代入すると

$Hx_{i}=\begin{pmatrix} y'_{i}h^{3T}x_{i}-w'_{i}h^{2T}x_{i} \\ w'_{i}h^{1T}x_{i}-x'_{i}h^{3T}x_{i} \\ x'_{i}h^{2T}x_{i}-y'_{i}h^{1T}x_{i} \end{pmatrix}$

hの表現を変え、Hxi を使って x'i × Hxi = 0 を書き換えると、

$\begin{pmatrix} 0^{T} & -w'_{i}x_{i}^{T} & y'_{i}x_{i}^{T} \\ w'_{i}x_{i}^{T} & 0^{T} & -x'_{i}x_{i}^{T} \\ -y'_{i}x_{i}^{T} & x'_{i}x_{i}^{T} & 0^{T} \end{pmatrix} \begin{pmatrix} h^{1} \\ h^{2} \\ h^{3}\end{pmatrix} =0$
$h^{1}=\begin{pmatrix} h^{1} \\ h^{2} \\ h^{3} \end{pmatrix} H=\begin{pmatrix} h_{1} & h_{2} & h_{3} \\ h_{4} & h_{5} & h_{6} \\ h_{7} & h_{8} & h_{9} \end{pmatrix}$

Hは等倍しても等価なので最後のrowは削除でき、Aih=0の形にできる。このときのAiはこの式で表現できる。

$A_{i}=\begin{pmatrix} 0^{T} & -w'_{i}x_{i}^{T} & y'x_{i}^{T} \\ w'_{i}x_{i}^{T} & 0^{T} & -x'_{i}x_{i}^{T} \end{pmatrix}$

(x'i, y'i, w'i)は同次座標なのでw'i = 1に取れる。このときx'i, y'iは実際の座標系の値をそのまま使用できる。次にxiも同次座標系(xi, yi, wi)を使って表現し、行列の要素を書き出すと、2x9の行列が求まる。

$A_{i}=\begin{pmatrix} 0^{T} & -w'_{i}x_{i}^{T} & y'x_{i}^{T} \\ w'_{i}x_{i}^{T} & 0^{T} & -x'_{i}x_{i}^{T} \end{pmatrix} =\begin{pmatrix} 0 & 0 & 0 & x_{i} & y_{i} & w_{i} & y_{i}'x_{i} & y'_{i}y_{i} & y'_{i}w_{i} \\ x_{i} & y_{i} & w_{i} & 0 & 0 & 0 & -x'_{i}x_{i} & -x'_{i}y_{i} & -x'_{i}w_{i} \end{pmatrix}$

Aih=0を同じように表現すると、
$$\begin{pmatrix} 0 & 0 & 0 & x_{i} & y_{i} & w_{i} & y_{i}'x_{i} & y'_{i}y_{i} & y'_{i}w_{i} \\ x_{i} & y_{i} & w_{i} & 0 & 0 & 0 & -x'_{i}x_{i} & -x'_{i}y_{i} & -x'_{i}w_{i} \end{pmatrix} \begin{pmatrix} h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ h_{8} \\ h_{9} \end{pmatrix} =0$

STEP2 n個の2x9行列を一つの2nx9の行列Aにする。

STEP1で計算したAiをi=1,2,3,4として4つ並べる。
$A=\begin{pmatrix} 0 & 0 & 0 & x_{1} & y_{1} & w_{1} & y'_{1}x_{1} & y'_{1}y_{1} & y'_{1}w_{1} \\ x_{1} & y_{1} & w_{1} & 0 & 0 & 0 & -x'_{1}x_{1} & -x'_{1}y_{1} & -x'_{1}w_{1} \\ 0 & 0 & 0 & x_{2} & y_{2} & w_{2} & y'_{2}x_{2} & y'_{2}y_{2} & y'_{2}w_{2} \\ x_{2} & y_{2} & w_{2} & 0 & 0 & 0 & -x'_{2}x_{2} & -x'_{2}y_{2} & -x'_{2}w_{2} \\ 0 & 0 & 0 & x_{3} & y_{3} & w_{3} & y'_{3}x_{3} & y'_{3}y_{3} & y'_{3}w_{3} \\ x_{3} & y_{3} & w_{3} & 0 & 0 & 0 & -x'_{3}x_{3} & -x'_{3}y_{3} & -x'_{3}w_{3} \\ 0 & 0 & 0 & x_{4} & y_{4} & w_{4} & y'_{4}x_{4} & y'_{4}y_{4} & y'_{4}w_{4} \\ x_{4} & y_{4} & w_{4} & 0 & 0 & 0 & -x'_{4}x_{4} & -x'_{4}y_{1} & -x'_{4}w_{4} \end{pmatrix}$

wi=1として、実際のx'i 等の値を代入する。

$A=\begin{pmatrix} 0 & 0 & 0 & 5 & 5 & 1 & 75 & 75 & 15 \\ 5 & 5 & 1 & 0 & 0 & 0 & -100 & -100 & -20 \\ 0 & 0 & 0 & 15 & 5 & 1 & 300 & 100 & 20 \\ 15 & 5 & 1 & 0 & 0 & 0 & -375 & -125 & -25 \\ 0 & 0 & 0 & 15 & 15 & 1 & 375 & 375 & 25 \\ 15 & 15 & 1 & 0 & 0 & 0 & -375 & -375 & -25 \\ 0 & 0 & 0 & 5 & 15 & 1 & 100 & 300 & 20 \\ 5 & 15 & 1 & 0 & 0 & 0 & -75 & -225 & -15 \end{pmatrix}$

STEP3 AのSVD (singular value decomposition, 特異値分解) を行う

STEP2で求めたAから、M=AA^T N=A^TAを求める。ここから数値計算ソフト使ってます。

$M=A^{t}A=\begin{pmatrix} 11526 & -15300 & 30401 & -37875 & 56776 & -56625 & 30401 & -22725 \\ -15300 & 20451 & -40400 & 50601 & -75500 & 75651 & -40400 & 30401 \\ 30401 & -40400 & 100651 & -125500 & 150801 & -150500 & 60551 & -45300 \\ -37875 & 50601 & -125500 & 157126 & -188125 & 188426 & -75500 & 56776 \\ 56776 & -75500 & 150801 & -188125 & 282326 & -281875 & 150801 & -112875 \\ -56625 & 75651 & -150500 & 188426 & -281875 & 282326 & -150500 & 113176 \\ 30401 & -40400 & 60551 & -75500 & 150801 & -150500 & 100651 & -75300 \\ -22725 & 30401 & -45300 & 56776 & -112875 & 113176 & -75300 & 56726 \end{pmatrix}$
$N=AA^{t}=\begin{pmatrix} 500 & 400 & 40 & 0 & 0 & 0 & -12125 & -9125 & -925 \\ 400 & 500 & 40 & 0 & 0 & 0 & -9125 & -10125 & -825 \\ 40 & 40 & 4 & 0 & 0 & 0 & -925 & -825 & -85 \\ 0 & 0 & 0 & 500 & 400 & 40 & 11000 & 9000 & 850 \\ 0 & 0 & 0 & 400 & 500 & 40 & 9000 & 11000 & 850 \\ 0 & 0 & 0 & 40 & 40 & 4 & 850 & 850 & 80 \\ -12125 & -9125 & -925 & 11000 & 9000 & 850 & 543125 & 420625 & 40375 \\ -9125 & -10125 & -825 & 9000 & 11000 & 850 & 420625 & 463125 & 36375 \\ -925 & -825 & -85 & 850 & 850 & 80 & 40375 & 36375 & 3525 \end{pmatrix}$

Nの固有値と固有ベクトルを計算する。ここも数値計算ソフトで計算しています。Nの固有値は、大きい順に λ＝929719.210541631, 80711.1190674933, 898.252012306525, 342.925445638261, 105.025927146904, 5.46731883242974, 0.91432069525331, 0.0853662562132915, 3.25254551524502E-13 となる。
それぞれの固有値に対応する固有ベクトルを並べるとこうなる。
$\begin{pmatrix} -0.0162896302429004 & 0.0174867033962668 & 0.485183442195986 & -0.0717208456407652 & -0.372990540295905 & 0.294326030956782 & -0.572600818755466 & 0.444825806571808 & 0.0859683627249945 \\ -0.0146232278814179 & -0.0167112389688832 & 0.512538485354574 & -0.0762946298017509 & 0.609715972594949 & 0.541180340040285 & 0.210455292133631 & -0.135271855740573 & -0.0614059733758339 \\ -0.00133685827734837 & 0.000152647903864893 & 0.0444860860149908 & -0.0336955645005135 & 0.010928031537983 & -0.0350396855621747 & 0.393532972688537 & 0.318767071026167 & 0.859683627267237 \\ \0.0153010896863559 & -0.00915944456408727 & 0.512595585916595 & 0.0720038604714055 & -0.597762363743942 & -0.0843855262482788 & 0.445640053137146 & -0.40407619870698 & -0.0736871680499479 \\ 0.0151572329038795 & 0.0258648577795804 & 0.484553400420063 & 0.115237168157642 & 0.361927320397211 & -0.760083307135421 & -0.20177976158078 & 0.0315160665229358 & 0.0245623893504672 \\ 0.00129502538142896 & 0.000710045126191261 & 0.0441873142750883 & 0.0311627569388701 & -0.0123566028885868 & -0.0868179424720417 & 0.481398915694935 & 0.718550511470746 & -0.49124778701298 \\ 0.738196676268601 & -0.67237267244566 & 0.00627063733708718 & -0.0470345721016824 & 0.00254607788395389 & -0.00377302891436078 & -0.0209191350155855 & 0.0162703119079957 & 0.0024562389349873 \\ 0.671330661990201 & 0.739301430794981 & -0.00492172816395286 & -0.0417833368336346 & 0.00495774310201827 & 0.0281647325287589 & 0.0111558810527539 & -0.00549070891304084 & -0.00245623893503336 \\ 0.0586028091048636 & -0.00314955728418357 & -0.0189742192713949 & 0.982092909158512 & 0.0225852079073412 & 0.161480195764102 & -0.0367636058166893 & 0.0365860082854019 & 0.0491247787008782 \end{pmatrix}$

行列全体はこちら。

一番小さい固有値 3.25254551524502E-13 に対応する固有ベクトルは、上の行列の一番右のcolumnで表されている。
右端のcolumn { 0.0859683627249945 -0.0614059733758339 0.859683627267237 -0.0736871680499479 0.0245623893504672 -0.49124778701298 0.0024562389349873 -0.00245623893503336 0.0491247787008782} が求める行列Hの要素である。

STEP4 SVDの結果からHを決定する。

{ 0.0859683627249945 -0.0614059733758339 0.859683627267237 -0.0736871680499479 0.0245623893504672 -0.49124778701298 0.0024562389349873 -0.00245623893503336 0.0491247787008782} を行列の形式に並び替える。
$H=\begin{pmatrix} 0.0859683627249945 & -0.0614059733758339 & 0.859683627267237 \\ -0.0736871680499479 & 0.0245623893504672 & -0.49124778701298 \\ 0.0024562389349873 & -0.00245623893503336 & 0.0491247787008782 \end{pmatrix}$

投影変換を表す行列Hはスカラー倍に対して等価なので、右下の要素が1になるように、 0.0491247787008782の逆数を全要素にかけると、Gauce-CVで求めた結果とほぼ同じになる。

$H=\begin{pmatrix} 1.7499999999686036 & -1.2499999999946292 & 17.50000000003804 \\ -1.4999999999721243 & 0.5000000000005721 & -10.000000000085457 \\ 0.04999999999884764 & -0.04999999999978525 & 1.0 \end{pmatrix}$