def __init__(self, points = np.array([]), dimensions = (None, None), granularity = None):
"""
Creates the Voronoi object
:param points: predefined points
:type points: numpy array of shape (w, 2) where w is the number of points [x, y] style, default None
:param dimensions: dimensions of the points, from [w, 2] where w is the highest value, this *cannot* be None if points is None
:type dimensions: tuple of ints, maximum (x,y) dimensions, default None
:param granularity: how many points to create, must be given if dimensions are given
:type granularity: int, default None
"""
if len(points) == 0 and dimensions == (None, None):
print('You can\'t have both points and dimensions be empty, try passing in some points or dimensions and granularity.')
return
if len(points) == 0 and dimensions != None and granularity == None:
print('Granularity can\'t be none if dimensions are passed in, try passing in a granularity.')
return
if len(points) != 0:
self.points = points
else:
points = np.random.random((granularity, 2))
points = list(map(lambda x: np.array([x[0]*dimensions[0], x[1]*dimensions[1]]), points))
self.points = np.array(points)
self.bounding_region = [min(self.points[:, 0]), max(self.points[:, 0]), min(self.points[:, 1]), max(self.points[:, 1])]
python类Voronoi()的实例源码
def _region_centroid(self, vertices):
"""
Finds the centroid of the voronoi region bounded by given vertices
See: https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
:param vertices: list of vertices that bound the region
:type vertices: numpy array of vertices from the scipy.spatial.Voronoi.regions (e.g. vor.vertices[region + [region[0]], :])
:return: list of centroids
:rtype: np.array of centroids
"""
signed_area = 0
C_x = 0
C_y = 0
for i in range(len(vertices)-1):
step = (vertices[i, 0]*vertices[i+1, 1])-(vertices[i+1, 0]*vertices[i, 1])
signed_area += step
C_x += (vertices[i, 0] + vertices[i+1, 0])*step
C_y += (vertices[i, 1] + vertices[i+1, 1])*step
signed_area = 1/2*signed_area
C_x = (1.0/(6.0*signed_area))*C_x
C_y = (1.0/(6.0*signed_area))*C_y
return np.array([[C_x, C_y]])
def relax_points(self, times=1):
"""
Relaxes the points after an initial Voronoi is created to refine the graph.
See: https://stackoverflow.com/questions/17637244/voronoi-and-lloyd-relaxation-using-python-scipy
:param times: Number of times to relax, default is 1
:type times: int
:return: the final voronoi diagrama
:rtype: scipy.spatial.Voronoi
"""
for i in range(times):
centroids = []
for region in self.filtered_regions:
vertices = self.vor.vertices[region + [region[0]], :]
centroid = self._region_centroid(vertices)
centroids.append(list(centroid[0, :]))
self.points = centroids
self.generate_voronoi()
return self.vor
def generate_voronoi(geoseries_polygons):
"""Generate Voronoi polygons from polygon edges
:param geoseries_polygons: GeoSeries of raw polygons
:return:
"""
edges = geoseries_polygons.unary_union.boundary
pnts = points_along_boundaries(edges, 0.75)
cent = pnts.centroid
tpnts = translate(pnts, -cent.x, -cent.y)
vor = Voronoi(tpnts)
polys = []
for region in vor.regions:
if len(region) > 0 and all([i > 0 for i in region]):
polys.append(Polygon([vor.vertices[i] for i in region]))
gs_vor = geopandas.GeoSeries(polys)
t_gs_vor = gs_vor.translate(cent.x, cent.y)
t_gs_vor.crs = geoseries_polygons.crs
return t_gs_vor, pnts
def vorEdges(vor, far):
"""
Given a voronoi tesselation, retuns the set of voronoi edges.
far is the length of the "infinity" edges
"""
edges = []
for simplex in vor.ridge_vertices:
simplex = numpy.asarray(simplex)
if numpy.all(simplex >= 0):
edge = {}
edge['p1'], edge['p2'] = vor.vertices[simplex, 0], vor.vertices[simplex, 1]
edge['p1'] = numpy.array([vor.vertices[simplex, 0][0], vor.vertices[simplex, 1][0]])
edge['p2'] = numpy.array([vor.vertices[simplex, 0][1], vor.vertices[simplex, 1][1]])
edge['t'] = (edge['p2'] - edge['p1']) / numpy.linalg.norm(edge['p2'] - edge['p1'])
edges.append(edge)
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 = numpy.asarray(simplex)
if numpy.any(simplex < 0):
i = simplex[simplex >= 0][0] # finite end Voronoi vertex
t = vor.points[pointidx[1]] - vor.points[pointidx[0]] # tangent
t /= numpy.linalg.norm(t)
n = numpy.array([-t[1], t[0]]) # normal
midpoint = vor.points[pointidx].mean(axis=0)
direction = numpy.sign(numpy.dot(midpoint - center, n)) * n
far_point = vor.vertices[i] + direction * ptp_bound.max() * far
edge = {}
edge['p1'], edge['p2'] = numpy.array([vor.vertices[i, 0], far_point[0]]), numpy.array(
[vor.vertices[i, 1], far_point[1]])
edge['p1'], edge['p2'] = vor.vertices[i, :], far_point
edge['t'] = (edge['p2'] - edge['p1']) / numpy.linalg.norm(edge['p2'] - edge['p1'])
edges.append(edge)
return edges
def phase_diagram_mesh(points, values,
title="Phase diagram",
xlabel="Pulse Duration (s)",
ylabel="Pulse Amplitude (V)",
shading="flat",
voronoi=False, **kwargs):
# fig = plt.figure()
if voronoi:
from scipy.spatial import Voronoi, voronoi_plot_2d
points[:,0] *= 1e9
vor = Voronoi(points)
cmap = mpl.cm.get_cmap('RdGy')
# colorize
for pr, v in zip(vor.point_region, values):
region = vor.regions[pr]
if not -1 in region:
polygon = [vor.vertices[i] for i in region]
plt.fill(*zip(*polygon), color=cmap(v))
else:
mesh = scaled_Delaunay(points)
xs = mesh.points[:,0]
ys = mesh.points[:,1]
plt.tripcolor(xs,ys,mesh.simplices.copy(),values, cmap="RdGy",shading=shading,**kwargs)
plt.xlim(min(xs),max(xs))
plt.ylim(min(ys),max(ys))
plt.title(title, size=18)
plt.xlabel(xlabel, size=16)
plt.ylabel(ylabel, size=16)
cb = plt.colorbar()
cb.set_label("Probability",size=16)
return mesh
def polyhedron(self, coords, j, L):
"""find the polyhedron for center molecule"""
vor = Voronoi(coords)
#get the vertices
points = [vor.vertices[x] for x in vor.regions[vor.point_region[j]] if x != -1]
return points
def compute_vc(self, points):
"""compute the Voronoi cell"""
#compute S and V
S = ConvexHull(points).area
V = ConvexHull(points).volume
#voronoi cell
eta = S ** 3 / (36 * np.pi * V ** 2)
return eta
def asphericity(self, freq = 1):
"""compute asphericity of the Voronoi cell"""
#progress bar
frames = int(self.traj.n_frames / freq)
bar = ChargingBar('Processing', max=frames,
suffix='%(percent).1f%% - %(eta)ds')
for i in range(0, self.traj.n_frames, freq):
for j in range(self.traj.n_atoms):
if self.traj.atom_names[i][j] == self.center:
#center coordinate
c = self.traj.coords[i][j]
#coordinates
cs = self.traj.coords[i]
#box_size
L = self.traj.box_size[i]
#new coordinates after wrapping
nc = self.wrap_box(c, cs, L)
points = self.polyhedron(nc, j, L)
e = self.compute_vc(points)
self.raw.append(e)
bar.next()
bar.finish()
def density_voronoi_2(points, cell_size):
r"""Density definition base on Voronoi-diagram. [Steffen2010]_
Compute Voronoi-diagram for each position :math:`\mathbf{p}_i` giving
cells :math:`A_i` to each agent :math:`i`. Area of Voronoi cell is
denoted :math:`|A_i|`. Area :math:`A` is the area inside which we
measure the density.
:math:`S = \{i : \mathbf{p}_i \in A\}` is the set of agents and
:math:`N = |S|` the number of agents that is inside area :math:`A`.
.. math::
D_{V^\prime} = \frac{N}{\sum_{i \in S} |A_i|}
Args:
points (numpy.ndarray):
Two dimensional points :math:`\mathbf{p}_{i \in \mathbb{N}}`
cell_size (float):
Cell size of the meshgrid. Each cell represents an area :math:`A`
inside which density is measured.
Returns:
numpy.ndarray: Grid of values density values for each cell.
"""
# Voronoi tesselation
vor = Voronoi(points)
new_regions, new_vertices = voronoi_finite_polygons_2d(vor)
# Compute areas for all of the Voronoi regions
area = np.array([polygon_area(new_vertices[i]) for i in new_regions])
return _core_2(points, cell_size, area, vor.point_region)
def generate_voronoi(self):
"""
Uses scipy.spatial.Voronoi to generate a voronoi diagram.
Filters viable regions and stashes them in filtered_regions, see https://stackoverflow.com/questions/28665491/getting-a-bounded-polygon-coordinates-from-voronoi-cells
:return: A voronoi diagram based on the points
:rtype: scipy.spatial.Voronoi
"""
eps = sys.float_info.epsilon
self.vor = Voronoi(self.points)
self.filtered_regions = []
for region in self.vor.regions:
flag = True
for index in region:
if index == -1:
flag = False
break
else:
x = self.vor.vertices[index, 0]
y = self.vor.vertices[index, 1]
if not (self.bounding_region[0] - eps <= x and x <= self.bounding_region[1] + eps and
self.bounding_region[2] - eps <= y and y <= self.bounding_region[3] + eps):
flag = False
break
if region != [] and flag:
self.filtered_regions.append(region)
return self.vor
def medial_axis(samples, contains):
'''
Given a set of samples on a boundary, find the approximate medial axis based
on a voronoi diagram and a containment function which can assess whether
a point is inside or outside of the closed geometry.
Arguments
----------
samples: (n,d) set of points on the boundary of the geometry
contains: function which takes (m,d) points and returns an (m) bool array
Returns
----------
lines: (n,2,2) set of line segments
'''
from scipy.spatial import Voronoi
from .path.io.load import load_path
# create the voronoi diagram, after vertically stacking the points
# deque from a sequnce into a clean (m,2) array
voronoi = Voronoi(samples)
# which voronoi vertices are contained inside the original polygon
contained = contains(voronoi.vertices)
# ridge vertices of -1 are outside, make sure they are False
contained = np.append(contained, False)
inside = [i for i in voronoi.ridge_vertices if contained[i].all()]
line_indices = np.vstack([stack_lines(i) for i in inside if len(i) >=2])
lines = voronoi.vertices[line_indices]
return load_path(lines)
