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

조회 수: 28 (최근 30일)
Egoshi
Egoshi 2021년 9월 28일
댓글: Egoshi 2021년 10월 5일
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個の座標をプロットしたもので,おおよそ円形に配置しています.この散布図上に円をフィッティングしたいのですがどうすればよいでしょうか.

채택된 답변

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

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Annotations에 대해 자세히 알아보기

제품


릴리스

R2019a

Community Treasure Hunt

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

Start Hunting!