はじめに

この記事では、ある点がポリゴン内に存在するかどうかを判定する方法について説明します。具体的な実装コードを基に、それぞれの関数やステップを説明していきます。

実行例は、google colab でも確認できます。

Point-in-Polygon (PIP) アルゴリズム

Point-in-Polygon (PIP) アルゴリズムは、与えられた点が多角形の内部にあるか、外部にあるか、境界上にあるかを判定するためのアルゴリズムです。このアルゴリズムは、GIS (地理情報システム) やコンピュータグラフィックスなどの分野で一般的に使用されます。

以下は、PIP アルゴリズムの一つである「Ray Casting（光線投射）アルゴリズム」の基本的な考え方です：

与えられた点から任意の方向に無限に伸びる半直線（光線）を投射します。
この光線が多角形の辺を何回交差するかをカウントします。
もし交差回数が奇数回であれば、点は多角形の内部にあります。
もし交差回数が偶数回であれば、点は多角形の外部にあります。

この方法が有効である理由は、点が多角形の外部にある場合、任意の方向に投射される光線は偶数回交差することになります（多角形の内部に入るたびと外に出るたびの2回でカウントされる）。しかし、点が多角形の内部にある場合、光線は奇数回交差します（最後の交差で多角形を出るため）。

ただし、光線が多角形の頂点や辺上に厳密に一致する場合などの特殊なケースを処理するための追加の手順やチェックが必要です。

Ray Casting アルゴリズム以外にも、Winding Number アルゴリズムなど、他の方法で PIP を解決するアルゴリズムもあります。

実装コードの概要

ポリゴンの頂点を生成
頂点をポリゴンの重心を基に反時計回りにソート
ランダムに生成された点がポリゴン内に存在するかを判定
結果をmatplotlibで可視化

1. ポリゴンの重心を計算する

def calculate_centroid(x, y):
    """Calculates the centroid (average) of a set of points."""
    return sum(x) / len(x), sum(y) / len(y)

この関数は、ポリゴンの頂点集合の重心（平均）を計算します。

2. 重心を基にポリゴンの頂点をソートする

def sort_points_based_on_centroid(x, y):
    """Sort points in counter-clockwise order based on angles from centroid."""
    centroid_x, centroid_y = calculate_centroid(x, y)

    def angle_from_centroid(point):
        """Calculate the angle between a point and the centroid."""
        return (math.atan2(point[1] - centroid_y, point[0] - centroid_x) + 2 * math.pi) % (2*math.pi)

    sorted_points = sorted(zip(x, y), key=angle_from_centroid)
    sorted_x, sorted_y = zip(*sorted_points)
    return list(sorted_x), list(sorted_y)

この関数は、ポリゴンの頂点をその重心に対する角度を基に反時計回りにソートします。これは、ポリゴンを描画する際に正しい順番で頂点を結ぶために必要です。

3. 点がポリゴン内に存在するかを判定する

def are_points_inside_polygon(x, y, poly_x, poly_y):
    """Determine if points lie inside the given polygon."""
    num = len(poly_x)
    inside = np.full(x.shape, False, dtype=bool)

    p1x, p1y = poly_x[0], poly_y[0]
    for i in range(1, num + 1):
        p2x, p2y = poly_x[i % num], poly_y[i % num]
        mask = ((y > min(p1y, p2y)) & (y <= max(p1y, p2y)) & (x <= max(p1x, p2x)) & (p1y != p2y))
        xinters = (y[mask] - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
        inside[mask] = np.logical_xor(inside[mask], x[mask] < xinters)
        p1x, p1y = p2x, p2y

    return inside

この関数は、与えられた点がポリゴン内に存在するかどうかを判定します。この判定は、線分交差テストに基づいて行われます。

4. 結果を可視化する

matplotlibを使用して、ポリゴンとランダムに生成された点をプロットします。そして、点がポリゴン内に存在するかどうかに基づいて色を付けることで、結果を可視化します。

plt.fill(poly_x, poly_y, "b-", alpha=0.4)
plt.scatter(test_x[inside], test_y[inside], c="green", marker=".", label="Inside")

実行結果

上記のコードを実行すると、ポリゴンとその内外の点を色分けして表示するグラフが生成されます。これにより、点がポリゴン内部に存在するか外部に存在するかを簡単に確認することができます。

コードの全文

qiita_point_in_polygon_fromrandom.py


#!/usr/bin/env python

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

def calculate_centroid(x, y):
    """Calculate the centroid of a set of points."""
    return sum(x) / len(x), sum(y) / len(y)

def sort_points_based_on_centroid(x, y):
    """Sort points in counter-clockwise order based on angles from centroid."""
    centroid_x, centroid_y = calculate_centroid(x, y)

    def angle_from_centroid(point):
        """Calculate the angle between a point and the centroid."""
        return (math.atan2(point[1] - centroid_y, point[0] - centroid_x) + 2 * math.pi) % (2*math.pi)

    sorted_points = sorted(zip(x, y), key=angle_from_centroid)
    sorted_x, sorted_y = zip(*sorted_points)
    return list(sorted_x), list(sorted_y)

def are_points_inside_polygon(x, y, poly_x, poly_y):
    """Determine if points lie inside the given polygon."""
    num = len(poly_x)
    inside = np.full(x.shape, False, dtype=bool)

    p1x, p1y = poly_x[0], poly_y[0]
    for i in range(1, num + 1):
        p2x, p2y = poly_x[i % num], poly_y[i % num]
        mask = ((y > min(p1y, p2y)) & (y <= max(p1y, p2y)) & (x <= max(p1x, p2x)) & (p1y != p2y))
        xinters = (y[mask] - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
        inside[mask] = np.logical_xor(inside[mask], x[mask] < xinters)
        p1x, p1y = p2x, p2y

    return inside

# Define random polygon
np.random.seed(42)  # Set seed for reproducibility
poly_x = np.random.uniform(-3, 3, 30)
poly_y = np.random.uniform(-3, 3, 30)
poly_x, poly_y = sort_points_based_on_centroid(poly_x, poly_y)

# Plot boundaries and points
plt.fill(poly_x, poly_y, "b-", alpha=0.4)  # ポリゴンの塗りつぶし
plt.plot(poly_x, poly_y, "co")  # ポリゴンの頂点

# Generate Random points
np.random.seed(0)
test_x = np.random.uniform(-6, 6, 1000)
test_y = np.random.uniform(-6, 6, 1000)

# Test points in or outside polygon
inside = are_points_inside_polygon(test_x, test_y, poly_x, poly_y)

# Plot results
plt.scatter(test_x[inside], test_y[inside], c="green", marker=".", label="Inside")
plt.scatter(test_x[~inside], test_y[~inside], c="red", marker=".", label="Outside")

plt.legend()
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Point in Polygon Test")
plt.grid(True)
plt.show()

実行結果

高速化

基本的に境界も判定する点も numpy を使うことで高速化を試みています。

まとめ

以上が、ポリゴン内の点判定の解説です。さまざまな地理的な問題やグラフィックスの点と領域問題などに役立てば幸いです。