MATLAB で計算幾何学 ~ polyshape() ~ 多角形オブジェクトで一括座標変換
まえがき
またTwitterで面白そうなのを発見!
各ボールは直線移動なのに、放物線を描いて移動しているように見えます。
最初に、MATLAB でのシミュレーション結果を示しておきます。
最初見たときは、
Hyperboloid Shower ~ 直線が曲面を作る ~ MATLABで双曲面アニメーション
と同じなのかな?と思いましたが、前のは実際に双曲面を成しているのに対し、こちらは後続につられてそう見えるだけの「錯覚」ですかね?
これを描くには、ボールの描画位置を、円に内接する、回転角度を変えた正方形の辺上を移動させれば良さそうです。
辺上の各座標を回転行列で回転させれば求められるでしょうが、面倒くさそう。(¬_¬)
計算幾何学
そこで、多分私は実際に使ったことのない、polyshape() を使ってみます。
pgon = polyshape([0 0 1 1],[1 0 0 1]); で、 4 点 (0,1)、(0,0)、(1,0)、(1,1) を頂点とする正方形オブジェクトが生成されます。
今回はガイドボックス表示に使用していますが、plot(pgon) だけで描画も行えます。
pgon = polyshape([0 0 1 1],[1 0 0 1]);
plot(pgon)
polyshape() は多角形を描画するだけでなく、頂点、固体領域、穴について記述されたプロパティをもつ polyshape オブジェクトを返します。
>> pgon = polyshape([0 0 1 1],[1 0 0 1])
pgon =
polyshape のプロパティ:
Vertices: [4×2 double]
NumRegions: 1
NumHoles: 0
これを用いて、後から境界を追加・削除したり、回転・移動・スケーリング等が行えます。
さらに、面積・重心・周囲長の計算、座標がその多角形内かどうかの判断、境界の縮小・拡大、一番近い頂点座標を求めたりもできます。矛盾があるような点を与えても、自動的に調整されます。
pgon = polyshape([0 0 1 1]*0.5+0.25,[1 0 0 1]*0.5+0.25);
plot(pgon)
xlim([0 1])
ylim([0 1])
hold on
d = rand(100,2);
dx = d(:,1); dy = d(:,2);
TFin = isinterior(pgon,dx,dy); % 領域内判定
scatter(dx(TFin),dy(TFin),'filled') % 領域内はマーカーfill
scatter(dx(~TFin),dy(~TFin),'x') % 領域外はx
hold off
> pgon = polyshape(rand(5,1), rand(5,1))
警告: 多角形状に、不正確な結果や予期しない結果を生じる可能性がある重複頂点、交差部分またはその他の矛盾があります。入力データは、適切に定義された多角形状を作成するように修正されています。
pgon =
polyshape のプロパティ:
Vertices: [13×2 double]
NumRegions: 3
NumHoles: 0
> pgon.Vertices
ans = 13×2
0.7060 0.8235
0.0971 0.0344
0.0631 0.6467
0.2769 0.3171
0.1323 0.7140
NaN NaN
0.0318 0.6948
0.0601 0.7002
0.0631 0.6467
NaN NaN
> plot(pgon)
> pgon.area
ans = 0.2066
> pgon.perimeter
ans = 3.7244
> [Cx Cy] = centroid(pgon)
Cx = 0.2913
Cy = 0.5267
th=0:pi/3:2*pi-pi/3;
pgon6 = polyshape({cos(th),0.5*cos(th)},{sin(th),0.5*sin(th)})
figure
tiledlayout flow
nexttile
plot(pgon6)
axis equal
xlim([-1 1])
ylim([-1 1])
pgon6.area
polyout = rmholes(pgon6) % 穴の削除
nexttile
plot(polyout)
axis equal
xlim([-1 1])
ylim([-1 1])
polyout.area
pgon6 =
polyshape のプロパティ:
Vertices: [13×2 double]
NumRegions: 1
NumHoles: 1
ans =
1.9486
polyout =
polyshape のプロパティ:
Vertices: [6×2 double]
NumRegions: 1
NumHoles: 0
ans =
2.5981
より詳しくは、公式ドキュメント「計算幾何学」をご覧ください。
多角形オブジェクト
今回はこの「多角形オブジェクト」を利用して、まず辺上にボール移動用の全ての点を持つ正方形オブジェクトを生成し、そのオブジェクトを回転させることで、回転後の全座標を一度に求めてしまおう、というわけです。
回転後の座標を使う場合は、Vertices プロパティにアクセスします。
polyshape() ではデフォルトだと不要な座標は削除されてしまうため、
'KeepCollinearPoints', true
にします。
これで、回転後の全座標が Vertices に入ります。
vq = -1:1; % 辺上の点を定義
vqL = length(vq);
pgon = polyshape([vq ones(1,vqL-1) -vq(2:end) -ones(1,vqL-1)], ...
[ones(1,vqL) -vq(2:end) -ones(1,vqL-1) vq(2:end)],'KeepCollinearPoints', true);
plot(pgon)
hold on
axis equal
idx=1:length(pgon.Vertices);
plot(pgon.Vertices(idx,1), pgon.Vertices(idx,2), 'o');
r = rotate(pgon,45);
plot(r)
plot(r.Vertices(idx,1), r.Vertices(idx,2), '*');
hold off
>> pgon.Vertices
ans =
-1 1
0 1
1 1
1 0
1 -1
0 -1
-1 -1
-1 0
>> r.Vertices
ans =
-1.4142 0
-0.7071 0.7071
0 1.4142
0.7071 0.7071
1.4142 0
0.7071 -0.7071
0 -1.4142
-0.7071 -0.7071
オブジェクト回転前後の各座標
>> pgon0 = polyshape([vq ones(1,vqL-1) -vq(2:end) -ones(1,vqL-1)], ...
[ones(1,vqL) -vq(2:end) -ones(1,vqL-1) vq(2:end)])
警告: 多角形状に、不正確な結果や予期しない結果を生じる可能性がある重複頂点、交差部分またはその他の矛盾があります。入力
データは、適切に定義された多角形状を作成するように修正されています。
> polyshape/checkAndSimplify (行 517) 内
polyshape (行 175) 内
pgon0 =
polyshape のプロパティ:
Vertices: [4×2 double]
NumRegions: 1
NumHoles: 0
>> pgon0.Vertices
ans =
-1 1
1 1
1 -1
-1 -1
デフォルト設定では不要点が削除される
これでもう、簡単に書けそう!
シミュレーション結果
各ボールは正方形の辺上を直線移動しているだけなのですが、放物線状に見えますね。(;¬_¬)
直線移動であることを確認しやすいよう、ガイドボックス表示をしてみましょう。
先頭のボールだけを追うと、直線移動であることが確認できると思います。
ほんと、人の認識はそのまま信用できませんね。(u_u)
最後に MATLAB スクリプトを貼っておきます。
boxf = true; にするとガイドボックスが表示され、we = true; にすると gif ファイルが書かれます。ライブスクリプトでは、これらを ”チェックボックス” にすると良いかと思います。
ビートルズが解散した頃にはまだ存在しなかった "computational geometry(計算幾何学)" 、楽しんでみてください。
では。
MATLAB スクリプト
boxf = false;
we = false;
filename = sprintf('BouncingBall%d.gif',boxf)
f = figure(1);
f.Color = [1 1 1];
vq = -1:0.2:1;
vqL = length(vq);
pgon = polyshape([vq ones(1,vqL-1) -vq(2:end) -ones(1,vqL-1)], ...
[ones(1,vqL) -vq(2:end) -ones(1,vqL-1) vq(2:end)],'KeepCollinearPoints', true);
len = length(pgon.Vertices);
col = [0.3020 0.7451 0.9333];
t=0:pi/50:2*pi;
a = sqrt(2)+0.05;
dth = 5;
st = true; % 1st frame
for l = 1:len*(1+boxf)
plot(a*cos(t), a*sin(t), 'r', LineWidth=2)
hold on
alpha = max(0.7 - abs(l/len-1), 0);
for th = 0:dth:45
r = rotate(pgon,th);
if boxf
if th==45
wid = 1.5;
boxCol = 'r';
else
wid = 0.7;
boxCol = 'b';
end
plot(r, EdgeColor=boxCol, EdgeAlpha=alpha, FaceColor='w', FaceAlpha=0, LineWidth=wid);
end
idx = mod(l + th/dth*2, len) + 1;
if (boxf && th==45)
p = plot(r.Vertices(idx,1), r.Vertices(idx,2), 'o', MarkerSize=10, ...
MarkerFaceColor= [col(1)+(1-col(1))*(alpha/0.7) col(2)*(1-alpha/0.7) col(3)*(1-alpha/0.7)]);
else
p = plot(r.Vertices(idx,1), r.Vertices(idx,2), 'o', MarkerSize=10, ...
MarkerFaceColor=col);
end
end
axis off
axis(gca,'equal')
hold off
drawnow
if we
frame = getframe(f); % capture with border frame
% frame = getframe(gca); % take image only
gr = frame2im(frame);
[bA,map] = rgb2ind(gr,256);
if (st)
imwrite(bA,map,filename,'gif','LoopCount',Inf,'DelayTime',1/15);
st = false;
else
imwrite(bA,map,filename,'gif','WriteMode','append','DelayTime',1/15);
end
end
end
この記事が気に入ったらサポートをしてみませんか?