散布図上で円のフィッティング

x1 = 172.0974;
y1 = -386.6972;

x2 = 203.1669;
y2 = -236.628;

x3 = 305.764;
y3 = -120.6805;

x4 = 456.0738;
y4 = -72.9798;

x5 = 614.1931;
y5 = -105.4311;

x6 = 736.6825;
y6 = -215.3598;

x7 = 789.8805;
y7 = -372.3829;

x8 = 757.1624;
y8 = -536.2989;

x9 = 650.8597;
y9 = -665.3943;

x10 = 492.0848;
y10 = -716.8889;

x11 = 326.1053;
y11 = -691.1757;

x12 = 204.3204;
y12 = -576.2446;
%%
x = [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12];
y = [y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12];
scatter(x,y,'filled');

axis([100 900 -800 0]);  %軸の範囲設定
grid on;  %軸の目盛り線
pbaspect([1 1 1]);  %データの縦横比

別のコマンドで入手した１２個の座標をプロットしたもので，おおよそ円形に配置しています．この散布図上に円をフィッティングしたいのですがどうすればよいでしょうか．

Answer 1

[Akira Agata]

いろいろな解決法があるかと思いますが、fminsearch を使って円をあらわす式の係数を推定するのはいかがでしょうか？

以下はその一例です。

x = [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12];
y = [y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12];

% 最小化したい関数 (円と(x,y)の二乗誤差の総和)
fcn = @(a) sum(((x-a(1)).^2 + (y-a(2)).^2 - a(3).^2).^2);

% 係数 a(1)~a(3) の初期値 (図から目視で見積もったおおよその値)
a0 = [500, -400, 400];

% 関数の出力が最小となる係数を計算
aOpt = fminsearch(fcn, a0);

% フィッティングした円の (x,y) 座標
theta = linspace(0,2*pi);
xFit = aOpt(3)*sin(theta) + aOpt(1);
yFit = aOpt(3)*cos(theta) + aOpt(2);

% 可視化
figure
scatter(x,y,'filled');
hold on
plot(xFit,yFit)
axis([100 900 -800 0]);  %軸の範囲設定
grid on;  %軸の目盛り線
pbaspect([1 1 1]);  %データの縦横比
box on

[Egoshi]

所望していたような円のフィッティングをすることができました，ありがとうございます．

余談ですが，Agata様がこのフィッティングの方法を思いつかれたのは経験からですか？

今後の研究活動に生かせればと思っているので，よろしければお答え頂けると幸いです．

[Akira Agata]

@Yosjomune Kuwajima 様

ご連絡ありがとうございます。

うまく円の fitting ができたとのことで安心しました。

>Agata様がこのフィッティングの方法を思いつかれたのは経験からですか？

簡単なようで、なかなか難しいご質問ですね。

経験も多少あるかもしれませんが、どちらかというと「ダメもとでも使えそうなアイデアをいくつ思いつくか」、だと思います。

たとえば今回のケースでの私の思考プロセス（...というほどのものでもないですが）は、

1. すべての点の (x,y) 座標の平均を円の中心として、そこから各点までの距離の平均を半径とすれば良いのでは？　➡　点の数が少ないので正確な結果にならなさそうなのでボツ。
2. 各点からの距離の分散を目的関数として、fminsearch で最小化すると良いのでは？　➡　目的関数を書くのに手間がかかりそう。もう少し楽に求めたいのでボツ。
3. 円の式をあらわす係数を推定すれば良いのでは？　➡　今回の回答

といった流れでした。

玉石混交でも良いので、こうした「使えるかもしれない」方法をなるべく多く出すことが課題解決につながるように思います。

[Egoshi]

Agata様，詳しく説明していただきありがとうございます．

まずはいろいろな可能性を考えて結果を導いていきたいと思います．

非常に勉強になりました．