複数円を囲む凸包
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
複数の円を囲む凸包を求める方法についての質問です。
複数の異なる半径を持つ円があるときに、それを囲む凸包を求めたいです。
半径が同一の円であれば、中心を囲む凸包をconvhullで求めた後に、半径をpolybufferでバッファーとすれば求められましたが、
半径が異なる場合はその方法が使えません。
解決策やアドバイスをお持ちの方がいましたら、ご教授いただきたいです。
よろしくお願いいたします。
6 commentaires
Atsushi Ueno
le 13 Avr 2023
Déplacé(e) : Atsushi Ueno
le 13 Avr 2023
>接線を求めて地道に求めていくしかなさそうですね
- 凸包を求める
- 接線を結ぶ点群に絞る(円周上の点群を削除)
- 接線同士の交点を求め点群に追加・凸包の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
Réponse acceptée
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
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!