Coverage for src/thunderlab/voronoi.py: 98%
460 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-20 15:01 +0000
1"""Analyze Voronoi diagrams based on scipy.spatial.
3## Classes
5- `class Voronoi`: Compute and analyse Voronoi diagrams.
6"""
9import numpy as np
10import scipy.spatial as sp
11try:
12 import matplotlib.pyplot as plt
13except ImportError:
14 pass
17def cross2d(x, y):
18 # numpy version 2.5 up require 3-dim vectors for np.cross()...
19 return x[..., 0] * y[..., 1] - x[..., 1] * y[..., 0]
22class Voronoi(object):
23 """
24 Parameters
25 ----------
26 points: ndarray of floats, shape (npoints, 2)
27 List of point coordinates.
28 radius: float or None
29 Radius for computing far points of infinite ridges.
30 If None twice the maximum extent of the input points is used.
31 qhull_options: string or None
32 Options to be passed on to QHull. From the manual:
34 - Qbb - scale last coordinate to [0,m] for Delaunay triangulations
35 - Qc - keep coplanar points with nearest facet
36 - Qx - exact pre-merges (skips coplanar and anglomaniacs facets)
37 - Qz - add point-at-infinity to Delaunay triangulation
38 - QJn - randomly joggle input in range [-n,n]
39 - Qs - search all points for the initial simplex
40 - Qz - add point-at-infinity to Voronoi diagram
41 - QGn - Voronoi vertices if visible from point n, -n if not
42 - QVn - Voronoi vertices for input point n, -n if not
44 Default is: "Qbb Qc Qz Qx" for ndim > 4 and "Qbb Qc Qz" otherwise.
46 Input points
47 ------------
48 The points from which the Voronoi diagram is constructed.
50 - `ndim`: int
51 The dimension of the input points, i.e. number of coordinates.
52 - `npoints`: int
53 The number of input points.
54 - `points`: 2-d ndarray
55 List of input point coordinates.
56 - `center`: ndarray of floats
57 Center of mass of the input points (i.e. mean coordinates).
58 - `min_bound`: ndarray of floats
59 Lower left corner of the bounding-box of the input points.
60 - `max_bound`: ndarray of floats
61 Upper right corner of the bounding-box of the input points.
63 Distances between input points
64 ------------------------------
65 - `ridge_points`: 2-d ndarray of ints
66 List of pairs of indices to `points` enclosing a Voronoi ridge.
67 - `ridge_distances`: ndarray of floats
68 For each ridge in `ridge_points` the distance between the enclosing points.
69 - `neighbor_points`: list of ndarrays of ints
70 For each point in `points` a list of indices of the Voronoi-cell's neighboring points.
71 - `neighbor_distances`: list of ndarrays of floats
72 For each point in `points` a list of distances to the Voronoi-cell's neighboring points,
73 matching `neighbor_points`.
74 - `nearest_points`: list of ints
75 For each point in `points` the index of its nearest neighbor.
76 - `nearest_distances`: ndarray of floats
77 For each point in `points` the distance to its nearest neighbor.
78 - `mean_nearest_distance`: float
79 Unbiased estimate of the mean nearest-neighbor distance. This is
80 the average nearest-neighbor distance corrected by one s.e.m,
81 since points close to the border of the region have too large
82 nearest-neighbor distances.
84 Voronoi diagram
85 ---------------
86 - `vertices`: 2-d ndarray of floats
87 List of vertex coordinates enclosing the Voronoi regions.
88 - `regions`: list of list of ints
89 List of lists of indices to vertices in `vertices` making up each Voronoi region.
90 - `ridge_vertices`: list of list of ints
91 List of lists of indices to `vertices` making up a Voronoi ridge.
92 The pairs of `points` to poth sides of the ridge are listed in `ridge_points`.
93 - `ridge_lengths()`: Length of Voronoi ridges between nearest neighbors.
94 - `areas()`: The areas of the Voronoi regions for each input point in `points`.
95 - `point_types()`: The type of Voronoi area (infinite, finite, inside) for each input point.
96 - `inside_vertices`: list of ints
97 Indices of `vertices` that are inside the convex hull of the input points.
98 - `infinite_vertices`: list of ndarrays of floats
99 For each ridge in `ridge_vertices` the coordinates of the far-point, if it is an infinite ridge.
100 - `infinite_regions`: list of list of ints
101 List of Voronoi regions with infinite ridges. If positive, the indices are indices to `vertices`.
102 If negative they are indices into `infinite_vertices` (`-index-1`).
104 Convex hull
105 -----------
106 - `hull_points`: list of ints
107 List of indices of the points in `points` making up the convex hull.
108 - `hull_center`: ndarray of floats
109 Center of mass of the points making up the convex hull.
110 - `inside_vertices`: ndarray of boolean
111 Indicates for each vertex in `vertices` whether it is inside the convex hull.
112 - `in_hull()`: Test if points are within the convex hull of the input points.
113 - `hull_area()`: The area contained in the convex hull.
115 Outer hull
116 ----------
117 The outer hull is constructed from the convex hull by expanding each
118 point of the convex hull away from the center of mass of the convex
119 hull such that randomly placing the same number of points in this
120 area results on average in the same mean nearest-neighbor
121 distance. This enlarged hull is needed for bootstrapping. Note that
122 in some cases the outer hull can be *smaller* than the convex hull.
124 - `outer_hull_points`: 2-d ndarray of floats
125 Array of coordinates of the points of the outer hull.
126 - `outer_min_bound`: ndarray of floats
127 Lower left corner of the bounding-box of the outer hull.
128 - `outer_max_bound`: ndarray of floats
129 Upper right corner of the bounding-box of the outer hull.
130 - `in_outer_hull()`: Test if points are within the outer hull.
131 - `outer_hull_area()`: The area contained in the outer hull.
133 Bootstrap Voronoi diagrams
134 --------------------------
135 - `random_points()`: Generate random points within the area defined by the input points.
137 Plotting the Voronoi diagram
138 ----------------------------
139 - `plot_points()`: Plot and optionally annotate the input points of the Voronoi diagram.
140 - `plot_center()`: Plot the center of mass of the input points.
141 - `plot_vertices()`: Plot and optionally annotate the vertices of the Voronoi diagram.
142 - `plot_distances()`: Plot lines connecting the neighbors in the Voronoi diagram.
143 - `plot_ridges()`: Plot the finite ridges of the Voronoi diagram.
144 - `plot_infinite_ridges(): Plot the infinite ridges of the Voronoi diagram.
145 - `fill_regions()`: Fill each finite region of the Voronoi diagram with a color.
146 - `fill_infinite_regions()`: Fill each infinite region of the Voronoi diagram with a color.
148 Plotting the convex hull
149 ------------------------
150 - `plot_hull()`: Plot the convex hull line containing the input points.
151 - `fill_hull()`: Fill the convex hull.
152 - `plot_hull_center()`: Plot the center of mass of the convex hull.
153 - `plot_outer_hull()`: Plot the outer hull edge.
154 - `fill_outer_hull()`: Fill the outer hull.
156 Usage
157 -----
158 Generate 20 random points in 2-D:
159 ```
160 import numpy as np
161 points = np.random.rand(20, 2)
162 ```
164 Calculate the Voronoi diagram:
165 ```
166 from thunderlab.voronoi import Voronoi
167 vor = Voronoi(points)
168 ```
170 Retrieve nearest-neighbor distances and compute Voronoi areas:
171 ```
172 dists = vor.nearest_distances
173 areas = vor.areas()
174 ```
176 Generate random points drawn from a Poisson process with the same
177 intensity and mean nearest neighbor distance as the data:
178 ```
179 new_points = vor.random_points(poisson=True, mode='outer')
180 ```
182 Plot Voronoi regions, distances, and input points:
183 ```
184 import matplotlib.pyplot as plt
185 vor.fill_regions(colors=['red', 'green', 'blue', 'orange', 'cyan'], alpha=0.3)
186 vor.plot_distances(color='red')
187 vor.plot_points(text='p%d', c='c', s=100)
188 ```
189 """
191 def __init__(self, points, radius=None, qhull_options=None):
192 self.vor = sp.Voronoi(points, furthest_site=False, incremental=False,
193 qhull_options=qhull_options)
194 if self.vor.ndim != 2:
195 raise ValueError("Only 2D input points are supported.")
196 self.hull = sp.Delaunay(points, furthest_site=False, incremental=False,
197 qhull_options=qhull_options)
198 self._compute_distances()
199 self._compute_infinite_vertices()
200 self._compute_hull(qhull_options)
201 self.ndim = self.vor.ndim
202 self.npoints = self.vor.npoints
203 self.points = self.vor.points
204 self.vertices = self.vor.vertices
205 self.regions = self.vor.regions
206 self.ridge_points = self.vor.ridge_points
207 self.ridge_vertices = self.vor.ridge_vertices
208 self.min_bound = self.vor.min_bound
209 self.max_bound = self.vor.max_bound
210 self.center = np.mean(self.points, axis=0)
212 def _compute_distances(self):
213 """Compute distances between points.
214 """
215 # For each ridge the distance of the points enclosing the ridge:
216 p1 = self.vor.points[self.vor.ridge_points[:,0]]
217 p2 = self.vor.points[self.vor.ridge_points[:,1]]
218 self.ridge_distances = sp.minkowski_distance(p1, p2)
220 # For each point all its Voronoi distances:
221 self.neighbor_points = [[] for k in range(len(self.vor.points))]
222 self.neighbor_distances = [[] for k in range(len(self.vor.points))]
223 for dist, points in zip(self.ridge_distances, self.vor.ridge_points):
224 self.neighbor_points[points[0]].append(points[1])
225 self.neighbor_points[points[1]].append(points[0])
226 self.neighbor_distances[points[0]].append(dist)
227 self.neighbor_distances[points[1]].append(dist)
228 for k in range(len(self.neighbor_points)):
229 inx = np.argsort(self.neighbor_distances[k])
230 self.neighbor_points[k] = np.array(self.neighbor_points[k])[inx]
231 self.neighbor_distances[k] = np.array(self.neighbor_distances[k])[inx]
233 # For each point the distance to its nearest neighbor:
234 self.nearest_points = []
235 self.nearest_distances = np.zeros(len(self.neighbor_distances))
236 for k in range(len(self.neighbor_distances)):
237 self.nearest_points.append(self.neighbor_points[k][0])
238 self.nearest_distances[k] = self.neighbor_distances[k][0]
239 self.mean_nearest_distance = np.mean(self.nearest_distances)
240 self.mean_nearest_distance *= 1.0-2.0*np.sqrt(1.0/np.pi-0.25)/np.sqrt(self.vor.npoints)
242 def _compute_infinite_vertices(self, radius=None):
243 """Compute far points of infinite ridges.
245 Parameters
246 ----------
247 radius: float or None
248 Radius for computing far points of infinite ridges.
249 If None twice the maximum extent of the input points is used.
251 Note
252 ----
253 Code inspired by http://stackoverflow.com/questions/20515554/colorize-voronoi-diagram
254 """
255 # For each ridge, compute far point:
256 center = self.vor.points.mean(axis=0)
257 if radius is None:
258 radius = np.max(np.ptp(2.0*self.vor.points, axis=0))
259 self.infinite_vertices = []
260 for points, vertices in zip(self.vor.ridge_points, self.vor.ridge_vertices):
261 vertices = np.asarray(vertices)
262 if np.all(vertices >= 0):
263 self.infinite_vertices.append(np.array([]))
264 else:
265 i = vertices[vertices >= 0][0] # finite end Voronoi vertex
266 t = self.vor.points[points[1]] - self.vor.points[points[0]] # tangent
267 t /= np.linalg.norm(t)
268 n = np.array([-t[1], t[0]]) # normal
269 midpoint = self.vor.points[points].mean(axis=0)
270 direction = np.sign(np.dot(midpoint - center, n)) * n
271 far_point = self.vor.vertices[i] + direction * radius
272 self.infinite_vertices.append(far_point)
273 # Assemble list of infinite regions:
274 # Indices to self.infinite_vertices are negative minus one.
275 self.infinite_regions = []
276 for rvertices in self.vor.regions:
277 if -1 in rvertices:
278 new_rvertices = []
279 prev_vertex = rvertices[-1]
280 # find index of data point enclosed by the region:
281 ridge_points = []
282 for p, v in zip(self.vor.ridge_points, self.vor.ridge_vertices):
283 if not -1 in v and set(v) <= set(rvertices):
284 ridge_points.extend(p)
285 region_point = None
286 for rp in ridge_points:
287 if ridge_points.count(rp) > 1:
288 region_point = rp
289 break
290 # fill in far points for each region:
291 for v_inx, v in enumerate(rvertices):
292 if v >= 0:
293 new_rvertices.append(v)
294 else:
295 for v1_inx, (points, vertices) in enumerate(zip(self.vor.ridge_points,
296 self.vor.ridge_vertices)):
297 if prev_vertex in vertices and -1 in vertices and \
298 (region_point is None or region_point in points):
299 new_rvertices.append(-v1_inx-1)
300 break
301 next_vertex = rvertices[0]
302 if v_inx+1 < len(rvertices):
303 next_vertex = rvertices[v_inx+1]
304 for v2_inx, (points, vertices) in enumerate(zip(self.vor.ridge_points,
305 self.vor.ridge_vertices)):
306 if next_vertex in vertices and -1 in vertices and \
307 (region_point is None or region_point in points) and \
308 new_rvertices[-1] != -v2_inx-1:
309 new_rvertices.append(-v2_inx-1)
310 break
311 prev_vertex = v
312 self.infinite_regions.append(new_rvertices)
314 def _flatten_simplices(self, simplices):
315 """Transforms list of simplex indices to list of vertex indices.
317 In particular, transforms the Delaunay.convex_hull to a list of points
318 of the hull that can then be easily plotted.
320 Parameters
321 ----------
322 simplices: 2-D array of ints
323 List of pairs of indices of points forming each ridge of a polygon.
325 Returns
326 -------
327 indices: list of ints
328 Indices of vertices of the polygon.
329 """
330 if len(simplices) == 0:
331 return []
332 indices = list(simplices[0])
333 simplices = np.delete(simplices, 0, 0)
334 while len(simplices) > 0:
335 for i, s in enumerate(simplices):
336 if indices[-1] in s:
337 if s[0] == indices[-1]:
338 indices.append(s[1])
339 else:
340 indices.append(s[0])
341 simplices = np.delete(simplices, i, 0)
342 break
343 return indices
345 def _compute_hull(self, qhull_options):
346 """Compute properties of the convex hull and set up the outer hull.
347 """
348 self.inside_vertices = self.in_hull(self.vor.vertices)
349 self.hull_points = self._flatten_simplices(self.hull.convex_hull)
350 self.hull_center = np.mean(self.hull.points[self.hull_points], axis=0)
352 # enlarge hull:
353 # estimate area needed for a poisson process with same mean distance:
354 new_area = 4.0 * self.mean_nearest_distance**2 * self.vor.npoints
355 fac = np.sqrt(new_area/self.hull_area())
356 self.outer_hull_points = np.zeros((len(self.hull_points), self.vor.ndim))
357 for k, point in enumerate(self.hull.points[self.hull_points]):
358 point -= self.hull_center
359 point *= fac
360 point += self.hull_center
361 self.outer_hull_points[k] = point
362 # compute outer hull:
363 self.outer_hull = sp.Delaunay(self.outer_hull_points,
364 furthest_site=False, incremental=False,
365 qhull_options=qhull_options)
366 self.outer_min_bound = np.min(self.outer_hull_points, axis=0)
367 self.outer_max_bound = np.max(self.outer_hull_points, axis=0)
369 def in_hull(self, p):
370 """Test if points `p` are within convex hull of the input points.
372 Parameters
373 ----------
374 p: 2-d ndarray
375 Array of points to be tested.
377 Returns
378 -------
379 inside: ndarray of booleans
380 For each point in `p` whether it is inside the hull.
381 """
382 inside = self.hull.find_simplex(p) >= 0
383 return inside
385 def in_outer_hull(self, p):
386 """Test if points `p` are within the outer hull.
388 Parameters
389 ----------
390 p: 2-d ndarray
391 Array of points to be tested.
393 Returns
394 -------
395 inside: ndarray of booleans
396 For each point in `p` whether it is inside the outer hull.
397 """
398 inside = self.outer_hull.find_simplex(p) >= 0
399 return inside
401 def point_types(self):
402 """The type of the Voronoi regions for each input point.
404 Returns
405 -------
406 points: ndarray of ints
407 Indicates the type of Voronoi region associated
408 with each point in `points`:
410 - 2: finite region with all vertices inside hull,
411 - 1: finite region,
412 - 0: infinite region.
413 """
414 points = np.zeros(len(self.vor.points), dtype=int) + 2
415 for i in range(len(self.vor.points)):
416 for j, rp in enumerate(self.vor.ridge_points):
417 if i in rp:
418 if not np.all(self.inside_vertices[self.vor.ridge_vertices[j]]) and\
419 points[i] > 0:
420 points[i] = 1
421 if not np.all(np.array(self.vor.ridge_vertices[j]) >= 0) and\
422 points[i] > -1:
423 points[i] = 0
424 return points
426 def ridge_lengths(self):
427 """Length of Voronoi ridges between nearest neighbors.
429 May be used, for example, as a weigth for `ridge_distances`.
431 Returns
432 -------
433 distances: ndarray of floats
434 The length of each ridge in `ridge_vertices`.
435 np.inf if one vertex is at infinity.
436 """
437 ridges = np.zeros(len(self.vor.ridge_vertices))
438 for k, p in enumerate(self.vor.ridge_vertices):
439 if np.all(np.array(p)>=0):
440 p1 = self.vor.vertices[p[0]]
441 p2 = self.vor.vertices[p[1]]
442 ridges[k] = sp.minkowski_distance(p1, p2)
443 else:
444 ridges[k] = np.inf
445 return ridges
447 def ridge_areas(self):
448 """For each ridge the triangular area of the Voronoi region
449 spanned by the center point and the ridge.
451 Returns
452 -------
453 areas: ndarray of floats
454 For each ridge in `ridge_points` or `ridge_vertices`
455 its corresponding triangular area.
456 np.inf for infinite ridges.
457 """
458 ridges = self.ridge_lengths()
459 heights = 0.5*self.ridge_distances
460 # area of a triangle:
461 areas = 0.5*ridges*heights
462 return areas
464 def areas(self, mode='finite'):
465 """The areas of the Voronoi regions for each input point.
467 Parameters
468 ----------
469 mode: string
470 - 'inside': Calculate area of finite Voronoi regions
471 whose vertices are all inside the hull,
472 set all other to zero.
473 - 'finite_inside': Calculate area of all Voronoi regions.
474 Consider only areas of finite ridges
475 whose vertices are all inside the hull.
476 - 'full': Calculate area of finite Voronoi regions only,
477 set all other to zero.
478 - 'finite': Calculate area of all Voronoi regions.
479 From infinite regions
480 only areas contributed by finite ridges are considered.
482 Returns
483 -------
484 areas: array of floats
485 For each point in `points` its corresponding area.
486 """
487 ridge_areas = self.ridge_areas()
488 areas = np.zeros(len(self.vor.points))
489 if mode == 'inside':
490 for i in range(len(self.vor.points)):
491 a = 0.0
492 for j, rp in enumerate(self.vor.ridge_points):
493 if i in rp:
494 if ridge_areas[j] != np.inf and \
495 np.all(self.inside_vertices[self.vor.ridge_vertices[j]]):
496 a += ridge_areas[j]
497 else:
498 a = 0.0
499 break
500 areas[i] = a
501 elif mode == 'full':
502 for i in range(len(self.vor.points)):
503 a = 0.0
504 for j, rp in enumerate(self.vor.ridge_points):
505 if i in rp:
506 if ridge_areas[j] != np.inf:
507 a += ridge_areas[j]
508 else:
509 a = 0.0
510 break
511 areas[i] = a
512 elif mode == 'finite_inside':
513 for i in range(len(self.vor.points)):
514 a = 0.0
515 for j, rp in enumerate(self.vor.ridge_points):
516 if i in rp and ridge_areas[j] != np.inf and \
517 np.all(self.inside_vertices[self.vor.ridge_vertices[j]]):
518 a += ridge_areas[j]
519 areas[i] = a
520 elif mode == 'finite':
521 for i in range(len(self.vor.points)):
522 a = 0.0
523 for j, rp in enumerate(self.vor.ridge_points):
524 if i in rp and ridge_areas[j] != np.inf:
525 a += ridge_areas[j]
526 areas[i] = a
527 else:
528 print('')
529 print(f'Voronoi.areas(): unknown value "{mode}" for the mode parameter:')
530 print('Use one of the following values:')
531 print(' inside: Finite Voronoi regions whose vertices are all inside the hull.')
532 print(' finite_inside: Use all areas corresponding to finite ridges whose vertices are all inside the hull.')
533 print(' full: Finite Voronoi regions only.')
534 print(' finite: Use all areas corresponding to finite ridges.')
535 print('')
536 return areas
538 def hull_area(self):
539 """The area of the convex hull of the input points.
541 Returns
542 -------
543 area: float
544 The area of the convex hull.
545 """
546 # two sides of the simplex triangles:
547 ab = self.hull.points[self.hull.simplices[:,0],:] - self.hull.points[self.hull.simplices[:,1],:]
548 cb = self.hull.points[self.hull.simplices[:,2],:] - self.hull.points[self.hull.simplices[:,1],:]
549 # area of each simplex is half of the absolute value of the cross product:
550 area = 0.5*np.sum(np.abs(cross2d(ab, cb)))
551 return area
553 def outer_hull_area(self):
554 """The area of the outer hull.
556 Returns
557 -------
558 area: float
559 The area of the outer hull.
560 """
561 # two sides of the simplex triangles:
562 ab = self.outer_hull.points[self.outer_hull.simplices[:,0],:] - self.outer_hull.points[self.outer_hull.simplices[:,1],:]
563 cb = self.outer_hull.points[self.outer_hull.simplices[:,2],:] - self.outer_hull.points[self.outer_hull.simplices[:,1],:]
564 # area of each simplex is half of the absolute value of the cross product:
565 area = 0.5*np.sum(np.abs(cross2d(ab, cb)))
566 return area
568 def random_points(self, n=None, poisson=True, mode='outer'):
569 """Generate random points.
571 Parameters
572 ----------
573 n: int or None
574 Number of random points to be generated.
575 If None `n` is set to the number of input points.
576 poisson: boolean
577 If `True` then draw the number of points from a Poisson distribution
578 with mean number of points given by `n`.
579 mode: string
580 - 'bbox': place points randomly in rectangular bounding box
581 of the Voronoi diagram.
582 - 'hull': place points randomly within convex hull of input data.
583 - 'outer' place points randomly within outer hull.
584 The mean nearest-neighbor distance
585 between the generated points will be close
586 to the observed one.
588 Returns
589 -------
590 points: 2-D array
591 List of randomly generated points.
592 """
593 # number of points:
594 if n is None:
595 n = self.npoints
596 nn = n
597 if poisson:
598 nn = np.random.poisson(n)
599 m = nn//2
600 if m < 5:
601 m = 5
602 # get bounding box:
603 if mode == 'outer':
604 min_bound = self.outer_min_bound
605 max_bound = self.outer_max_bound
606 else:
607 min_bound = self.min_bound
608 max_bound = self.max_bound
609 delta = np.max(max_bound - min_bound)
610 points = np.zeros((0, self.ndim))
611 while len(points) < nn:
612 # random points within bounding box:
613 newpoints = np.random.rand(m, self.ndim)
614 newpoints *= delta
615 newpoints += min_bound
616 if mode == 'outer':
617 # only take the ones within outer hull:
618 inside = self.in_outer_hull(newpoints)
619 points = np.vstack((points, newpoints[inside]))
620 elif mode == 'hull':
621 # only take the ones within hull:
622 inside = self.in_hull(newpoints)
623 points = np.vstack((points, newpoints[inside]))
624 elif mode == 'bbox':
625 points = np.vstack((points, newpoints[np.all(newpoints<max_bound, axis=1),:]))
626 else:
627 print('')
628 print(f'Voronoi.random_points(): unknown value "{mode}" for the mode parameter:')
629 print('Use one of the following values:')
630 print(' bbox: Place points within rectangular bounding box.')
631 print(' hull: Place points inside the convex hull.')
632 print(' outer: Place points inside the outer hull.')
633 print('')
634 return
635 return points[:nn]
637 def plot_points(self, ax, text=None, text_offs=(0, 0.05),
638 text_align='center', **kwargs):
639 """Plot and optionally annotate the input points of the Voronoi diagram.
641 Parameters
642 ----------
643 ax: matplotlib.Axes
644 The axes to be used for plotting.
645 text: string or None
646 If not None the string that is placed at each point.
647 A '%d' is replaced by the index of the point.
648 text_offset: tuple of numbers
649 The offset of the point labels.
650 text_align: string
651 The horizontal alignment of the point labels.
652 **kwargs: dict
653 Key-word arguments that are passed on to the matplotlib.scatter() function.
654 """
655 ax.scatter(self.points[:,0], self.points[:,1], **kwargs)
656 if text is not None:
657 for i, p in enumerate(self.points):
658 s = text
659 if '%' in text:
660 s = text % i
661 ax.text(p[0]+text_offs[0], p[1]+text_offs[1], s, ha=text_align)
663 def plot_center(self, ax, **kwargs):
664 """Plot the center of mass of the input points.
666 Parameters
667 ----------
668 ax: matplotlib.Axes
669 The axes to be used for plotting.
670 **kwargs: dict
671 Key-word arguments that are passed on to the `matplotlib.plot()` function.
672 """
673 ax.plot(self.center[0], self.center[1], 'o', **kwargs)
675 def plot_vertices(self, ax, text=None, text_offs=(0, 0.05),
676 text_align='center', **kwargs):
677 """Plot and optionally annotate the vertices of the Voronoi diagram.
679 Parameters
680 ----------
681 ax: matplotlib.Axes
682 The axes to be used for plotting.
683 text: string or None
684 If not None the string that is placed at each vertex.
685 A '%d' is replaced by the index of the vertex.
686 text_offset: tuple of numbers
687 The offset of the vertex labels.
688 text_align: string
689 The horizontal alignment of the vertex labels.
690 **kwargs: dict
691 Key-word arguments that are passed on to the matplotlib.scatter() function.
692 """
693 ax.scatter(self.vertices[:,0], self.vertices[:,1], **kwargs)
694 if text is not None:
695 for i, p in enumerate(self.vertices):
696 s = text
697 if '%' in text:
698 s = text % i
699 ax.text(p[0]+text_offs[0], p[1]+text_offs[1], s, ha=text_align)
701 def plot_distances(self, ax, **kwargs):
702 """Plot lines connecting the nearest neighbors in the Voronoi diagram.
704 Parameters
705 ----------
706 ax: matplotlib.Axes
707 The axes to be used for plotting.
708 **kwargs: dict
709 Key-word arguments that are passed on to the matplotlib.plot() function.
710 """
711 for i, p in enumerate(self.ridge_points):
712 ax.plot(self.points[p, 0], self.points[p, 1], **kwargs)
714 def plot_ridges(self, ax, **kwargs):
715 """Plot the finite ridges of the Voronoi diagram.
717 Parameters
718 ----------
719 ax: matplotlib.Axes or None
720 The axes to be used for plotting.
721 **kwargs: dict
722 Key-word arguments that are passed on to the matplotlib.plot() function.
723 """
724 for i, p in enumerate(self.ridge_vertices):
725 if np.all(np.array(p)>=0):
726 ax.plot(self.vertices[p, 0], self.vertices[p, 1], **kwargs)
728 def plot_infinite_ridges(self, ax, **kwargs):
729 """Plot the infinite ridges of the Voronoi diagram.
731 Parameters
732 ----------
733 ax: matplotlib.Axes
734 The axes to be used for plotting.
735 **kwargs: dict
736 Key-word arguments that are passed on to the matplotlib.plot() function.
737 """
738 for far_point, vertices in zip(self.infinite_vertices, self.vor.ridge_vertices):
739 vertices = np.asarray(vertices)
740 if not np.all(vertices >= 0):
741 i = vertices[vertices >= 0][0] # finite end Voronoi vertex
742 ax.plot([self.vor.vertices[i][0], far_point[0]],
743 [self.vor.vertices[i][1], far_point[1]], **kwargs)
745 def fill_regions(self, ax, inside=None, colors=None, **kwargs):
746 """Fill each finite region of the Voronoi diagram with a color.
748 Parameters
749 ----------
750 ax: matplotlib.Axes
751 The axes to be used for plotting.
752 inside: boolean or None
753 True: plot only finite regions with all vertices inside the hull
754 False: plot only finite regions with at least one vertex outside the hull
755 None: plot all finite regions
756 colors: list of colors or None
757 If not None then these colors are used in turn to fill the regions.
758 **kwargs: dict
759 Key-word arguments that are passed on to the matplotlib.fill() function.
760 """
761 c = 0
762 for region in self.regions:
763 if not -1 in region:
764 polygon = self.vertices[region]
765 if len(polygon) > 0:
766 inside_hull = self.in_hull(polygon)
767 if inside is None or (inside and all(inside_hull)) or (not inside and any(inside_hull)):
768 c += 1
769 if colors is None:
770 ax.fill(polygon[:, 0], polygon[:, 1], lw=0, **kwargs)
771 else:
772 ax.fill(polygon[:, 0], polygon[:, 1],
773 color=colors[c % len(colors)], lw=0, **kwargs)
775 def fill_infinite_regions(self, ax, colors=None, **kwargs):
776 """Fill each infinite region of the Voronoi diagram with a color.
778 Parameters
779 ----------
780 ax: matplotlib.Axes
781 The axes to be used for plotting.
782 colors: list of colors or None
783 If not None then these colors are used in turn to fill the regions.
784 **kwargs: dict
785 Key-word arguments that are passed on to the matplotlib.fill() function.
786 """
787 c = 0
788 for region in self.infinite_regions:
789 polygon = []
790 for p in region:
791 if p >= 0:
792 polygon.append(self.vertices[p])
793 else:
794 polygon.append(self.infinite_vertices[-p-1])
795 if len(polygon) > 0:
796 polygon = np.asarray(polygon)
797 c += 1
798 if colors is None:
799 ax.fill(polygon[:, 0], polygon[:, 1], lw=0, **kwargs)
800 else:
801 ax.fill(polygon[:, 0], polygon[:, 1],
802 color=colors[c % len(colors)], lw=0, **kwargs)
804 def plot_hull(self, ax, **kwargs):
805 """Plot the hull line containing the input points.
807 Parameters
808 ----------
809 ax: matplotlib.Axes
810 The axes to be used for plotting.
811 **kwargs: dict
812 Key-word arguments that are passed on to the matplotlib.plot() function.
813 """
814 ax.plot(self.hull.points[self.hull_points, 0],
815 self.hull.points[self.hull_points, 1], **kwargs)
817 def fill_hull(self, ax, **kwargs):
818 """Fill the hull containing the input points with a color.
820 Parameters
821 ----------
822 ax: matplotlib.Axes
823 The axes to be used for plotting.
824 **kwargs: dict
825 Key-word arguments that are passed on to the matplotlib.fill() function.
826 """
827 ax.fill(self.hull.points[self.hull_points, 0],
828 self.hull.points[self.hull_points, 1], lw=0, **kwargs)
830 def plot_hull_center(self, ax, **kwargs):
831 """Plot the center of mass of the convex hull of the input points.
833 Parameters
834 ----------
835 ax: matplotlib.Axes
836 The axes to be used for plotting.
837 **kwargs: dict
838 Key-word arguments that are passed on to the matplotlib.plot() function.
839 """
840 ax.plot(self.hull_center[0], self.hull_center[1], 'o', **kwargs)
842 def plot_outer_hull(self, ax, **kwargs):
843 """Plot the hull line containing the input points and the vertices of the Voronoi diagram.
845 Parameters
846 ----------
847 ax: matplotlib.Axes or None
848 The axes to be used for plotting.
849 **kwargs: dict
850 Key-word arguments that are passed on to the matplotlib.plot() function.
851 """
852 ax.plot(self.outer_hull_points[:, 0],
853 self.outer_hull_points[:, 1], **kwargs)
855 def fill_outer_hull(self, ax, **kwargs):
856 """Fill the hull containing the input points and the vertices of the Voronoi diagram.
858 Parameters
859 ----------
860 ax: matplotlib.Axes
861 The axes to be used for plotting.
862 **kwargs: dict
863 Key-word arguments that are passed on to the matplotlib.fill() function.
864 """
865 ax.fill(self.outer_hull_points[:, 0],
866 self.outer_hull_points[:, 1], lw=0, **kwargs)
869def main():
870 import scipy.stats as st
871 import matplotlib.pyplot as plt
873 # generate random points:
874 rs = np.random.randint(0xffffffff)
875 #rs = 1718315706 # n=10
876 #rs = 803645102 # n=10
877 #rs = 2751318392 # double infinite ridges at vertex 0 (n=10)
878 #rs = 4226093154 out-hull fac < 1.0
880 print(f'random seed: {rs}')
881 np.random.seed(rs)
882 n = 20 # number of points
883 ar = 1.0 # aspect ratio of area in which the points should be placed
884 points = np.random.rand(int(n//ar), 2)
885 points = points[points[:, 1] < ar, :]
887 # calculate Voronoi diagram:
888 vor = Voronoi(points)
890 # what we get is:
891 print(f'dimension: {vor.ndim}')
892 print(f'number of points: {vor.npoints}')
893 print('')
894 print('distances of nearest neighbors:')
895 print(vor.nearest_distances)
896 print('for each point its nearest neighbor:')
897 print(vor.nearest_points)
898 print('for each point all Voronoi distances:')
899 print(vor.neighbor_distances)
900 print('for each point all its neighbors:')
901 print(vor.neighbor_points)
902 print('for each ridge distances of neighbors:')
903 print(vor.ridge_distances)
904 print('corresponding neighbors enclosing ridges:')
905 print(vor.ridge_points)
906 print('')
907 print('length of corresponding ridges:')
908 print(vor.ridge_lengths())
909 print('area of corresponding triangles:')
910 print(vor.ridge_areas())
911 for mode in ['inside', 'finite_inside', 'full', 'finite']:
912 print(f'Voronoi area of each point ({mode}):')
913 print(vor.areas(mode))
914 print('Comparison of sum of Voronoi areas (inside < finite_inside < full < finit):')
915 a = [f'{np.sum(vor.areas(mode)):.2f}' for mode in ['inside', 'finite_inside', 'full', 'finite']]
916 print(' < '.join(a))
917 print('Type of Voronoi area of each point:')
918 print(vor.point_types())
920 # plot Voronoi diagram:
921 fig, ax = plt.subplots()
922 ax.set_title('Voronoi')
923 vor.fill_regions(ax, inside=True, colors=['red', 'green', 'blue', 'orange', 'cyan', 'magenta'], alpha=1.0, zorder=0)
924 vor.fill_regions(ax, inside=False, colors=['red', 'green', 'blue', 'orange', 'cyan', 'magenta'], alpha=0.4, zorder=0)
925 vor.fill_infinite_regions(ax, colors=['red', 'green', 'blue', 'orange', 'cyan', 'magenta'], alpha=0.1, zorder=0)
926 vor.plot_distances(ax, color='red')
927 #vor.plot_ridges(ax, c='k', lw=2)
928 #vor.plot_infinite_ridges(ax, c='k', lw=2, linestyle='dashed')
929 vor.plot_points(ax, text='p%d', c='c', s=100, zorder=10)
930 vor.plot_vertices(ax, text='v%d', c='r', s=60)
931 ex = 0.3
932 ax.set_xlim(vor.min_bound[0]-ex, vor.max_bound[0]+ex)
933 ax.set_ylim(vor.min_bound[1]-ex, vor.max_bound[1]+ex)
934 ax.set_aspect('equal')
936 # Convex hull:
937 print('Convex hull:')
938 print(f'Area of convex hull: {vor.hull_area():g}')
939 print(f'Area of outer hull: {vor.outer_hull_area():g}')
940 new_points = vor.random_points(poisson=True, mode='outer')
942 # plot convex hull:
943 fig, ax = plt.subplots()
944 vor.fill_outer_hull(ax, color='black', alpha=0.1)
945 vor.plot_outer_hull(ax, color='m', lw=2, label='outer hull')
946 vor.fill_hull(ax, color='black', alpha=0.2)
947 vor.plot_hull(ax, color='r', lw=2, label='convex hull')
948 vor.plot_hull_center(ax, color='r', ms=16)
949 vor.plot_hull_center(ax, color='k', ms=4)
950 vor.plot_center(ax, color='c', ms=16)
951 vor.plot_center(ax, color='k', ms=4)
952 vor.plot_points(ax, text='p%d', c='c', s=100, zorder=10,
953 label='input points')
954 ax.plot(new_points[:, 0], new_points[:, 1], 'ob', ms=10,
955 label='random points')
956 ax.set_xlim(vor.min_bound[0]-ex, vor.max_bound[0]+ex)
957 ax.set_ylim(vor.min_bound[1]-ex, vor.max_bound[1]+ex)
958 ax.set_aspect('equal')
959 ax.legend()
961 # bootstrap:
962 def bootstrapped_nearest_distances(vor, n, poisson, mode, ax1, ax2, ax3, bins):
963 ps = 'fixed'
964 if poisson:
965 ps = 'poisson'
966 print(f'bootstrap {mode} {ps} ...')
967 hist_distances = np.zeros((n, len(bins)-1))
968 mean_distances = []
969 cvs = []
970 for j in range(n):
971 points = vor.random_points(poisson=poisson, mode=mode)
972 if len(points) < 4:
973 continue
974 try:
975 bvor = Voronoi(points)
976 except:
977 continue
978 hist_distances[j,:] = np.histogram(bvor.nearest_distances,
979 bins=bins, density=True)[0]
980 mean_distances.append(np.mean(bvor.nearest_distances))
981 cvs.append(np.std(bvor.nearest_distances)/np.mean(bvor.nearest_distances))
982 ax1.set_title(f'bootstrap {mode} {ps} ...')
983 ax1.hist(mean_distances, 30, density=True, label='bootstrap')
984 de = np.mean(vor.nearest_distances)
985 sem = 2.1*0.26136*np.sqrt(vor.outer_hull_area())/vor.npoints # note factor 2.1!
986 h = 0.4/sem
987 x = np.linspace(de-4.0*sem, de+4.0*sem, 200)
988 p = st.norm.pdf(x, loc=de, scale=sem)
989 ax1.plot(x, p, 'r', lw=2, label='Gaussian')
990 ax1.plot([de, de], [0.0, h], 'r', lw=4, label='observed a.n.n.')
991 bmd = np.mean(mean_distances)
992 ax1.plot([bmd, bmd], [0.0, h], 'g', lw=4, label='bootstrapped a.n.n.')
993 ax1.set_xlabel('mean nearest-neighbor distance')
994 ax1.legend()
995 ax2.set_title(f'bootstrap {mode} {ps} ...')
996 ax2.hist(cvs, 30, density=True, label='bootstrap')
997 # observed cv:
998 cvo = np.std(vor.nearest_distances)/np.mean(vor.nearest_distances)
999 # significance from bootstrap:
1000 pb = 0.01*st.percentileofscore(cvs, cvo)
1001 if pb > 0.5:
1002 pb = 1.0-pb
1003 pb *= 2.0
1004 # significance from z-score:
1005 secv = 2.0*np.sqrt(1.0/np.pi-0.25)/np.sqrt(vor.npoints)
1006 zcv = (cvo - 2.0*0.26136)/secv
1007 pz = 2.0*(1.0 - st.norm.cdf(np.abs(zcv)))
1008 h = 0.4/secv
1009 ax2.plot([cvo, cvo], [0.0, h], 'r', lw=4, label=rf'observed CV $\alpha=${100.0*pz:.0f}%')
1010 cvb = np.mean(cvs)
1011 x = np.linspace(cvb-4.0*secv, cvb+4.0*secv, 200)
1012 p = st.norm.pdf(x, loc=cvb, scale=secv)
1013 ax2.plot(x, p, 'g', lw=2, label='Gaussian')
1014 ax2.plot([cvb, cvb], [0.0, h], 'g', lw=4, label=rf'bootstrapped CV $\alpha=${100.0*pb:.0f}%')
1015 ax2.set_xlabel('CV')
1016 ax2.legend()
1017 # nearest-neighbor distance histogram:
1018 ax3.set_title(f'bootstrap {mode} {ps} ...')
1019 ax3.fill_between(bins[:-1], *np.percentile(hist_distances, [2.5, 97.5], axis=0), facecolor='blue', alpha=0.5)
1020 ax3.plot(bins[:-1], np.median(hist_distances, axis=0), 'b', lw=2)
1021 hd = np.histogram(vor.nearest_distances, bins=bins, density=True)[0]
1022 ax3.plot(bins[:-1], hd, 'r', lw=2)
1023 ax3.set_xlim(0.0, 0.5)
1024 ax3.set_xlabel('nearest-neighbor distance')
1026 print(f'Mean distance: {np.mean(vor.nearest_distances):g}')
1027 print(f'CV: {np.std(vor.nearest_distances)/np.mean(vor.nearest_distances):g}')
1028 bins = np.linspace(0.0, 1.0, 30)
1029 nb = 300
1030 fig1, axs1 = plt.subplots(2, 2, sharex=True)
1031 fig2, axs2 = plt.subplots(2, 2, sharex=True)
1032 fig3, axs3 = plt.subplots(2, 2, sharex=True)
1033 for k, poisson in enumerate([False, True]):
1034 for j, mode in enumerate(['hull', 'outer']):
1035 bootstrapped_nearest_distances(vor, nb, poisson, mode,
1036 axs1[k,j], axs2[k,j], axs3[k,j],
1037 bins)
1038 for f in [fig1, fig2, fig3]:
1039 f.tight_layout()
1040 print('... done.')
1041 plt.tight_layout()
1042 plt.show()
1045if __name__ == "__main__":
1046 main()