テンソルによるアフィン変換の応用～基本から物体検出まで～

テンソルによるアフィン変換を使って前処理を簡単にしようという記事です。アフィン変換は画像に対しても適用できますが、同一の変換をBounding Boxなどアノテーションにも適用することができます。またテンソルに拡張することで、1つのオブジェクトに対して異なる変換を同時にかけることができます。

アフィン変換の基本

以下の数式で点の移動を表したのがアフィン変換。

$$\begin{bmatrix}x' \\ y' \\ 1 \end{bmatrix} = \begin{bmatrix}a & b & c \\ d & e & f \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x \\ y \\ 1 \end{bmatrix} \tag{1}$$

アフィン変換の漠然とした理解を見たい方はこちらの記事もどうぞ。

3x3行列と3x1行列の積を取っていますが、3x3行列が変換を定義したもの、3x1行列が移動前の点です。

具体例1：平行移動

点(10, 20)を$x$方向に100、$y$方向に50移動したとき、移動後の座標は(110, 70)ですが、次のように表されます。

$$\begin{bmatrix}110 \\ 70 \\ 1 \end{bmatrix} = \begin{bmatrix}1 & 0 & 100 \\ 0 & 1 & 50 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}10 \\ 20 \\ 1 \end{bmatrix} \tag{2}$$

「こんなの行列使うまでもないじゃないか」と思うかもしれませんが、行列で表せるのがポイントです。コードで書けば次の通りです。

affine = np.array([[1, 0, 100], [0, 1, 50], [0, 0, 1]])
source = np.array([10, 20, 1])[:, None]
dest = np.dot(affine, source)
print(dest)
# [[110]
# [ 70]
# [  1]]

$x$方向の移動を$t_x$, $y$方向の移動を$t_y$とすると、平行移動の変換行列は以下の通りです。

$$\begin{bmatrix}1 & 0 & t_x \\ 0 & 1 & t_y \\ 0 & 0 & 1 \end{bmatrix} \tag{3}$$

なお、各行列の3行目は特に意味のない数値です。行列積の計算の都合上1を入れています。

具体例2：拡大縮小

点(50, 100)を$x$方向に2倍、$y$方向に0.8倍したとき、移動後の座標は(100, 80)ですが、次のように表されます。

$$\begin{bmatrix}100 \\ 80 \\ 1 \end{bmatrix} = \begin{bmatrix}2 & 0 & 0 \\ 0 & 0.8 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}50 \\ 100 \\ 1 \end{bmatrix} \tag{4}$$

affine = np.array([[2, 0, 0], [0, 0.8, 0], [0, 0, 1]])
source = np.array([50, 100, 1])[:, None]
dest = np.dot(affine, source)
print(dest)
# [[100.]
#  [ 80.]
#  [  1.]]

$x$方向の拡大を$s_x$, $y$方向の拡大を$s_y$とすると、拡大縮小の変換行列は以下の通りです。

$$\begin{bmatrix}s_x & 0 & 0 \\ 0 & s_y & 0 \\ 0 & 0 & 1 \end{bmatrix} \tag{5}$$

複数点に対するアフィン変換

アフィン変換は行列計算なので、点の数を任意に拡張できます。$N$個の点に対してアフィン変換を求めるときは、3x3行列と3xN行列の積を取ります。結果は3xN行列になります。

$$\begin{bmatrix} x_1' & \cdots & x_N' \\ y_1' & \cdots & y_N' \\ 1 & \cdots & 1 \end{bmatrix} = \begin{bmatrix} a & b & c \\ d & e & f \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 & \cdots & x_N \\ y_1 & \cdots & y_N \\ 1 & \cdots & 1 \end{bmatrix} \tag{6}$$

具体例3：四角形のアフィン変換

$N=4$個の点をアフィン変換すると、四角形のアフィン変換になります。$(x, y)$方向に$(2, 3)$倍拡大し、$(5, -1)$平行移動するアフィン変換を、四角形の頂点に適用すると次のようになります。

w, h = 2, 1
points = np.array([[0, w, w, 0], [0, 0, h, h], [1, 1, 1, 1]], np.float32) # (3, 4)
affine = np.array([[2, 0, 5], [0, 3, -1], [0, 0, 1]])  # (3, 3)
dest = np.dot(affine, points)

plt.scatter(points[0,:], points[1,:], color="cyan")
plt.scatter(dest[0,:], dest[1,:], color="magenta")
plt.show()

具体例4:回転する四角形とBoudning Box

これは物体検出の前処理で使える例です。物体検出のData Augmentationで画像を回転させたとき、Bounding Boxも回転させる必要があります。Bounding Boxは左上と右下の2点で定義できますが、頂点の4点を取ることで、回転後のBounding Boxを簡単に計算できます。回転後のx, yについて最小・最大値を取ると、回転後のBounding Boxの左上、右下の座標が求められます。

回転もアフィン変換の1つで、原点を中心に反時計回りに$\theta$回転させたときの変換行列は、

$$\begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \tag{7}$$

であり、回転前のBounding Box（四角形）の頂点が回転後に$(x_1', y_1'), \cdots, (x_4', y_4')$に移動したとき、傾いた四角形に外接する新たなBounding Boxの座標は、

• 左上 $(\min(x_1', \cdots, x_4'), \min(y_1', \cdots, y_4'))$
• 右下 $(\max(x_1', \cdots, x_4'), \max(y_1', \cdots, y_4'))$

で求めることができます。なぜこの計算をする必要があるのかというと、回転後の元の頂点はBounding Boxとして不適切なため（$xy$軸に平行な長方形ではないから）、調整する必要があるからです。詳しくは下の動画を見てください。

プロットすると次のようになります。プロットの関係上5点を（原点を重複させて）移動させていますが、計算だけなら4点の移動でOKです。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches

def rotate_box():
    w, h = 2, 1
    max_wh = max(w, h)
    points = np.array([[0, w, w, 0, 0], [0, 0, h, h, 0], [1, 1, 1, 1, 0]], np.float32) # 元のBounding Box
    for theta in range(0, 360, 10):
        rad = np.radians(theta)
        rotate_matrix = np.array([
            [np.cos(rad), -np.sin(rad), 0],
            [np.sin(rad), np.cos(rad), 0],
            [0, 0, 1]], np.float32)
        dest_points = np.dot(rotate_matrix, points)[:2, :] # 回転後の四角形
        rectangle = np.concatenate([np.min(dest_points, axis=-1),
                                    np.max(dest_points, axis=-1)]) # 新しいBounding Box

        plt.clf()
        plt.plot(dest_points[0,:], dest_points[1,:], linewidth=2, marker="o")

        ax = plt.gca()
        rect = patches.Rectangle(rectangle[:2], *(rectangle[2:] - rectangle[:2]),
                                 linewidth=1, edgecolor="magenta", fill=False)
        ax.add_patch(rect)
        plt.ylim(-max_wh*2, max_wh*2)
        plt.xlim(-max_wh * 2, max_wh * 2)
        plt.title("degree = " + str(theta))
        plt.show()

アフィン変換の合成

アフィン変換同士は行列積を取ることで合成できます。$A_1->A_2$という変換をするとき、$A_2A_1P$（$P$は点の行列）という積を取ります。順序が逆になるのに注意してください。また交換法則が成り立たなく、順序を入れ替えると結果が異なります。

例えば、$A_1$が$x,y$を2倍拡大、$A_2$が$x$方向に50、$y$方向に100移動する変換とします。このとき、$A_2A_1(A_1 \to A_2)$は、

$$A_2A_1 = \begin{bmatrix}1 & 0 & 50 \\ 0 & 1 & 100 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix}2 & 0 & 50 \\ 0 & 2 & 100 \\ 0 & 0 & 1 \end{bmatrix} \tag{8}$$

ですが、$A_1A_2(A_2 \to A_1)$は、

$$A_1A_2 = \begin{bmatrix}2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}1 & 0 & 50 \\ 0 & 1 & 100 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix}2 & 0 & 100 \\ 0 & 2 & 200 \\ 0 & 0 & 1 \end{bmatrix} \tag{9}$$

となり、平行移動の大きさが変わります。これは「平行移動してから拡大するか、拡大してから平行移動するか」の差を表します。順序がわからなくなったら、このような単純な例で計算してみるのがおすすめです。

具体例5：回転しながら迫りくる頭の悪い人

先程アフィン変換は任意の数の点に適用できるので、点の数が数千あっても大丈夫です。ここでは「頭の悪い人」（フリー素材より）

を点群データに変換し、拡大と回転のアフィン変換の合成を行いながらプロットしてみます。例4とほとんど同じです。

from PIL import Image

def atamanowaruihito():
    with Image.open("atamanowaruihito.png") as img:
        img = img.resize((img.width // 2, img.height // 2))        
        img = img.convert("L").point(lambda x: 255 if x >= 128 else 0)  # グレー化、二値化
        points = np.stack(np.where(np.array(img) == 0)[::-1], axis=0)  # yx -> xyして点群の行列化
        points[1,:] = img.height - points[1,:]  # y軸を下向きが正から上向きが正に修正 (2, 5912)
        points = np.concatenate([points, np.ones_like(points[0:1, :])], axis=0) # 3行目に1をつける (3, 5912)

    for theta in range(0, 360, 1):
        rad = np.radians(theta)
        rotate_matrix = np.array([
            [np.cos(rad), -np.sin(rad), 0],
            [np.sin(rad), np.cos(rad), 0],
            [0, 0, 1]], np.float32) # 回転行列
        scale_matrix = np.eye(3, dtype=np.float32) * (1 + theta / 180)
        # アフィン変換の合成は変換行列の行列積
        # A1->A2ならA2A1の順に積を取る（順序に注意）
        affine = np.dot(scale_matrix, rotate_matrix)  # 回転して拡大
                
        dest_points = np.dot(affine, points)[:2, :] # 回転後の点群
        rectangle = np.concatenate([np.min(dest_points, axis=-1),
                                    np.max(dest_points, axis=-1)]) # 点群に対するBounding Box

        plt.clf()
        plt.scatter(dest_points[0,:], dest_points[1,:], s=1)

        ax = plt.gca()
        rect = patches.Rectangle(rectangle[:2], *(rectangle[2:] - rectangle[:2]),
                                 linewidth=1, edgecolor="magenta", fill=False)
        ax.add_patch(rect)
        plt.ylim(-750, 750)
        plt.xlim(-750, 750)
        plt.show()

Bounding Boxは四角形の例と同様にして計算できます。

複数のアフィン変換を同時に適用する

ここからがテンソル計算です。アフィン変換の合成ではなく、同一の点に対し、複数のアフィン変換を同時に適用する方法を考えます。式では表しづらいのでコードで考えます。

もしアフィン変換が1つだったら、変換行列は(3, 3)というshapeでした。では変換が2個つまり、(2, 3, 3)というshapeならどうでしょうか？　つまり、

# points(4点): (3, 4), affines: (2, 3, 3)
output = np.zeros((2, 3, 4)
for i in range(affines.shape[0]):
        output[i] = np.dot(affines[i], points)

という計算をしたいのです。実はこれはforループを使わずにワンライナーで、

output = np.matmul(affines, points)

で計算できます。

具体例6：四角形の複数のアフィン変換

例3のアフィン変換を3つにしてみます。

1. 元の四角形そのまま（恒等変換）
2. $(x, y)$方向に$(2, 3)$倍拡大し、$(5, -1)$平行移動
3. $(x, y)$方向に$(3, 1)$倍拡大し、$(10, 2)$平行移動

from matplotlib import cm

def rectangle_multi():
    w, h = 2, 1
    points = np.array([[0, w, w, 0], [0, 0, h, h], [1, 1, 1, 1]], np.float32)  # (3, 4)
    a1 = np.eye(3) # 恒等変換
    a2 = np.array([[2, 0, 5], [0, 3, -1], [0, 0, 1]]) # (3, 3)
    a3 = np.array([[3, 0, 10], [0, 1, 2], [0, 0, 1]])  # (3, 3)
    affine = np.stack([a1, a2, a3], axis=0) # (3, 3, 3)
    dest = np.matmul(affine, points)  # (3, 3, 4)

    cmap = cm.get_cmap("tab10")
    for i in range(3):
        plt.scatter(dest[i,0,:], dest[i, 1,:], color=cmap(i))
    plt.show()

たった1回の計算で3つの異なる変換を実行することができました。このように一対多の変換もテンソルによるアフィン変換を使えばできます。

具体例7：物体検出のAnchor Box

物体検出ではニューラルネットワークの出力の各点（Anchor）から見たBounding Boxの予測を行います。Bounding Boxは生の画像の（座標）に対応しますが、物体検出では多くの座標が登場します。例えば、

• 画像をニューラルネットワークの入力サイズに合わせて拡大縮小したため、軸のスケールが変わる
• アンカーから見たBounding Boxの値を予測するため、各アンカーに対して平行移動と同様の操作が入る
• RoI Poolingを使う場合は、RoIの中の座標の概念がある（クロップ分のスケール・平行移動がある）

つまり「座標から座標へ」の柔軟な点の変換が必要になります。この座標の変換をテンソルによるアフィン変換で統一することができます。

ここで、次のことを考えます。

• ニューラルネットワークの出力解像度は4x4で、Anchor Boxも4x4個ある
• ニューラルネットワークの入力から出力で、解像度は1/16になる。つまり、Anchor Boxの1点は16x16の解像度に対応する。
• 前処理で画像を横方向に1.5倍、縦方向に0.8倍リサイズしており、Ground TruthのBounding Boxのデータは変形前のピクセル値に基づく
• 画像に対して回転や左右反転は行わない（したがってBounding Boxの左上と右下だけ考えれば良い）
• Bounding Boxが(ymin, xmin, ymax, xmax)の形式のとき、Ground Truthで(10, 20, 30, 40), (30, 30, 50, 50)の2つある。このとき、各Anchor Boxから見たGround TruthのBounding Boxの座標を求めたい。ただし、各Anchorは1x1のピクセルの中点にあり、$x,y$ともに0.5分のオフセットを取った地点を各アンカーの原点とする。

方針としては2つのアフィン変換を考えます。

1. 前処理（元の画像→入力画像）のアフィン変換。$x$方向に1.5倍、$y$方向に0.8倍。(3, 3)のshape。
2. 入力画像→anchorのアフィン変換。(4, 4, 3, 3)のshapeで、(3, 3)のアフィン変換を縦横4x4に重ねたもの。(i, j, :, :)$(0\leq i,j \leq 3)$の位置のアフィン変換は以下のとおりです。Anchor Boxから見た平行移動なので符号がマイナスになります。

$$\begin{bmatrix}1/16 & 0 & -0.5+i \\ 0 & 1/16 & -0.5+j \\ 0 & 0 & 1 \end{bmatrix} \tag{10}$$

1, 2の2つを合成すればいいわけです。これが美味しいのは逆変換（anchor→元の画像）をしたいときは、1,2の合成テンソルの逆行列を取る、つまり

inv_transform = np.linalg.inv(combined_affine)

をするだけで逆変換のテンソルが取れることです。np.linalg.invの入力がテンソルの場合は、最後の2軸について逆行列を取ったものを積み重ねます。

物体検出では$(x, y)$の表記より、$(y, x)$の表記のほうが都合がいいので（画像のテンソルが$(B, H, W, C)$なので）、平行移動$t_x, t_y$、拡大縮小$s_x, s_y$のアフィン変換は、

$$\begin{bmatrix}s_y & 0 & t_y \\ 0 & s_x & t_X \\ 0 & 0 & 1 \end{bmatrix} \tag{11}$$

で表されます。軸を入れ替えてもアフィン変換として機能します。

def anchor_box():
    bounding_boxes = np.array([[10, 20, 30, 40], [30, 30, 50, 50]])  # (2, 4)
    points = bounding_boxes.reshape(-1, 2).T  # (4, 2) -> (2, 4)
    points = np.concatenate([points, np.ones_like(points[0:1,:])], axis=0)  # (3, 4)
    # yx座標で考える
    a1 = np.array([[0.8, 0, 0], [0, 1.5, 0], [0, 0, 1]])  # y方向に0.8倍、x方向に1.5倍
    # アンカー
    offset_x, offset_y = np.meshgrid(-(np.arange(4) + 0.5), -(np.arange(4) + 0.5)) # アンカーから入力画像の原点を見るのでマイナス平行移動
    a2 = np.zeros((4, 4, 3, 3)) + np.eye(3).reshape(1, 1, 3, 3) / 16.0 # (4, 4, 3, 3)にブロードキャスト
    a2[:,:,0,2] = offset_y
    a2[:,:, 1, 2] = offset_x
    a2[:,:, 2, 2] = 1.0
    # アフィン合成
    affine = np.matmul(a2, a1)
    # アフィン変換
    raw_dest = np.matmul(affine, points)  # (4, 4, 3, 4)
    dest = raw_dest.swapaxes(-1, -2)[:,:,:,:2]  # (4, 4, 4, 3) 転置のテンソル版 -> (4, 4, 4, 2)
    dest = dest.reshape(4, 4, 2, 4)
    print("元画像のBounding Boxを、各アンカーでの見たときのアフィン変換後の座標")
    print(dest)

    # 逆変換して確かめる
    raw_inv = np.matmul(np.linalg.inv(affine), raw_dest)
    inv = raw_inv.swapaxes(-1, -2)[:,:,:,:2]
    inv = inv.reshape(4, 4, 2, 4)  # bounding boxesと一致
    print("変換後の座標を逆変換すると元の値に戻る")
    print(inv)

    # 逆変換アフィン（確認）
    inv_transoform = np.linalg.inv(affine)
    print("逆変換のアフィン（デバッグ用）")
    print(inv_transoform)

クリックして出力を表示

元画像のBounding Boxを、各アンカーでの見たときのアフィン変換後の座標
[[[[ 0.      1.375   1.      3.25  ]
   [ 1.      2.3125  2.      4.1875]]

  [[ 0.      0.375   1.      2.25  ]
   [ 1.      1.3125  2.      3.1875]]

  [[ 0.     -0.625   1.      1.25  ]
   [ 1.      0.3125  2.      2.1875]]

  [[ 0.     -1.625   1.      0.25  ]
   [ 1.     -0.6875  2.      1.1875]]]


 [[[-1.      1.375   0.      3.25  ]
   [ 0.      2.3125  1.      4.1875]]

  [[-1.      0.375   0.      2.25  ]
   [ 0.      1.3125  1.      3.1875]]

  [[-1.     -0.625   0.      1.25  ]
   [ 0.      0.3125  1.      2.1875]]

  [[-1.     -1.625   0.      0.25  ]
   [ 0.     -0.6875  1.      1.1875]]]


 [[[-2.      1.375  -1.      3.25  ]
   [-1.      2.3125  0.      4.1875]]

  [[-2.      0.375  -1.      2.25  ]
   [-1.      1.3125  0.      3.1875]]

  [[-2.     -0.625  -1.      1.25  ]
   [-1.      0.3125  0.      2.1875]]

  [[-2.     -1.625  -1.      0.25  ]
   [-1.     -0.6875  0.      1.1875]]]


 [[[-3.      1.375  -2.      3.25  ]
   [-2.      2.3125 -1.      4.1875]]

  [[-3.      0.375  -2.      2.25  ]
   [-2.      1.3125 -1.      3.1875]]

  [[-3.     -0.625  -2.      1.25  ]
   [-2.      0.3125 -1.      2.1875]]

  [[-3.     -1.625  -2.      0.25  ]
   [-2.     -0.6875 -1.      1.1875]]]]
変換後の座標を逆変換すると元の値に戻る
[[[[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]]


 [[[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]]


 [[[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]]


 [[[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]

  [[10. 20. 30. 40.]
   [30. 30. 50. 50.]]]]
逆変換のアフィン（デバッグ用）
[[[[20.          0.         10.        ]
   [ 0.         10.66666667  5.33333333]
   [ 0.          0.          1.        ]]

  [[20.          0.         10.        ]
   [ 0.         10.66666667 16.        ]
   [ 0.          0.          1.        ]]

  [[20.          0.         10.        ]
   [ 0.         10.66666667 26.66666667]
   [ 0.          0.          1.        ]]

  [[20.          0.         10.        ]
   [ 0.         10.66666667 37.33333333]
   [ 0.          0.          1.        ]]]


 [[[20.          0.         30.        ]
   [ 0.         10.66666667  5.33333333]
   [ 0.          0.          1.        ]]

  [[20.          0.         30.        ]
   [ 0.         10.66666667 16.        ]
   [ 0.          0.          1.        ]]

  [[20.          0.         30.        ]
   [ 0.         10.66666667 26.66666667]
   [ 0.          0.          1.        ]]

  [[20.          0.         30.        ]
   [ 0.         10.66666667 37.33333333]
   [ 0.          0.          1.        ]]]


 [[[20.          0.         50.        ]
   [ 0.         10.66666667  5.33333333]
   [ 0.          0.          1.        ]]

  [[20.          0.         50.        ]
   [ 0.         10.66666667 16.        ]
   [ 0.          0.          1.        ]]

  [[20.          0.         50.        ]
   [ 0.         10.66666667 26.66666667]
   [ 0.          0.          1.        ]]

  [[20.          0.         50.        ]
   [ 0.         10.66666667 37.33333333]
   [ 0.          0.          1.        ]]]


 [[[20.          0.         70.        ]
   [ 0.         10.66666667  5.33333333]
   [ 0.          0.          1.        ]]

  [[20.          0.         70.        ]
   [ 0.         10.66666667 16.        ]
   [ 0.          0.          1.        ]]

  [[20.          0.         70.        ]
   [ 0.         10.66666667 26.66666667]
   [ 0.          0.          1.        ]]

  [[20.          0.         70.        ]
   [ 0.         10.66666667 37.33333333]
   [ 0.          0.          1.        ]]]]

元画像→アンカー、アンカー→元画像の逆変換をしても元のBounding Boxの座標に戻っているのが確認できます。逆変換が行列計算で1本なのは強いです。アンカーの変換がうまく行ってるか見たいときは、逆変換のアフィン（逆行列）を確認するのがわかりやすいです。

具体例8：とにかく頭の悪そうな”頭の悪い人”

テンソルによるアフィン変換を活用するとこのようなこともできます。以下の動画を見てください。

解説は割愛しますが気になる方はコードを見てみてください。

クリックしてコードを表示

def atamanowaruihito2():
    with Image.open("atamanowaruihito.png") as img:
        img = img.resize((img.width // 2, img.height // 2))        
        img = img.convert("L").point(lambda x: 255 if x >= 128 else 0)  # グレー化、二値化
        points = np.stack(np.where(np.array(img) == 0)[::-1], axis=0)  # yx -> xyして点群の行列化
        points[1,:] = img.height - points[1,:]  # y軸を下向きが正から上向きが正に修正 (2, 5912)
        points = np.concatenate([points, np.ones_like(points[0:1, :])], axis=0) # 3行目に1をつける (3, 5912)

    # 回転行列用乱数
    a = np.random.uniform(1.0, 5.0, size=(4, 4))
    b = np.random.uniform(-180, 180, size=a.shape)
    # 拡大行列用乱数
    c = np.random.uniform(0.5, 1.5, size=a.shape)
    d = np.random.uniform(1.0, 5.0, size=a.shape)
    e = np.random.uniform(-180, 180, size=a.shape)
    f = np.random.uniform(0.5, 1.0, size=a.shape) + c
    # 平行移動用乱数
    e = np.random.uniform(0, 200, size=a.shape)
    g = np.random.uniform(1.0, 5.0, size=a.shape)
    h = np.random.uniform(-180, 180, size=a.shape)
    
    for theta in range(0, 360, 10):
        # 回転変換
        rad = np.radians(a * theta + b)
        rotate_tensor = np.broadcast_to(np.eye(3)[None, None,:], (4, 4, 3, 3)).copy()
        rotate_tensor[:,:, 0, 0] = np.cos(rad)
        rotate_tensor[:,:, 0, 1] = -np.sin(rad)
        rotate_tensor[:,:, 1, 0] = np.sin(rad)
        rotate_tensor[:,:, 1, 1] = np.cos(rad)
        # 拡大変換
        rad = np.radians(d * theta + e)
        scale_tensor = np.broadcast_to(np.eye(3)[None, None,:], (4, 4, 3, 3)).copy()
        scale_tensor[:,:, 0, 0] = c * np.sin(rad) + f
        scale_tensor[:,:, 1, 1] = c * np.sin(rad) + f
        # 個別の平行移動
        rad = np.radians(g * theta + h)
        transform_tensor = np.broadcast_to(np.eye(3)[None, None,:], (4, 4, 3, 3)).copy()
        transform_tensor[:,:, 0, 2] = e * np.cos(rad)        
        transform_tensor[:,:, 1, 2] = e * np.sin(rad)
        # 全体の平行移動
        shift_x, shift_y = np.meshgrid(np.arange(4), np.arange(4))
        anchor_tensor = np.broadcast_to(np.eye(3)[None, None,:], (4, 4, 3, 3)).copy()
        anchor_tensor[:,:, 0, 2] = shift_x * 500
        anchor_tensor[:,:, 1, 2] = shift_y * 500

        # 回転→拡大→平行移動→全体の平行移動
        affine = np.matmul(anchor_tensor, transform_tensor)
        affine = np.matmul(affine, scale_tensor)
        affine = np.matmul(affine, rotate_tensor)
        dest_points = np.matmul(affine, points)[:, :, :2, :] # 回転後の点群
        rectangle = np.concatenate([np.min(dest_points, axis=-1),
                                    np.max(dest_points, axis=-1)], axis=-1)  # 点群に対するBounding Box
                                    
        dest_points = dest_points.swapaxes(-1, -2).reshape(-1, 2)
        rectangle = rectangle.reshape(-1, 4)

        plt.clf()
        plt.scatter(dest_points[:, 0], dest_points[:, 1], s=1)

        ax = plt.gca()
        for i in range(rectangle.shape[0]):
            rect = patches.Rectangle(rectangle[i, :2], *(rectangle[i, 2:] - rectangle[i, :2]),
                                    linewidth=1, edgecolor="magenta", fill=False)
            ax.add_patch(rect)
        plt.ylim(-750, 2500)
        plt.xlim(-750, 2500)
        plt.title("theta = " + str(theta))
        plt.show()