Shapelyを用いた幾何学的オブジェクトのアフィン変換の例

# Shapelyを使用した幾何学的オブジェクトの変換

[Shapelyライブラリ](https://shapely.readthedocs.io/en/stable/manual.html#shapely.affinity.rotate)はPythonで幾何学的オブジェクトを操作するための強力なツールです。この記事では、ShapelyのAffine変換機能を用いて、幾何学的オブジェクトの様々な変換（平行移動、回転、拡大/縮小、せん断）を行う方法を解説します。

Shapelyのマニュアルは、

https://shapely.readthedocs.io/en/stable/manual.html#shapely.affinity.rotate

を参考にしてください。

この記事は、

https://colab.research.google.com/drive/1Uyk2G-Vg91qlyVbArpY282W3LUb1YOp1?usp=sharing

から実行結果が閲覧できます。

## Affine変換の基本

Affine変換には、オブジェクトの形状や位置を変更するためのさまざまな操作が含まれます。これには、平行移動、回転、拡大/縮小、せん断などがあります。ShapelyとMatplotlibを組み合わせることで、これらの変換を視覚的に確認することができます。

### Affine変換の基本的な使用方法

ShapelyでAffine変換を行うには、`shapely.affinity`モジュールの関数を使用します。主な関数は以下のとおりです：

- `affine_transform(geometry, matrix)`: 一般的なAffine変換を適用します。
- `rotate(geometry, angle, origin='center', use_radians=False)`: オブジェクトを回転させます。
- `scale(geometry, xfact, yfact, origin='center')`: オブジェクトを拡大縮小します。
- `translate(geometry, xoff, yoff)`: オブジェクトを平行移動させます。

使い方の詳細は、github にコード `shapely/affinity.py` が公開されてるので、

https://github.com/shapely/shapely/blob/main/shapely/affinity.py

これを読むのが早いかもしれません。

### affine_transform

`affine_transform` 関数は、一般的なアフィン変換（平行移動、回転、拡大/縮小、せん断）を適用します。関数の形式は以下の通りです：

``` python:
affine_transform(geometry, matrix)
```

- geometry: 変換するShapelyの幾何学的オブジェクトです。
- matrix: 6要素のタプルまたはリストで、アフィン変換行列の値を含みます。この行列は `[a, b, d, e, xoff, yoff]` の形をとり、2x3の変換行列 `[a, b, xoff], [d, e, yoff]` に相当します。

普通は、２次元のアフィン変換は、2x3の変換行列ですが、Shapelyでは一行6列の配列で渡すことだけ注意が必要です。

### サンプルコード

以下のサンプルコードでは、四角形オブジェクトに対して、様々なAffine変換を適用し、その結果をプロットしています。

```python
import matplotlib.pyplot as plt
import numpy as np
import shapely.geometry as geometry
from shapely.affinity import affine_transform, rotate, scale, skew, translate

# 四角形の幅と高さを設定
width = 5
height = 10

# 四角形を作成
pixel = geometry.box(0, 0, width, height)

# Affine変換を適用
angle = 45.0  # 回転角度
matrix = [
    np.cos(np.radians(angle)), -np.sin(np.radians(angle)),
    np.sin(np.radians(angle)), np.cos(np.radians(angle)),
    5.0, 10.0
]
affine_pixel = affine_transform(pixel, matrix)

# 他の変換を適用
rotate_pixel = rotate(pixel, angle, origin=(10, 20))
scale_pixel = scale(pixel, xfact=0.8, yfact=2.0, origin=(-10, -20))
skew_pixel = skew(pixel, xs=0, ys=30, origin=(-50, 200))
translate_pixel = translate(pixel, xoff=-10.0, yoff=30.0)

# Matplotlibでプロット
fig, ax = plt.subplots()
ax.fill(*pixel.exterior.xy, alpha=0.5, label="original")
ax.fill(*affine_pixel.exterior.xy, alpha=0.5, label="affine")
ax.fill(*rotate_pixel.exterior.xy, alpha=0.5, label="rotate")
ax.fill(*scale_pixel.exterior.xy, alpha=0.5, label="scale")
ax.fill(*skew_pixel.exterior.xy, alpha=0.5, label="skew")
ax.fill(*translate_pixel.exterior.xy, alpha=0.5, label="translate")

ax.set_aspect('equal', adjustable='box')
plt.grid(alpha=0.3)
plt.legend()
plt.show()
``` 

このコードにより得られた結果がこちらです。

![スクリーンショット 2023-12-17 1.14.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/1470a68c-353b-f09c-02d8-a2ccba69f9b1.png)

`skew` の平行移動の量は、`origin=(-50, 200)` のはずですが、その通りではないような気がします。。
`rotate` は、この定義を用いると、反時計回りに正の角度で指定されています。

# アフィン変換についての補足

## 同次座標系とは？

同次座標系は、通常の座標系に1つの追加次元を加えることで、点を表現します。2次元空間の点 (x, y) は同次座標系では (x, y, 1) となります。この追加された次元は、主に変換の計算を容易にするために用いられます。

### なぜ3x3行列が必要なのか？

2次元空間の点に対するアフィン変換（平行移動、回転、拡大/縮小、せん断）を行うためには、同次座標を使用してこれらの変換を一つの行列で表現することが望ましいです。

アフィン変換行列は以下のように表されます：

```
[ a  b  tx ]
[ c  d  ty ]
[ 0  0  1  ]
```

ここで、`(a, b, c, d)`は回転、拡大/縮小、せん断などの線形変換を表し、`(tx, ty)`は平行移動を表します。最後の行 `[0 0 1]` は同次座標の形式を維持します。

２次元平面上の点 `(x, y, 1)` にこの変換行列を適用すると、以下のようになります：

```
[ a  b  tx ]   [ x ]
[ c  d  ty ] * [ y ]
[ 0  0  1  ]   [ 1 ]
```

この計算を行うと、点 `(x, y)` に対して回転、拡大/縮小、せん断、平行移動を一度に適用することができます。3x3の行列を使用することで、これらの変換を単一の行列演算で行うことが可能になるのです。これにより、複数の変換を効率的に組み合わせることができます。

# 応用例

これを用いて、四角の画素をタイル状に並べる例を紹介します。

## 基本ジオメトリの設計

まずは、基本ジオメトリを作成します。この例では、四角の画素に、配線の足が２本付いているような形を基本とします。

``` python:
import matplotlib.pyplot as plt
import numpy as np
from shapely.affinity import affine_transform, rotate, scale, skew, translate
from shapely.geometry import Polygon

# 使用例
width = 8
height = 10
wire_width = 1
leg1_offset = 2

# make global geometry, and place its left bottom edge at (0,0)
g_pixel = Polygon([(0, 0), (width, 0), (width, height), (0, height)])
g_leg1 = Polygon([(0, 0), (width + leg1_offset, 0), (width + leg1_offset, wire_width), (0, wire_width)])
g_leg1side = Polygon([(width + leg1_offset - wire_width, 0), (width + leg1_offset, 0),\
                    (width + leg1_offset, height - 2*wire_width), (width + leg1_offset - wire_width, height - 2*wire_width)])
g_leg2 = Polygon([(0, height - wire_width), (width + leg1_offset, height - wire_width), (width + leg1_offset, height), (0, height)])
g_leg2side = Polygon([(width + leg1_offset - wire_width, height), (width + leg1_offset - wire_width, 2*wire_width),\
                    (width + leg1_offset, 2*wire_width), (width + leg1_offset, height)])

# move the center of the geometry to (0,0)
g_pixel = translate(g_pixel, xoff=-width/2, yoff=-height/2)
g_leg1 = translate(g_leg1, xoff=-width/2, yoff=-height/2)
g_leg1side = translate(g_leg1side, xoff=-width/2, yoff=-height/2)
g_leg2 = translate(g_leg2, xoff=-width/2, yoff=-height/2)
g_leg2side = translate(g_leg2side, xoff=-width/2, yoff=-height/2)

fig, ax = plt.subplots()

def plot_geo(fig, ax, polys, labels="original"):
    # polysがリストでない場合は、リストに変換
    if not isinstance(polys, list):
        polys = [polys]
    if not isinstance(labels, list):
        labels = [labels]

    for poly, label in zip(polys, labels):
        ax.fill(*poly.exterior.xy, alpha=0.5, label=label)

plot_geo(fig, ax, [g_pixel, g_leg1, g_leg1side, g_leg2, g_leg2side], ["pixel", "leg1", "leg1side", "leg2", "leg2side"])

ax.set_aspect('equal', adjustable='box')
plt.grid(alpha=0.3)
plt.legend()
plt.savefig("shapely_check_affine_2.png")
plt.show()
```

![スクリーンショット 2023-12-25 11.14.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/41f02ff9-6441-acd1-258e-1adef15fe53f.png)


## 画素をタイル張りする方法

タイル状に画素を配置する例を示します。現実的な理由によって対称性が崩れる場合を想定して、４象限ごとにジオメトリーを微調整して生成する方法を紹介します。

``` python:
def gen_pixel(matrix, xoff = 0, yoff = 0.0, legside="left"):
    trans_pixel = affine_transform(g_pixel, matrix)
    trans_leg1 = affine_transform(g_leg1, matrix)
    trans_leg2 = affine_transform(g_leg2, matrix)

    if legside == "left":
#        print("left")
        trans_legside = affine_transform(g_leg2side, matrix)
    elif legside == "right":
#        print("right")
        trans_legside = affine_transform(g_leg1side, matrix)
    else:
        print("ledside is either left or right.")

    return [trans_pixel, trans_leg1, trans_leg2, trans_legside]

pixel_lists = []

legside = "left" # leg side

nrows = 3
ncols = 3

angle = 0.0  # 回転角度
for i in range(nrows):
    for j in range(ncols):
        xoff = j * (height + leg1_offset*2) + width/2 + wire_width
        yoff = i * (height + leg1_offset*2) + 3 * wire_width + height/2
        matrix = [np.cos(np.radians(angle)), -np.sin(np.radians(angle)), np.sin(np.radians(angle)), np.cos(np.radians(angle)), xoff, yoff]
        pixel_lists.append(gen_pixel(matrix, xoff = xoff, yoff = yoff, legside=legside))


angle = 180.0  # 回転角度
legside = "right" # leg side
for i in range(nrows):
    for j in range(ncols):
        xoff = - (j * (height + leg1_offset*2) + width/2 + wire_width)
        yoff = i * (height + leg1_offset*2)    + 3 * wire_width + height/2
        matrix = [np.cos(np.radians(angle)), -np.sin(np.radians(angle)), np.sin(np.radians(angle)), np.cos(np.radians(angle)), xoff, yoff]
        pixel_lists.append(gen_pixel(matrix, xoff = xoff, yoff = yoff, legside=legside))

angle = 0.0  # 回転角度
legside = "right" # leg side
for i in range(nrows):
    for j in range(ncols):
        xoff = j * (height + leg1_offset*2) + width/2 + wire_width
        yoff = - (i * (height + leg1_offset*2)    + 3 * wire_width + height/2)
        matrix = [np.cos(np.radians(angle)), -np.sin(np.radians(angle)), np.sin(np.radians(angle)), np.cos(np.radians(angle)), xoff, yoff]
        pixel_lists.append(gen_pixel(matrix, xoff = xoff, yoff = yoff, legside=legside))

angle = 180.0  # 回転角度
legside = "right" # leg side
for i in range(nrows):
    for j in range(ncols):
        xoff = - (j * (height + leg1_offset*2) + width/2 + wire_width)
        yoff = - (i * (height + leg1_offset*2)    + 3 * wire_width + height/2)
        matrix = [np.cos(np.radians(angle)), -np.sin(np.radians(angle)), np.sin(np.radians(angle)), np.cos(np.radians(angle)), xoff, yoff]
        pixel_lists.append(gen_pixel(matrix, xoff = xoff, yoff = yoff, legside=legside))

# Matplotlibを使用したプロット
fig, ax = plt.subplots()

def plot_geo(fig, ax, polys, labels="original"):
    # polysがリストでない場合は、リストに変換
    if not isinstance(polys, list):
        polys = [polys]
    if not isinstance(labels, list):
        labels = [labels]

    for poly, label in zip(polys, labels):
        ax.fill(*poly.exterior.xy, alpha=0.5, label=label)

for onepixel in pixel_lists:
    plot_geo(fig, ax, onepixel, ["pixel", "leg1", "leg1side", "leg2", "leg2side"])

ax.set_aspect('equal', adjustable='box')
plt.grid(alpha=0.3)
#plt.legend()
plt.savefig("shapely_check_affine_2.png")
plt.show()
```

![スクリーンショット 2023-12-25 11.18.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/44fe09d7-b27b-add3-7097-d4b39a3c691a.png)

このような、少し形が異なる形状をタイル状に配置する例になります。

# 最後に

このように、Shapelyのアフィン変換を上手に使いことで、

https://qiita.com/yamadasuzaku/items/2b22f22d46e5fa855f50

など、gdstk を用いた GDSファイルの生成も効率的にできるはずです。