凹型領域内の一様乱数を用いた点のばら撒き方

項番	目次
1	目的
2	結論
3	手順
4	最後に

目的

ゲームフィールド内のランダムな地点にアクションを起こしたかったため, 問題解決を一般化したアルゴリズムを作成しました.
また、凹型に対応した一様分布アルゴリズムが記載された日本語ドキュメントを見つけられなかったため、この記事を公開しようと考えました.

結論

多角形を三角形分割するアルゴリズムと, 三角形内の一様分布アルゴリズムを用いることで, 凹凸対応の多角形一様分布を実現しました.

コード全文

※テストコードのため雑に記述されています

import numpy as np
from math import pi, atan2
import math
import matplotlib.pyplot as plt
from math import atan2, pi
    
def sort_vertices(vertices):
    center = tuple(map(lambda x: sum(x) / len(vertices), zip(*vertices)))
    vertex_angles = []
    for vertex in vertices:
        x, y = vertex
        dx, dy = x - center[0], y - center[1]
        angle = atan2(dy, dx)
        vertex_angles.append((angle, vertex))
    vertex_angles.sort()
    return [vertex for angle, vertex in vertex_angles]

def ear_clipping(p):
    polygon = p.copy()
    triangles = []
    trialCount = 0
    errLimit = 1000

    if len(polygon) < 3:
        return triangles

    while len(polygon) > 3:
        for i in range(len(polygon)):
           prev_idx = (i - 1 + len(polygon)) % len(polygon)
           next_idx = (i + 1) % len(polygon)
           prev_point = polygon[prev_idx]
           point = polygon[i]
           next_point = polygon[next_idx]

           # イアーの条件を満たさない場合はスキップ
           if not is_convex_vertex(prev_point, point, next_point, polygon): continue
           triangle = [prev_point, point, next_point]
           triangles.append(triangle)

           # ポリゴンからイアーを削除
           polygon.pop(i)
           break

        trialCount += 1
        if trialCount > errLimit:
            print("--- Error: ear_clipping ---")
            break

    # 残りの三角形を追加
    triangles.append(polygon)

    return triangles

def is_convex_vertex(prev_point, point, next_point, polygon):
    # 引数が存在しない場合はスキップ
    if point is None or prev_point is None or next_point is None:
        return False
    
    # 三角形の頂点が同じ場合はスキップ
    if point == prev_point or point == next_point or prev_point == next_point:
        return False
    # 三角形の内側にある頂点がないかをチェック
    for p in polygon:
        if p != prev_point and p != point and p != next_point:
            if is_point_in_triangle(p, prev_point, point, next_point):
                return False

    # 頂点の回転方向を確認
    prev_vec = (prev_point[0] - point[0], prev_point[1] - point[1])
    next_vec = (next_point[0] - point[0], next_point[1] - point[1])
    cross_product = prev_vec[0] * next_vec[1] - prev_vec[1] * next_vec[0]

    return cross_product < 0

def is_point_in_triangle(p, a, b, c):
    area_abc = triangle_area(a, b, c)
    area_pbc = triangle_area(p, b, c)
    area_apc = triangle_area(a, p, c)
    area_apb = triangle_area(a, b, p)

    return abs(area_abc - area_pbc - area_apc - area_apb) < 1e-9

def triangle_area(a, b, c):
    return 0.5 * abs(a[0]*b[1] + b[0]*c[1] + c[0]*a[1] - a[0]*c[1] - c[0]*b[1] - b[0]*a[1])

def plot(vertices, triangles, points):
    # 頂点座標のx座標とy座標をそれぞれリストに分ける
    xs, ys = zip(*vertices)
    
    # 多角形を描写
    plt.figure(figsize=(6, 6))  # 図のサイズを指定
    plt.plot(xs + (xs[0],), ys + (ys[0],), 'r-')  # 最初の点と最後の点を結ぶ
    
    # 三角形を描写
    for triangle in triangles:
        triangle_xs, triangle_ys = zip(*triangle)
        plt.fill(triangle_xs + (triangle_xs[0],), triangle_ys + (triangle_ys[0],), alpha=0.3)

    # 頂点に円を描画
    for x, y in vertices:
        plt.plot(x, y, 'ro')  # 赤い円を描画
    
    # 軸の範囲を設定
    plt.xlim(min(xs) - 1, max(xs) + 1)
    plt.ylim(min(ys) - 1, max(ys) + 1)

    # ランダムな点を描写
    random_xs, random_ys = zip(*points)
    plt.plot(random_xs, random_ys, 'bo', markersize=1)

    plt.show()

def get_polygon_area(vertices):
    area = 0
    n = len(vertices)
    for i in range(n):
        j = (i + 1) % n
        area += vertices[i][0] * vertices[j][1]
        area -= vertices[j][0] * vertices[i][1]
    area = abs(area) / 2
    return area

def get_random_triangle(triangles, total_area):
    random_area = np.random.uniform(0, total_area)
    current_area = 0
    for triangle in triangles:
        area = triangle_area(triangle[0], triangle[1], triangle[2])
        current_area += area
        if current_area >= random_area:
            return triangle

def get_random_point_in_triangle(triangle):
    a, b, c = triangle
    s, t = np.random.uniform(0, 1, 2)
    if s + t >= 1:
        s, t = 1 - s, 1 - t
    x = a[0] + s * (b[0] - a[0]) + t * (c[0] - a[0])
    y = a[1] + s * (b[1] - a[1]) + t * (c[1] - a[1])
    return x, y

# ランダムな多角形を生成
num_vertices = np.random.randint(10, 40)  # 3~10個の頂点数
vertices = [(np.random.randint(-10, 11), np.random.randint(-10, 11)) for _ in range(num_vertices)]
sorted_vertices = sort_vertices(vertices)
polygon_area = get_polygon_area(sorted_vertices)

triangles = ear_clipping(sorted_vertices)

points = []
trial_number = 10000

for _ in range(trial_number):
    # ランダムな三角形を取得
    triangle = get_random_triangle(triangles, polygon_area)
    try:
        # 三角形内のランダムな点を取得
        point = get_random_point_in_triangle(triangle)
        points.append(point)
    except:
        print("--- Error: get_random_point_in_triangle ---")
        continue

# 多角形を描写
plot(sorted_vertices, triangles, points)


# ランダムな多角形を生成
num_vertices = np.random.randint(10, 40)  # 3~10個の頂点数
vertices = [(np.random.randint(-10, 11), np.random.randint(-10, 11)) for _ in range(num_vertices)]
sorted_vertices = sort_vertices(vertices)
polygon_area = get_polygon_area(sorted_vertices)

triangles = ear_clipping(sorted_vertices)

points = []
trial_number = 10000

for _ in range(trial_number):
    # ランダムな三角形を取得
    triangle = get_random_triangle(triangles, polygon_area)
    point = get_random_point_in_triangle(triangle)
    points.append(point)

# 多角形を描写
plot(sorted_vertices, triangles, points)

手順

1. 三角形分割

耳刈り取り法を用いました.
wikiによると Ο(nlogn) のアルゴリズムがあるようですが, 処理が複雑であったため, Ο(n^2) の耳刈り取り法を用いました.
また私のシステム上 Ο(n^2) でも問題がなかったためです.
各自のシステムに応じて, 編集してください.

コード

def ear_clipping(p):
    polygon = p.copy()
    triangles = []
    ears = []

    if len(polygon) < 3:
        return triangles

    for i in range(len(polygon)):
        prev_idx = (i - 1) % len(polygon)
        next_idx = (i + 1) % len(polygon)
        prev_point = polygon[prev_idx]
        point = polygon[i]
        next_point = polygon[next_idx]

        is_ear = is_convex_vertex(prev_point, point, next_point, polygon)
        if is_ear:
            ears.append(i)

    while len(polygon) > 3:
        if len(ears) == 0:
            break

        ear_idx = ears.pop(0)
        prev_idx = (ear_idx - 1 + len(polygon)) % len(polygon)
        next_idx = (ear_idx + 1) % len(polygon)
        triangle = [polygon[prev_idx], polygon[ear_idx], polygon[next_idx]]
        triangles.append(triangle)
        polygon.pop(ear_idx)

        # 新しい耳の検出
        new_ears = []
        for i in range(len(polygon)):
            prev_idx = (i - 1) % len(polygon)
            next_idx = (i + 1) % len(polygon)
            prev_point = polygon[prev_idx]
            point = polygon[i]
            next_point = polygon[next_idx]

            if is_convex_vertex(prev_point, point, next_point, polygon):
                new_ears.append(i)
        ears = new_ears

    triangles.append(polygon)

    return triangles

def is_convex_vertex(prev_point, point, next_point, polygon):
    # 三角形の内側にある頂点がないかをチェック
    for p in polygon:
        if p != prev_point and p != point and p != next_point:
            if is_point_in_triangle(p, prev_point, point, next_point):
                return False

    # 頂点の回転方向を確認
    prev_vec = (prev_point[0] - point[0], prev_point[1] - point[1])
    next_vec = (next_point[0] - point[0], next_point[1] - point[1])
    cross_product = prev_vec[0] * next_vec[1] - prev_vec[1] * next_vec[0]

    return cross_product < 0

def is_point_in_triangle(p, a, b, c):
    area_abc = triangle_area(a, b, c)
    area_pbc = triangle_area(p, b, c)
    area_apc = triangle_area(a, p, c)
    area_apb = triangle_area(a, b, p)

    return abs(area_abc - area_pbc - area_apc - area_apb) < 1e-9

def triangle_area(a, b, c):
    return 0.5 * abs(a[0]*b[1] + b[0]*c[1] + c[0]*a[1] - a[0]*c[1] - c[0]*b[1] - b[0]*a[1])

2. 面積比に応じた三角形の取得

耳刈り取り法で取得したいくつかの三角形からランダムに一つを取得します.
面積に応じて取得確率を調整しています.

コード

def get_polygon_area(vertices):
    area = 0
    n = len(vertices)
    for i in range(n):
        j = (i + 1) % n
        area += vertices[i][0] * vertices[j][1]
        area -= vertices[j][0] * vertices[i][1]
    area = abs(area) / 2
    return area

def get_random_triangle(triangles, total_area):
    random_area = np.random.uniform(0, total_area)
    current_area = 0
    for triangle in triangles:
        area = triangle_area(triangle[0], triangle[1], triangle[2])
        current_area += area
        if current_area >= random_area:
            return triangle

3. 三角形内の一様乱数

お好きな方法をお使いください.

コード

def get_random_point_in_triangle(triangle):
    a, b, c = triangle
    s, t = np.random.uniform(0, 1, 2)
    if s + t >= 1:
        s, t = 1 - s, 1 - t
    x = a[0] + s * (b[0] - a[0]) + t * (c[0] - a[0])
    y = a[1] + s * (b[1] - a[1]) + t * (c[1] - a[1])
    return x, y

使用例

points = [] # 打点の座標
trial_number = 10000 # 打点の数
num_vertices = np.random.randint(10, 40)  # 3~10個の頂点数
vertices = [(np.random.randint(-10, 11), np.random.randint(-10, 11)) for _ in range(num_vertices)]
sorted_vertices = sort_vertices(vertices) # 辺が交わらないようにするため整列させる
polygon_area = get_polygon_area(sorted_vertices) # 多角形の面積
triangles = ear_clipping(sorted_vertices) # 分割された三角形

for _ in range(trial_number):
    # ランダムな三角形を取得
    triangle = get_random_triangle(triangles, polygon_area)
    point = get_random_point_in_triangle(triangle)
    points.append(point)

# 結果を描写
plot(sorted_vertices, triangles, points)

最後に

変数名にgetがあると蕁麻疹が出る方は病院に行ってお薬をもらってきてください. もちろん精神科です.
最初に与えられている頂点が並び替えられている点からわかるように, 交点が存在する多角形や穴開き多角形には対応していません.