takminの書きっぱなし備忘録 @はてなブログ

主にコンピュータビジョンなど技術について、たまに自分自身のことや思いついたことなど

【改訂版】三次元空間上の楕円の公式から回転行列と並進ベクトルを計算する

三次元空間上の楕円は以下の式で表されます。
(1)


例えば三次元空間上の点群があり、この点群に上記の楕円をあてはめたい時などは、RANSACや最小二乗法などを用いてaからjまでの各係数を計算できます。


で、せっかくこういう係数がもとまったは良いけど、そこから楕円の中心座標や長軸/短軸の長さ、楕円がある平面の法線ベクトルなどを算出する方法が意外とややこしく、色々と調べたけどやり方が見つからなかったため自分で導出してみました。


---------- (2015/11/10追記) ----------

XY平面における楕円の式は(1)'で現れます。
(1)'

(1)式が10個の変数を持つのに対し、(1)式は6個の変数を持ちます。したがって(1)式を直接求めるよりも3次元上の点を一度最小二乗法などで求めた2次元平面上に射影した上で、(1)'式から楕円の公式を求めて、3次元空間へ戻した方が正しく求まります。

以下は3次元空間上における公式(1)から楕円パラメータを求める方法ですが、手順は2次元空間上の公式(1)'に対しても使用可能です。

--------------------------------------

なお、二次元空間上の楕円の公式から中心座標や長軸/短軸の長さ、傾きなどを求める方法は以下のサイトに解説がありました。


画像処理ソリューション「最小二乗法による楕円近似」
http://imagingsolution.blog107.fc2.com/blog-entry-20.html



まずXY平面上に、原点を中心とする傾きなしの標準的な楕円を考えます。



この楕円上の点の三次元空間上における座標は以下で表されます。

(2)


次にこの楕円を三次元空間上の任意の場所へ移動させることを考えます。
任意の三次元空間上の移動は回転行列Rと並進ベクトルTで表すことができます。

(3)


ここで三次元空間上の楕円座標X'、回転行列R、並進ベクトルTは以下のように表します。




Rは回転行列なので、それぞれの列ベクトルは直交しており、ノルムは1です。


(3)式を展開すると
(4)
(5)


となり、sinとcosを消去すると、以下の式が成り立ちます。

(6)


(6)式は以下のように変形することができます。

(7)

ここでQは以下の通りです。

(8)


(8)式を展開します。

(9)


(9)式のX'の係数と(1)式の係数を比較することで、行列Qと並進ベクトルTを求めることができます。
Qは対象行列なので、以下のように求まります。
(10)


これを用いてTは以下のように求まります。
(11)
(12)


---------- (2015/11/10追記) ----------

また定数項では以下の式が成り立ちます。
(A)

さて、(1)式は左辺をs倍しても成り立ちます。その時

のように、変形されるため(A)式は、
(B)

のように変形され、(B)式を満たすsは
(C)

から求まります。

このsを用いてQを以下のように変形します。

--------------------------------------


さて、Qはそもそも(8)式で表せるので、(8)式の両辺に回転行列の列ベクトルr1, r2, r3をかけてやると、
(13)
(14)
(15)


と表せます。
すなわち回転行列RはQの固有ベクトルから求まり、またその第1/第2固有値から短軸/長軸の長さが求まります。


手順をまとめると、
1. RANSACや最小二乗法等で(1)式の係数a-jを求める。
2. (10)式から行列Qを求める。
3. (12)式から並進ベクトルTを求める。
4. (C)式からスケールsを求める。
5. 行列Qをs倍する。
5. Qの第1固有値λ1と固有ベクトルr1, 第2固有値λ2と固有ベクトルr2を求める。
6. r1とr2の外積を計算しr3(楕円がある平面の法線ベクトル)とする。
7. 回転行列 R=(r1 r2 r3)
8. 短軸の長さ:A=sqrt(1/λ1)、長軸の長さ:B=sqrt(1/λ2)


以上、我流で求めたので使用するときは自己責任でお願いします。
また、もし何か誤りがあればご指摘ください。