python，matplotlibにおけるボロノイ図の描画について

Question

特定の図形内に指定した点によるボロノイ図を作成したい

pythonのライブラリであるmatplotlibにて，ボロノイ図の描画を行おうと思っております．

用意している情報はポリゴンの座標リストとボロノイ図作成の際に必要な母点の座標リストです．

発生している問題・エラーメッセージ

shapely.errors.TopologicalError: The operation 'GEOSIntersection_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x0000029A93E85940>

該当のソースコード

from scipy.spatial import Delaunay, delaunay_plot_2d, Voronoi, voronoi_plot_2d
import matplotlib.pyplot as plt
import matplotlib.patches as patch
from matplotlib.collections import PolyCollection
from scipy.spatial import Voronoi, voronoi_plot_2d
from shapely.geometry import Polygon
import numpy as np
import math

data = np.array([[11.322,5.82],
                [64.808,6.401],
                [50.481,35.431],
                [53.738,50.02],
                [52.727,76.305],
                [29.215,72.951],
                [24.786,55.766],
                [10.013,45.077]
                ])

 pt = np.array([[15.22142778, 38.74503721],
                 [26.01477792, 38.86749356],
                 [30.86884446, 48.17364433],
                 [40.64079104, 61.38215327],
                 [36.53142362, 10.81861143],
                 [20.89040778, 19.60481585],
                 [52.0941051, 20.63233987],
                 [46.06085561, 70.15568748],
                 [31.48504392, 19.81353478],
                 [15.90721556, 10.57598896],
                 [26.01360206, 28.55506018],
                 [48.96167152, 55.27385493],
                 [41.17787351, 48.50723961],
                 [46.83406199, 11.28722635],
                 [33.1668619, 68.509277  ],
                 [57.14529894, 11.43505146],
                 [41.78480506, 20.3275005 ],
                 [30.74029729, 58.48498758],
                 [26.21789971, 10.77473159],
                 [46.62100914, 39.72495837],
                 [46.67382691, 29.41005014],
                 [36.31321427, 39.41182202],
                 [20.54161943, 47.6073278 ],
                 [15.56658469, 28.43798703],
                 [36.3599739, 29.09789437]])

rudius = 5.196

def bounded_voronoi(bnd, pnts):
    """
    有界なボロノイ図を計算・描画する関数．
    """

    # すべての母点のボロノイ領域を有界にするために，ダミー母点を3個追加
    gn_pnts = np.concatenate([pnts, np.array([[100, 100], [100, -100], [-100, 0]])])

    # ボロノイ図の計算
    vor = Voronoi(gn_pnts)

    # 分割する領域をPolygonに
    bnd_poly = Polygon(bnd)

    # 各ボロノイ領域をしまうリスト
    vor_polys = []

    # ダミー以外の母点についての繰り返し
    for i in range(len(gn_pnts) - 3):

        # 閉空間を考慮しないボロノイ領域
        vor_poly = [vor.vertices[v] for v in vor.regions[vor.point_region[i]]]
        # 分割する領域をボロノイ領域の共通部分を計算
        i_cell = bnd_poly.intersection(Polygon(vor_poly))

        # 閉空間を考慮したボロノイ領域の頂点座標を格納
        vor_polys.append(list(i_cell.exterior.coords[:-1]))


    # ボロノイ図の描画
    fig = plt.figure(figsize=(7, 6))
    ax = fig.add_subplot(111)

    # 母点
    ax.scatter(pnts[:,0], pnts[:,1])

    # ボロノイ領域
    poly_vor = PolyCollection(vor_polys, edgecolor="black",
                              facecolors="None", linewidth = 1.0)
    ax.add_collection(poly_vor)

    xmin = np.min(bnd[:,0])
    xmax = np.max(bnd[:,0])
    ymin = np.min(bnd[:,1])
    ymax = np.max(bnd[:,1])

    ax.set_xlim(xmin-0.1, xmax+0.1)
    ax.set_ylim(ymin-0.1, ymax+0.1)
    ax.set_aspect('equal')

    fig.savefig('voronoi.pdf')

    return vor_polys

# ボロノイ図の計算・描画
vor_polys = bounded_voronoi(data, pt)
print(vor_polys)

試したこと

具体的なコードはキータ(URL：https://qiita.com/supbon2/items/30e0cb49c9338e721b8c)を参考にさせていただきました．
基本的には描画も可能ではあるのですが，母点の座標データを15個以上に設定した場合にのみこのようなエラーが発生するようです．

上記のコードは25個の座標が指定してありますが，この数を14個以下に変更した際は想定通り動いていたので間違いないと思われます．

補足情報

他にも必要な情報等ございましたらコメントにてよろしくお願いします。

Answer 1

• すべてのボロノイ領域を有界にするために，ダミーの母点を3個追加する．

これが実現できていないからです。
全ての点がダミーの母点でできる三角形の中に入らなければならないのですが、
[33.1668619, 68.509277 ]は直線y = 0.5x + 50よりも上にありますね。

Pythonで有界な（閉じた）ボロノイ図を計算・描画するの場合は、全ての点の座標が0と１の間にあるので、ダミーの母点でできる三角形の中に入っています。

    gn_pnts = np.concatenate([pnts, np.array([[100, 100], [100, -100], [-100, 0]])])

を```python
gn_pnts = np.concatenate([pnts, np.array([[1000, 1000], [1000, -1000], [-1000, 0]])])

に変更しましょう。

Comments

[narimi]

なるほど，，，
最初に設定した母点から成る三角形内に収まっていないのがたまたま15個目だったってことですね．

とても分かりやすかったです．
ありがとうございました！