質問: 3次元の近似直線 こんにちは。2次元で実験データなどの点列から近似直線を求めるのは、最小二乗法の基本問題ですが、3次元の点群から直線の方程式(x-x0)/a=(y-y0)/b=(z-z0)/cを求めるにはどんなアルゴリズムを使いますか?スマートな方法があれば教えていただけたら幸いです。よろしくお願いします。
回答: 3次元空間の曲面ではなく、直線に乗ると仰るのだから、 (1) x, y, zのどれかを与えて、残りの2つを推定する問題。 (2) と直線との距離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) 点と直線との最短距離を求める。 直線上の任意の点と点の距離を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 = r x + p y y q z z r と表せます。rは3×3の行列で rの転置をr'とすると rr' = r' r = 単位行列 となる行列です。各点をこの変換で<x[i],y[i],z[i]>に写したとすると、 直線、すなわちz軸との最短距離はx[i]^2 + y[i]^2ですから、他のどんな回転・平行移動の仕方に比べても u=σ(x[i]^2 + y[i]^2) (i=1,2,....,n) が最小になっている筈で、しかも s=u です。 さて、uはz[i]の値とは無関係ですからz[i]を求める必要はない。さらに座標系をz軸の周りで回転してもuは変化しません。従って、 x = r x + p y y q z r =p(α)q(β) p(α)=cosα 0 -sinα 0 1 0 q(β)= 1 0 0 0 cosβ -sinβ 0 sinβ cosβ とすれば良いのです。展開すれば x[i] = x[i]cosα-y[i]sinαsinβ-z[i]sinαcosβ+p y[i] = y[i]cosβ-z[i]sinβ+q ですね。 ここでα、β、p、qを決めたい訳です。 始めに(1)の問題を解けば、α、β、p、qの大体の値を求めることができます。これを使ってu(α,β,p,q)を計算します。 それから、u(α,β,p,q)が小さくなるようにα、β、p、qをちょっとずつ改良して行けば良いでしょう。これには微小な角度δα、δβを使って、 p(δα)=cosδα 0 -sinδα 0 1 0 sinδα 0 cosδα q(δβ)= 1 0 0 0 cosδβ -sinδβ 0 sinδβ cosδβ を作り、p(α)、q(α)にそれぞれ掛け算すれば良い。 p(α+δα)=p(δα)p(α) q(β+δβ)=q(δβ)q(β) だからです。さらにここで、δα、δβは微小だから、 cosδα≒1、cosδβ≒1、sinδα≒δα、sinδβ≒δβ (δα)^2≒0、(δβ)^2≒0、δαδβ≒0 という近似をしても構わないでしょう。 この近似を利用すると計算は一層簡単になり、uを最小にするようにδα、δβ、p、qを求める問題は線形最小二乗法(一次式の最小二乗法)になってしまい、簡単に解けます。 それを解いてから、真面目にp(α)、q(α)を計算しなおし、また線形最小二乗法を解く。これを収束するまで繰り返せば良いのです。 なお、stomachmanは計算間違いの常習犯ですから、チェックは慎重に。