вырезание диаграммы вороного python

Я вычисляю диаграмму Вороного из набора точек следующим образом:

from scipy.spatial import Voronoi
import numpy as np


np.random.seed(0)
points = np.random.uniform(-0.5, 0.5, (100, 2))
# Compute Voronoi
v = Voronoi(points)
voronoi_plot_2d(v)
plt.show()

Это создает изображение следующим образом:

Вороной

Как видно, это создает вершины, которые уходят в бесконечность (пунктирные линии), а также за пределы исходной ограничивающей рамки для точек, которая:

 bbox = np.array([[-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5]])

Что я хотел бы сделать, так это прикрепить диаграмму Вороного к этой ограничивающей рамке, т.е. спроецировать за пределы и бесконечные вершины на соответствующие места в этой ограничивающей рамке. Таким образом, вершины должны быть переставлены и спроецированы обратно в правильные точки пересечения из бесконечности или конечных вершин, но которые находятся за пределами моей области отсечения.


person Luca    schedule 17.03.2016    source источник
comment
@unutbu Хотя кажется правдой, что можно получить ответ по ссылке на раскраску диаграммы Вороного, с моей точки зрения, связь не является прямой или очевидной. Лука задавал аналогичный вопрос ранее, и в то время я интерпретировал, что ему нужны только вершины внутри ограничивающей рамки (и мой ответ в то время был для этого). На этот раз ему, кажется, нужны точки пересечения между ограничивающей рамкой и бесконечными областями Вороного. Не обязательно для сюжета. Я не считаю это дубликатом.   -  person armatita    schedule 17.03.2016
comment
Да, это как-то пересечь плоскость Вороного с выходящими за границы и бесконечными вершинами и ограничительной рамкой, заданной [[-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5]]. Я не уверен, делает ли это другой код... Я смотрю на него сейчас.   -  person Luca    schedule 17.03.2016
comment
Я не уверен, есть ли более быстрый способ, но проверьте код для пересечения двух строк (пример: stackoverflow.com/questions/20677795/). Очевидно, вам придется сделать это для каждой из вершин, которые входят в бесконечные области.   -  person armatita    schedule 17.03.2016
comment
Просто упомянем, что я написал версию C++ того, что вам нужно, она доступна здесь . Он также содержит пользовательский пакет для cgal-swig-bindings, поэтому версия Python может не быть слишком сложным (хотя я только приложил усилия, чтобы заставить его работать для java).   -  person sloriot    schedule 18.03.2016
comment
Вот аналогичный (повторяющийся?) вопрос с ответом   -  person jorgeh    schedule 11.01.2017
comment
Отвечает ли это на ваш вопрос? Получение координат ограниченного многоугольника из ячеек Вороного   -  person Gabriel    schedule 06.04.2020


Ответы (1)


Это легко сделать с помощью Shapely. Вы можете установить его из Conda Forge: conda install shapely -c conda-forge

Код, который вам нужен, находится на github.gist, на основе ответа @Габриэль и @pv.:

# coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
from shapely.geometry import Polygon

def voronoi_finite_polygons_2d(vor, radius=None):
    """
    Reconstruct infinite voronoi regions in a 2D diagram to finite
    regions.
    Parameters
    ----------
    vor : Voronoi
        Input diagram
    radius : float, optional
        Distance to 'points at infinity'.
    Returns
    -------
    regions : list of tuples
        Indices of vertices in each revised Voronoi regions.
    vertices : list of tuples
        Coordinates for revised Voronoi vertices. Same as coordinates
        of input vertices, with 'points at infinity' appended to the
        end.
    """

    if vor.points.shape[1] != 2:
        raise ValueError("Requires 2D input")

    new_regions = []
    new_vertices = vor.vertices.tolist()

    center = vor.points.mean(axis=0)
    if radius is None:
        radius = vor.points.ptp().max()*2

    # Construct a map containing all ridges for a given point
    all_ridges = {}
    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))

    # Reconstruct infinite regions
    for p1, region in enumerate(vor.point_region):
        vertices = vor.regions[region]

        if all(v >= 0 for v in vertices):
            # finite region
            new_regions.append(vertices)
            continue

        # reconstruct a non-finite region
        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]

        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1
            if v1 >= 0:
                # finite ridge: already in the region
                continue

            # Compute the missing endpoint of an infinite ridge

            t = vor.points[p2] - vor.points[p1] # tangent
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # normal

            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[v2] + direction * radius

            new_region.append(len(new_vertices))
            new_vertices.append(far_point.tolist())

        # sort region counterclockwise
        vs = np.asarray([new_vertices[v] for v in new_region])
        c = vs.mean(axis=0)
        angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
        new_region = np.array(new_region)[np.argsort(angles)]

        # finish
        new_regions.append(new_region.tolist())

    return new_regions, np.asarray(new_vertices)

# make up data points
np.random.seed(1234)
points = np.random.rand(15, 2)

# compute Voronoi tesselation
vor = Voronoi(points)

# plot
regions, vertices = voronoi_finite_polygons_2d(vor)

min_x = vor.min_bound[0] - 0.1
max_x = vor.max_bound[0] + 0.1
min_y = vor.min_bound[1] - 0.1
max_y = vor.max_bound[1] + 0.1

mins = np.tile((min_x, min_y), (vertices.shape[0], 1))
bounded_vertices = np.max((vertices, mins), axis=0)
maxs = np.tile((max_x, max_y), (vertices.shape[0], 1))
bounded_vertices = np.min((bounded_vertices, maxs), axis=0)



box = Polygon([[min_x, min_y], [min_x, max_y], [max_x, max_y], [max_x, min_y]])

# colorize
for region in regions:
    polygon = vertices[region]
    # Clipping polygon
    poly = Polygon(polygon)
    poly = poly.intersection(box)
    polygon = [p for p in poly.exterior.coords]

    plt.fill(*zip(*polygon), alpha=0.4)

plt.plot(points[:, 0], points[:, 1], 'ko')
plt.axis('equal')
plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1)
plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1)

plt.savefig('voro.png')
plt.show()
person Sklavit    schedule 26.03.2017