出典: Wikimura

ホモグラフィについて書かれた資料参照。ホモグラフィ(平面射影)は、空間内の平面と、画像(平面)の投影を調べるということ。

4つ以上の対応点が分かれば求められるらしい。

計算方法

「How to compute a homography」および「Computing the plane to plane homography」参照。 2次元の点対応が与えられると、ホモグラフィが定義できる。ホモグラフィは3x3の行列。スケールが関係するので、PとP'は数値的には等しくないが、同じものを表している。

$\begin{align} &\mathbf{p}, \mathbf{P} \in \mathcal{P}^2 \\ &\mathbf{P}' = \mathbf{H} \mathbf{p}\\ &\mathbf{P} = \frac{1}{W'}\mathbf{P}' \end{align}$

ここで鍵となるのは...9個の未知数を持つ行列についての式を、3つの3次元ベクトルについての式にする方法。行列を行ベクトルで表示すると、行ベクトルを3つ重ねた9次元のベクトルと、3行9列の行列との積の形になる。

$\mathbf{A} \mathbf{v} = \mathbf{0} \ \Rightarrow\ \begin{bmatrix} \mathbf{a}^{1 \mathrm{T}} \\ \mathbf{a}^{2 \mathrm{T}} \\ \mathbf{a}^{3 \mathrm{T}} \end{bmatrix} \ \mathbf{v}=\mathbf{0} \ \Rightarrow\ \left\{ \begin{align} \mathbf{a}^{1 \mathrm{T}} \mathbf{v} = 0\\ \mathbf{a}^{2 \mathrm{T}} \mathbf{v} = 0\\ \mathbf{a}^{3 \mathrm{T}} \mathbf{v} = 0 \end{align} \right. \ \Rightarrow\ \begin{bmatrix} \mathbf{v}^\mathrm{T} & \mathbf{0}^\mathrm{T} & \mathbf{0}^\mathrm{T}\\ \mathbf{0}^\mathrm{T} & \mathbf{v}^\mathrm{T} & \mathbf{0}^\mathrm{T}\\ \mathbf{0}^\mathrm{T} & \mathbf{0}^\mathrm{T} & \mathbf{v}^\mathrm{T} \end{bmatrix} \begin{bmatrix} \mathbf{a}^1 \\ \mathbf{a}^2 \\ \mathbf{a}^3 \end{bmatrix} = \mathbf{0}$

与えられた対応点からホモグラフィ行列を求めるには、行列をベクトルの形にする。既知なのは二次元の同次座標だけ。大文字で表したものは、小文字のものと同じ点を指す同次座標であるが、スケールは未知。

$\mathbf{p} = \begin{bmatrix} x\\ y\\ 1 \end{bmatrix} \ ,\ \mathbf{P} = \begin{bmatrix} X\\ Y\\ 1 \end{bmatrix} \ ,\ \mathbf{P}' = \begin{bmatrix} X'\\ Y'\\ W' \end{bmatrix} \ ,\ \mathbf{p} = \frac{1}{W'}\mathbf{P}$

ここで、スケールが次数を1つ減らす性質を利用すると、同次方程式からスケールについて解くことができる。 Wをh3とpで表すことができる。

$\mathbf{P}' = \mathbf{H}\mathbf{p} \ \Rightarrow \ \left\{ \begin{align} X' = \mathbf{h}^{1\mathrm{T}} \mathbf{p} \\ Y' = \mathbf{h}^{2\mathrm{T}} \mathbf{p} \\ W' = \mathbf{h}^{3\mathrm{T}} \mathbf{p} \end{align} \right. \ \Rightarrow \ \left\{ \begin{align} \mathbf{p}^\mathrm{T} \mathbf{h}^1 - X \mathbf{p}^\mathrm{T} \mathbf{h}^3 = 0\\ \mathbf{p}^\mathrm{T} \mathbf{h}^2 - Y \mathbf{p}^\mathrm{T} \mathbf{h}^3 = 0 \end{align} \right.$

最終的に次式で表すことができる。

$\begin{bmatrix} \mathbf{p}^\mathrm{T} & \mathbf{0}^\mathrm{T} & -X \mathbf{p}^\mathrm{T} \\ \mathbf{0}^\mathrm{T} & \mathbf{p}^\mathrm{T} & -Y \mathbf{p}^\mathrm{T} \end{bmatrix} \begin{bmatrix} \mathbf{h}^1 \\ \mathbf{h}^2 \\ \mathbf{h}^3 \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix}$

1つの点対応で同次方程式が2つが立てられるので、8自由度全てを決定するには4つの点対応が必要となる。もちろん5つ以上の点対応があっても良いが、その場合は固有値分解や特異値分解を使うらしい。

同次方程式の解

ホモグラフィを求める式は同次方程式となった。同次方程式(Ax = 0)の解は行列Aの性質で変わってくる。行列が特異値(singular value)を持たなければ、自明な解(x=0)しか持たない。特異値を持つなら、そうでない解が(どれくらいの自由度かはものによるが)得られるらしい。「コンピュータビジョン3章:同時系と固有値問題」参照。

「同次連立1次方程式の理論」によると、方程式と未知数の数が一致するとき方程式が解をもつ条件は、係数行列の行列式が0となることとある。

過制約の場合、誤差の絶対値の二乗和が最小となるようにする。同次方程式の場合、以下のようになる。

$\begin{align} \mathbf{e} &= \mathbf{A}\mathbf{x} \\ E &= | \mathbf{A}\mathbf{x} |^2 \\ &= \mathbf{x}^\mathrm{T}\mathbf{A}^\mathrm{T}\mathbf{A}\mathbf{x} \end{align}$

コンピュータビジョンの本によれば、誤差を最小化するベクトル(求めたいもの)は、対称行列 $A T A$ の最小の固有値に対応する固有ベクトルらしい。この時のEの最小値は、最小固有値の2乗となるらしい。

Matlabで計算してみる

解となる固有ベクトルを得るなら、ヤコビ変換、三重対角行列からのQR分解、特異値分解があるらしい。

|x|=1の制約のもと、||Ax||を最小化するベクトルは、 $A T A$ の最小固有値に対する固有ベクトルで与えられる。この固有値は行列AのSVD(固有値分解)から直接得ることができる。

試しにMatlab以下の4点対応についてホモグラフィ行列を求めてみる。これらの点を同次座標で表現し、行ベクトルにして重ね合わせた行列をP,Qとしている。

(1,1) - (4,4)
(-1,1) - (-4,4)
(-1,-1) - (-4,-4)
(1,-1) - (4,-4)

Matlabでテスト(固有値を使う場合)

Matlabのeig関数の戻り値Vは固有ベクトル(列)の集まり、Eは対応する固有値の行列となる。 Vのi列目にある固有ベクトルに対する固有値は、Eのi列目にある。

P = [
    1   1   1;
    -1  1   1;
    -1  -1  1;
    1   -1  1;
    ];

Q = [
    4   4   1;
    -4  4   1;
    -4  -4  1;
    4   -4  1;
    ];

z = zeros(1,3);

A = [
    P(1,:)  z       -Q(1,1) * P(1,:);
    z       P(1,:)  -Q(1,2) * P(1,:);
    P(2,:)  z       -Q(2,1) * P(2,:);
    z       P(2,:)  -Q(2,2) * P(2,:);
    P(3,:)  z       -Q(3,1) * P(3,:);
    z       P(3,:)  -Q(3,2) * P(3,:);
    P(4,:)  z       -Q(4,1) * P(4,:);
    z       P(4,:)  -Q(4,2) * P(4,:);
    ];

[V E] = eig( A'*A)

------------------------------------------------------------------------------------------
V =

   -0.6963    0.0000   -0.0000    0.6887   -0.1602         0   -0.0000   -0.0000   -0.1231
         0         0         0         0         0    1.0000         0         0         0
    0.0000    0.9920   -0.0066   -0.0000    0.0000         0   -0.0000   -0.1259   -0.0000
         0   -0.0000    0.0000    0.2265    0.9740         0         0   -0.0000         0
   -0.6963   -0.0000   -0.0000   -0.6887    0.1602         0   -0.0000    0.0000   -0.1231
    0.0000   -0.0066   -0.9920    0.0000   -0.0000         0    0.1259   -0.0000    0.0000
   -0.0000    0.1259   -0.0008    0.0000   -0.0000         0    0.0000    0.9920    0.0000
    0.0000   -0.0008   -0.1259    0.0000   -0.0000         0   -0.9920   -0.0000   -0.0000
   -0.1741   -0.0000   -0.0000    0.0000   -0.0000         0   -0.0000   -0.0000    0.9847


E =

   -0.0000         0         0         0         0         0         0         0         0
         0    1.9688         0         0         0         0         0         0         0
         0         0    1.9688         0         0         0         0         0         0
         0         0         0    4.0000         0         0         0         0         0
         0         0         0         0    4.0000         0         0         0         0
         0         0         0         0         0    4.0000         0         0         0
         0         0         0         0         0         0  130.0312         0         0
         0         0         0         0         0         0         0  130.0312         0
         0         0         0         0         0         0         0         0  132.0000

これで $\mathbf{v}^\mathrm{T}\mathbf{A}^\mathrm{T}\mathbf{A}\mathbf{v}$ を求めると、対応する固有値になる。実は以下のような関係がある。

$\mathbf{E} = \mathbf{V}^\mathrm{T}\mathbf{A}^\mathrm{T}\mathbf{A}\mathbf{V}$

V' * A' * A * V
ans =

    0.0000   -0.0000   -0.0000    0.0000   -0.0000         0   -0.0000   -0.0000   -0.0000
   -0.0000    1.9688   -0.0000    0.0000   -0.0000         0    0.0000   -0.0000   -0.0000
   -0.0000    0.0000    1.9688   -0.0000    0.0000         0    0.0000    0.0000    0.0000
    0.0000    0.0000   -0.0000    4.0000   -0.0000         0   -0.0000    0.0000   -0.0000
    0.0000   -0.0000    0.0000   -0.0000    4.0000         0    0.0000   -0.0000    0.0000
         0         0         0         0         0    4.0000         0         0         0
   -0.0000    0.0000    0.0000   -0.0000    0.0000         0  130.0312    0.0000   -0.0000
   -0.0000   -0.0000    0.0000    0.0000   -0.0000         0    0.0000  130.0312    0.0000
   -0.0000    0.0000   -0.0000   -0.0000    0.0000         0   -0.0000    0.0000  132.0000

最終的にホモグラフィ行列が求まった。ちゃんとスケーリングになった。

H = [ V(1,1) V(2,1) V(3,1); V(4,1) V(5,1) V(6,1);V(7,1) V(8,1) V(9,1)] / V(9,1);
------------------------------------------------------------------------------------------
H =
    4.0000         0   -0.0000
         0    4.0000   -0.0000
    0.0000   -0.0000    1.0000

Matlabでテスト(特異値分解を使う場合)

「行列の特異値分解」に、 $A T A$ との固有ベクトルが $A$ の特異値分解から求まることが示されている。

ランクrの行列Aは3つの行列の積に分解できるという。 UとVはそれぞれ左特異ベクトル、右特異ベクトルと呼ぶ。中央の $Σ$ は特異値(Singular value)からなる対角行列となる。

$\mathop{A}_{n \times m}=\mathop{U}_{n \times r} \mathop{\Sigma}_{r \times r} \mathop{V'}_{r \times m}$

UとVは正規直行ベクトルを列に持つため、以下のような性質がある。

$U'U = V'V = I_{r \times r}$

この性質を利用すると、固有値・固有ベクトル分解の形になる。

$\begin{align} A'A &= V \Sigma'\Sigma V'\\ AA' &= U \Sigma'\Sigma U' \end{align}$

この関係から、以下のことが言える。

Vの列ベクトルはA'Aの固有ベクトル
Uの列ベクトルはAA'の固有ベクトル

先ほどのプログラムに以下を加えてみた。固有値を求める関数を使った時と同じ結果になった。ただ、固有値を求めるeig関数は固有値が小さい順に固有ベクトルを並べるが、特異値分解をするsvd関数は固有値が大きい順に並べるらしい。そのため、ホモグラフィ行列の計算の仕方が変わってくる。

[U S V] = svd( A);
H_svd = [ V(1,9) V(2,9) V(3,9); V(4,9) V(5,9) V(6,9);V(7,9) V(8,9) V(9,9)] / V(9,9)
------------------------------------------------------------------------------------------
H_svd =

    4.0000    0.0000    0.0000
   -0.0000    4.0000    0.0000
    0.0000   -0.0000    1.0000

Matlabでテスト(まとめ)

P = [
    1   1   1;
    -1  1   1;
    -1  -1  1;
    1   -1  1;
    ];

Q = [
    4   4   1;
    -4  4   1;
    -4  -4  1;
    4   -4  1;
    ];

z = zeros(1,3);

A = [
    P(1,:)  z       -Q(1,1) * P(1,:);
    z       P(1,:)  -Q(1,2) * P(1,:);
    P(2,:)  z       -Q(2,1) * P(2,:);
    z       P(2,:)  -Q(2,2) * P(2,:);
    P(3,:)  z       -Q(3,1) * P(3,:);
    z       P(3,:)  -Q(3,2) * P(3,:);
    P(4,:)  z       -Q(4,1) * P(4,:);
    z       P(4,:)  -Q(4,2) * P(4,:);
    ];


[V E] = eig( A'*A);
H_eig = [ V(1,1) V(2,1) V(3,1); V(4,1) V(5,1) V(6,1);V(7,1) V(8,1) V(9,1)] / V(9,1)

[U S V] = svd( A);
H_svd = [ V(1,9) V(2,9) V(3,9); V(4,9) V(5,9) V(6,9);V(7,9) V(8,9) V(9,9)] / V(9,9)
------------------------------------------------------------------------------------------
H_eig =

    4.0000         0   -0.0000
         0    4.0000   -0.0000
    0.0000   -0.0000    1.0000

H_svd =

    4.0000    0.0000    0.0000
   -0.0000    4.0000    0.0000
    0.0000   -0.0000    1.0000

参考文献

ホモグラフィ行列を求める