質問: 平面の計算方法 3つ以上の座標点(n個)から最小二乗法を用いて,平面の中心座標,xy・xz・yz方向の傾きなどを計算したいのですが,どのように計算したらよろしいでしょうか?
回答: 考え方としてはano.2の参考urlにあるとおりですが大変面倒な計算になります。 うまく収束するかどうか分かりませんが、excelのソルバーを使って、残差2乗和が最小となるような、5個の未知パラメータを求めるのはどうでしょうか。 結果だけ書きますが、n個のデータ点と平面との距離の2乗の和 s が最小となるような次式の a, b, x0, y0, z0 を求めれば、平面が決定されます。 s = σ [ i = 1, 2, ... n ] [ a*( xi - x0 ) + b*( yi - y0 ) + √{ 1 - ( a^2 + b^2 ) } *( zi - z0 ) ]^2 ( xi, yi, zi ) はn個のデータ点(i = 1,2,.. n )です。a、b は平面の法線ベクトル(平面に垂直なベクトル)のx成分とy成分です。( x0, y0 , z0 ) この平面が通る点の座標です。法線ベクトルのz成分がないのは、a^2 + b^2 + c^2 = 1 と規格化したので、c = √{ 1 - ( a^2 + b^2 ) } と自動的に決まるからです。なぜこのような式になるかについては【詳細】に書きました。 未知パラメータ a, b, x0, y0, z0 のうち、x0, y0, z0 は、ano.1さんのコメントにある「中心座標の計算はn個の座標を足してnで割れば」というのは、常に成り立つかどうか分かりません。ただ、初期値として、x0 = σ(xi)、y0 = σ(yi)、z0 = σ(zi) とするのは良いと思います(excelのソルバーは初期値の選択を誤ると収束しないので、もっともらしい値を初期値とすれば収束しやすくなります)。 【詳細】 法線ベクトル(平面に垂直なベクトル)が (a,b,c) で、点(x0,y0,z0) を通る平面の方程式は次式で表されます。 a*( x - x0 )+b*( y - y0 )+c*( z - z0 ) = 0 --- (1) 一方、n 個の点 ( xi, yi, zi ) [ i = 1, 2, ... n ] と平面との距離を li [ i = 1, 2, ... n ] とすれば li^2 = {a*( xi - x0 ) + b*( yi - y0 ) + c*( zi - z0 ) }^2 / ( a^2 + b^2 + c^2 ) ---- (2) となります。したがって、n個の点と平面(1)との距離の2乗の総和を s とすれば、 s = σ [ i = 1, 2, ... n ] ( li^2 ) = σ [ i = 1, 2, ... n ] [ {a*( xi - x0 ) + b*( yi - y0 ) + c*( zi - z0 ) }^2 / ( a^2 + b^2 + c^2 ) ] --- (3) となります。これが最小 となるような、a, b, c, x0, y0, z0 を求めるというのが最終目的です。 s が最小 であるとき、∂s/∂a = ∂s/∂b = ∂s/∂c = ∂s/∂x0 = ∂s/∂y0 = ∂s/∂z0 = 0 が成り立ちます。sを偏微分すると、次の式が得られます(和の範囲はいつも [ i = 1, 2, ... n ] なので省略します)。 a*( b^2+c^2 )*σ( xi - x0 )^2 + b*σ( xi - x0 )*( yi - y0 ) + c*σ( zi - z0 )*( xi - x0 ) = 0 --- (4) b*( c^2+a^2 )*σ( yi - y0 )^2 + a*σ( xi - x0 )*( yi - y0 ) + c*σ( yi - y0 )*( zi - z0 ) = 0 --- (5) c*( a^2+b^2 )*σ( zi - z0 )^2 + b*σ( yi - y0 )*( zi - z0 ) + a*σ( zi - z0 )*( xi - x0 ) = 0 --- (6) a*σ( xi - x0 ) + b*σ( xi - x0 ) + c*σ( zi - z0 ) = 0 --- (7) ()などを展開などして、この計算を続ければano.2の参考urlにある式(4)に行く着くでしょうが、大変複雑な結果となります。 手っ取り早くには、式(3)の s が最小になるような、a, b, c, x0, y0, z0 を求めればいいわけです。 未知パラメータは6個あるように思えますが、実は5個です。a,b,c は比が分かればいいので、 a^2 + b^2 + c^2 = 1 と規格化してしまいます(とすれば |a|≦1、 |b|≦1、 |c|≦1 です)。すると、c^2 = 1 - ( a^2 + b^2 ) ≧ 0 ですから、a,bが決まれば c が決まりますので、c は求めなくていいのです。c の符号は±の2通りありますが、法線ベクトルというのは 、平面の表と裏を区別しないなら、(a,b,c) も (-a,-b,-c) も同じことですから、cの符号は+か-かどちらに決めてしまえばいいのです。したがって c > 0 とすることにすれば、c = √{ 1 - ( a^2 + b^2 ) } です。すると式(3)は s = σ [ i = 1, 2, ... n ] [ a*( xi - x0 ) + b*( yi - y0 ) + √{ 1 - ( a^2 + b^2 ) } *( zi - z0 ) ]^2 ---- (2') となります。 なお、ano1さんの「中心座標の計算はn個の座標を足してnで割れば」というのは、式(7)で、σ( xi - x0 ) =0、 σ( xi - x0 ) = 0 ς( zi - z0 ) = 0 としてx0, y0, z0 を求めるものですが、σ() <> 0でも式(7)が成り立つ場合もあるかと思います。