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