Python で凸包を計算して表示する

凸包とは、点の集合を含む最小の凸多角形です。PythonではSciPyのConvexHullで凸包を計算し、Matplotlibで境界を描画できます。アルゴリズムを学びたい場合は、小規模なQuickhullの自作実装も役立ちます。

凸包とは？

平面上の点に対する凸包は、すべての点を含む最小の凸多角形です。点の集合をぴんと張った輪ゴムで囲む様子を想像してください。輪ゴムが触れる点が凸包の頂点で、輪ゴムの内側にある点は境界には含まれません。多角形の辺は、外側の頂点を順番につなぎます。

ここでいう凸包は点の集合に対するもので、すでに凸である物体の境界を指すものではありません。凸包は幾何学、衝突判定の境界、空間解析などで役立ちますが、点の内部での配置は意図的に捨てられます。

SciPyで凸包を計算して描画する

スクリプトを実行するPython環境にNumPy、SciPy、Matplotlibをインストールします。たとえば、python -m pip install numpy scipy matplotlibを実行してください。ConvexHullには、各行が1点を表し、2列にx座標とy座標を格納した配列を渡します。vertices属性には境界のインデックスが順番に格納されます。この順序で点を結び、先頭の点を末尾にもう一度加えると、多角形を閉じられます。

次の例を実行します。

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([
    [0.0, 0.0], [1.0, 0.0], [2.0, 0.0],
    [2.0, 1.0], [2.0, 2.0], [1.0, 2.0],
    [0.0, 2.0], [0.0, 1.0], [1.0, 1.0],
])
hull = ConvexHull(points)
boundary = points[hull.vertices]
closed_boundary = np.vstack((boundary, boundary[0]))

plt.scatter(points[:, 0], points[:, 1], label="Points")
plt.plot(closed_boundary[:, 0], closed_boundary[:, 1], "r-", label="Hull")
plt.gca().set_aspect("equal", adjustable="box")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

内側の点は入力に含まれていますが、凸包には含まれません。ここではhull.simplicesを使って辺を1本ずつ描くより、hull.verticesを使う方が適しています。順序付きの頂点を使えば、閉じた輪郭を1本の線として描けます。この例を実行するとMatplotlibのウィンドウが開きます。ノートブックでは、図が出力領域に表示されます。図中のラベル「Points」と「Hull」は英語のまま表示されます。

PythonでQuickhullを実装する

Quickhullは、2次元の凸包を求める再帰的な方法です。まず、x座標が最小と最大の点を選びます。この2点が候補となる境界を2つの側に分けます。向きを持つ各辺について、アルゴリズムは正の外積が最大となる点を選びます（この値は垂直距離に比例します）。次に、その点の両側にある2つの小さな点集合を調べます。線分の左側に候補点が残っていなければ、その線分は境界の一部です。厳密な側判定によって、選んだ点は後続の再帰呼び出しから除かれるため、有限個の相異なる点であれば再帰は終了します。

次の実装は、2次元座標ペアの反復可能オブジェクトを受け取り、最初の頂点を重複させずに、境界に沿った順序で多角形の頂点を返します。再帰探索の前に重複する座標を取り除きます。相異なる点が3つ未満の場合や、すべての点が一直線上にある場合は、利用可能な両端の点を返します。この例はPython標準ライブラリだけを使うため、単一のスクリプトとして実行できます。返されたリストは、後でプロット用ライブラリに渡せます。

次の例を実行します。

def cross(a, b, point):
    """Return the signed 2D cross product of AB and A-to-point."""
    return (b[0] - a[0]) * (point[1] - a[1]) - (b[1] - a[1]) * (point[0] - a[0])


def quickhull(points):
    """Return the vertices of a 2D convex hull in boundary order."""
    unique = sorted(set(tuple(point) for point in points))
    if len(unique) <= 2:
        return unique

    left = min(unique, key=lambda point: (point[0], point[1]))
    right = max(unique, key=lambda point: (point[0], point[1]))

    def chain(start, end, candidates):
        candidates = [point for point in candidates if cross(start, end, point) > 0]
        if not candidates:
            return []
        farthest = max(candidates, key=lambda point: cross(start, end, point))
        first = [point for point in candidates if cross(start, farthest, point) > 0]
        second = [point for point in candidates if cross(farthest, end, point) > 0]
        return chain(start, farthest, first) + [farthest] + chain(farthest, end, second)

    upper = chain(left, right, unique)
    lower = chain(right, left, unique)
    return [left] + upper + [right] + lower


points = [
    (0, 0), (1, 0), (2, 0), (2, 1), (2, 2),
    (1, 2), (0, 2), (0, 1), (1, 1),
]
hull_points = quickhull(points)
closed_hull = hull_points + hull_points[:1]
print(hull_points)

外積の符号は、点が向き付き直線のどちら側にあるかを示し、その大きさは同じ側にある候補点の順位付けに使われます。再帰呼び出しでは、新しい2つの線分の外側にある点だけを残します。返されるリストは多角形の境界に沿った順序になっているため、閉じた線を描くには先頭の点を末尾に追加します。この学習用実装は2次元の点を対象としており、本番環境で検証済みの幾何学ライブラリの代わりにはなりません。また、単純な再帰処理のため、点の分布によっては計算量が大きくなる場合があります。

方法の選び方と特殊な入力への対処

保守されたライブラリの処理、凸包のインデックス、または2次元を超える入力への対応が必要なら、SciPyを使います。再帰的な分割の仕組みを理解したい場合や、小規模な2次元処理を応用したい場合は、自作のQuickhullを使います。速度が優れていると決めつけて選ぶのは避けてください。処理速度はデータや実行環境によって異なり、この記事ではベンチマーク結果を示していません。

入力の形状と幾何学的な配置が重要です。SciPyで2次元の多角形を作るには、少なくとも3つの非共線な点が必要です。点が少なすぎる場合やすべての点が一直線上にある場合は、そのような多角形を定義できず、Qhullは幾何学的なエラーを報告します。アプリケーションでこの状態に対処する必要がある場合は、scipy.spatial.QhullErrorを捕捉するか、ConvexHullを呼び出す前に入力を確認してください。重複した点は境界の頂点を増やしません。上の自作関数は重複を取り除き、共線な入力では利用可能な両端の点を返しますが、有限の数値座標ペアを前提としています。信頼できないデータを渡す前に検証してください。3次元以上の凸包には、この2次元の自作アルゴリズムではなく、SciPyのライブラリ実装を使います。

よくある質問

Quickhullはどのように境界を見つけますか？

両端の点から始め、各側で最も遠い候補点を選び、新しい線分の外側にある点集合を再帰的に調べます。候補点が残っていない側は、凸包の辺になります。

3Dでも凸包を計算できますか？

はい。SciPyのConvexHullは、座標列が2つを超える点配列に対応しています。ここで紹介した自作実装は2次元の外積を使うため、3次元の入力には対応していません。

点が一直線上にある場合はどうなりますか？

2次元の多角形にはなりません。SciPyはQhullの幾何学的なエラーを報告します。自作関数は両端の2点を返します（相異なる点がそれより少なければ、返される点も少なくなります）。

どの方法を選べばよいですか？

アプリケーションのコードや高次元の入力にはSciPyを選びます。2次元アルゴリズムを学ぶには自作版を選び、より大きなシステムに応用する前に入力検証とテストを追加してください。