見出し画像

MATLAB で計算幾何学 ～ polyshape() ～ 多角形オブジェクトで一括座標変換

leftbank

まえがき

またTwitterで面白そうなのを発見！

各ボールは直線移動なのに、放物線を描いて移動しているように見えます。

The balls are going straight.


最初に、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() は多角形を描画するだけでなく、頂点、固体領域、穴について記述されたプロパティをもつ 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') % 領域外はｘ
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

この記事が気に入ったらサポートをしてみませんか？