６点からの楕円の近似計算（最小二乗法、xyz次元）

みなさん、こんにちわ。
楕円についてプログラムしないといけなく、わからない点が多いので質問させてください。
ある６点（座標はxyz）からxyz楕円の最小二乗法をプログラムすることは可能でしょうか？この掲示板の「楕円のプログラム」を参考にするとある５点以上の測定点からxy座標（２次元）のグラフにすることは可能であると書いてあり、納得しました。
分からない点はxyzの楕円の基本方程式です。
なお、使用する環境はｃ言語です。参考になるサイトなどありましたら教えてください。

No.2 ベストアンサー 回答者： noocyte 回答日時： 2006/11/04 13:28

ちょっと訂正＆補足です．

> 合計８自由度ですから，８点以上ないと決定できないと思います (いまいち自信不足)．

直観的に考えてみると，５点以上あれば決定できるはずですね．m(_ _)m
ただし，この５点は同一平面上にないといけませんが．
そのことが自由度にどう影響しているのか，
また「３次元空間内の平面図形」ということがどう影響しているのか，
別途考え直してみます．

自由度 (Wikipedia)
http://ja.wikipedia.org/wiki/%E8%87%AA%E7%94%B1% …


３次元空間内の楕円に直接最小自乗法を適用すると，計算がとても大変そうなので，
次のようにした方がよいかもしれません．

(1) 与えられた点から，楕円が存在する平面を最小自乗法で決定する．
(2) 与えられた点の座標を上記の平面の座標に変換する．
(3) 平面座標から最小自乗法で楕円を決定する．

No.1 回答者： noocyte 回答日時： 2006/11/04 04:54

xyzということは３次元ですか？

それじゃ楕円ではなくて楕円体ですね．
それとも３次元空間内での楕円ということでしょうか？
多分後者だと思いますので，その前提で書きます．


> この掲示板の「楕円のプログラム」を参考にするとある５点以上の測定点から
> xy座標（２次元）のグラフにすることは可能であると書いてあり、納得しました。

２次元では，楕円の自由度は次のようになります．

中心座標：２
長軸の長さ：１
短軸の長さ：１
長軸の方向 (Ｘ軸に対する角度)：１

合計５自由度ですから，５点以上あれば楕円を決定することができます．

３次元での楕円の自由度は次のようになります．

中心座標：３
長軸の長さ：１
短軸の長さ：１
長軸の方向：２
長軸まわりの回転角：１
(あるいは最後の２つをまとめて，３次元での回転＝３自由度としてもよい)

合計８自由度ですから，８点以上ないと決定できないと思います (いまいち自信不足)．


> 分からない点はxyzの楕円の基本方程式です。

普通の楕円の方程式は探せばいくらでもあると思いますが，
３次元内での楕円の方程式はちょっとなさそうですね．

まず，楕円を含む平面をＵＶ平面とし，それに垂直なＷ軸を考えます．
楕円の中心をＵＶ平面の原点，長軸，短軸方向をＵ軸，Ｖ軸方向とし，
それぞれの長さを 2 * a，2 * b とすると，

(u/a)^2 ＋ (v/b)^2 ＝ 1．
w ＝ 0．

ＵＶＷ系とＸＹＺ系の原点が一致している (つまり楕円の原点のＸＹＺ座標が (0, 0, 0))
の場合には，ＸＹＺ座標をＵＶＷ座標に変換する回転行列をＲとすると，

U ＝ R X (ただし U＝(u, v, w)，X＝(x, y, z))．

楕円の原点が X0＝(x0, y0, z0) にある場合は，

U ＝ R (X - X0)．

Ｒは３×３の行列だが，これを２×３の行列Ｒuv (ＸＹＺをＵＶに変換する) と，
１×３の行列Ｒw (ＸＹＺをＷに変換する) に分けて，

R ＝ (
Ruv
Rw
)

とすると，

(u, v) ＝ Ruv (X - X0)，
w = Rw (X - X0)．

また

A ＝ (
1/(a^2), 0
0, 1/(b^2)
)

とすると，楕円の方程式は
t(u, v) A (u, v) ＝ 1．(t は転置)
これに (u, v) ＝ Ruv (X - X0) を代入して
t(X - X0) tRuv A Ruv (X - X0) ＝ 1 … (1)
w ＝ Rw (X - X0) ＝ 0 … (2)

(1)，(2) が楕円の方程式になります．



> 最小二乗法をプログラムすることは可能でしょうか？

最小二乗法だと，何を誤差とするかによりますね．
素直に考えれば，点と楕円の (最短) 距離だと思いますが，
かなり面倒な式になりそうですね．




楕円 (Wikipedia)
http://ja.wikipedia.org/wiki/%E6%A5%95%E5%86%86

Ellipsoid (楕円体，Wikipedia (英語))
http://en.wikipedia.org/wiki/Ellipsoid

二次曲面 (Wikipedia)
http://ja.wikipedia.org/wiki/%E4%BA%8C%E6%AC%A1% …