Shapelyで曲線に沿って幅が変化するポリゴンを生成する方法

# 可変幅のポリゴンを生成するためのPythonの小技

この記事では、[`Shapely`](https://shapely.readthedocs.io/en/stable/manual.html)を用いて曲線に沿って幅が変化するポリゴンを生成する方法を紹介します。
`NumPy`、`Shapely`、`Matplotlib`を使って、複雑な幾何学的形状を作成し、可視化する方法なので、いろいろな応用も効くかなと思います。

https://colab.research.google.com/drive/1r6qs_2QQcPYYuogP5EtHjev1dozeivWe?usp=sharing

に、この記事の実行結果を掲載しています。

下記のような[gdstk](https://heitzmann.github.io/gdstk/) と組み合わせると、GDSファイルの生成にも使えます。

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


## 可変幅のオフセットポイント生成

### ポイントオフセットの基本

中心となるのは`create_offset_points`関数です。この関数の目的は、与えられたベースラインに沿って左右対称なオフセット点を生成し、それらを結ぶことで幅が変化するポリゴンを形成することです。ここでのキーは、各ポイントにおけるベースラインの方向ベクトルに垂直なベクトル（`perpendicular_vector`関数で計算）を使用して、ポイントを左右にシフトすることです。

``` python:
def perpendicular_vector(v):
    """与えられたベクトルに垂直なベクトルを返す。"""
    return np.array([-v[1], v[0]])

def create_offset_points(base_line, start_width, stop_width, num_points=1000, xoffset=0, yoffset=0):
    """ベースラインに沿ってオフセットされた点のペアを生成し、ラインの太さを徐々に変化させる。"""
    offset_points_left = []
    offset_points_right = []
    points = [base_line.interpolate(i / num_points * base_line.length) for i in range(num_points + 1)]
    widths = [(i / num_points) * (stop_width - start_width) + start_width for i in range(num_points + 1)]
    offset_vec = np.array([xoffset, yoffset])

    for i in range(num_points):
        direction_vector = np.array((points[i + 1].x, points[i + 1].y)) - np.array((points[i].x, points[i].y))
        direction_vector /= np.linalg.norm(direction_vector)
        perp_vector = perpendicular_vector(direction_vector) * widths[i] / 2
        offset_points_left.append((np.array((points[i].x, points[i].y)) + offset_vec - perp_vector).tolist())
        offset_points_right.append((np.array((points[i].x, points[i].y)) + offset_vec + perp_vector).tolist())

    return offset_points_left, offset_points_right
```

### 幅の変化

幅は、`start_width`から`stop_width`まで線形に変化します。これは、各ポイントにおいて、ベースラインの長さに対する相対的な位置に基づいて計算されます。具体的には、ベースライン上の各ポイントにおいて、そのポイントの位置に応じた幅でオフセットされた点を計算し、これらの点を結んでポリゴンを形成します。

## 応用: ジグザグと円弧の組み合わせ

このスクリプトでは、ジグザグ形状と円弧を組み合わせてユニークなベースラインを作成し、それに沿って可変幅のポリゴンを生成しています。

### ジグザグポイント

`create_zigzag_points`関数では、ジグザグの頂点を定義し、それらを指定されたスケールで変換しています。

``` python:
def create_zigzag_points(xfactor, yfactor):
    """ジグザグポイントを生成し、スケーリングする。"""
    points = [(0, 0), (5, -40), (10, -10), (20, 30), (40, 30), (50, 60), (60, 60)]
    return list(map(lambda x: (x[0] * xfactor, x[1] * yfactor), points))
```

### 円弧ポイント

`create_circle_points`関数では、指定された半径と角度に基づいて円弧のポイントを生成しています。

``` python:
def create_circle_points(cx, cy, radius, start_angle, end_angle, num_points):
    """円弧のポイントを生成する。"""
    angles = np.linspace(start_angle, end_angle, num_points)
    x_coords = radius * np.cos(angles) + cx
    y_coords = radius * np.sin(angles) + cy
    return list(zip(x_coords, y_coords))
```

## 可視化: Matplotlibでの描画

最後に、Matplotlibを使用して生成されたポリゴンと元のベースラインを描画します。このビジュアライゼーションは、生成された形状の理解を深め、さらなる分析や調整に役立ちます。

# コード全容

## 全コード

``` python:
import numpy as np
from shapely.geometry import Polygon, LineString
import matplotlib.pyplot as plt

def perpendicular_vector(v):
    """与えられたベクトルに垂直なベクトルを返す。"""
    return np.array([-v[1], v[0]])

def create_offset_points(base_line, start_width, stop_width, num_points=1000, xoffset=0, yoffset=0):
    """ベースラインに沿ってオフセットされた点のペアを生成し、ラインの太さを徐々に変化させる。"""
    offset_points_left = []
    offset_points_right = []
    points = [base_line.interpolate(i / num_points * base_line.length) for i in range(num_points + 1)]
    widths = [(i / num_points) * (stop_width - start_width) + start_width for i in range(num_points + 1)]
    offset_vec = np.array([xoffset, yoffset])

    for i in range(num_points):
        direction_vector = np.array((points[i + 1].x, points[i + 1].y)) - np.array((points[i].x, points[i].y))
        direction_vector /= np.linalg.norm(direction_vector)
        perp_vector = perpendicular_vector(direction_vector) * widths[i] / 2
        offset_points_left.append((np.array((points[i].x, points[i].y)) + offset_vec - perp_vector).tolist())
        offset_points_right.append((np.array((points[i].x, points[i].y)) + offset_vec + perp_vector).tolist())

    return offset_points_left, offset_points_right

def create_zigzag_points(xfactor, yfactor):
    """ジグザグポイントを生成し、スケーリングする。"""
    points = [(0, 0), (5, -40), (10, -10), (20, 30), (40, 30), (50, 60), (60, 60)]
    return list(map(lambda x: (x[0] * xfactor, x[1] * yfactor), points))

def create_circle_points(cx, cy, radius, start_angle, end_angle, num_points):
    """円弧のポイントを生成する。"""
    angles = np.linspace(start_angle, end_angle, num_points)
    x_coords = radius * np.cos(angles) + cx
    y_coords = radius * np.sin(angles) + cy
    return list(zip(x_coords, y_coords))

# ジグザグと円弧のポイントを生成
zigzag_points = create_zigzag_points(1, 0.3)
circle_points = create_circle_points(60, 60, 30, np.pi, 0.5 * np.pi, 30)

# ベースラインを生成
base_line = LineString(zigzag_points + circle_points)

# オフセットポリゴンを生成
offset_polygon_1 = Polygon(create_offset_points(base_line, 1, 5)[0] + create_offset_points(base_line, 1, 5)[1][::-1])
offset_polygon_2 = Polygon(create_offset_points(base_line, 5, 1, xoffset=100, yoffset=10)[0] + create_offset_points(base_line, 5, 1, xoffset=100, yoffset=10)[1][::-1])

# 描画
fig, ax = plt.subplots()
x, y = offset_polygon_1.exterior.xy
ax.fill(x, y, alpha=0.5, fc='cyan', ec='none', label="Poly1")
x, y = offset_polygon_2.exterior.xy
ax.fill(x, y, alpha=0.5, fc='magenta', ec='none', label="Poly2")
x, y = zip(*circle_points)
ax.plot(x, y, color='red', marker='.', label="Circle")
x, y = zip(*zigzag_points)
ax.plot(x, y, color='blue', marker='.', label="Zigzag")
plt.legend()

ax.set_aspect('equal', adjustable='box')
plt.show()
```

## 実行結果

![スクリーンショット 2023-12-15 16.28.11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/d8da1599-ba74-37c3-7f34-fc211b4cccc8.png)


# まとめ

このコードは、幾何学的な問題に対する一つの解決策です。より高度なPythonのライブラリを組み合わせることで、複雑な問題にも応用できます。

## 可変幅ポリゴン生成の高度なアプローチとアルゴリズム

カーブフィッティングとスプライン補間、など補完方法は様々あります。

- **ベジエ曲線**: ベジエ曲線は、複雑な形状をスムーズに表現するのに適しています。制御点を用いて曲線を生成し、これをベースラインとして使用することで、より自然な曲線に沿ったポリゴンを生成できます。
- **B-スプライン**: B-スプライン補間を使用すると、与えられたポイントセットを通る滑らかな曲線を生成できます。これは特に、データポイントが多く、複雑な形状を生成する必要がある場合に有効です。

https://docs.scipy.org/doc/scipy/reference/interpolate.html

など、scipy にも多様な補完方法が実装されています。

その他、特定の制約条件下で最も効率的なポリゴンの生成(例えば、面積、周囲長、あるいは物理的特性を最適化)を問題に応じて考える必要があります。