如何为 Voronoi 镶嵌设置固定的外边界?
How to set a fixed outer boundary to Voronoi tessellations?
我画了一个 Voronoi 镶嵌(采矿业的爆炸图案)。我必须绘制 Voronoi 镶嵌的外边界,但我不想要盒子的边界;我想设置固定的外部单元格边界。
- 我得到这个结果:
- 我想要的结果是:
代码:
import xlrd
import operator
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d
wb = xlrd.open_workbook('C:/Users/s.gaur/desktop/Collar Coordinates 2620 S3C 5007 P2.xls')
sh1 = wb.sheet_by_name(u'2620-s3c-5007')
x = sh1.col_values(0)
y = sh1.col_values(1)
L = sorted(zip(x,y), key = operator.itemgetter(0))
Point = (L)
vor = Voronoi(Point)
voronoi_plot_2d(vor)
plt.show()
如何将外边缘边界固定到 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()*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)
# compute Voronoi tesselation
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),"brown", alpha = 0.4, edgecolor = 'black')
plt.plot(x, y, 'ko')
plt.plot(Dx,Dy, 'ko',markerfacecolor = 'red', markersize = 10)
plt.title("Blast 2620 S3C 5009 P1")
plt.show()
我画了一个 Voronoi 镶嵌(采矿业的爆炸图案)。我必须绘制 Voronoi 镶嵌的外边界,但我不想要盒子的边界;我想设置固定的外部单元格边界。
- 我得到这个结果:
- 我想要的结果是:
代码:
import xlrd
import operator
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d
wb = xlrd.open_workbook('C:/Users/s.gaur/desktop/Collar Coordinates 2620 S3C 5007 P2.xls')
sh1 = wb.sheet_by_name(u'2620-s3c-5007')
x = sh1.col_values(0)
y = sh1.col_values(1)
L = sorted(zip(x,y), key = operator.itemgetter(0))
Point = (L)
vor = Voronoi(Point)
voronoi_plot_2d(vor)
plt.show()
如何将外边缘边界固定到 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()*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)
# compute Voronoi tesselation
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),"brown", alpha = 0.4, edgecolor = 'black')
plt.plot(x, y, 'ko')
plt.plot(Dx,Dy, 'ko',markerfacecolor = 'red', markersize = 10)
plt.title("Blast 2620 S3C 5009 P1")
plt.show()