Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/spatial/_spherical_voronoi.py: 16%
79 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2Spherical Voronoi Code
4.. versionadded:: 0.18.0
6"""
7#
8# Copyright (C) Tyler Reddy, Ross Hemsley, Edd Edmondson,
9# Nikolai Nowaczyk, Joe Pitt-Francis, 2015.
10#
11# Distributed under the same BSD license as SciPy.
12#
14import numpy as np
15import scipy
16from . import _voronoi
17from scipy.spatial import cKDTree
19__all__ = ['SphericalVoronoi']
22def calculate_solid_angles(R):
23 """Calculates the solid angles of plane triangles. Implements the method of
24 Van Oosterom and Strackee [VanOosterom]_ with some modifications. Assumes
25 that input points have unit norm."""
26 # Original method uses a triple product `R1 . (R2 x R3)` for the numerator.
27 # This is equal to the determinant of the matrix [R1 R2 R3], which can be
28 # computed with better stability.
29 numerator = np.linalg.det(R)
30 denominator = 1 + (np.einsum('ij,ij->i', R[:, 0], R[:, 1]) +
31 np.einsum('ij,ij->i', R[:, 1], R[:, 2]) +
32 np.einsum('ij,ij->i', R[:, 2], R[:, 0]))
33 return np.abs(2 * np.arctan2(numerator, denominator))
36class SphericalVoronoi:
37 """ Voronoi diagrams on the surface of a sphere.
39 .. versionadded:: 0.18.0
41 Parameters
42 ----------
43 points : ndarray of floats, shape (npoints, ndim)
44 Coordinates of points from which to construct a spherical
45 Voronoi diagram.
46 radius : float, optional
47 Radius of the sphere (Default: 1)
48 center : ndarray of floats, shape (ndim,)
49 Center of sphere (Default: origin)
50 threshold : float
51 Threshold for detecting duplicate points and
52 mismatches between points and sphere parameters.
53 (Default: 1e-06)
55 Attributes
56 ----------
57 points : double array of shape (npoints, ndim)
58 the points in `ndim` dimensions to generate the Voronoi diagram from
59 radius : double
60 radius of the sphere
61 center : double array of shape (ndim,)
62 center of the sphere
63 vertices : double array of shape (nvertices, ndim)
64 Voronoi vertices corresponding to points
65 regions : list of list of integers of shape (npoints, _ )
66 the n-th entry is a list consisting of the indices
67 of the vertices belonging to the n-th point in points
69 Methods
70 -------
71 calculate_areas
72 Calculates the areas of the Voronoi regions. For 2D point sets, the
73 regions are circular arcs. The sum of the areas is `2 * pi * radius`.
74 For 3D point sets, the regions are spherical polygons. The sum of the
75 areas is `4 * pi * radius**2`.
77 Raises
78 ------
79 ValueError
80 If there are duplicates in `points`.
81 If the provided `radius` is not consistent with `points`.
83 Notes
84 -----
85 The spherical Voronoi diagram algorithm proceeds as follows. The Convex
86 Hull of the input points (generators) is calculated, and is equivalent to
87 their Delaunay triangulation on the surface of the sphere [Caroli]_.
88 The Convex Hull neighbour information is then used to
89 order the Voronoi region vertices around each generator. The latter
90 approach is substantially less sensitive to floating point issues than
91 angle-based methods of Voronoi region vertex sorting.
93 Empirical assessment of spherical Voronoi algorithm performance suggests
94 quadratic time complexity (loglinear is optimal, but algorithms are more
95 challenging to implement).
97 References
98 ----------
99 .. [Caroli] Caroli et al. Robust and Efficient Delaunay triangulations of
100 points on or close to a sphere. Research Report RR-7004, 2009.
102 .. [VanOosterom] Van Oosterom and Strackee. The solid angle of a plane
103 triangle. IEEE Transactions on Biomedical Engineering,
104 2, 1983, pp 125--126.
106 See Also
107 --------
108 Voronoi : Conventional Voronoi diagrams in N dimensions.
110 Examples
111 --------
112 Do some imports and take some points on a cube:
114 >>> import numpy as np
115 >>> import matplotlib.pyplot as plt
116 >>> from scipy.spatial import SphericalVoronoi, geometric_slerp
117 >>> from mpl_toolkits.mplot3d import proj3d
118 >>> # set input data
119 >>> points = np.array([[0, 0, 1], [0, 0, -1], [1, 0, 0],
120 ... [0, 1, 0], [0, -1, 0], [-1, 0, 0], ])
122 Calculate the spherical Voronoi diagram:
124 >>> radius = 1
125 >>> center = np.array([0, 0, 0])
126 >>> sv = SphericalVoronoi(points, radius, center)
128 Generate plot:
130 >>> # sort vertices (optional, helpful for plotting)
131 >>> sv.sort_vertices_of_regions()
132 >>> t_vals = np.linspace(0, 1, 2000)
133 >>> fig = plt.figure()
134 >>> ax = fig.add_subplot(111, projection='3d')
135 >>> # plot the unit sphere for reference (optional)
136 >>> u = np.linspace(0, 2 * np.pi, 100)
137 >>> v = np.linspace(0, np.pi, 100)
138 >>> x = np.outer(np.cos(u), np.sin(v))
139 >>> y = np.outer(np.sin(u), np.sin(v))
140 >>> z = np.outer(np.ones(np.size(u)), np.cos(v))
141 >>> ax.plot_surface(x, y, z, color='y', alpha=0.1)
142 >>> # plot generator points
143 >>> ax.scatter(points[:, 0], points[:, 1], points[:, 2], c='b')
144 >>> # plot Voronoi vertices
145 >>> ax.scatter(sv.vertices[:, 0], sv.vertices[:, 1], sv.vertices[:, 2],
146 ... c='g')
147 >>> # indicate Voronoi regions (as Euclidean polygons)
148 >>> for region in sv.regions:
149 ... n = len(region)
150 ... for i in range(n):
151 ... start = sv.vertices[region][i]
152 ... end = sv.vertices[region][(i + 1) % n]
153 ... result = geometric_slerp(start, end, t_vals)
154 ... ax.plot(result[..., 0],
155 ... result[..., 1],
156 ... result[..., 2],
157 ... c='k')
158 >>> ax.azim = 10
159 >>> ax.elev = 40
160 >>> _ = ax.set_xticks([])
161 >>> _ = ax.set_yticks([])
162 >>> _ = ax.set_zticks([])
163 >>> fig.set_size_inches(4, 4)
164 >>> plt.show()
166 """
167 def __init__(self, points, radius=1, center=None, threshold=1e-06):
169 if radius is None:
170 raise ValueError('`radius` is `None`. '
171 'Please provide a floating point number '
172 '(i.e. `radius=1`).')
174 self.radius = float(radius)
175 self.points = np.array(points).astype(np.double)
176 self._dim = self.points.shape[1]
177 if center is None:
178 self.center = np.zeros(self._dim)
179 else:
180 self.center = np.array(center, dtype=float)
182 # test degenerate input
183 self._rank = np.linalg.matrix_rank(self.points - self.points[0],
184 tol=threshold * self.radius)
185 if self._rank < self._dim:
186 raise ValueError("Rank of input points must be at least {0}".format(self._dim))
188 if cKDTree(self.points).query_pairs(threshold * self.radius):
189 raise ValueError("Duplicate generators present.")
191 radii = np.linalg.norm(self.points - self.center, axis=1)
192 max_discrepancy = np.abs(radii - self.radius).max()
193 if max_discrepancy >= threshold * self.radius:
194 raise ValueError("Radius inconsistent with generators.")
196 self._calc_vertices_regions()
198 def _calc_vertices_regions(self):
199 """
200 Calculates the Voronoi vertices and regions of the generators stored
201 in self.points. The vertices will be stored in self.vertices and the
202 regions in self.regions.
204 This algorithm was discussed at PyData London 2015 by
205 Tyler Reddy, Ross Hemsley and Nikolai Nowaczyk
206 """
207 # get Convex Hull
208 conv = scipy.spatial.ConvexHull(self.points)
209 # get circumcenters of Convex Hull triangles from facet equations
210 # for 3D input circumcenters will have shape: (2N-4, 3)
211 self.vertices = self.radius * conv.equations[:, :-1] + self.center
212 self._simplices = conv.simplices
213 # calculate regions from triangulation
214 # for 3D input simplex_indices will have shape: (2N-4,)
215 simplex_indices = np.arange(len(self._simplices))
216 # for 3D input tri_indices will have shape: (6N-12,)
217 tri_indices = np.column_stack([simplex_indices] * self._dim).ravel()
218 # for 3D input point_indices will have shape: (6N-12,)
219 point_indices = self._simplices.ravel()
220 # for 3D input indices will have shape: (6N-12,)
221 indices = np.argsort(point_indices, kind='mergesort')
222 # for 3D input flattened_groups will have shape: (6N-12,)
223 flattened_groups = tri_indices[indices].astype(np.intp)
224 # intervals will have shape: (N+1,)
225 intervals = np.cumsum(np.bincount(point_indices + 1))
226 # split flattened groups to get nested list of unsorted regions
227 groups = [list(flattened_groups[intervals[i]:intervals[i + 1]])
228 for i in range(len(intervals) - 1)]
229 self.regions = groups
231 def sort_vertices_of_regions(self):
232 """Sort indices of the vertices to be (counter-)clockwise ordered.
234 Raises
235 ------
236 TypeError
237 If the points are not three-dimensional.
239 Notes
240 -----
241 For each region in regions, it sorts the indices of the Voronoi
242 vertices such that the resulting points are in a clockwise or
243 counterclockwise order around the generator point.
245 This is done as follows: Recall that the n-th region in regions
246 surrounds the n-th generator in points and that the k-th
247 Voronoi vertex in vertices is the circumcenter of the k-th triangle
248 in self._simplices. For each region n, we choose the first triangle
249 (=Voronoi vertex) in self._simplices and a vertex of that triangle
250 not equal to the center n. These determine a unique neighbor of that
251 triangle, which is then chosen as the second triangle. The second
252 triangle will have a unique vertex not equal to the current vertex or
253 the center. This determines a unique neighbor of the second triangle,
254 which is then chosen as the third triangle and so forth. We proceed
255 through all the triangles (=Voronoi vertices) belonging to the
256 generator in points and obtain a sorted version of the vertices
257 of its surrounding region.
258 """
259 if self._dim != 3:
260 raise TypeError("Only supported for three-dimensional point sets")
261 _voronoi.sort_vertices_of_regions(self._simplices, self.regions)
263 def _calculate_areas_3d(self):
264 self.sort_vertices_of_regions()
265 sizes = [len(region) for region in self.regions]
266 csizes = np.cumsum(sizes)
267 num_regions = csizes[-1]
269 # We create a set of triangles consisting of one point and two Voronoi
270 # vertices. The vertices of each triangle are adjacent in the sorted
271 # regions list.
272 point_indices = [i for i, size in enumerate(sizes)
273 for j in range(size)]
275 nbrs1 = np.array([r for region in self.regions for r in region])
277 # The calculation of nbrs2 is a vectorized version of:
278 # np.array([r for region in self.regions for r in np.roll(region, 1)])
279 nbrs2 = np.roll(nbrs1, 1)
280 indices = np.roll(csizes, 1)
281 indices[0] = 0
282 nbrs2[indices] = nbrs1[csizes - 1]
284 # Normalize points and vertices.
285 pnormalized = (self.points - self.center) / self.radius
286 vnormalized = (self.vertices - self.center) / self.radius
288 # Create the complete set of triangles and calculate their solid angles
289 triangles = np.hstack([pnormalized[point_indices],
290 vnormalized[nbrs1],
291 vnormalized[nbrs2]
292 ]).reshape((num_regions, 3, 3))
293 triangle_solid_angles = calculate_solid_angles(triangles)
295 # Sum the solid angles of the triangles in each region
296 solid_angles = np.cumsum(triangle_solid_angles)[csizes - 1]
297 solid_angles[1:] -= solid_angles[:-1]
299 # Get polygon areas using A = omega * r**2
300 return solid_angles * self.radius**2
302 def _calculate_areas_2d(self):
303 # Find start and end points of arcs
304 arcs = self.points[self._simplices] - self.center
306 # Calculate the angle subtended by arcs
307 cosine = np.einsum('ij,ij->i', arcs[:, 0], arcs[:, 1])
308 sine = np.abs(np.linalg.det(arcs))
309 theta = np.arctan2(sine, cosine)
311 # Get areas using A = r * theta
312 areas = self.radius * theta
314 # Correct arcs which go the wrong way (single-hemisphere inputs)
315 signs = np.sign(np.einsum('ij,ij->i', arcs[:, 0],
316 self.vertices - self.center))
317 indices = np.where(signs < 0)
318 areas[indices] = 2 * np.pi * self.radius - areas[indices]
319 return areas
321 def calculate_areas(self):
322 """Calculates the areas of the Voronoi regions.
324 For 2D point sets, the regions are circular arcs. The sum of the areas
325 is `2 * pi * radius`.
327 For 3D point sets, the regions are spherical polygons. The sum of the
328 areas is `4 * pi * radius**2`.
330 .. versionadded:: 1.5.0
332 Returns
333 -------
334 areas : double array of shape (npoints,)
335 The areas of the Voronoi regions.
336 """
337 if self._dim == 2:
338 return self._calculate_areas_2d()
339 elif self._dim == 3:
340 return self._calculate_areas_3d()
341 else:
342 raise TypeError("Only supported for 2D and 3D point sets")