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

1"""Analyze Voronoi diagrams based on scipy.spatial. 

2 

3## Classes 

4 

5- `class Voronoi`: Compute and analyse Voronoi diagrams. 

6""" 

7 

8 

9import numpy as np 

10import scipy.spatial as sp 

11try: 

12 import matplotlib.pyplot as plt 

13except ImportError: 

14 pass 

15 

16 

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: 

28  

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 

38  

39 Default is: "Qbb Qc Qz Qx" for ndim > 4 and "Qbb Qc Qz" otherwise. 

40  

41 Input points 

42 ------------ 

43 The points from which the Voronoi diagram is constructed. 

44  

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. 

57  

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. 

78  

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`). 

98  

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. 

109  

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. 

118  

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. 

127  

128 Bootstrap Voronoi diagrams 

129 -------------------------- 

130 - `random_points()`: Generate random points within the area defined by the input points. 

131 

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. 

142  

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. 

150  

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 ``` 

158  

159 Calculate the Voronoi diagram: 

160 ``` 

161 from thunderlab.voronoi import Voronoi 

162 vor = Voronoi(points) 

163 ``` 

164  

165 Retrieve nearest-neighbor distances and compute Voronoi areas: 

166 ``` 

167 dists = vor.nearest_distances 

168 areas = vor.areas() 

169 ``` 

170  

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 ``` 

176  

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 """ 

185 

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) 

206 

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) 

214 

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] 

227 

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) 

236 

237 def _compute_infinite_vertices(self, radius=None): 

238 """Compute far points of infinite ridges. 

239 

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. 

245 

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) 

308 

309 def _flatten_simplices(self, simplices): 

310 """Transforms list of simplex indices to list of vertex indices. 

311 

312 In particular, transforms the Delaunay.convex_hull to a list of points 

313 of the hull that can then be easily plotted. 

314 

315 Parameters 

316 ---------- 

317 simplices: 2-D array of ints 

318 List of pairs of indices of points forming each ridge of a polygon. 

319 

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 

339 

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) 

346 

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) 

363 

364 def in_hull(self, p): 

365 """Test if points `p` are within convex hull of the input points. 

366 

367 Parameters 

368 ---------- 

369 p: 2-d ndarray 

370 Array of points to be tested. 

371 

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 

379 

380 def in_outer_hull(self, p): 

381 """Test if points `p` are within the outer hull. 

382 

383 Parameters 

384 ---------- 

385 p: 2-d ndarray 

386 Array of points to be tested. 

387 

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 

395 

396 def point_types(self): 

397 """The type of the Voronoi regions for each input point. 

398 

399 Returns 

400 ------- 

401 points: ndarray of ints 

402 Indicates the type of Voronoi region associated 

403 with each point in `points`: 

404 

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 

420 

421 def ridge_lengths(self): 

422 """Length of Voronoi ridges between nearest neighbors. 

423 

424 May be used, for example, as a weigth for `ridge_distances`. 

425 

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 

441 

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. 

445 

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 

458 

459 def areas(self, mode='finite'): 

460 """The areas of the Voronoi regions for each input point. 

461 

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. 

476 

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 

532 

533 def hull_area(self): 

534 """The area of the convex hull of the input points. 

535  

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 

547 

548 def outer_hull_area(self): 

549 """The area of the outer hull. 

550  

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 

562 

563 def random_points(self, n=None, poisson=True, mode='outer'): 

564 """Generate random points. 

565 

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. 

582 

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] 

631 

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. 

635 

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) 

657 

658 def plot_center(self, ax, **kwargs): 

659 """Plot the center of mass of the input points. 

660 

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) 

669 

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. 

673 

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) 

695 

696 def plot_distances(self, ax, **kwargs): 

697 """Plot lines connecting the nearest neighbors in the Voronoi diagram. 

698 

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) 

708 

709 def plot_ridges(self, ax, **kwargs): 

710 """Plot the finite ridges of the Voronoi diagram. 

711 

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) 

722 

723 def plot_infinite_ridges(self, ax, **kwargs): 

724 """Plot the infinite ridges of the Voronoi diagram. 

725 

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) 

739 

740 def fill_regions(self, ax, inside=None, colors=None, **kwargs): 

741 """Fill each finite region of the Voronoi diagram with a color. 

742 

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) 

769 

770 def fill_infinite_regions(self, ax, colors=None, **kwargs): 

771 """Fill each infinite region of the Voronoi diagram with a color. 

772 

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) 

798 

799 def plot_hull(self, ax, **kwargs): 

800 """Plot the hull line containing the input points. 

801 

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) 

811 

812 def fill_hull(self, ax, **kwargs): 

813 """Fill the hull containing the input points with a color. 

814 

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) 

824 

825 def plot_hull_center(self, ax, **kwargs): 

826 """Plot the center of mass of the convex hull of the input points. 

827 

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) 

836 

837 def plot_outer_hull(self, ax, **kwargs): 

838 """Plot the hull line containing the input points and the vertices of the Voronoi diagram. 

839 

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) 

849 

850 def fill_outer_hull(self, ax, **kwargs): 

851 """Fill the hull containing the input points and the vertices of the Voronoi diagram. 

852 

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) 

862 

863 

864def main(): 

865 import scipy.stats as st 

866 import matplotlib.pyplot as plt 

867 

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 

874 

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, :] 

881 

882 # calculate Voronoi diagram: 

883 vor = Voronoi(points) 

884 

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()) 

914 

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') 

930 

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') 

936 

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() 

955 

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') 

1020 

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() 

1038 

1039 

1040if __name__ == "__main__": 

1041 main()