複数円を囲む凸包

複数の円を囲む凸包を求める方法についての質問です。

複数の異なる半径を持つ円があるときに、それを囲む凸包を求めたいです。

半径が同一の円であれば、中心を囲む凸包をconvhullで求めた後に、半径をpolybufferでバッファーとすれば求められましたが、

半径が異なる場合はその方法が使えません。

解決策やアドバイスをお持ちの方がいましたら、ご教授いただきたいです。

よろしくお願いいたします。

[YA]

説明が不足していました。

凸包の角は円周を使うのではなく直線で角を作るものとします。

[Atsushi Ueno]

>凸包の角は円周を使うのではなく直線で角を作るものとします

つまり円の接線で構成され面積最小の多角形を求める事になりますね。難しいですね。

下記は円周を覆う凸包になりました。

x = []; y = []; th = (0:pi/5:2*pi)';
for n = 1:10
    r = randi(50,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]);
plot(x,y,'*')
hold on
plot(x(K),y(K),'r')

[YA]

ご回答ありがとうございました。

円を点の集合としているので，convhullを使うと角が滑らかになってしまいますね。

接線を求めて地道に求めていくしかなさそうですね。

ありがとうございました。

[Atsushi Ueno]

>接線を求めて地道に求めていくしかなさそうですね

もう作ってしまったかもしれませんが...

中心を囲む凸包をconvhullで求めた後に、半径分離れた直線(まだ繋がっていない)を引いてみました。あとは直線同士の交点を求めれば良さそうですが、円の中心を囲む凸包に含まれない円の半径が大きいと凸包の外にはみ出してしまうので、もう少し検討すべきかなと思いました。

th = (0:0.2:2*pi)';
r = [randi(800,10,2) + 100, randi(200,10,1)+50]; % x座標、y座標、半径
K = convhull(r(:,1:2)); % (x,y)座標の組から凸包のインデックスを得る
plot(r(K,1),r(K,2),'o-'); hold on % 円の中心を囲む凸包を結ぶ
plot(r(:,1)'+cos(th)*r(:,3)',r(:,2)'+sin(th)*r(:,3)','.'); % 円を点線で描画

for k = [K(1:end-1)'; circshift(K(1:end-1),-1)']
    theta = atan2(diff(r(k,2)),diff(r(k,1))); % 凸包の一辺の角度を求める
    x = r(k,1) + cos(theta - pi/2)*r(k,3); % 円を囲む凸包の座標
    y = r(k,2) + sin(theta - pi/2)*r(k,3); % 円を囲む凸包の座標
    plot(x,y,'LineWidth',2); % 円を囲む凸包の描画（まだ角が繋がっていない)
end

[Atsushi Ueno]

>接線を求めて地道に求めていくしかなさそうですね

もう作ってしまったかもしれませんが...あくまで 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]

上記のプログラムを若干整理して回答にしました。

Answer 1

[Atsushi Ueno]

>半径が同一の円であれば、中心を囲む凸包を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

[YA]

ご回答ありがとうございます。

点群を囲う凸包から、円の接線のみ抽出するというアイデアはとても参考になりました。

ありがとうございました。