Effacer les filtres
Effacer les filtres

複数円を囲む凸包

4 vues (au cours des 30 derniers jours)
YA
YA le 10 Avr 2023
Commenté : YA le 16 Avr 2023
複数の円を囲む凸包を求める方法についての質問です。
複数の異なる半径を持つ円があるときに、それを囲む凸包を求めたいです。
半径が同一の円であれば、中心を囲む凸包をconvhullで求めた後に、半径をpolybufferでバッファーとすれば求められましたが、
半径が異なる場合はその方法が使えません。
解決策やアドバイスをお持ちの方がいましたら、ご教授いただきたいです。
よろしくお願いいたします。
  6 commentaires
Atsushi Ueno
Atsushi Ueno le 13 Avr 2023
Déplacé(e) : Atsushi Ueno le 13 Avr 2023
>接線を求めて地道に求めていくしかなさそうですね
もう作ってしまったかもしれませんが...あくまで 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(:); % 3つおきに交点の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 le 14 Avr 2023
上記のプログラムを若干整理して回答にしました。

Connectez-vous pour commenter.

Réponse acceptée

Atsushi Ueno
Atsushi Ueno le 14 Avr 2023
Modifié(e) : Atsushi Ueno le 14 Avr 2023
>半径が同一の円であれば、中心を囲む凸包を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 commentaire
YA
YA le 16 Avr 2023
ご回答ありがとうございます。
点群を囲う凸包から、円の接線のみ抽出するというアイデアはとても参考になりました。
ありがとうございました。

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur 境界領域 dans Help Center et File Exchange

Produits


Version

R2022b

Community Treasure Hunt

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

Start Hunting!