最小二乗円

[noel]

円の輪郭を示す、ランダムな座標群(Xi,Yi)がある時に、その座標を用いて、最小二乗法で、円の半径や中心を求める方法はありませんか？
よろしくお願いいたします。

[box]

座標群(X<sub>i</sub>,Y<sub>i</sub>)は、必ず円周上にあるのですか？
また、どうしても最小二乗法を使わないとダメですか？

もし、3点以上から成る座標群が必ず円周上にあって、
最小二乗法を使わなくてもよいのであれば、
座標群から任意の2点を2組取り出し、各2点の垂直二等分線の交点から
中心座標が求まり、座標群の任意の点と中心座標との距離が円の半径となります。
ご存じだとは思いますけれども。

最小二乗法を今回の問題にどう適用すればよいかは、よくわかりません。
他の方からの回答をお待ちください。

[noel]

比類無きプログラマーさん、ありがとうございます。

必ず円周上にあるわけではなく、円周周辺にランダムにあります。
しかも、一部角度の座標は、欠落している場合もあります。
そこでそれらの点群からベストフィットする円を求めたいのです。
なので最小二乗法が良いのかと？思っています。
特に最小二乗法にこだわるわけではありません、ベストフィットする円を求めたいだけなのです。
よろしくお願いいたします。

[noel]

すみません、お名前を間違えてしましました。
boxさん、ありがとうございます。

[box]

> 必ず円周上にあるわけではなく、円周周辺にランダムにあります。

それはそうですよね。
円の周辺にある、ばらつきを持った点の集まりをもとに、
いちばんフィットする円を見つけ出すのが最小二乗法の目的ですものね。

今回の場合はXもYも2次式ですからね。
どういう風にすれば最小二乗法（あるいはその他の方法）を使って
答えが求まるか、ちょっと見当がつかないです。
申し訳ありません。

[組木紙織]

携帯から失礼します。
詳しく書けませんが、
円を与えて、その円との差を考えてみてください。
多変量解析の辺りが参考になると思います。

[管理人]

全ての組み合わせで２点の直線を作り、その出てきた交点の中心を取ったのではだめなんでしょうかね・・。
すみません、よくわかりませんが。

[noel]

回答ありがとうございます。
円周上の座標が正確な場所と微妙に不正確な座標があり、正確な場所の数の方が多いのですが、不正確な座標に足を引っ張られるケースが多いのです。
ですので、最小二乗法などを使うとベストフィットするのでは？？？と考えたしだいであります。
きっと組木紙織さまの言われる「円を与えて、その円との差を考えてみてください。」と考え方は同じになるのかな？なんて思っています。
でも 円の最小二乗法のやり方を思いつきません…

[組木紙織]

少し考えてみたのですが、かなり難しいです。
いや考え方は難しくはないのですが、式がかなり複雑になるので計算が大変でした。

一番単純な
「中心が原点にある場合」
をやってみます。

与えられた座標を(x_i,y_i)とします。
中心からその点までの距離をd',円の半径をr,円から点までの距離をdとします。

d' = sqrt((x_i)^2 + (y_i)^2)
d =|d' - r|

よって、
d^2 = (x_i)^2 + (y_i)^2 - 2r*sqrt((x_i)^2 + (y_i)^2) + r^2

総和を考えて
sig(d^2) = sig((x_i)^2 + (y_i)^2 - 2r*sqrt((x_i)^2 + (y_i)^2) + r^2)

この式の最小値を取るrの値が、この場合の答えです。

微分して考えるのが一番楽なのかな。
これは原点に中心がある場合なのでrで微分すれば話を進められるのですが、
中心の場所が分からない場合はそれぞれの変数で偏微分して考えないと解けないと思う。

大学(数学科？)レベルの問題ですね。

今はこのくらいまでですが、自分が気になるのでもう少し考えてみます。

[しっぽ]

こんなのでどうでしょうか。

１．座標を配列にし、解析しやすいように変形する。
　座標の配列をX座標でソートする
　座標をX座標で３領域（A,B,C）に分割する

２．３つの座標から円の中心座標を求める。
　座標の３領域から１点づつ選び出す（A[j], B[j], C[j]）

３．おおまかな円の中心座標を探す。
　２の操作を（座標数 / 3）点行い、中心座標の候補を配列に保存する
　上記配列のX座標のメディアン値を円の中心座標の探索初期値とする
　Y座標にについても同様の処理をする

４．中心座標の探索初期値に対して、２乗誤差が最も小さくなる半径を
　半径の探索初期値とする

５．あとは
　探索初期値の中心座標を初期の否探査領域とし、非探査領域を囲むく領域と
　比較する



言葉で書くと何のこっちゃ！言葉は難しいです。コードで語りたいのですが

▼【エラー】 申し訳ありません。メッセージが長すぎます。

だそうです。