球の中心と半径の求め方

カテゴリーでは数学部門だと思うのですが…

現在、ある球体を測定してその物の球の中心と半径を求めようと思っています。

球体の表面をある機械で測定しています。それぞれの機械から計測された３次元の形状データとしては１００００点ほど有ります。
そこで、４点をピックアップして球の公式に入れて求めると、全てがかなりばらつきのあるデータが出てきます。

これではどのあたりに球の中心があり半径がどうなのか分かりません。

できれば、統計的、数学的に中心と半径は信頼度が高く求める方法があれば教えて下さい。

ちなみに、４点からのデータから中心と半径を１００回求め、平均してもばらつきがひどくてこの方法は使わないようにしています。

そこで、１０点をピックアップして、最小二乗法から中心と半径を１００回求め、平均する方法をしたりしました。
まだ、この方法が有用かなと思いましたが…
データのばらつきが少ないからいいのかなと思いました。

大変理解していただくには難しい内容かもしれません。随時、応えさせていただきますので、色々なご意見宜しくお願いいたします。

No.4 ベストアンサー 回答者： boobee0125 回答日時： 2008/11/29 09:43

#2です。

本件は球の中心位置と半径を正確に求めるのが目的かと思っていましたが、機械の測定誤差を評価・検証する場合には別の考慮も必要になります。

一般に計測装置の誤差にはノンリニアリティー、スパン誤差、オフセット誤差、バラツキなどがありますが、バラツキ以外は装置によっては補正をかけることも可能です。

球を使った場合はノンリニアリティー以外の評価が可能で、半径の誤差がスパン誤差に、中心位置の誤差がオフセット誤差に、また #2手順(2)の rms がバラツキの標準偏差 σ に対応します。また半径の異なる複数の球を使えばノンリニアリティーの評価も可能になります。測定のバラツキが他の誤差評価に影響を与えないようにするには測定点数 n を大きく取る必要があります。n の目安は

rms/√n ＜＜ 要求されるスパン誤差、オフセット誤差評価の標準偏差

で求められます。＜＜ は３～１０程度の比と考えてください。
　
>この教えていただいたニュートン法でのループを減らし、または誤差二乗和（二乗平均誤差）を利用すれば今回の私が求めたいと考えていることに近づけるのでしょうか？

n = 10000程度のデータでしたらパソコンでの計算ではループ回数＝15で異常データ排除ループを入れてもおそらく数秒で終る程度ですからループ回数を減らす必要はなさそうですし、またバラツキが不明の場合は#2手順(2)の終了判定は使わないでループ回数を固定するほうがよいでしょう。装置の性格によっては異常データは排除しない方が良い場合もありますし、異常データを排除する場合は異常データを含めた rms も評価の対象にすべきです。

この回答へのお礼

お返事遅くなり大変申し訳ございません。

大変参考になるご意見をいただけ、感謝をしております。
今までは違う事を研究していたため、専門ではなかった測定器の精度検定で、上司からもう少し違う側面から専門的にと言われ、どのようにしたら良いのか迷っていました。

今、自分なりにプログラムを組んでシュミレーションをしたり、色々やっています。

しかし、今の所それほど上手く言っている感じ（あくまでも自分の感覚ですが…）がしないため、もう少し頑張ってみたいと思います。

今後も多分、分からないことが出ると思うので、そのときはご教授いただければ幸いと考えております。
返事が大変遅くなったうえに、今後のお願いまでさせていただき申し訳有りませんでしたが、今後も宜しくお願い申し上げます。

お礼日時：2008/12/04 00:08

No.3 回答者： arrysthmia 回答日時： 2008/11/28 23:43

No.1 の者です。

基本方針はアレでよいと思うのですが、

＞ 結局、１００００点の重心を求めるだけですから、計算も簡単です。

の部分には、重大な思い違いがありました。スミマセン。
No.2 さんの手順を参考にして下さい。

No.2 回答者： boobee0125 回答日時： 2008/11/27 23:23

基本的に #1 の方と同じ意見です。

追加質問のアルゴリズムに関して原理をスキップして手順のみ示します。以下ベクトルを ~ 、転置を ' で表します。

[定義]
計測データ数を n とし、下記のベクトル、行列を定義する。

●計測データベクトル ui~ = (xi,yi,zi)'
●パラメータベクトル p~ = (a,b,c,R)'
●式誤差ベクトル e~ = (E0,E1,...,En-1)
ここに Ei = F(ui~,p~)
= (a-xi)^2 + (b-yi)^2 + (c-zi)^2 - R^2

●偏微分行列 G = {gij} （サイズは n x 4）

gi0 = ∂F(ui~,p~)/∂a = 2*(a-xi)
gi1 = ∂F(ui~,p~)/∂b = 2*(b-yi)
gi2 = ∂F(ui~,p~)/∂c = 2*(c-zi)
gi3 = ∂F(ui~,p~)/∂R = -2*R


[ニュートン法により最良パラメータを求める手順]

(1) パラメータベクトル p~ に初期値として（何らかの方法で求めた）近似値を与える。

(2) 誤差ベクトル e~ を作り、誤差の自乗平均誤差

rms = √(ΣEi^2/n)

が指定値以下であれば現パラメータベクトル p~ を最終解として処理を終了する。

(3) 偏微分行列 G を作る。

(4) G から擬似逆行列 inv(G'*G)*G' を作る。

(5) パラメータベクトル p~ を下記のように改良する。

p~ = p~ - k*inv(G'*G)*G'*e~

ここに k はいわゆる加速係数で 0＜k≦1 の範囲に選ぶ。この k が小さすぎると収束が遅くなり、大きすぎると収束が不安定になったり発散することがある。

(6) 手順(2)にジャンプする。　

[補足１]
実際のプログラムでは手順(3)～(5)の inv(G'G)*G'*e~ の計算を下記のようにしてメモリ量を節約できる。Σは for ループによる i = 0～n-1の積算を意味する。G'*e~、G'*G に対応する配列をそれぞれ GTE[4]、GTG[4][4]とするとき

GTE[0] = Σ{gi0*Ei} = 2*Σ{(a-xi)*Ei}
GTE[1] = Σ{gi1*Ei} = 2*Σ{(b-yi)*Ei}
GTE[2] = Σ{gi2*Ei} = 2*Σ{(c-zi)*Ei}
GTE[3] = Σ{gi3*Ei} = -2*R*Σ{Ei}

GTG[j][k]=Σ{gij*gik}

[補足２]
手順(2)の終了判断を行わずにループ回数を固定する方が高い精度が得られる場合がある。今回の問題の場合では加速係数=0.5、ループ回数=15回程度で十分な収束が可能と思われる。

[補足３]
計測データの中に異常データが含まれる場合には、上記の手順完了後に

|Ei|＞k1・rms

となるデータを全て排除し、残りのデータのみを用いて再度(1)～(6)を実行する。この k1 は異常データの発生状況にもよるが 2～3 程度にとる。この操作を３～４回繰り返せば異常データを排除した精度の良い計測が可能となる。

この回答への補足

ありがとうございます。
僕の場合は単純な数学で物事を解決しようとしていました。
ですから、このような方法を教えていただけると、目からうろこがと言う新鮮な気持ちに慣れます。

少しずつですが、プログラムを組んでいます。

さて、今回の皆様から頂いた意見を参考にデータ処理について考えていましたが、このような方法を使うと言うことは、つまり１００００点のデータからより良いデータを抽出し、中心と半径が一定の値に収束していくと思われるのですが、如何でしょうか？

もう少しこの質問の背景を書かせていただきますと、
私は複数の機械の精度判定をしています。
そこで、考え付いたのが同じ球を用いて、測定をしその結果によって、機械の精度を測りたい事です。

基準となる機械を一台所有していて、新しく作成した機械がこの機械の性能まで近づいているかを判断したいと思っています。

考え付いたのが、単純な物を測定して相互比較をするということです。
球の半径と中心を求めて、それを統計的に比較するだけの事です。

問題となったのが、球の中心と半径を求める最良の方法です。
単純に高校で習った４点から球の半径を求めると、かなりのばらつきが生じてきます。
基準となる測定機でもばらつきが生じたので、そこで今までこの機械でやってきた方法を用いました。
つまり、１０点を任意に選択して、球の中心と半径を最小二乗法にて求める方法です。
しかし、あまり応用数学の知能を持ち合わせていないためこれ以上の解決方法を見出すことが出来ませんでした。


私としては、得られたデータを全てを活用してそこから精度判定できたらと考えています。

教えていただいた。ニュートン法になると精度良い中心と半径が得られるのはわかりますが、測定点を判別し、抽出した結果で得られるものだと解釈しました。

これまでの説明不足のために皆様に多大なるご迷惑と多くの時間を私のために割いていただいたことに深く感謝し、そしてお詫び申し上げます。

この教えていただいたニュートン法でのループを減らし、または誤差二乗和（二乗平均誤差）を利用すれば今回の私が求めたいと考えていることに近づけるのでしょうか？

それとも、上記の説明を踏まえてそれ以外の方法があれば教えていただければ幸いと考えております。

補足日時：2008/11/29 01:08

No.1 回答者： arrysthmia 回答日時： 2008/11/27 01:37

なぜ、１００００点のデータから、１０点×１００回だけ取り出すのでしょう？

最適性の定義が最小二乗法で良いのなら、１００００点全部使って
半径の誤差二乗和が最小になるような中心を選んでみては、どうでしょう。
点数が多い程、中心極限定理が効いてくるし、
結局、１００００点の重心を求めるだけですから、計算も簡単です。

１００００点の中に、とんでもない外れ値が含まれていて、値がズレる
ということならば、(かなり場当たり的ですが、)
上の方法で、仮の中心と半径を一旦求め、仮の半径を参考に
仮の中心からの距離に閾値を設けてデータを選別し、残ったデータでもう一度
最小二乗法を行う　という手もあるでしょう。

この回答への補足

ありがとうございます。
とても、興味深いお返事をいただき、感謝しております。

まず、１０点を１００回取り出すのは、昔使用していた機械で球の中心と半径を求める時に、表面の形状データ１０点を任意に選択して、中心と半径を求めていたためです。
この機械の中での計算も最小二乗法を用いていました。

１００００点あるのは新しい機械で球の表面を測定したら、この点数だけデータを採取できました。

つまり、今回は計測して得られた１００００点のデータからPC上で任意に１０点を選択して、最小二乗法をかけていました。

さて、ここで質問が…
”半径の誤差二乗和が最小になるような中心を選んでみては…”
とありますが、誤差２乗和とはどのように計算式を組んだらよろしいのでしょうか。
プログラムとしてはVisual　C++　6.0をWindows環境下で組んでいます。
できれば、アルゴリズム的な事を教えていただければ幸いです。

また、この誤差２乗和と中心極限定理との兼ね合いと言うか関係がどのようになっているのかも、少々分かり辛いので、もう少し解説していただければ幸いです。

補足日時：2008/11/27 02:21