如何使用 matplotlib 在 Python 中创建 Voronoi 镶嵌多边形的轮廓?
How do you create outlines to the polygons of a Voronoi tessellation in Python using matplotlib?
我正在尝试从 trajectory.xyz 文件创建一个 voronoi 镶嵌,其中包含 305 个粒子和多帧数据。我在这里使用了在 GitHub 上找到的这个脚本:
https://gist.github.com/pv/8036995
并重新设计了一些东西以更准确地满足我绘制单帧曲面细分的目标(下面的代码)。然而,我的问题是,每当我 运行 代码时,生成的任何多边形都没有 border/outline。
这让整个事情看起来有点混乱和混乱。我相当有信心它与曲面细分的图形有关,但我不确定如何修复它。即使我 运行 GitHub 脚本与随机点完全一样,我的多边形也没有轮廓。有什么想法吗?
代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
from curved_analysis import read_xyz, read_nfo, project_frame
def voronoi(traj):
voronoi_frames = []
coords = np.array(read_xyz("traj0.xyz"))
for k in range(coords.shape[0]):
points = project_frame(coords[k])
vor = Voronoi(points[:, :2], qhull_options=('Qz'))
voronoi_frames.append(vor)
return points, voronoi_frames
def voronoi_finite_polygons_2d(vor, radius=None):
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
def voronoi_plot(vor, points):
regions, vertices = voronoi_finite_polygons_2d(vor)
for region in regions:
polygon = vertices[region]
plt.fill(*zip(*polygon), alpha=0.4)
plt.plot(points[:,0], points[:,1], 'ko', marker=".", markersize=2)
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()
if __name__ == "__main__":
frame_num = 166
points, vor = voronoi("traj0.xyz")
voronoi_plot(vor[frame_num], points)
替换行
plt.fill(*zip(*polygon), alpha=0.4)
与:
plt.fill(*zip(*polygon), alpha=0.4, lw=1, ec='k')
lw
设置多边形边的宽度,ec
颜色。
我正在尝试从 trajectory.xyz 文件创建一个 voronoi 镶嵌,其中包含 305 个粒子和多帧数据。我在这里使用了在 GitHub 上找到的这个脚本: https://gist.github.com/pv/8036995 并重新设计了一些东西以更准确地满足我绘制单帧曲面细分的目标(下面的代码)。然而,我的问题是,每当我 运行 代码时,生成的任何多边形都没有 border/outline。
这让整个事情看起来有点混乱和混乱。我相当有信心它与曲面细分的图形有关,但我不确定如何修复它。即使我 运行 GitHub 脚本与随机点完全一样,我的多边形也没有轮廓。有什么想法吗?
代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
from curved_analysis import read_xyz, read_nfo, project_frame
def voronoi(traj):
voronoi_frames = []
coords = np.array(read_xyz("traj0.xyz"))
for k in range(coords.shape[0]):
points = project_frame(coords[k])
vor = Voronoi(points[:, :2], qhull_options=('Qz'))
voronoi_frames.append(vor)
return points, voronoi_frames
def voronoi_finite_polygons_2d(vor, radius=None):
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
def voronoi_plot(vor, points):
regions, vertices = voronoi_finite_polygons_2d(vor)
for region in regions:
polygon = vertices[region]
plt.fill(*zip(*polygon), alpha=0.4)
plt.plot(points[:,0], points[:,1], 'ko', marker=".", markersize=2)
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()
if __name__ == "__main__":
frame_num = 166
points, vor = voronoi("traj0.xyz")
voronoi_plot(vor[frame_num], points)
替换行
plt.fill(*zip(*polygon), alpha=0.4)
与:
plt.fill(*zip(*polygon), alpha=0.4, lw=1, ec='k')
lw
设置多边形边的宽度,ec
颜色。