ぱたへね

はてなダイアリーはrustの色分けができないのでこっちに来た

オイラー角

cv

写真から作る3次元CGから

任意の回転運動は、Z軸まわりの回転α、新しいY軸まわりの回転β、新しいZ軸まわりの回転γの3段階の回転で実現できる。
このα、β、γをオイラー角といい、回転行列Rは次のように表される。

R=R(Z,\alpha )R(Y,\beta )R(Z,\gamma )=\begin{pmatrix} \cos \alpha & -\sin \alpha & 0 \\ \sin \alpha & \cos \alpha & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \cos \beta & 0 & \sin \beta \\ 0 & 1 & 0 \\ -\sin \beta & 0 & \cos \beta \end{pmatrix} \begin{pmatrix} \cos \gamma & -\sin \gamma & 0 \\ \sin \gamma & \cos \gamma & 0 \\ 0 & 0 & 1 \end{pmatrix}

回転のオイラー角による表現は唯一では無いが、例えばαの範囲を -π/2 < α < π/2 に制限すれば、唯一になる。

Iterative minimization method

cv

Multiple View Geometry in Computer Vision Second Edition4章から。
直接線形方程式を解くDLTアルゴリズムに比べた、イテレーティブに最小値を求める方法の特徴。

  • 遅い
  • イテレーションをはじめる前に初期値が必要
  • 収束しない可能性と、最小値でなくローカルな極小値に収束しない可能性がある。
  • イテレーションの終了方法が、トリッキーになることがある。

仕事で使っている感覚から行くと、適切な初期値と適切な終了条件があればDLTよりは速い。終了条件はどうしても現物あわせになってしまう。例えばGoogleの最小二乗法ライブラリCERESでもイテレーションの終了条件を細かく指定できる。裏返すと、一般的な終了条件というのが設定しにくいと言うこと。

イテレーションは次のステップで行われる。

  • Cost function コスト関数を選択する。
  • Parametrization 計算する変換を、有限のパラメータで表現する。一般的に、パラメータの数が最小である必要はなく、実際にはオーバーパラメータの方が良い。
  • Function specification 関数はパラメータの組で、コストが表現される。
  • Initialization 適切な初期パラメータが計算される。一般的にはDLTアルゴリズムのような線形のアルゴリズムが使われる。
  • Iteration 初期値からはじめて、コスト関数を最小にするGoalまで、反復計算を行う。

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

cv

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}

リファレンス

Gauche-CVで投影変換Hを求める

cv scheme

投影変換Hを求めるのに必要な条件の続き。
つまり4つの点の対応が分かれば、投影変換を表す行列Hは計算できます。

この図のような変換を使って実験してみます。（変換は青→黄）
座標の対応はこのようになります。

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

この4点の対応を引数にcv-get-perspective-transformを呼ぶと、投影変換を表すHが計算できます。

(use srfi-1)
(use cv)

(require "../proj2d/proj2d")
(import proj2d)
(load "../cvutil/cvutil.scm")

;; 変換元の4角形
(define src-rect '((5 5) (15 5) (15 15) (5 15)))

;; 行列Hの計算
(define pts-src (map (lambda (p) (make-cv-point-2d32f (first p) (second p))) src-rect))
(define pts-dst (map (lambda (p) (make-cv-point-2d32f (first p) (second p))) '((20 15) (25 20) (25 25) (15 20))))
(define H (make-cv-mat 3 3 CV_32F)) 
(cv-get-perspective-transform pts-src pts-dst H)

;; 結果の確認と表示
(format #t "Projective matrix H~%") 
(print-2d-matrix H)

(define display-transform
  (lambda (H)
    (lambda (x)
      (let ((xh (make-hm-point-2d (first x) (second x))))
        (let ((xdash (transform H xh)))
          (format #t "~A -> ~A~%" x xdash))))))
(map (display-transform H) src-rect)

計算した行列Hの中身はこうなります。
[ [ 1.75 -1.25 17.5 ] [ 1.5 -0.5 10.0 ] [ 0.05000000074505806 -0.05000000074505806 1.0 ]]

この投影変換Hを使って、もとの点を変換してみると、全ての点が正しく変換されていることがわかります。
(5 5) -> (20.0, 15.0)
(15 5) -> (25.0, 20.0)
(15 15) -> (25.0, 25.0)
(5 15) -> (15.0, 20.0)

ソースはこちら
https://github.com/natsutan/computervision/blob/master/mvg/chap4/perspective-transform.scm

投影変換Hを求めるのに必要な条件

cv

Multiple View Geometry in Computer Vision Second Edition4章から

  • 投影変換Hは3x3の行列だが、定数倍に関して等価なので、自由度(dof)は8
  • 点 xi -> x'iへの対応一つにつき、2次元空間で(x, y) -> (x', y')の対応が取れるので、制約が2つつく
  • 最低限4点の対応が取れれば、投影変換Hは求められる。
  • 点は正確に測定できないため、4点だけでは上手く行かない
  • より多くの点を対応づけたとき、ベストなHを求める問題に直面する。
  • Golden Standard アルゴリズムを使う。

2章 練習問題

cv

写真から作る3次元CGの2章練習問題。

問題1

D^-1を求めよ。D^-1はDの逆運動である。

答
Dは剛体変換(rigid transformation)で、回転行列Rと並進行列tを使って以下のように表現できる。Dの定義は投射行列、もしくは2次元変換の階層のClass I Isometries transformations に。

D=\begin{pmatrix} R & t \\ 0^{T}_{3} & 1 \end{pmatrix}

逆運動であるD'=D^-1も同じく剛体変換なので同じように表現できる。

D'=\begin{pmatrix} R' & t' \\ 0^{T}_{3} & 1 \end{pmatrix}

ここでDD'は何もしない（変換を元に戻す）ので、
DD'=I
ここから
RR'=I
Rt'+t=0

Rは回転行列なので以下の式を満たす。
RR^{T}=I

よって、次の2式からD'が求められる。
R'=R^{T}
t'=-R^{T}t

問題2

θ=π/2、αu=αvの時、A^-1を求めよ。

答
行列Aに問題の条件を代入する。（ただし、αu=αv=α）
A=\begin{pmatrix} \alpha _{u} & -\alpha _{u}\cot \theta & u_{0} \\ 0 & \frac{\alpha _{v}}{\sin \theta } & v_{0} \\ 0 & 0 & 1 \end{pmatrix} =\begin{pmatrix} \alpha & 0 & u_{0} \\ 0 & \alpha & v_{0} \\ 0 & 0 & 1 \end{pmatrix}
逆行列は
A^{-1} = \frac{ 1 }{ \alpha } \begin{pmatrix} 1 & 0 & -u_{0} \\ 0 & 1 & -v_{0} \\ 0 & 0 & \alpha \end{pmatrix}

問題3

Paを任意のアフィン投影行列とするとき、APaDもアフィン投影行列であることを確かめよ。

答
カメラ内部行列A、剛体変換行列Dは、共にアフィン投射行列である。アフィン投射行列の積は、またアフィン投射行列である。

問題4

Mi, i=1, ... n は3次元空間のn点であり、mi, i=1,... nは、それぞれのアフィン投影である。sm~ = PaM~ を照明せよ。m~ は、mi, i=1 ... nの重心であり、M~は、Mi, i=1 ... nの重心である。

答
3次元アフィン変換では、面の平行関係、体積比、重心は不変量。よって、Mの重心のアフィン投影は、mの重心に等しい。ただし、定数倍は等価。

問題5

任意の3x1ベクトルについて以下の式を示せ
r^{2}_{x}=rr^{T}-r^{T}r
（左辺はrの外積）

答
左辺の定義がよく分からないので、答え丸写しです。

r^{2}_{x}=\begin{pmatrix} -r_{3}^{2}-r_{2}^{2} & r_{1}r_{2} & r_{1}r_{3} \\ r_{1}r_{2} & -r_{1}^{2}-r_{3}^{2} & r_{2}r_{3} \\ r_{1}r_{3} & r_{2}r_{3} & -r_{1}^{2}^r_{2}^{2} \end{pmatrix} =rr^{T}-(r^{T}r)