最近2台のカメラのピクセル座標(2次元)から、ある一点のワールド座標(3次元)を求める際にDLTと特異値分解について学習しました。今回はその備忘録として書いておきたいと思います。
※ 私は独学で勉強している身なので、この記事には間違いがあると予想されます。ご了承ください。
DLTとは
DLT(Direct Linear Transformation)とは、変数の集合を比例関係の集合を用いて解くアルゴリズムです。
式に表すと以下の通りになります。
\mathbf{x_k} \propto \mathbf{A}\mathbf{y_k} for k=1,…,N – ①
このうち\mathbf{x_k}と\mathbf{y_k}はベクトル、Aは行列です。
ここで一般的な式\mathbf{x_k} = \mathbf{A}\mathbf{y_k} for k=1,…,Nを考えます。
この式を行列で表した式\mathbf{X}=\mathbf{A}\mathbf{Y}の\mathbf{A}を求めようとすると以下の操作を行います。
- 両辺に右から\mathbf(Y)^Tを掛ける
- \mathbf{Y}\mathbf{Y}^Tの逆行列を両辺に右から掛ける
※正方行列\mathbf{Y}\mathbf{Y}^Tが正則行列である場合にしか成り立ちませんが、ここでは簡単化のために正則行列であるとさせてください
\mathbf{A}=\mathbf{X}\mathbf{Y}^T(\mathbf{Y}\mathbf{Y}^T)^{-1}
上の操作で行列\mathbf{A}を求めることが出来ました。これがスタンダードな求め方です。
スタンダードな求め方ではスケールの値は当然既知である必要があります。
要は\alpha\mathbf{X}=\mathbf{A}\mathbf{Y}から解を求めることは出来ません。
一方でDLTは先ほどの方程式\alpha\mathbf{X}=\mathbf{A}\mathbf{Y}に不明なスケールの値\alphaが含まれていても計算できるアルゴリズムになります。
このスケールの値が不明でも計算が出来るという部分を表現したのが最初に紹介した①式となります。
DLTの仕組み
ではDLTはどのようにして解を求めるのでしょうか。DLTは方程式を操作し、左辺(右辺)に現れる項=0という状態を作り出します。そして左辺(右辺)が0になると、その逆の右辺(左辺)も0になります。すると左辺(右辺)が0なので、不明なスケールの値(スカラー値)を割ることができます。
以上より、スケールの値が不明であっても解を求めることが可能になる.. という訳です。
例えば0=\alpha\mathbf{A}\mathbf{x}という場合、スカラー値\alphaを取り除くことが出来るため\alphaが不明でも解を求めることが出来ます
DLTと特異値分解
コンピュータビジョン等ではDLTと特異値分解を組み合わせて使うケースがみられます。それは何故かというと、ノイズによって解を求められないケースが多いからです。
画像の座標から3次元の座標を求める際などに\mathbf{A}\mathbf{x}=0という式がよく出てきますが、実際の画像から得られる座標にはノイズがあるため式が成り立たないことがあります。
※\mathbf{A}\mathbf{x}=0は平行なベクトルの外積が0になる特性を用いて導き出したりします
そのため\mathbf{A}\mathbf{x}=0の解を求めるのではなく\mathbf{A}\mathbf{x}=\mathbf{w}の\mathbf{w}の大きさを最小化する解を求めるということを行います。この際に特異値分解が出現します。
実際に解いてみる
先ほどの式のベクトル\mathbf{x}が求めたい解であり、行列\mathbf{A}(ここではm\times{n}行列とする)が既に求まっているとします。
今回の式\mathbf{A}\mathbf{x}=0は右辺が0であるため、解は以下の2パターンに分類できます
- \mathbf{x}=0
- \mathbf{x}\neq0
1の解は自明なので、ここでは飛ばします。2の場合の解は右辺が0であるため解をスカラー倍したものも解になります。しがたって一意に求めることは出来ません。
なので\mathbf{x}=\alpha\mathbf{y} where \alpha\in\mathbf{R}, \alpha\neq{0}, ||\mathbf{y}||=1とします。
さて\mathbf{A}を特異値分解すると\mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{V}^Tが得られます。すると\mathbf{w}=\mathbf{A}\mathbf{x}=\mathbf{U}\mathbf{S}\mathbf{V}^T\mathbf{x}という式が得られます。
ここで\mathbf{w}の大きさを考えます。
||\mathbf{w}|| = ||\mathbf{U}\mathbf{S}\mathbf{V}^T\mathbf{x}|| – ②
特異値分解の定義より\mathbf{U},\mathbf{V}は直行行列、\mathbf{S}は対角行列です。直行行列はノルムを不変に保つ性質(||\mathbf{A}\mathbf{x}|| = ||\mathbf{x}||)があるので、②式から以下の式が得られます。
||\mathbf{w}|| = ||\mathbf{S}\mathbf{V}^T\mathbf{x}|| – ③
ここで\mathbf{x}=\alpha\mathbf{y}としたので、③式は以下のように変形できます
||\mathbf{w}||=\alpha||\mathbf{S}\mathbf{V}^T\mathbf{y}|| – ③’
さらにここで\mathbf{z}=\mathbf{V}^T\mathbf{y}と置きます。このとき||\mathbf{z}||=||\mathbf{V}^T\mathbf{y}|| = ||\mathbf{y}|| = 1です。
すると③’は以下のようになります
||\mathbf{w}||=\alpha||\mathbf{S}\mathbf{z}|| – ③”
③”より||\mathbf{w}||の最小値は||\mathbf{S}\mathbf{z}||の最小値により決まります。(\alphaはスケール値で不定)
\mathbf{S}は特異値で対角行列です。特異値は取り扱いを簡単にするため降順に並べることが多く、ここでも降順に並んでいると仮定したとき||\mathbf{w}||が最小となるのは\mathbf{z}=[0,0,…,1]^Tのときになります。
以上より
\mathbf{y}=(V^T)^{-1}\mathbf{z}=\mathbf{v}_n
\mathbf{x}=\alpha\mathbf{v}_n
となります。
※ コンピュータビジョンの場合だと求めるベクトルの最後の値を1にすることが多いので\alpha = \frac{1}{v_{nn}}となっているケースを良く見かける気がします。(気がするだけかもしれない)
参考ページ
最後に参考ページを記載します。これらのページには非常に参考になりました。これらのページの著者らに、この場を借りて御礼申し上げます。
- Direct linear transformation – Wikipedia
https://en.wikipedia.org/wiki/Direct_linear_transformation - 特異値分解 – Wikipedia
https://ja.wikipedia.org/wiki/%E7%89%B9%E7%95%B0%E5%80%A4%E5%88%86%E8%A7%A3 - Direct Linear Transform (DLT)
https://temugeb.github.io/computer_vision/2021/02/06/direct-linear-transorms.html - 同次式の最小二乗問題 – Daily Tech blog
https://daily-tech.hatenablog.com/entry/2018/03/18/180439
コメント