Colorize Voronoi Diagram

Question

Я пытаюсь раскрасить диаграмму Вороного, созданную с помощью scipy.spatial.Voronoi. Вот мой код:

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.spatial import Voronoi, voronoi_plot_2d 

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

# compute Voronoi tesselation 
vor = Voronoi(points) 

# plot 
voronoi_plot_2d(vor) 

# colorize 
for region in vor.regions: 
    if not -1 in region: 
     polygon = [vor.vertices[i] for i in region] 
     plt.fill(*zip(*polygon)) 

plt.show() 

Полученное изображение:

Как вы можете увидеть некоторые из областей Вороного на границе изображения не окрашиваются. Это связано с тем, что некоторые индексы к вершинам Вороного для этих областей установлены на -1, то есть для тех вершин за пределами диаграммы Вороного. Согласно документации:

регионов: (список списка Интс, формы (nregions, *)) Индексы вершин Вороным, образующим каждой Вороной области. -1 указывает вершину вне диаграммы Вороного.

Для того, чтобы раскрасить эти регионы, а также, я пытался просто удалить эту «внешнюю» вершину из многоугольника, но это не сработало. Я думаю, мне нужно заполнить некоторые пункты на границе области изображения, но я не могу понять, как достичь этого разумно.

Может ли кто-нибудь помочь?

Answer 1

Структура данных Voronoi содержит всю необходимую информацию для построения позиций для «точек на бесконечности». Qhull также сообщает их просто как -1 индексы, поэтому Scipy не вычисляет их для вас.

https://gist.github.com/pv/8036995

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.spatial import Voronoi 

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() 

    # 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) 
print "--" 
print regions 
print "--" 
print vertices 

# colorize 
for region in regions: 
    polygon = vertices[region] 
    plt.fill(*zip(*polygon), alpha=0.4) 

plt.plot(points[:,0], points[:,1], 'ko') 
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.show() 

Answer 2

Я не думаю, что из данных, доступных в структуре vor, достаточно информации, чтобы понять это, не выполняя, по крайней мере, некоторые вычисления voronoi. Так как это так, вот соответствующие части исходной функции voronoi_plot_2d, которую вы должны использовать, чтобы извлекать точки, которые пересекаются с vor.max_bound или vor.min_bound, которые являются нижним левым и верхним правыми углами диаграммы в укажите другие координаты для ваших полигонов.

for simplex in vor.ridge_vertices: 
    simplex = np.asarray(simplex) 
    if np.all(simplex >= 0): 
     ax.plot(vor.vertices[simplex,0], vor.vertices[simplex,1], 'k-') 

ptp_bound = vor.points.ptp(axis=0) 
center = vor.points.mean(axis=0) 
for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices): 
    simplex = np.asarray(simplex) 
    if np.any(simplex < 0): 
     i = simplex[simplex >= 0][0] # finite end Voronoi vertex 

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

     midpoint = vor.points[pointidx].mean(axis=0) 
     direction = np.sign(np.dot(midpoint - center, n)) * n 
     far_point = vor.vertices[i] + direction * ptp_bound.max() 

     ax.plot([vor.vertices[i,0], far_point[0]], 
       [vor.vertices[i,1], far_point[1]], 'k--')