３次元の近似直線

質問

こんにちは。２次元で実験データなどの点列から近似直線を求めるのは、最小二乗法の基本問題ですが、３次元の点群から直線の方程式(x-x0)/a=(y-y0)/b=(z-z0)/cを求めるにはどんなアルゴリズムを使いますか？スマートな方法があれば教えていただけたら幸いです。よろしくお願いします。

回答No.2

stomachman · ベストアンサー

３次元空間の曲面ではなく、直線に乗ると仰るのだから、
(1) x, y, zのどれかを与えて、残りの２つを推定する問題。
(2) <x[i],y[i],z[i]>と直線との距離d[i]の二乗和が最小になる直線を求める問題。
と分類すべきでしょう。

(1)の場合は、たとえばzを与えてx,yを求めたいのであれば、
・zからxを求める問題。
・zからyを求める問題。
の二つを別々に解けばおしまいです。
それぞれの解は（x=Az+B, yは任意）という平面と、（y=Cz+D, xは任意）という平面を定めますから、この二つの平面の交線が、求める直線ということですね。

(2)の場合はやっかいです。
[1]ちょっと手抜きしながらも、まともにやってみましょう。
(i) 直線をどう表すか。
ご質問の式を見ると、この直線はx軸、y軸、z軸のどれとも平行でも垂直でもないことが仮定されています。
ですから、zをパラメータとして
x=az+c
y=bz+d
と書いても良いでしょう。a,b,c,dが決められれば良い訳です。
(ii) 点<p,q,r>と直線との最短距離を求める。
直線上の任意の点<az+c,bz+d, z>と点<p,q,r>の距離をdとすると
d^2 = (az+c-p)^2+(bz+d-q)^2+(z-r)^2
= (az)^2+2az(c-p)+(c-p)^2+(bz)^2+2bz(d-q)+(d-q)^2+z^2-2rz+r^2
です。これが最小になるzを求めると、
0=∂(d^2)/∂z = 2(az+c-p)a+2(bz+d-q)b+2(z-r)
ゆえに
z=(ap+bq+r-ac-bd)/(a^2+b^2+1)
であって、このときの最短距離h(p,q,r)は
h(p,q,r)^2 = (ap+bq+r-ac-bd)^2/(a^2+b^2+1)+2(ap+bq+r-ac-bd)(ac-ap+bd-bq-r)/(a^2+b^2+1)+(c-p)^2+(d-q)^2+r^2
わあ、とんでもないですね。
(iii) じゃあ、直線を求めるには？
S=Σ(h(x[i],y[i],z[i]))^2　　（i=1,2,...,N)
を最小化するには
∂S/∂a = 0
∂S/∂b = 0
∂S/∂c = 0
∂S/∂d = 0
を解く必要があります。言い換えれば
∂(h(x[i],y[i],z[i]))/∂a
∂(h(x[i],y[i],z[i]))/∂b
∂(h(x[i],y[i],z[i]))/∂c
∂(h(x[i],y[i],z[i]))/∂d
を求めておいて
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂a)=0
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂b)=0
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂c)=0
Σh(x[i],y[i],z[i]) (∂(h(x[i],y[i],z[i]))/∂d)=0
という連立方程式を解くことになります。
これがa,b,c,dについて非線形である（一次式でない）ことは言うまでもありません。一筋縄では行かず、反復計算で徐々に収束させていくしかありません。

[2]手抜き
もうすこし手抜きの方法を考えてみましょう。
この座標系を回転・平行移動した座標系をX-Y-Zとします。そして、求めたい直線がZ軸と一致するようにしたとします。回転と平行移動は行列を使って
Ｘ　＝　Ｒ　ｘ　＋　ｐ
Ｙ　　　　　ｙ　　　ｑ
Ｚ　　　　　ｚ　　　ｒ
と表せます。Ｒは３×３の行列で Ｒの転置をＲ'とすると ＲＲ' = Ｒ' Ｒ = 単位行列
となる行列です。各点<x[i],y[i],z[i]>をこの変換で<Ｘ[i],Ｙ[i],Ｚ[i]>に写したとすると、
直線、すなわちＺ軸との最短距離はＸ[i]^2 + Ｙ[i]^2ですから、他のどんな回転・平行移動の仕方に比べても
U=Σ(Ｘ[i]^2 + Ｙ[i]^2)　　（i=1,2,....,N)
が最小になっている筈で、しかも
S=U
です。
　さて、UはＺ[i]の値とは無関係ですからＺ[i]を求める必要はない。さらに座標系をＺ軸の周りで回転してもUは変化しません。従って、
Ｘ　＝　Ｒ　ｘ　＋　ｐ
Ｙ　　　　　ｙ　　　ｑ
　　　　　　ｚ
Ｒ ＝Ｐ（α）Ｑ（β）
Ｐ（α）＝ｃｏｓα　　０　　－ｓｉｎα
　　　　　　　０　　　１　　　　０
Ｑ（β）＝　１　　０　　　　　０
　　　　　　０　ｃｏｓβ　－ｓｉｎβ
　　　　　　０　ｓｉｎβ　　ｃｏｓβ
とすれば良いのです。展開すれば
Ｘ[i] = x[i]ｃｏｓα－y[i]ｓｉｎαｓｉｎβ－z[i]ｓｉｎαｃｏｓβ＋ｐ
Ｙ[i] = y[i]ｃｏｓβ－z[i]ｓｉｎβ＋ｑ
ですね。
ここでα、β、ｐ、ｑを決めたい訳です。

　始めに（１）の問題を解けば、α、β、ｐ、ｑの大体の値を求めることができます。これを使ってU(α,β,ｐ,ｑ)を計算します。
　それから、U(α,β,ｐ,ｑ)が小さくなるようにα、β、ｐ、ｑをちょっとずつ改良して行けば良いでしょう。これには微小な角度Δα、Δβを使って、
Ｐ（Δα）＝ｃｏｓΔα　　０　　－ｓｉｎΔα
　　　　　　　　０　　　　１　　　　０
　　　　　　ｓｉｎΔα　　０　　　ｃｏｓΔα
Ｑ（Δβ）＝　１　　０　　　　　　０
　　　　　　　０　ｃｏｓΔβ　－ｓｉｎΔβ
　　　　　　　０　ｓｉｎΔβ　　ｃｏｓΔβ
を作り、Ｐ（α）、Ｑ（α）にそれぞれ掛け算すれば良い。
Ｐ（α＋Δα）＝Ｐ（Δα）Ｐ（α）
Ｑ（β＋Δβ）＝Ｑ（Δβ）Ｑ（β）
だからです。さらにここで、Δα、Δβは微小だから、
ｃｏｓΔα≒１、ｃｏｓΔβ≒１、ｓｉｎΔα≒Δα、ｓｉｎΔβ≒Δβ
(Δα)^2≒０、(Δβ)^2≒０、ΔαΔβ≒０
という近似をしても構わないでしょう。
この近似を利用すると計算は一層簡単になり、Uを最小にするようにΔα、Δβ、ｐ、ｑを求める問題は線形最小二乗法（一次式の最小二乗法）になってしまい、簡単に解けます。
それを解いてから、真面目にＰ（α）、Ｑ（α）を計算しなおし、また線形最小二乗法を解く。これを収束するまで繰り返せば良いのです。

なお、stomachmanは計算間違いの常習犯ですから、チェックは慎重に。

お礼

ありがとうございました。非常に理解しやすく、活用することができました。また、内容も大変興味深いので、仕事の合間を見つけて非線形最小二乗法も勉強しようと思います。今後とも、よろしくお願いいたします。よいお年をお迎えください。