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

27 vues (au cours des 30 derniers jours)
Egoshi
Egoshi le 28 Sep 2021
Commenté : Egoshi le 5 Oct 2021
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]); %データの縦横比
別のコマンドで入手した12個の座標をプロットしたもので,おおよそ円形に配置しています.この散布図上に円をフィッティングしたいのですがどうすればよいでしょうか.

Réponse acceptée

Akira Agata
Akira Agata le 28 Sep 2021
いろいろな解決法があるかと思いますが、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
  3 commentaires
Akira Agata
Akira Agata le 1 Oct 2021
ご連絡ありがとうございます。
うまく円の fitting ができたとのことで安心しました。
>Agata様がこのフィッティングの方法を思いつかれたのは経験からですか?
簡単なようで、なかなか難しいご質問ですね。
経験も多少あるかもしれませんが、どちらかというと「ダメもとでも使えそうなアイデアをいくつ思いつくか」、だと思います。
たとえば今回のケースでの私の思考プロセス(...というほどのものでもないですが)は、
  1. すべての点の (x,y) 座標の平均を円の中心として、そこから各点までの距離の平均を半径とすれば良いのでは? ➡ 点の数が少ないので正確な結果にならなさそうなのでボツ。
  2. 各点からの距離の分散を目的関数として、fminsearch で最小化すると良いのでは? ➡ 目的関数を書くのに手間がかかりそう。もう少し楽に求めたいのでボツ。
  3. 円の式をあらわす係数を推定すれば良いのでは? ➡ 今回の回答
といった流れでした。
玉石混交でも良いので、こうした「使えるかもしれない」方法をなるべく多く出すことが課題解決につながるように思います。
Egoshi
Egoshi le 5 Oct 2021
Agata様,詳しく説明していただきありがとうございます.
まずはいろいろな可能性を考えて結果を導いていきたいと思います.
非常に勉強になりました.

Connectez-vous pour commenter.

Plus de réponses (0)

Produits


Version

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!