Interactively test the spherical Voronoi algorithm implementation proposed by Ross Hemsley at PyData London 2015
In [62]:
%load_ext autoreload
%autoreload 2
import matplotlib.pyplot as plt
%matplotlib inline
import voronoi_utility
import numpy as np
import scipy as sp
from mpl_toolkits.mplot3d import Axes3D
In [63]:
#pre-requisite step: generate the random coordinates on the sphere
prng = np.random.RandomState(117)
random_coordinate_array = voronoi_utility.generate_random_array_spherical_generators(1000,1.0,prng)
In [64]:
#confirm reasonable distribution of points on the unit sphere
fig_initial_sphere = plt.figure()
ax = fig_initial_sphere.add_subplot('111', projection = '3d')
ax.scatter(random_coordinate_array[...,0], random_coordinate_array[...,1], random_coordinate_array[...,2], c = 'blue', edgecolor = 'none')
fig_initial_sphere.set_size_inches(6,6)
In [65]:
#step 1: place an additional generator at the centre of the sphere (the origin of the coordinate system)
random_coordinate_array = np.concatenate((random_coordinate_array, np.zeros((1,3))))
random_coordinate_array[-5:,...]
Out[65]:
In [66]:
#confirm appropriate coordinates for the center (in red)
fig_additional_generator = plt.figure()
ax = fig_additional_generator.add_subplot('111', projection = '3d')
ax.scatter(random_coordinate_array[:-1,0], random_coordinate_array[:-1,1], random_coordinate_array[:-1,2], c = 'blue', edgecolor = 'none')
ax.scatter(random_coordinate_array[-1,0], random_coordinate_array[-1,1], random_coordinate_array[-1,2], c = 'red', edgecolor = 'none')
fig_additional_generator.set_size_inches(6,6)
In [67]:
#step 2: perform 3D Delaunay triangulation on data set that includes the extra generator
tri = sp.spatial.Delaunay(random_coordinate_array)
In [68]:
#confirm reasonable-looking tetrahedralization of sphere
fig_initial_tetrahedralization = plt.figure()
ax = fig_initial_tetrahedralization.add_subplot('111', projection = '3d')
for simplex in tri.points[tri.simplices]:
ax.scatter(simplex[...,0], simplex[...,1], simplex[...,2], c = 'k', edgecolor = 'none')
ax.scatter(random_coordinate_array[...,0], random_coordinate_array[...,1], random_coordinate_array[...,2], c = 'blue', edgecolor = 'none', s = 70, alpha = 0.25)
fig_initial_tetrahedralization.set_size_inches(20,20)
Although connecting the tetrahedra in a non-intersecting fashion would probably be a bit more work, the superposition of generators on the tetrahedral vertices is sensible given that each generator should indeed be a vertex in the 3D tetrahedralization.
In [116]:
#step 3a: produce circumspheres / circumcenters of tetrahedra from 3D Delaunay
#based on some rough testing, it turns out that ALL simplices (tetrahedra) have one vertex at the origin [I think Ross was trying to explain this to me at PyData London 2015]
#the unit normal to the outer Delaunay facet (on sphere surface) is a Voronoi vertex (as in my previous algorithm)
#because each tetrahedron contains the origin as a vertex, the circumcenter is located halfway between the origin and the Voronoi vertex [othwerise the other vertices would illegally be located inside the circumsphere]
#furthermore, if we exclude the tetrahedral vertex at the origin, the remaining three vertices form a face (triangle) on the sphere surface, the circumcenter of which may be projected to the sphere surface to produce the coordinates of the corresponding Voronoi vertex
import circumcircle
list_circumcenter_coordinates = []
for simplex in tri.points[tri.simplices]: #iterate over tetrahedra
#need to find circumcenter of tetrahedron circumsphere
tetrahedron_circumsphere_circumcenter = circumcircle.calc_circumcenter_circumsphere_tetrahedron_2(simplex)
list_circumcenter_coordinates.append(tetrahedron_circumsphere_circumcenter)
array_circumcenter_coords = np.array(list_circumcenter_coordinates)
In [117]:
array_circumcenter_coords.shape
Out[117]:
In [118]:
array_circumcenter_coords
Out[118]:
In [122]:
#I'm pretty sure the circumcenter coords should fall within the unit sphere because the circumsphere of the simplex should not contain any generators
#do some isolated plotting to test this out
In [123]:
first_tetrahedron = tri.points[tri.simplices][0]
first_tetrahedron
first_circumcenter = circumcircle.calc_circumcenter_circumsphere_tetrahedron_2(first_tetrahedron)
fig_test_circumcenter = plt.figure()
ax = fig_test_circumcenter.add_subplot('111', projection = '3d')
ax.scatter(first_tetrahedron[...,0], first_tetrahedron[...,1], first_tetrahedron[...,2], c = 'k')
ax.scatter(first_circumcenter[...,0], first_circumcenter[...,1], first_circumcenter[...,2], c= 'blue', edgecolor = 'none')
Out[123]:
In [124]:
#ok, the circumcenter of circumsphere calculation is now starting to look more reasonable with the new implementation
In [127]:
#plot the full set of circumcenter coordinates to confirm that they are all looking reasonable relative to the generators
fig_test_all_circumcenters = plt.figure()
ax = fig_test_all_circumcenters.add_subplot('111', projection = '3d')
ax.scatter(array_circumcenter_coords[...,0], array_circumcenter_coords[...,1], array_circumcenter_coords[...,2], c = 'blue', edgecolor = 'none', label = 'tetrahedron circumcenters')
ax.scatter(random_coordinate_array[...,0], random_coordinate_array[...,1], random_coordinate_array[...,2], c = 'black', edgecolor = 'none', label = 'generators')
ax.legend()
fig_test_all_circumcenters.set_size_inches(6,6)
In [137]:
#attempt to project the tetrahedron circumcenters up to the surface of the sphere, to produce the Voronoi vertices
# based on http://stackoverflow.com/a/9604279
radius = 1.0
array_vector_lengths = sp.spatial.distance.cdist(array_circumcenter_coords,np.zeros((1,3)))
array_Voronoi_vertices = (radius / np.abs(array_vector_lengths)) * array_circumcenter_coords
In [138]:
array_Voronoi_vertices.shape
Out[138]:
In [139]:
#confirm reasonable positions for Voronoi vertices relative to generators by visual inspection:
fig_assess_Voronoi_vertices = plt.figure()
ax = fig_assess_Voronoi_vertices.add_subplot('111', projection = '3d')
ax.scatter(array_Voronoi_vertices[...,0], array_Voronoi_vertices[...,1], array_Voronoi_vertices[...,2], c = 'blue', edgecolor = 'none', label = 'Voronoi vertices')
ax.scatter(random_coordinate_array[...,0], random_coordinate_array[...,1], random_coordinate_array[...,2], c = 'black', edgecolor = 'none', label = 'generators')
ax.legend()
fig_assess_Voronoi_vertices.set_size_inches(14,14)