空間上の測定された点群から最小二乗法により円（中心、直径）を算出する方法

xyz空間で測定された点群（６点）から最小二乗法により円（球ではなく）を検出したく。
円の中心点と直径が分かれば結構です。この点群が既知のax+by+cZ+1の平面上にあることは
分かっているものとします。
よろしくお願いいたします。

質問者からの補足コメント

• horahukisann2017様　回答ありがとうございます。
  
  ６つの点は測定により得た点なので、同一の円周上に完全には乗らないと思われますが、最小二乗法により一番近い円を算出したいと思っています（どうやって最小二乗法を使うのかもわからないのですが）。ただし、６点は同一平面上にあることは分かっているとします。
  horahukisann2017様の暫定的なm、rというのはどのように確定的？なm、rに持っていくのでしょうか。
  よろしくお願いいたします。

  No.1の回答に寄せられた補足コメントです。 補足日時：2017/09/01 08:48

• 実際に測定したわけではないので６つの点は（x1,y1,z1）〜（x6,y6,z6）として、一般的な解法が
  知りたいのですが。

  No.2の回答に寄せられた補足コメントです。 補足日時：2017/09/01 09:28

• 平面ax+by+cz+1=0のa,b,cは既知であるものとします。当然mもこの平面上にあるものとします。
  添付画像のような点群の分布から最小二乗法で円を算出したいのですが。

  
  No.3の回答に寄せられた補足コメントです。 補足日時：2017/09/01 11:15

• そうです。とにかく同一平面上にあることは分かっているものとします。

  No.4の回答に寄せられた補足コメントです。 補足日時：2017/09/01 17:18

• 申し訳ありません。xy平面での話であれば私にもわかるのですが、既知のax+by+cz+1＝0の平面上の点をxy平面へ変換する方法を教えて頂けないでしょうか。

  No.6の回答に寄せられた補足コメントです。 補足日時：2017/09/02 19:19

No.1 回答者： horahukisann2017 回答日時： 2017/09/01 08:15

六つの点は、ある点(たとえば m)を中心とした円周上の点なのでしょうか。

それなら、六つの点の座標を単純に足し合わせた平均が、まず、暫定的な円の中心値 (mx, my) になりそうな気がします。

つぎに、六点個々の、この暫定的な中心点(m)からの距離の値を求めて、その平均値が(r)になるのではないでしょうか。

直径は 2r として求めます。

この暫定的な m, r を用いて、偏差を最小にするように考えるのでしょうか。

No.2 回答者： horahukisann2017 回答日時： 2017/09/01 09:11

測定された六つの点を教えてください。

No.3 回答者： horahukisann2017 回答日時： 2017/09/01 10:17

仮に測定された六つの点は、平面 ax + by + cz + 1 = 0 上の点であり、その平面上に円を描くならば、ｍ もこの平面上にあると考えてよいでしょうか。

平面を決定するには、三点があれば決定できそうなので、六つの点のうち、三点を選べば、a, b, c の候補は計算できそうです。

六つのうちの三点を選ぶ方法は 20 通りあるので、(a, b, c) の候補は 20 組できます。この 20 組のデータの平均の組 (am, bm, cm) は暫定的な平面の式とできるでしょうか。

もしあなたが六つの点を測定して、自分の技量に信頼を与えてみたとき、どの六つの点も同じくらいに確からしいデータとして扱いませんか？

m も r も、二次数値かもしれませんが、一次データ全部をある意味信頼して計算されているものと思います。

できれば、模擬データがあって演習出来たほうが、いろいろ考えることがあって面白いと思うのですが。

No.4 回答者： horahukisann2017 回答日時： 2017/09/01 11:35

これらの 1 から 6 の点は、平面 ax + by + cz + 1 = 0 に投影した点なのでしょうか。

No.5 回答者： horahukisann2017 回答日時： 2017/09/02 05:57

L = Σ{√((xi - a)^2 + (yi - b)^2) - r}^2 として、a, b, r を求めることを考える。

∂L/∂a = 0
∂L/∂b = 0
∂L/∂r = 0

となる条件を考えてみてください。

No.6 回答者： diff3356 回答日時： 2017/09/02 08:19

この定平面を回転によりｘｙ平面と重なるようにします。

６点が同一平面上にのっていることから、こうすることにより「２次元」の問題となります。以下、変換後の話です。
与えられた６点からの距離の平方和(の平均)が最小となるように未知の定点P(x, y)を決定します。
(1/6)*Σ[i=1~6]{(x-x[i])^2+(y-y[i])^2}=f(x, y) とすると、
(x, y)=((1/6)*Σx[i], (1/6)*Σy[i]) のときf(x, y)が最小であることがすぐにわかります。
このとき、「円」の半径Rは、
R^2=(1/6)*{(-1/6)*{Σx[i]}^2+Σ{x[i]}^2 - (1/6)*{Σy[i]}^2+Σ{y[i]}^2}.
です。

No.7 回答者： stomachman 回答日時： 2017/09/07 17:14

平面の話のようなので、「x-yの直交座標で点の座標データ(X[i], Y[i]) (i=1,2,…,N)が与えられている。

これらに最もフィットするような円周の中心座標(P,Q)と半径rが欲しい。ただし、点はどれも概ね円周上にあって、ばらついている（円周の1〜2箇所の付近に集まってはいない）」という前提で考えてみましょう。

[1] もし点がどれも円周上にあるのなら
(1)　(X[i]-P)^2 + (Y[i]-Q)^2 = r^2
を満たしているはずですが、実測値には誤差がある、という話です。
　まずは、中心の座標(P,Q)の大雑把な近似値を出します。たとえば、与えられた点全部の重心だとか、x,y座標それぞれの最大値と最小値の中点とか。ま、かなり大雑把でいいんです。
　それを(P0, Q0)としまして、
(2)　x[i] = X[i]-P0, y[i] = Y[i]-Q0, p = P-P0, q = Q-Q0
と変換しておきます。そして
(3)　z = r^2 - p^2 - q^2
とおいて、(1)を
(1')　(x[i]^2 + y[i]^2) - 2 px[i] - 2qy[i] - z = ε[i]
（ε[i]は誤差を表す）と書き換える。
　すると、左辺第一項はデータから決まり、これは定数。左辺の残りの項は未知数p, q, zを持つ一次式ですから、尺度
(4)　E = ∑(ε[i])^2
を最小化するp, q, zを計算するのは線形最小二乗法の問題です。従って鶴亀算で簡単に解けて、これでp, q, rの（かなり良い）近似値が得られます。
　なお、得られた近似値を改めて(P0,Q0)にして[1]を繰り返せば、（(P0,Q0)が正解に近いお蔭で）もっと良い近似が得られます。

[2] さてここからの話は、「最もフィットする」をどう定義するかに依存します。「どう定義するか」とは、すなわち、
(1)　(X[i]-P)^2 + (Y[i]-Q)^2 = r^2
のどこに誤差を表す項を組み込むかです。たとえば、「円周と点との最短距離の２乗和を最小にしたい」というのであれば、誤差δ[i]は
(5)　(X[i]-P)^2 + (Y[i]-Q)^2 = (r+δ[i])^2
と表されることになり、尺度
(6)　D = ∑(δ[i])^2
を最小化するというのは非線形最小二乗問題です。非線形最小二乗問題は数値的に解くしかなく、一般にその計算方法は

　　まず (P,Q,r)の近似値（出発値）を（なんとか）出して、
　　Dが小さくなるように（なんとか）P, Q, rを改良することを繰り返す。
　　Dがほとんど改良されなくなったら、収束したとみなして終わる。

ということです。
　さて、もし最初に与える出発値の近似精度が良ければ、（Dにかぎらず、「最もフィットする」というココロに沿うような尺度ならどんな尺度であっても、）かなり速く収束する。そして、良い出発値を決めるのには[1]の近似計算法が利用できます。
　また「（なんとか）P, Q, rを改良する」の部分のやり方は様々工夫されているのですが、理論的に素直な「ニュートン法」もしくは「ガウス・ニュートン法」を用いると速いでしょう。というのは、これらの方法は出発値の近似精度が悪い場合には発散してしまうおそれがあって、一般には発散を抑える工夫が必要なのですが、出発値の近似精度が良いときには発散の心配は要りません。だからシンプルなプログラムで計算できるというわけです。

　なお（線形・非線形の）最小二乗法の実践的な教科書としては、毎度
「最小二乗法による実験データ解析」 (UP応用数学選書 7)
をお薦めしています。