Конечные граничные клетки Питона Вороного

Я пытаюсь адаптировать код, который я нашел в stackoverflow, чтобы создать ячейку voronoi с конечными границами. Я нашел код ниже наhttps://stackoverflow.com/a/20678647/2443944 однако моя проблема в том, что, хотя клетки вороной не уходят в бесконечность на границах, они все еще слишком далеко. Даже с радиусом = 0 вершины гребня находятся слишком далеко. В идеале я хочу, чтобы граничные вершины вороной располагались примерно на том же уровне, что и остальные ячейки вороной в центре, т.е. я хочу, чтобы размеры ячеек вороной на границах были такими же, как и в центре.

Точки данных, которые я использую:

points = [[-30.0, 30.370371], [-27.777777, 35.925926], [-34.444443, 58.51852], [-2.9629631, 57.777779], [-17.777779, 75.185181], [-29.25926, 58.148151], [-11.111112, 33.703705], [-11.481482, 40.0], [-27.037037, 40.0], [-7.7777777, 94.444443], [-2.2222223, 122.22222], [-20.370371, 106.66667], [1.1111112, 125.18518], [-6.2962961, 128.88889], [6.666667, 133.7037], [11.851852, 136.2963], [8.5185184, 140.74074], [20.370371, 92.962959], [17.777779, 114.81482], [12.962962, 97.037041], [13.333334, 127.77778], [22.592592, 120.37037], [16.296295, 127.77778], [11.851852, 50.740742], [20.370371, 54.814816], [19.25926, 47.40741], [32.59259, 122.96296], [20.74074, 130.0], [24.814816, 84.814819], [26.296295, 91.111107], [56.296295, 131.48149], [60.0, 141.85185], [32.222221, 136.66667], [53.703705, 147.03703], [87.40741, 196.2963], [34.074074, 159.62964], [34.444443, -2.5925925], [36.666668, -1.8518518], [34.074074, -7.4074073], [35.555557, -18.888889], [76.666664, -39.629627], [35.185184, -37.777779], [25.185184, 14.074074], [42.962959, 32.962963], [35.925926, 9.2592592], [52.222221, 77.777779], [57.777779, 92.222221], [47.037041, 92.59259], [82.222221, 54.074074], [48.888889, 24.444445], [35.925926, 47.777779], [50.740742, 69.259254], [51.111111, 51.851849], [56.666664, -12.222222], [117.40741, -4.4444447], [59.629631, -5.9259262], [66.666664, 134.07408], [91.481483, 127.40741], [66.666664, 141.48149], [53.703705, 4.0740738], [85.185181, 11.851852], [69.629631, 0.37037039], [68.518517, 99.259262], [75.185181, 100.0], [70.370369, 113.7037], [74.444443, 82.59259], [82.222221, 93.703697], [72.222221, 84.444443], [77.777779, 167.03703], [88.888893, 168.88889], [73.703705, 178.88889], [87.037041, 123.7037], [78.518517, 97.037041], [95.555557, 52.962959], [85.555557, 57.037041], [90.370369, 23.333332], [100.0, 28.51852], [88.888893, 37.037037], [87.037041, -42.962959], [89.259262, -24.814816], [93.333328, 7.4074073], [98.518517, 5.185185], [92.59259, 1.4814816], [85.925919, 153.7037], [95.555557, 154.44444], [92.962959, 150.0], [97.037041, 95.925919], [106.66667, 115.55556], [92.962959, 114.81482], [108.88889, 56.296295], [97.777779, 50.740742], [94.074081, 89.259262], [96.666672, 91.851852], [102.22222, 77.777779], [107.40741, 40.370369], [105.92592, 29.629629], [105.55556, -46.296295], [118.51852, -47.777779], [112.22222, -43.333336], [112.59259, 25.185184], [115.92592, 27.777777], [112.59259, 31.851852], [107.03704, -36.666668], [118.88889, -32.59259], [114.07408, -25.555555], [115.92592, 85.185181], [105.92592, 18.888889], [121.11111, 14.444445], [129.25926, -28.51852], [127.03704, -18.518518], [139.25926, -12.222222], [141.48149, 3.7037036], [137.03703, -4.814815], [153.7037, -26.666668], [-2.2222223, 5.5555558], [0.0, 9.6296301], [10.74074, 20.74074], [2.2222223, 54.074074], [4.0740738, 50.740742], [34.444443, 46.296295], [11.481482, 1.4814816], [24.074076, -2.9629631], [74.814819, 79.259254], [67.777779, 152.22223], [57.037041, 127.03704], [89.259262, 12.222222]]
points = np.array(points)

Картины, которые я возвращаю, ниже для радиуса радиуса = 0.

 piccolo23 янв. 2016 г., 22:18
@ MGC Да это было бы хорошо
 mgc23 янв. 2016 г., 22:17
Было бы хорошо обрезать этот результат выпуклой оболочкой вашего набора точек? (или слегка забуференный выпуклый корпус)

Ответы на вопрос(2)

и снова используяvoronoi_finite_polygons_2d отhttps://stackoverflow.com/a/43023639/855617Вот решение для обрезания тесселяции Вороного до произвольной формы (здесь из бинарной маски). Единственная дополнительная работа здесь - создание полигона из вашей маски. Я уверен, что есть другие (и, вероятно, лучшие) способы полигонизации маски, подобные этой, но это сработало для моих целей.

import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage.morphology import binary_erosion
from scipy.spatial import Voronoi
from shapely.geometry import Point, Polygon
from skimage import draw
from sklearn.neighbors import KDTree

def get_circular_se(radius=2):

    N = (radius * 2) + 1
    se = np.zeros(shape=[N,N])
    for i in range(N):
        for j in range(N):
                se[i,j] = (i - N / 2)**2 + (j - N / 2)**2 <= radius**2
    se = np.array(se, dtype="uint8")
    return se

def polygonize_by_nearest_neighbor(pp):
    """Takes a set of xy coordinates pp Numpy array(n,2) and reorders the array to make
    a polygon using a nearest neighbor approach.

    """

    # start with first index
    pp_new = np.zeros_like(pp)
    pp_new[0] = pp[0]
    p_current_idx = 0

    tree = KDTree(pp)

    for i in range(len(pp) - 1):

        nearest_dist, nearest_idx = tree.query([pp[p_current_idx]], k=4)  # k1 = identity
        nearest_idx = nearest_idx[0]

        # finds next nearest point along the contour and adds it
        for min_idx in nearest_idx[1:]:  # skip the first point (will be zero for same pixel)
            if not pp[min_idx].tolist() in pp_new.tolist():  # make sure it's not already in the list
                pp_new[i + 1] = pp[min_idx]
                p_current_idx = min_idx
                break

    pp_new[-1] = pp[0]
    return pp_new


#generates a circular mask
side_len = 512
rad = 100
mask = np.zeros(shape=(side_len, side_len))
rr, cc = draw.circle(side_len/2, side_len/2, radius=rad, shape=mask.shape)
mask[rr, cc] = 1

#makes a polygon from the mask perimeter
se = get_circular_se(radius=1)
contour = mask - binary_erosion(mask, structure=se)
pixels_mask = np.array(np.where(contour==1)[::-1]).T
polygon = polygonize_by_nearest_neighbor(pixels_mask)
polygon = Polygon(polygon)

#generates random seeds
points_x = np.random.random_integers(0,side_len,250)
points_y = np.random.random_integers(0,side_len,250)
points = (np.vstack((points_x,points_y))).T

# returns a list of the centroids that are contained within the polygon
new_points = []
for point in points:
    if polygon.contains(Point(point)):
        new_points.append(point)

#performs voronoi tesselation
if len(points) > 3: #otherwise the tesselation won't work
    vor = Voronoi(new_points)
    regions, vertices = voronoi_finite_polygons_2d(vor)

    #clips tesselation to the mask
    new_vertices = []
    for region in regions:
        poly_reg = vertices[region]
        shape = list(poly_reg.shape)
        shape[0] += 1
        p = Polygon(np.append(poly_reg, poly_reg[0]).reshape(*shape)).intersection(polygon)
        poly = (np.array(p.exterior.coords)).tolist()
        new_vertices.append(poly)

    #plots the results
    fig, ax = plt.subplots()
    ax.imshow(mask,cmap='Greys_r')
    for poly in new_vertices:
        ax.fill(*zip(*poly), alpha=0.7)
    ax.plot(points[:,0],points[:,1],'ro',ms=2)
    plt.show()

Решение Вопроса

вы могли бы добиться этого, обрезав свой результат выпуклой оболочкой ваших точек. Для этого я бы, вероятно, использовалshapely модуль. УчитываяТАК сообщение Вы связали, я предполагаю, что вы используетеvoronoi_finite_polygons_2d Функция написана в посте. Так что я думаю, что это может сделать работу:

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import MultiPoint, Point, Polygon
from scipy.spatial import Voronoi

points = [[-30.0, 30.370371], [-27.777777, 35.925926], [-34.444443, 58.51852], [-2.9629631, 57.777779], [-17.777779, 75.185181], [-29.25926, 58.148151], [-11.111112, 33.703705], [-11.481482, 40.0], [-27.037037, 40.0], [-7.7777777, 94.444443], [-2.2222223, 122.22222], [-20.370371, 106.66667], [1.1111112, 125.18518], [-6.2962961, 128.88889], [6.666667, 133.7037], [11.851852, 136.2963], [8.5185184, 140.74074], [20.370371, 92.962959], [17.777779, 114.81482], [12.962962, 97.037041], [13.333334, 127.77778], [22.592592, 120.37037], [16.296295, 127.77778], [11.851852, 50.740742], [20.370371, 54.814816], [19.25926, 47.40741], [32.59259, 122.96296], [20.74074, 130.0], [24.814816, 84.814819], [26.296295, 91.111107], [56.296295, 131.48149], [60.0, 141.85185], [32.222221, 136.66667], [53.703705, 147.03703], [87.40741, 196.2963], [34.074074, 159.62964], [34.444443, -2.5925925], [36.666668, -1.8518518], [34.074074, -7.4074073], [35.555557, -18.888889], [76.666664, -39.629627], [35.185184, -37.777779], [25.185184, 14.074074], [42.962959, 32.962963], [35.925926, 9.2592592], [52.222221, 77.777779], [57.777779, 92.222221], [47.037041, 92.59259], [82.222221, 54.074074], [48.888889, 24.444445], [35.925926, 47.777779], [50.740742, 69.259254], [51.111111, 51.851849], [56.666664, -12.222222], [117.40741, -4.4444447], [59.629631, -5.9259262], [66.666664, 134.07408], [91.481483, 127.40741], [66.666664, 141.48149], [53.703705, 4.0740738], [85.185181, 11.851852], [69.629631, 0.37037039], [68.518517, 99.259262], [75.185181, 100.0], [70.370369, 113.7037], [74.444443, 82.59259], [82.222221, 93.703697], [72.222221, 84.444443], [77.777779, 167.03703], [88.888893, 168.88889], [73.703705, 178.88889], [87.037041, 123.7037], [78.518517, 97.037041], [95.555557, 52.962959], [85.555557, 57.037041], [90.370369, 23.333332], [100.0, 28.51852], [88.888893, 37.037037], [87.037041, -42.962959], [89.259262, -24.814816], [93.333328, 7.4074073], [98.518517, 5.185185], [92.59259, 1.4814816], [85.925919, 153.7037], [95.555557, 154.44444], [92.962959, 150.0], [97.037041, 95.925919], [106.66667, 115.55556], [92.962959, 114.81482], [108.88889, 56.296295], [97.777779, 50.740742], [94.074081, 89.259262], [96.666672, 91.851852], [102.22222, 77.777779], [107.40741, 40.370369], [105.92592, 29.629629], [105.55556, -46.296295], [118.51852, -47.777779], [112.22222, -43.333336], [112.59259, 25.185184], [115.92592, 27.777777], [112.59259, 31.851852], [107.03704, -36.666668], [118.88889, -32.59259], [114.07408, -25.555555], [115.92592, 85.185181], [105.92592, 18.888889], [121.11111, 14.444445], [129.25926, -28.51852], [127.03704, -18.518518], [139.25926, -12.222222], [141.48149, 3.7037036], [137.03703, -4.814815], [153.7037, -26.666668], [-2.2222223, 5.5555558], [0.0, 9.6296301], [10.74074, 20.74074], [2.2222223, 54.074074], [4.0740738, 50.740742], [34.444443, 46.296295], [11.481482, 1.4814816], [24.074076, -2.9629631], [74.814819, 79.259254], [67.777779, 152.22223], [57.037041, 127.03704], [89.259262, 12.222222]]

points = np.array(points)

vor = Voronoi(points)

regions, vertices = voronoi_finite_polygons_2d(vor)

pts = MultiPoint([Point(i) for i in points])
mask = pts.convex_hull
new_vertices = []
for region in regions:
    polygon = vertices[region]
    shape = list(polygon.shape)
    shape[0] += 1
    p = Polygon(np.append(polygon, polygon[0]).reshape(*shape)).intersection(mask)
    poly = np.array(list(zip(p.boundary.coords.xy[0][:-1], p.boundary.coords.xy[1][:-1])))
    new_vertices.append(poly)
    plt.fill(*zip(*poly), alpha=0.4)
plt.plot(points[:,0], points[:,1], 'ko')
plt.title("Clipped Voronois")
plt.show()

Вообще говоря (т.е. без использованияvoronoi_finite_polygons_2d но используя непосредственно выводVoronoi если это соответствует моим потребностям), я бы сделал:

import numpy as np
import matplotlib.pyplot as plt
from shapely.ops import polygonize,unary_union
from shapely.geometry import LineString, MultiPolygon, MultiPoint, Point
from scipy.spatial import Voronoi
points = [[-30.0, 30.370371], [-27.777777, 35.925926], [-34.444443, 58.51852], [-2.9629631, 57.777779], [-17.777779, 75.185181], [-29.25926, 58.148151], [-11.111112, 33.703705], [-11.481482, 40.0], [-27.037037, 40.0], [-7.7777777, 94.444443], [-2.2222223, 122.22222], [-20.370371, 106.66667], [1.1111112, 125.18518], [-6.2962961, 128.88889], [6.666667, 133.7037], [11.851852, 136.2963], [8.5185184, 140.74074], [20.370371, 92.962959], [17.777779, 114.81482], [12.962962, 97.037041], [13.333334, 127.77778], [22.592592, 120.37037], [16.296295, 127.77778], [11.851852, 50.740742], [20.370371, 54.814816], [19.25926, 47.40741], [32.59259, 122.96296], [20.74074, 130.0], [24.814816, 84.814819], [26.296295, 91.111107], [56.296295, 131.48149], [60.0, 141.85185], [32.222221, 136.66667], [53.703705, 147.03703], [87.40741, 196.2963], [34.074074, 159.62964], [34.444443, -2.5925925], [36.666668, -1.8518518], [34.074074, -7.4074073], [35.555557, -18.888889], [76.666664, -39.629627], [35.185184, -37.777779], [25.185184, 14.074074], [42.962959, 32.962963], [35.925926, 9.2592592], [52.222221, 77.777779], [57.777779, 92.222221], [47.037041, 92.59259], [82.222221, 54.074074], [48.888889, 24.444445], [35.925926, 47.777779], [50.740742, 69.259254], [51.111111, 51.851849], [56.666664, -12.222222], [117.40741, -4.4444447], [59.629631, -5.9259262], [66.666664, 134.07408], [91.481483, 127.40741], [66.666664, 141.48149], [53.703705, 4.0740738], [85.185181, 11.851852], [69.629631, 0.37037039], [68.518517, 99.259262], [75.185181, 100.0], [70.370369, 113.7037], [74.444443, 82.59259], [82.222221, 93.703697], [72.222221, 84.444443], [77.777779, 167.03703], [88.888893, 168.88889], [73.703705, 178.88889], [87.037041, 123.7037], [78.518517, 97.037041], [95.555557, 52.962959], [85.555557, 57.037041], [90.370369, 23.333332], [100.0, 28.51852], [88.888893, 37.037037], [87.037041, -42.962959], [89.259262, -24.814816], [93.333328, 7.4074073], [98.518517, 5.185185], [92.59259, 1.4814816], [85.925919, 153.7037], [95.555557, 154.44444], [92.962959, 150.0], [97.037041, 95.925919], [106.66667, 115.55556], [92.962959, 114.81482], [108.88889, 56.296295], [97.777779, 50.740742], [94.074081, 89.259262], [96.666672, 91.851852], [102.22222, 77.777779], [107.40741, 40.370369], [105.92592, 29.629629], [105.55556, -46.296295], [118.51852, -47.777779], [112.22222, -43.333336], [112.59259, 25.185184], [115.92592, 27.777777], [112.59259, 31.851852], [107.03704, -36.666668], [118.88889, -32.59259], [114.07408, -25.555555], [115.92592, 85.185181], [105.92592, 18.888889], [121.11111, 14.444445], [129.25926, -28.51852], [127.03704, -18.518518], [139.25926, -12.222222], [141.48149, 3.7037036], [137.03703, -4.814815], [153.7037, -26.666668], [-2.2222223, 5.5555558], [0.0, 9.6296301], [10.74074, 20.74074], [2.2222223, 54.074074], [4.0740738, 50.740742], [34.444443, 46.296295], [11.481482, 1.4814816], [24.074076, -2.9629631], [74.814819, 79.259254], [67.777779, 152.22223], [57.037041, 127.03704], [89.259262, 12.222222]]
points = np.array(points)
vor = Voronoi(points)
lines = [
    LineString(vor.vertices[line])
    for line in vor.ridge_vertices if -1 not in line
]

convex_hull = MultiPoint([Point(i) for i in points]).convex_hull.buffer(2)
result = MultiPolygon(
    [poly.intersection(convex_hull) for poly in polygonize(lines)])
result = MultiPolygon(
    [p for p in result]
    + [p for p in convex_hull.difference(unary_union(result))])

plt.plot(points[:,0], points[:,1], 'ko')
for r in result:
    plt.fill(*zip(*np.array(list(
        zip(r.boundary.coords.xy[0][:-1], r.boundary.coords.xy[1][:-1])))),
        alpha=0.4)
plt.show()

За исключением небольшого буфера на выпуклой оболочке, результат должен выглядеть одинаково:

Или, если вы хотите получить результат, который немного менее «сырой», вы можете попробовать поиграть с буферным методом (и егоresolution/join_style/cap_style свойства) ваших точек (и / или буфера выпуклой оболочки):

pts = MultiPoint([Point(i) for i in points])
mask = pts.convex_hull.union(pts.buffer(10, resolution=5, cap_style=3))
result = MultiPolygon(
    [poly.intersection(mask) for poly in polygonize(lines)])

И получить что-то вроде (вы можете достичь лучшего ..!)

 DrBwts26 янв. 2016 г., 16:21
Ах, да, это было так. Все еще не получая вышеупомянутые графики все же.
 mgc25 янв. 2016 г., 19:10
@DrBwts это вашpoints переменная массив с такой же формой, как в примере? Я скопировал / вставил мой пример и получил тот же результат, что и тот, который виден в ответе.
 piccolo24 янв. 2016 г., 01:58
Привет @mgc, как бы вы получили формат polygon = vertice [region] для необработанного графа? большое спасибо
 mgc23 янв. 2016 г., 23:24
В моем случае я использовалgdf = geopandas.GeoDataFrame(geometry=[i for i in result]) ; gdf.plot() построить результат. Вы также можете добавить эти красивые многоугольники вmatplotlib.collections и построить сюжет с регулярным Matplotlib API.
 mgc22 янв. 2019 г., 12:54
@AaronBramson Я изменил второй фрагмент, чтобы включить в него многоугольники тезисов.
 piccolo23 янв. 2016 г., 23:32
В идеале я просто хочу обновить вершины в списке полигонов в исходном коде. Это будет возможно?
 piccolo24 янв. 2016 г., 01:28
Не могли бы вы переформатировать первый фрагмент кода для первого «сырого» графа, чтобы получить его также в том же формате polygon = vertice [region], пожалуйста?
 Aaron Bramson19 февр. 2019 г., 08:25
Хорошо, теперь они разделены там, где они раньше были неправильно объединены, но теперь они создают больше (а не меньше) многоугольников Вороного, чем количество точек (используя реальную геометрию ГИС в качестве ограничивающей области вместо выпуклой оболочки).
 mgc23 янв. 2019 г., 13:03
@Phonix мой плохой за упущения, которые вы исправили;)
 piccolo23 янв. 2016 г., 22:58
Это выглядит великолепно. Как бы я объединил это с оригинальным кодом, чтобы обновить список полигонов?
 mgc23 янв. 2016 г., 23:50
Хм, наверное, я постараюсь выяснить.
 Phonix23 янв. 2019 г., 08:38
@mgc Я имел в виду, что пишу полный код (включая часть кода вопроса, например voronoi_finite_polygons_2d) в каждом фрагменте, поэтому я просто копирую и запускаю sample.Right сейчас я должен найти каждую часть кода в разных местах. Спасибо
 piccolo24 янв. 2016 г., 00:58
Список полигонов не обновляется, чтобы включить новые вершины на границах. Есть ли способ, которым вы могли бы обновить polygon = vertices [region], чтобы включить обновления вершин?
 mgc24 янв. 2016 г., 00:43
@ user2443944 я отредактировал свое сообщение, но может быть действительно более элегантный способ сделать это ..!
 piccolo23 янв. 2016 г., 23:08
Я в основном использую этоstackoverflow.com/a/20678647/2443944
 piccolo24 янв. 2016 г., 01:20
Не обращайте внимания на мой предыдущий комментарий. Я получаю то же изображение. Спасибо!
 Aaron Bramson23 янв. 2019 г., 06:36
Спасибо! Теперь он включает в себя многоугольники, но он объединяет воедино отдельные многоугольники, которые разделяют бесконечную вершину ... или мне кажется, что в моем случае три разных, но близлежащих точки около края пространства отсечения помещаются в один и тот же многоугольник с использованием вашего метода , но они находятся в отдельных полигонах, используяvoronoi_finite_polygons_2d функция.
 Aaron Bramson22 янв. 2019 г., 11:27
Я изменил более низкие версии с маской, сделанной из области, определенной объединением кругов. Это работает, чтобы связать большинство полигонов Вороного, но некоторые точки в регионе вообще не имеют полигонов. Я думал, что это должно было обрезать многоугольники, которые уходят в бесконечность, а не исключать их. Можно ли это отрегулировать так, чтобы эти полигоны тоже были прикреплены к маске?
 mgc23 янв. 2019 г., 12:32
@Phonix thevoronoi_finite_polygons_2d довольно длинный, поэтому я не буду копировать его в своем посте. Но это нужно только для первого фрагмента. Второй фрагмент завершен и может быть запущен как есть. Третий фрагмент представляет собой замену, которая будет использоваться во втором фрагменте. Я не вижу необходимости повторять линии, связанные с построением или созданием точек.
 mgc24 янв. 2016 г., 01:05
Я отредактировал (и упростил) мой код, он создаетǹew_vertice список с новыми вершинами в том же формате, что и вершиныpolygon = vertice[region] объект (т. е. в виде пустого массива и без последней точки, закрывающей кольцо).
 Phonix22 янв. 2019 г., 12:14
пожалуйста, напишите полный код для каждого графика. Это меня очень смущает!
 mgc23 янв. 2016 г., 23:06
Я не знаю, как выглядит ваш оригинальный код?
 DrBwts25 янв. 2016 г., 17:34
Ммммм, я не получаю ничего похожего на диаграммы в ответе. Подумайте, может быть, мое масштабирование выключено? Такжеplt.plot(points[:,0], points[:,1], 'ko') выдает ошибку кортежа
 mgc25 янв. 2016 г., 19:03
@ user2443944 извините, но я не могу обновитьvertice массив (потому чтоregion это просто список индексов, ссылающихся на вершины региона), и из-за операции отсечения число вершин изменилось (вот почему я заполнил эти новые вершины вnew_vertices переменная, это список, где каждый элемент содержит вершины "области").
 mgc22 янв. 2019 г., 12:55
@Phonix Я изменил код, чтобы код для первого и второго фрагментов был кодом для полного графика.

Ваш ответ на вопрос