Starting with version 0.18.0, scipy provides the class scipy.spatial.SphericalVoronoi, that defines a spherical Voronoi diagram associated to a set of points on a sphere.
The example posted at the above link suggests to visualize the Voronoi diagram as a mpl_toolkits.mplot3d.art3d.Poly3DCollection. The collection of polygons has only the vertices on the sphere and it represents an approximation of the spherical Voronoi diagram.
In this Jupyter Notebook we project the sides of the polygons onto the sphere, in order to get the true boundaries of the spherical Voronoi regions. The spherical Voronoi diagram is visualized as an interactive Plotly plot.
A second Plotly plot displays a colored polyhedral approximation of a spherical Voronoi diagram.
In [2]:
import scipy
scipy. __version__
Out[2]:
In [3]:
from scipy.spatial import SphericalVoronoi
from scipy.interpolate import splev, splprep
import numpy as np
In [4]:
import plotly.plotly as py
from plotly.graph_objs import *
We are working with the unit sphere S(center=O, radius=1).
In [5]:
def sphere():
theta=np.linspace(0, 2*np.pi, 200)
phi=np.linspace(0, np.pi, 100)
theta, phi=np.meshgrid(theta, phi)
x=np.cos(theta)*np.sin(phi)
y=np.sin(theta)*np.sin(phi)
z=np.cos(phi)
return x, y ,z
Define the SphericalVoronoi object associated to a set of points
In [6]:
def s_Voronoi(points):
if not isinstance(points, np.ndarray):
raise ValueError('points must be a numpy array of shape (n,3)')
center = np.zeros(3)
radius = 1.0
sv = SphericalVoronoi(points, radius, center)
sv.sort_vertices_of_regions()# sort vertices
return sv
sv.vertices is the numpy array of all Voronoi vertices (3d points on the unit sphere).sv.regions is a list of lists. An inner list stores the indices of vertices associated to a region:N points uniformly distributed on the unit sphere are defined as unit vectors associated to N 3D-vectors, of independent coordinates, generated from the standard normal distribution:
In [7]:
def points_sphere(N=80):
points=np.random.randn(N,3)
return points/np.linalg.norm(points, axis=1)[:, None]
In order to get the spherical arcs as boundaries of the Voronoi regions we proceed as follows:
For each pair of consecutive points, $P_k, P_{k+1}$, defining a Voronoi region, compute 5 points on the segment having these points as ends. The five points are projected (via a central projection of center O) onto the sphere of radius $R$, and the corresponding points are spline interpolated to get a spherical arc:
In [8]:
def spheric_arcs(sv, R=1.005):
#sv is an instance of SphericalVoronoi class
#R is a radius slightly greater than 1; the arcs are plotted on the shere of radius R
#to avoid to be covered by the sphere surface
t=np.array([0, 0.25, 0.5, 0.75, 1.0])# five parameters for convex combination of points
p=(1-t)[:, None]
q=t[:, None]
#Xa, Ya, Za are the lists of point coordinates to be plotted as spherical arcs
Xa=[]
Ya=[]
Za=[]
for region in sv.regions:
P=sv.vertices[region]#P is an array whose rows are the vertices of the Voronoi points on the sphere
L=P.shape[0]
for k in range(L):
B=np.array([P[k,:]]*5)
C=np.array([P[(k+1)%L, :]]*5)
A=B*p+C*q#A is an array of 5 points on the segment of ends P[k,:], P[(k+1)%L, :]
A=R*A/np.linalg.norm(A, axis=1)[:, None]#central projection of the points in A onto the sphere
tck,u=splprep([A[:,0],A[:,1],A[:,2]],s=0)
xi,yi, zi= splev(np.linspace(0,1,20),tck)#spline interpolation of the five points on sphere
Xa+=xi.tolist()
Ya+=yi.tolist()
Za+=zi.tolist()
Xa+=[None] #after processing a region insert None in each list to avoid
Ya+=[None] #unwanted lines from one region to another
Za+=[None]
return Xa, Ya, Za
The planar approximations of the spherical Voronoi regions are polygons in the 3d space. The union of their sides is returned by the function:
In [9]:
def polygons(sv):
Xp=[]
Yp=[]
Zp=[]
for region in sv.regions:
V=sv.vertices[region]
Xp+=V[:,0].tolist()+[V[0,0], None]
Yp+=V[:,1].tolist()+[V[0,1], None]
Zp+=V[:,2].tolist()+[V[0,2], None]
return Xp, Yp, Zp
The collection of all polygons that represent the Voronoi regions defines a polyhedral approximation of the spherical Voronoi diagram. This approximation can be visualized as a Plotly Mesh3d object.
Namely, since each Voronoi region is a convex polygon, we can triangulate it adding diagonals from the vertex 0 to all non-adjacent vertices. The triangles of a region are colored with the same color chosen randomly from a list colors.
The following function triangulates each region and assigns to its triangles (simplices) a color:
In [10]:
def triangulate_regions(sv, colors):
simplices=[]# a list of 3-lists of integers giving the vertex indices of a triangle
facecolor=[]# the list of colors associated to each simplex in the triangulation
for k, region in enumerate(sv.regions):
color=colors[np.random.randint(0, len(colors))]# choose the color for the region region
#triangulate the region
simplices+=[[region[0], region[j], region[j+1]] for j in range(1, len(region)-1)]
facecolor+=[color for _ in range(1, len(region)-1)]
return simplices, facecolor
Depending on the type of plot, spherical arcs on sphere or polyhedral approximation of the Voronoi diagram,
the function get_data defines the Plotly objects involved in each plot (sphere, points, arcs, respectively,
points, line segments, colored Voronoi regions):
In [11]:
def get_data(points, R=1.005, arcs=True, colorscale=[], colors=[]):
#for arcs=True, and colorscale of length at least 2, the function set up data for plotting the sphere, data points and
#spherical arcs
# for arcs=False, and length of colors, non zero -> data for the polyhedral approximation: points, sides of polygons,
# colored planar Voronoi regions
sv=s_Voronoi(points)
if arcs:
if len(colorscale)<2:
raise ValueError('the colorscale must have at least length=2')
x, y, z=sphere()
sphere_surf=Surface(x=x, y=y, z=z, colorscale=colorscale, showscale=False, name='sphere')
data_pts=Scatter3d(x=R*points[:,0], y=R*points[:,1], z=R*points[:,2], name='points',
mode='markers', marker=dict(color='black', size=3))
Xa, Ya, Za=spheric_arcs(sv,R=R)
lines=Scatter3d(x=Xa, y=Ya, z=Za, name='spheric arc', mode='lines',
line=dict(width=2, color='rgb(10,10,10)'))
return Data([sphere_surf, data_pts, lines ])
else:
if len(colors)==0:
raise ValueError('the list of colors is empty')
simplices, facecolor=triangulate_regions(sv, colors=colors)
I, J, K=np.array(simplices).T
x,y,z=sv.vertices.T
triangles=Mesh3d(x=x,
y=y,
z=z,
facecolor=facecolor,
i=I,
j=J,
k=K,
name='',
)
data_pts=Scatter3d(x=points[:,0], y=points[:,1], z=points[:,2], name='points',
mode='markers', marker=dict(color='black', size=2))
Xp, Yp, Zp=polygons(sv)
lines=Scatter3d(x=Xp, y=Yp, z=Zp, name='spheric arc', mode='lines',
line=dict(width=2, color='rgb(10,10,10)'))
return Data([triangles, data_pts, lines])
Set the plot layout (with axes or not):
In [12]:
axis = dict(
showbackground=True,
backgroundcolor='rgb(40,40,40)', #"rgb(230, 230,230)",
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
)
noaxis=dict(showbackground=False,
showgrid=False,
showline=False,
showticklabels=False,
ticks='',
title='',
zeroline=False)
In [13]:
def plot_layout(ax=noaxis):
return Layout(title='Spherical Voronoi Diagram',
font=dict(family='Balto', size=16),
width=700,
height=700,
showlegend=False,
scene=Scene(xaxis=XAxis(ax),
yaxis=YAxis(ax),
zaxis=ZAxis(ax),
aspectratio=dict(x=1,
y=1,
z=1
),
)
)
In [15]:
points=points_sphere(N=80)
Define a custom Plotly colorscale to plot the sphere:
In [16]:
pl_col=[[0.0, 'rgb(230,230,230)'], [1.0, 'rgb(230,230,230)']]
In [17]:
data1=get_data(points, R=1.005, arcs=True, colorscale=pl_col, colors=[])
In [18]:
fig1 = Figure(data=data1, layout=plot_layout(ax=axis))
py.sign_in('empet', 'jhog40sb94')
py.iplot(fig1, filename='sph-voronoi-axes')
Out[18]:
List of colors for planar Voronoi regions:
In [19]:
colors=['rgb(53,195,176)',
'rgb(168,201,121)',
'rgb(255,210,181)',
'rgb(255,169,164)',
'rgb(255,140,148)']
In [20]:
points=points_sphere(N=100)
data2=get_data(points, R=1.005, arcs=False, colorscale=[], colors=colors)
In [21]:
fig2 = Figure(data=data2, layout=plot_layout(ax=noaxis))
fig2['layout'].update(title='Polyhedral approximation of a spherical Voronoi diagram')
py.iplot(fig2, filename='polyhedral-voronoi')
Out[21]:
In [22]:
from IPython.core.display import HTML
def css_styling():
styles = open("./custom.css", "r").read()
return HTML(styles)
css_styling()
Out[22]: