从 Voronoi 单元获取有界多边形坐标
Getting a bounded polygon coordinates from Voronoi cells
我有点(例如,基站位置的纬度、经度对),我需要获取它们形成的 Voronoi 单元的多边形。
from scipy.spatial import Voronoi
tower = [[ 24.686 , 46.7081],
[ 24.686 , 46.7081],
[ 24.686 , 46.7081]]
c = Voronoi(towers)
现在,我需要获取每个单元格的经纬度坐标中的多边形边界(以及此多边形围绕的质心是什么)。我也需要这个 Voronoi 有界。这意味着边界不会无限大,而是在边界框内。
给定一个矩形边界框,我的第一个想法是在这个边界框和 scipy.spatial.Voronoi
产生的 Voronoï 图之间定义一种交集操作。想法不一定好,因为这需要编写大量计算几何的基本函数。
但是,这是我想到的第二个想法(hack?):计算平面中一组 n
点的 Voronoï 图的算法的时间复杂度为 O(n ln(n))
.如何添加点以将初始点的 Voronoï 单元约束在边界框中?
有界 Voronoï 图的解法
一张图胜过一场精彩的演讲:
我在这里做了什么?这很简单!初始点(蓝色)位于 [0.0, 1.0] x [0.0, 1.0]
。然后我根据x = 0.0
(边界框的左边缘)通过反射对称性得到左侧(即[-1.0, 0.0] x [0.0, 1.0]
)的点(蓝色)。根据 x = 1.0
、y = 0.0
和 y = 1.0
(边界框的其他边缘)的反射对称性,我得到了完成这项工作所需的所有点(蓝色)。
那我运行scipy.spatial.Voronoi
。上图描绘了生成的 Voronoï 图(我使用 scipy.spatial.voronoi_plot_2d
)。
接下来要做什么?只需根据边界框过滤点、边或面。并根据well-known公式得到每个面的质心,计算centroid of polygon。这是结果的图像(质心为红色):
展示代码前的一些乐趣
太棒了!它似乎工作。如果在一次迭代后我尝试 re-run 质心(红色)而不是初始点(蓝色)上的算法怎么办?如果我一次又一次地尝试呢?
步骤 2
步骤 10
第 25 步
酷! Voronoï 细胞倾向于最小化它们的 能量 ...
这是代码
import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
eps = sys.float_info.epsilon
n_towers = 100
towers = np.random.rand(n_towers, 2)
bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max]
def in_box(towers, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
towers[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= towers[:, 1],
towers[:, 1] <= bounding_box[3]))
def voronoi(towers, bounding_box):
# Select towers inside the bounding box
i = in_box(towers, bounding_box)
# Mirror points
points_center = towers[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
# Filter regions
regions = []
for region in vor.regions:
flag = True
for index in region:
if index == -1:
flag = False
break
else:
x = vor.vertices[index, 0]
y = vor.vertices[index, 1]
if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
flag = False
break
if region != [] and flag:
regions.append(region)
vor.filtered_points = points_center
vor.filtered_regions = regions
return vor
def centroid_region(vertices):
# Polygon's signed area
A = 0
# Centroid's x
C_x = 0
# Centroid's y
C_y = 0
for i in range(0, len(vertices) - 1):
s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
A = A + s
C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
A = 0.5 * A
C_x = (1.0 / (6.0 * A)) * C_x
C_y = (1.0 / (6.0 * A)) * C_y
return np.array([[C_x, C_y]])
vor = voronoi(towers, bounding_box)
fig = pl.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
# Plot ridges points
for region in vor.filtered_regions:
vertices = vor.vertices[region, :]
ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
ax.plot(centroid[:, 0], centroid[:, 1], 'r.')
ax.set_xlim([-0.1, 1.1])
ax.set_ylim([-0.1, 1.1])
pl.savefig("bounded_voronoi.png")
sp.spatial.voronoi_plot_2d(vor)
pl.savefig("voronoi.png")
我在使用 scipy 的 voronoi 函数和创建 CVD 时遇到了很多麻烦,所以这些精彩的帖子和评论提供了极大的帮助。作为一名编程新手,我试图理解 Flabetvvibes 答案中的代码,我将分享我对它如何工作的解释以及 Energya 和我自己的修改。我还在这个答案的底部完整地发布了我的代码版本
import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
import copy
eps = sys.float_info.epsilon
# Returns a new np.array of towers that within the bounding_box
def in_box(towers, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
towers[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= towers[:, 1],
towers[:, 1] <= bounding_box[3]))
in_box 函数使用 numpy 的 logical_and 方法来 return 布尔数组,表示来自塔的哪些坐标在边界框中。
# Generates a bounded vornoi diagram with finite regions in the bounding box
def bounded_voronoi(towers, bounding_box):
# Select towers inside the bounding box
i = in_box(towers, bounding_box)
# Mirror points left, right, above, and under to provide finite regions for the
# edge regions of the bounding box
points_center = towers[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
Flabetvvibes 对点进行镜像,以允许沿着边界框内边缘的区域是有限的。 Scipy 的 voronoi 方法 returns -1 用于未定义的顶点,因此镜像点允许边界框内的所有区域都是有限的,而所有无限区域都在边界框外的镜像区域中稍后将被丢弃的边界框。
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
# creates a new attibute for points that form the diagram within the region
vor.filtered_points = points_center
# grabs the first fifth of the regions, which are the original regions
vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]]
return vor
bounded_voronoi 方法的最后一位调用 scipy 的 voronoi 函数并为边界框内的过滤点和区域添加新属性。 Energya 建议删除 Flabetvvibe 的代码,该代码手动找到边界框内的所有有限区域,其中一个衬里获得区域的前五分之一,这些区域是原始输入以及构成边界框的点。
def generate_CVD(points, iterations, bounding_box):
p = copy.copy(points)
for i in range(iterations):
vor = bounded_voronoi(p, bounding_box)
centroids = []
for region in vor.filtered_regions:
# grabs vertices for the region and adds a duplicate
# of the first one to the end
vertices = vor.vertices[region + [region[0]], :]
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
p = np.array(centroids)
return bounded_voronoi(p, bounding_box)
我采用了执行 loyd 算法迭代的 Flabetvvibe 代码,并将其形成为易于使用的方法。对于每次迭代,调用前一个 bounded_voronoi 函数,然后为每个单元找到质心,它们成为下一次迭代的新点集。 vertices = vor.vertices[region + [region[0]], :]
简单地抓取当前区域的所有顶点,并将第一个顶点复制到最后,这样第一个和最后一个顶点在计算质心时是相同的。
感谢 Flabetvvibes 和 Energya。您的 posts/answers 教会了我如何比其文档更好地使用 scipy 的 voronoi 方法。我还将代码作为一个整体发布到任何其他寻找 copy/paste.
的人下面
import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
import copy
eps = sys.float_info.epsilon
# Returns a new np.array of towers that within the bounding_box
def in_box(towers, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
towers[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= towers[:, 1],
towers[:, 1] <= bounding_box[3]))
# Generates a bounded vornoi diagram with finite regions
def bounded_voronoi(towers, bounding_box):
# Select towers inside the bounding box
i = in_box(towers, bounding_box)
# Mirror points left, right, above, and under to provide finite regions for the edge regions of the bounding box
points_center = towers[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
vor.filtered_points = points_center # creates a new attibute for points that form the diagram within the region
vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]] # grabs the first fifth of the regions, which are the original regions
return vor
# Finds the centroid of a region. First and last point should be the same.
def centroid_region(vertices):
# Polygon's signed area
A = 0
# Centroid's x
C_x = 0
# Centroid's y
C_y = 0
for i in range(0, len(vertices) - 1):
s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
A = A + s
C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
A = 0.5 * A
C_x = (1.0 / (6.0 * A)) * C_x
C_y = (1.0 / (6.0 * A)) * C_y
return np.array([[C_x, C_y]])
# Performs x iterations of loyd's algorithm to calculate a centroidal vornoi diagram
def generate_CVD(points, iterations, bounding_box):
p = copy.copy(points)
for i in range(iterations):
vor = bounded_voronoi(p, bounding_box)
centroids = []
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :] # grabs vertices for the region and adds a duplicate of the first one to the end
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
p = np.array(centroids)
return bounded_voronoi(p, bounding_box)
# returns a pyplot of given voronoi data
def plot_vornoi_diagram(vor, bounding_box, show_figure):
# Initializes pyplot stuff
fig = pl.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
# Plot ridges points
for region in vor.filtered_regions:
vertices = vor.vertices[region, :]
ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# stores references to numbers for setting axes limits
margin_percent = .1
width = bounding_box[1]-bounding_box[0]
height = bounding_box[3]-bounding_box[2]
ax.set_xlim([bounding_box[0]-width*margin_percent, bounding_box[1]+width*margin_percent])
ax.set_ylim([bounding_box[2]-height*margin_percent, bounding_box[3]+height*margin_percent])
if show_figure:
pl.show()
return fig
我有点(例如,基站位置的纬度、经度对),我需要获取它们形成的 Voronoi 单元的多边形。
from scipy.spatial import Voronoi
tower = [[ 24.686 , 46.7081],
[ 24.686 , 46.7081],
[ 24.686 , 46.7081]]
c = Voronoi(towers)
现在,我需要获取每个单元格的经纬度坐标中的多边形边界(以及此多边形围绕的质心是什么)。我也需要这个 Voronoi 有界。这意味着边界不会无限大,而是在边界框内。
给定一个矩形边界框,我的第一个想法是在这个边界框和 scipy.spatial.Voronoi
产生的 Voronoï 图之间定义一种交集操作。想法不一定好,因为这需要编写大量计算几何的基本函数。
但是,这是我想到的第二个想法(hack?):计算平面中一组 n
点的 Voronoï 图的算法的时间复杂度为 O(n ln(n))
.如何添加点以将初始点的 Voronoï 单元约束在边界框中?
有界 Voronoï 图的解法
一张图胜过一场精彩的演讲:
我在这里做了什么?这很简单!初始点(蓝色)位于 [0.0, 1.0] x [0.0, 1.0]
。然后我根据x = 0.0
(边界框的左边缘)通过反射对称性得到左侧(即[-1.0, 0.0] x [0.0, 1.0]
)的点(蓝色)。根据 x = 1.0
、y = 0.0
和 y = 1.0
(边界框的其他边缘)的反射对称性,我得到了完成这项工作所需的所有点(蓝色)。
那我运行scipy.spatial.Voronoi
。上图描绘了生成的 Voronoï 图(我使用 scipy.spatial.voronoi_plot_2d
)。
接下来要做什么?只需根据边界框过滤点、边或面。并根据well-known公式得到每个面的质心,计算centroid of polygon。这是结果的图像(质心为红色):
展示代码前的一些乐趣
太棒了!它似乎工作。如果在一次迭代后我尝试 re-run 质心(红色)而不是初始点(蓝色)上的算法怎么办?如果我一次又一次地尝试呢?
步骤 2
步骤 10
第 25 步
酷! Voronoï 细胞倾向于最小化它们的 能量 ...
这是代码
import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
eps = sys.float_info.epsilon
n_towers = 100
towers = np.random.rand(n_towers, 2)
bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max]
def in_box(towers, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
towers[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= towers[:, 1],
towers[:, 1] <= bounding_box[3]))
def voronoi(towers, bounding_box):
# Select towers inside the bounding box
i = in_box(towers, bounding_box)
# Mirror points
points_center = towers[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
# Filter regions
regions = []
for region in vor.regions:
flag = True
for index in region:
if index == -1:
flag = False
break
else:
x = vor.vertices[index, 0]
y = vor.vertices[index, 1]
if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
flag = False
break
if region != [] and flag:
regions.append(region)
vor.filtered_points = points_center
vor.filtered_regions = regions
return vor
def centroid_region(vertices):
# Polygon's signed area
A = 0
# Centroid's x
C_x = 0
# Centroid's y
C_y = 0
for i in range(0, len(vertices) - 1):
s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
A = A + s
C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
A = 0.5 * A
C_x = (1.0 / (6.0 * A)) * C_x
C_y = (1.0 / (6.0 * A)) * C_y
return np.array([[C_x, C_y]])
vor = voronoi(towers, bounding_box)
fig = pl.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
# Plot ridges points
for region in vor.filtered_regions:
vertices = vor.vertices[region, :]
ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
ax.plot(centroid[:, 0], centroid[:, 1], 'r.')
ax.set_xlim([-0.1, 1.1])
ax.set_ylim([-0.1, 1.1])
pl.savefig("bounded_voronoi.png")
sp.spatial.voronoi_plot_2d(vor)
pl.savefig("voronoi.png")
我在使用 scipy 的 voronoi 函数和创建 CVD 时遇到了很多麻烦,所以这些精彩的帖子和评论提供了极大的帮助。作为一名编程新手,我试图理解 Flabetvvibes 答案中的代码,我将分享我对它如何工作的解释以及 Energya 和我自己的修改。我还在这个答案的底部完整地发布了我的代码版本
import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
import copy
eps = sys.float_info.epsilon
# Returns a new np.array of towers that within the bounding_box
def in_box(towers, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
towers[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= towers[:, 1],
towers[:, 1] <= bounding_box[3]))
in_box 函数使用 numpy 的 logical_and 方法来 return 布尔数组,表示来自塔的哪些坐标在边界框中。
# Generates a bounded vornoi diagram with finite regions in the bounding box
def bounded_voronoi(towers, bounding_box):
# Select towers inside the bounding box
i = in_box(towers, bounding_box)
# Mirror points left, right, above, and under to provide finite regions for the
# edge regions of the bounding box
points_center = towers[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
Flabetvvibes 对点进行镜像,以允许沿着边界框内边缘的区域是有限的。 Scipy 的 voronoi 方法 returns -1 用于未定义的顶点,因此镜像点允许边界框内的所有区域都是有限的,而所有无限区域都在边界框外的镜像区域中稍后将被丢弃的边界框。
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
# creates a new attibute for points that form the diagram within the region
vor.filtered_points = points_center
# grabs the first fifth of the regions, which are the original regions
vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]]
return vor
bounded_voronoi 方法的最后一位调用 scipy 的 voronoi 函数并为边界框内的过滤点和区域添加新属性。 Energya 建议删除 Flabetvvibe 的代码,该代码手动找到边界框内的所有有限区域,其中一个衬里获得区域的前五分之一,这些区域是原始输入以及构成边界框的点。
def generate_CVD(points, iterations, bounding_box):
p = copy.copy(points)
for i in range(iterations):
vor = bounded_voronoi(p, bounding_box)
centroids = []
for region in vor.filtered_regions:
# grabs vertices for the region and adds a duplicate
# of the first one to the end
vertices = vor.vertices[region + [region[0]], :]
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
p = np.array(centroids)
return bounded_voronoi(p, bounding_box)
我采用了执行 loyd 算法迭代的 Flabetvvibe 代码,并将其形成为易于使用的方法。对于每次迭代,调用前一个 bounded_voronoi 函数,然后为每个单元找到质心,它们成为下一次迭代的新点集。 vertices = vor.vertices[region + [region[0]], :]
简单地抓取当前区域的所有顶点,并将第一个顶点复制到最后,这样第一个和最后一个顶点在计算质心时是相同的。
感谢 Flabetvvibes 和 Energya。您的 posts/answers 教会了我如何比其文档更好地使用 scipy 的 voronoi 方法。我还将代码作为一个整体发布到任何其他寻找 copy/paste.
的人下面import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
import copy
eps = sys.float_info.epsilon
# Returns a new np.array of towers that within the bounding_box
def in_box(towers, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
towers[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= towers[:, 1],
towers[:, 1] <= bounding_box[3]))
# Generates a bounded vornoi diagram with finite regions
def bounded_voronoi(towers, bounding_box):
# Select towers inside the bounding box
i = in_box(towers, bounding_box)
# Mirror points left, right, above, and under to provide finite regions for the edge regions of the bounding box
points_center = towers[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
vor.filtered_points = points_center # creates a new attibute for points that form the diagram within the region
vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]] # grabs the first fifth of the regions, which are the original regions
return vor
# Finds the centroid of a region. First and last point should be the same.
def centroid_region(vertices):
# Polygon's signed area
A = 0
# Centroid's x
C_x = 0
# Centroid's y
C_y = 0
for i in range(0, len(vertices) - 1):
s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
A = A + s
C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
A = 0.5 * A
C_x = (1.0 / (6.0 * A)) * C_x
C_y = (1.0 / (6.0 * A)) * C_y
return np.array([[C_x, C_y]])
# Performs x iterations of loyd's algorithm to calculate a centroidal vornoi diagram
def generate_CVD(points, iterations, bounding_box):
p = copy.copy(points)
for i in range(iterations):
vor = bounded_voronoi(p, bounding_box)
centroids = []
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :] # grabs vertices for the region and adds a duplicate of the first one to the end
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
p = np.array(centroids)
return bounded_voronoi(p, bounding_box)
# returns a pyplot of given voronoi data
def plot_vornoi_diagram(vor, bounding_box, show_figure):
# Initializes pyplot stuff
fig = pl.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
# Plot ridges points
for region in vor.filtered_regions:
vertices = vor.vertices[region, :]
ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# stores references to numbers for setting axes limits
margin_percent = .1
width = bounding_box[1]-bounding_box[0]
height = bounding_box[3]-bounding_box[2]
ax.set_xlim([bounding_box[0]-width*margin_percent, bounding_box[1]+width*margin_percent])
ax.set_ylim([bounding_box[2]-height*margin_percent, bounding_box[3]+height*margin_percent])
if show_figure:
pl.show()
return fig