複数円を囲む凸包

3 次查看（过去 30 天）
YA
YA 2023-4-10
评论： YA 2023-4-16
複数の円を囲む凸包を求める方法についての質問です。
複数の異なる半径を持つ円があるときに、それを囲む凸包を求めたいです。
半径が同一の円であれば、中心を囲む凸包をconvhullで求めた後に、半径をpolybufferでバッファーとすれば求められましたが、
半径が異なる場合はその方法が使えません。
解決策やアドバイスをお持ちの方がいましたら、ご教授いただきたいです。
よろしくお願いいたします。
  6 个评论
显示 4更早的评论
Atsushi Ueno
Atsushi Ueno 2023-4-13
移动：Atsushi Ueno 2023-4-13
Ran in:
>接線を求めて地道に求めていくしかなさそうですね
もう作ってしまったかもしれませんが...あくまで convhull 関数を生かす事に拘りました。一応要件を満たす動きになったと思います。プログラムは回りくどいです。
  • 凸包を求める
  • 接線を結ぶ点群に絞る(円周上の点群を削除)
  • 接線同士の交点を求め点群に追加・凸包のindexに挿入
th = (0:pi/1000:2*pi)';
x = []; y = [];
for n = 1:20
r = randi(100,1,1) + 50;
x = [x; r * cos(th) + randi(1000,1,1)]; % 複数の円状に並んだ点群
y = [y; r * sin(th) + randi(1000,1,1)];
end
K = convhull(x,y); % 複数の円状に並んだ点群を囲む凸包のindexを求める
dst = hypot(diff(x(K)),diff(y(K))); % 凸包の点間距離を求める
idx = dst > mean(dst); % 円周上の点群距離は平均値で、接線のみ比較的長い
K = K([idx; false] | [false; idx]); % 接線の始点と終点のみ選り分ける
temp = [K circshift(K,-1) circshift(K,-2) circshift(K,-3)]';
K = [reshape(K,2,[]); size(x,1)+1:size(x,1)+size(temp,2)/2];
K = K(:); % ３つおきに交点のindexを挿入
for k = temp(:,1:2:end) % 4点ずつ並べて2つおきに走査する
[x(end+1),y(end+1)] = intersection(x(k),y(k)); % 交点を求める
end
plot(x,y); hold on
plot(x(K),y(K),'o-','lineWidth',2)
function [xi,yi] = intersection(x,y) % 交点を求める
numx = (x(1)*y(2)-y(1)*x(2))*(x(3)-x(4))-(x(1)-x(2))*(x(3)*y(4)-y(3)*x(4));
numy = (x(1)*y(2)-y(1)*x(2))*(y(3)-y(4))-(y(1)-y(2))*(x(3)*y(4)-y(3)*x(4));
den = (x(1)-x(2))*(y(3)-y(4))-(y(1)-y(2))*(x(3)-x(4));
xi = numx / den;
yi = numy / den;
end
Atsushi Ueno
Atsushi Ueno 2023-4-14
上記のプログラムを若干整理して回答にしました。

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Atsushi Ueno
Atsushi Ueno 2023-4-14
编辑：Atsushi Ueno 2023-4-14
Ran in:
>半径が同一の円であれば、中心を囲む凸包をconvhullで求めた後に、半径をpolybufferでバッファーとすれば求められましたが、半径が異なる場合はその方法が使えません。解決策やアドバイスをお持ちの方がいましたら、ご教授いただきたいです
複数の円状に並んだ点群を囲う凸包から、円の接線のみ抽出して繋ぐアイデアです。
円周と接線の区別は場当たり的で、単純に線分の長さで判定しています。(角度刻みがやたら細かいのはその為)
%% 複数の円状に並んだ点群のサンプルデータを作成
N = 20; Rmax = 200; Rmin = 50; a = [];
th = (0:pi/1000:2*pi)'; % 角度刻み (2001行1列)
r = randi(Rmax - Rmin,1,N) + Rmin; % N円の半径( 1行N列)
x = cos(th) * r + randi(1000,1,N); % N円の点群(2001行N列)
y = sin(th) * r + randi(1000,1,N); % N円の点群(2001行N列)
%% 複数の円状に並んだ点群の「接線による凸包」を求める
K = convhull(x,y); % 点群に対する凸包を求める
dst = hypot(diff(x(K)),diff(y(K))); % 凸包の点間距離を求める
idx = dst > mean(dst); % 円周上の点間距離は短く接線は長い事を利用する
K = K([idx; false] | [false; idx]); % 接線の始点と終点のみ選り分ける
temp = [K circshift(K,-1) circshift(K,-2) circshift(K,-3)]';
for k = temp(:,1:2:end) % 4点ずつ循環シフトして並べ、2つおきに走査する
a(end+1,:) = intersection(x(k),y(k)); % 交点を求める
end
plot(x,y); hold on
plot(polyshape(a));
function ans = intersection(x,y) % 交点を求める
numx = (x(1)*y(2)-y(1)*x(2))*(x(3)-x(4))-(x(1)-x(2))*(x(3)*y(4)-y(3)*x(4));
numy = (x(1)*y(2)-y(1)*x(2))*(y(3)-y(4))-(y(1)-y(2))*(x(3)*y(4)-y(3)*x(4));
den = (x(1)-x(2))*(y(3)-y(4))-(y(1)-y(2))*(x(3)-x(4));
[numx / den, numy / den];
end
  1 个评论
YA
YA 2023-4-16
ご回答ありがとうございます。
点群を囲う凸包から、円の接線のみ抽出するというアイデアはとても参考になりました。
ありがとうございました。

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 境界領域 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!