In [1]:
%load_ext autoreload
%autoreload 2
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
import numpy
import scipy
import voronoi_utility
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import pandas
import math

In [2]:
#for testing, produce a mesh of points on a sphere:
circumdiameter = 4.0
u, v = numpy.mgrid[0.01:2*numpy.pi-0.01:20j, 0.01:numpy.pi-0.01:10j]
x=circumdiameter/2.0 * (numpy.cos(u)*numpy.sin(v))
y=circumdiameter/2.0 * (numpy.sin(u)*numpy.sin(v))
z=circumdiameter/2.0 * (numpy.cos(v))
#and plot the points to confirm spherical shape by visual inspection:
fig = plt.figure()
ax1 = fig.add_subplot(111, projection='3d')
#ax1.plot_wireframe(x, y, z, color='r',alpha=0.5)
ax1.scatter(x.ravel(),y.ravel(),z.ravel(),c='k')
fig.set_size_inches(9.0,9.0)



In [3]:
#now, let's see if the Delaunay Triangulation from my code looks sensible or even works at all:
input_sphere_coordinate_array = numpy.zeros((200,3))
input_sphere_coordinate_array[...,0] = x.ravel()
input_sphere_coordinate_array[...,1] = y.ravel()
input_sphere_coordinate_array[...,2] = z.ravel()
voronoi_instance = voronoi_utility.Voronoi_Sphere_Surface(input_sphere_coordinate_array)

In [4]:
Delaunay_point_array = voronoi_instance.delaunay_triangulation_spherical_surface()

In [5]:
Delaunay_point_array.shape,x.ravel().shape


Out[5]:
((396, 3, 3), (200,))

In [6]:
#try to plot the triangles of the Delaunay tesselation to see if they look reasonable by visual inspection relative to the data points:
fig2 = plt.figure()
fig2.set_size_inches(14,14)
ax = fig2.add_subplot(111, projection='3d')
for triangle_coordinate_array in Delaunay_point_array:
    #plot the Delaunay edges in red--using separate calls to close the triangles
    ax.plot(triangle_coordinate_array[...,0],triangle_coordinate_array[...,1],triangle_coordinate_array[...,2],c='r',alpha=0.6)
    connecting_array = numpy.delete(triangle_coordinate_array,1,0)
    ax.plot(connecting_array[...,0],connecting_array[...,1],connecting_array[...,2],c='r',alpha=0.6)
    ax.scatter(x.ravel(),y.ravel(),z.ravel(),c='k') #the original data points in black


The above Delaunay triangulation looks reasonable by visual inspection, although specific observations are rather confounded by the 3D perspective (triangle stacking).


In [7]:
#testing 
array_facet_normals = voronoi_instance.Voronoi_vertices_spherical_surface()


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-7-ec3f03ae9fad> in <module>()
      1 #testing
----> 2 array_facet_normals = voronoi_instance.Voronoi_vertices_spherical_surface()

AttributeError: Voronoi_Sphere_Surface instance has no attribute 'Voronoi_vertices_spherical_surface'

In [8]:
array_facet_normals


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-bc40c1e16c89> in <module>()
----> 1 array_facet_normals

NameError: name 'array_facet_normals' is not defined

In [9]:
generator_Voronoi_region_dictionary,dictionary_sorted_Voronoi_point_coordinates_for_each_generator = voronoi_instance.Voronoi_polygons_spherical_surface()


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-9-fb57b031b8dc> in <module>()
----> 1 generator_Voronoi_region_dictionary,dictionary_sorted_Voronoi_point_coordinates_for_each_generator = voronoi_instance.Voronoi_polygons_spherical_surface()

AttributeError: Voronoi_Sphere_Surface instance has no attribute 'Voronoi_polygons_spherical_surface'

In [10]:
dictionary_sorted_Voronoi_point_coordinates_for_each_generator


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-ca968110e971> in <module>()
----> 1 dictionary_sorted_Voronoi_point_coordinates_for_each_generator

NameError: name 'dictionary_sorted_Voronoi_point_coordinates_for_each_generator' is not defined

In [11]:
#plot the results
import matplotlib.colors as colors
fig3 = plt.figure()
fig3.set_size_inches(14,14)
ax = fig3.add_subplot(111, projection='3d')
#colour each Voronoi region and generator point randomly to increase confidence in the results (i.e., in the connection between generators and their Voronoi regions assigned by my code)
for generator_index,voronoi_polygon_vertices in dictionary_sorted_Voronoi_point_coordinates_for_each_generator.iteritems():    
    random_color = colors.rgb2hex(scipy.rand(3))
    generator_coordinate = input_sphere_coordinate_array[generator_index]
    ax.scatter(generator_coordinate[...,0],generator_coordinate[...,1],generator_coordinate[...,2],facecolor=random_color,lw=0,s=20,)
    ax.plot(voronoi_polygon_vertices[...,0],voronoi_polygon_vertices[...,1],voronoi_polygon_vertices[...,2],c='black')
    #want to close the polygons as well
    connecting_vertex_array = numpy.vstack((voronoi_polygon_vertices[0,...],voronoi_polygon_vertices[-1,...]))
    ax.plot(connecting_vertex_array[...,0],connecting_vertex_array[...,1],connecting_vertex_array[...,2],c='black')
    polygon = Poly3DCollection([voronoi_polygon_vertices],alpha=0.6)
    polygon.set_color(random_color)
    ax.add_collection3d(polygon)
ax.set_ylim(-2,2)     
ax.set_xlim(-2,2)     
ax.set_zlim(-2,2)     
#ax.scatter(array_facet_normals[...,0],array_facet_normals[...,1],array_facet_normals[...,2],c='blue') #should be Voronoi vertices when working properly and multiplied by correct factor


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-11-33d85e8a15ce> in <module>()
      5 ax = fig3.add_subplot(111, projection='3d')
      6 #colour each Voronoi region and generator point randomly to increase confidence in the results (i.e., in the connection between generators and their Voronoi regions assigned by my code)
----> 7 for generator_index,voronoi_polygon_vertices in dictionary_sorted_Voronoi_point_coordinates_for_each_generator.iteritems():
      8     random_color = colors.rgb2hex(scipy.rand(3))
      9     generator_coordinate = input_sphere_coordinate_array[generator_index]

NameError: name 'dictionary_sorted_Voronoi_point_coordinates_for_each_generator' is not defined

So, the above result for a regular point grid on a sphere actually looks reasonable now. It would be useful to perform a similar visual test using a random distribution of points on a sphere.


In [12]:
#generate a random distribution of points on the unit sphere (http://mathworld.wolfram.com/SpherePointPicking.html)
import numpy.random
#go for 1000 random points
u = numpy.random.random((1000,)) #200 points on interval [0,1); ideally want (0,1), but perhaps close enough?
v = numpy.random.random((1000,))
theta_array = 2 * math.pi * u
phi_array = numpy.arccos((2*v - 1.0))
r_array = numpy.ones((1000,))
spherical_polar_coord_array = numpy.column_stack((r_array,theta_array,phi_array))
#convert to Cartesian coordinates
cartesian_coord_array = voronoi_utility.convert_spherical_array_to_cartesian_array(spherical_polar_coord_array)

In [13]:
#test plot the random points on unit sphere
fig4 = plt.figure()
fig4.set_size_inches(14,14)
ax = fig4.add_subplot(111, projection='3d')
ax.scatter(cartesian_coord_array[...,0],cartesian_coord_array[...,1],cartesian_coord_array[...,2])
ax.set_ylim(-1,1)     
ax.set_xlim(-1,1)     
ax.set_zlim(-1,1)


Out[13]:
(-1, 1)

In [14]:
#Not entirely convinced I got a true random distribution there (still have pole bias?!), but should be random enough to test with my code
random_dist_voronoi_instance = voronoi_utility.Voronoi_Sphere_Surface(cartesian_coord_array)

In [15]:
random_dictionary_sorted_Voronoi_point_coordinates_for_each_generator = random_dist_voronoi_instance.voronoi_region_vertices_spherical_surface()


---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-15-c9f3b44d7a5d> in <module>()
----> 1 random_dictionary_sorted_Voronoi_point_coordinates_for_each_generator = random_dist_voronoi_instance.voronoi_region_vertices_spherical_surface()

/Users/edd/voronoi_edd/py_sphere_Voronoi/voronoi_utility.py in voronoi_region_vertices_spherical_surface(self)
    452             #print 'Voronoi_point_index:', Voronoi_point_index
    453             #print '5 closest distances:', numpy.sort(Voronoi_point_distance_array)[0:5]
--> 454             assert indices_of_generators_for_which_this_Voronoi_point_is_a_polygon_vertex.size >= 3, "By definition, a Voronoi vertex must be equidistant to at least 3 generators, but in this case only got {num_gens} generators for Voronoi vertex at {coords}, which has 5 closest distances: {distances}.".format(num_gens=indices_of_generators_for_which_this_Voronoi_point_is_a_polygon_vertex.size,coords=array_Voronoi_vertices[Voronoi_point_index],distances=numpy.sort(Voronoi_point_distance_array)[0:5])
    455             generator_Voronoi_region_dictionary[Voronoi_point_index] = indices_of_generators_for_which_this_Voronoi_point_is_a_polygon_vertex #so dictionary looks like 0: array(12,17,27), ...
    456 

AssertionError: By definition, a Voronoi vertex must be equidistant to at least 3 generators, but in this case only got 1 generators for Voronoi vertex at [ 0.95845165  0.2507226  -0.13604634], which has 5 closest distances: [ 0.063  0.103  0.148  0.149  0.155].

In [20]:
#plot the results
import matplotlib.colors as colors
fig5 = plt.figure()
fig5.set_size_inches(14,14)
ax = fig5.add_subplot(111, projection='3d')
for generator_index,voronoi_polygon_vertices in random_dictionary_sorted_Voronoi_point_coordinates_for_each_generator.iteritems():    
    random_color = colors.rgb2hex(scipy.rand(3))
    generator_coordinate = cartesian_coord_array[generator_index]
    #ax.scatter(generator_coordinate[...,0],generator_coordinate[...,1],generator_coordinate[...,2],facecolor=random_color,lw=0,s=50,)
    ax.plot(voronoi_polygon_vertices[...,0],voronoi_polygon_vertices[...,1],voronoi_polygon_vertices[...,2],c='black')
    #want to close the polygons as well
    connecting_vertex_array = numpy.vstack((voronoi_polygon_vertices[0,...],voronoi_polygon_vertices[-1,...]))
    ax.plot(connecting_vertex_array[...,0],connecting_vertex_array[...,1],connecting_vertex_array[...,2],c='black')
    polygon = Poly3DCollection([voronoi_polygon_vertices],alpha=1.0)
    polygon.set_color(random_color)
    ax.add_collection3d(polygon)
ax.set_ylim(-1,1)     
ax.set_xlim(-1,1)     
ax.set_zlim(-1,1)


Out[20]:
(-1, 1)

I played with the presence / absence of the generators and the fill alpha values in the above plot. Considering the risk of ultra-close generators with randomly generated data, I'd say these results look excellent for a less regular data set.


In [202]:
#some rough work to generate image files for upcoming Dengue Club presentation slide:
fig3.set_size_inches(4,4)
fig5.set_size_inches(4,4)
ax3 = fig3.gca()
ax3.set_xticks([-2,0,2])
ax3.set_yticks([-2,0,2])
ax3.set_zticks([-2,0,2])
ax5 = fig5.gca()
ax5.set_xticks([-1,0,1])
ax5.set_yticks([-1,0,1])
ax5.set_zticks([-1,0,1])
fig3.savefig('regular_points_test.png',dpi=300)
fig5.savefig('random_points_test.png',dpi=300)

In [56]:
#I now want to attempt to reconstitute (sum up) the surface areas of the Voronoi polygons to see if this matches the actual surface area of the sphere
#for the surface area of an arbitrary spherical polygon: http://mathworld.wolfram.com/SphericalPolygon.html
unit_circle_great_arc_array = numpy.array([[-1,0,0],[0,1,0]])    
unit_circle_great_arc_array_deriv_slope_one = numpy.array([[-1,0,0],[-math.sqrt(0.5),math.sqrt(0.5),0]])

In [59]:
derivative_vector_3D = voronoi_utility.calculate_derivative_great_circle_arc_specified_point(unit_circle_great_arc_array_deriv_slope_one,1.0)
derivative_vector_3D


Out[59]:
array([[ -7.07350915e-01,   7.06862563e-01,   6.12323400e-17],
       [ -7.07106781e-01,   7.07106781e-01,   0.00000000e+00]])

In [60]:
numpy.diff(derivative_vector_3D,axis=0)


Out[60]:
array([[  2.44133372e-04,   2.44217690e-04,  -6.12323400e-17]])

In [67]:
#problematic_polygon_array = numpy.array([[-0.70041688,  0.07294857, -0.70993082],
#[-0.78297057 ,-0.03756652 ,-0.62084854],
 #[-0.77464149 ,-0.07389853 ,-0.62799413],
 #[-0.65546009 ,-0.19355067 ,-0.72994334],
 #[-0.62186187 ,-0.14154995 ,-0.77016783]])
problematic_polygon_array = numpy.array([[-0.12278101, 0.38828208, 0.90397089],
 [-0.18533492 ,0.28384049, 0.9317119 ],
 [ 0.07210294 ,0.29806975, 0.94284522],
 [ 0.1316095  ,0.32464041, 0.92751769]])
 #[ 0.08298007 ,0.36485282 ,0.91814936]])

In [120]:
fig6 = plt.figure()
fig6.set_size_inches(19,19)
ax = fig6.add_subplot(111, projection='3d')
ax.scatter(problematic_polygon_array[...,0],problematic_polygon_array[...,1],problematic_polygon_array[...,2],color='None',s=60)
ax.plot(problematic_polygon_array[...,0],problematic_polygon_array[...,1],problematic_polygon_array[...,2],color='black')
ax.scatter(cartesian_coord_array[...,0],cartesian_coord_array[...,1],cartesian_coord_array[...,2],alpha=0.1) #original generators
ax.plot(final_subtriangle[...,0],final_subtriangle[...,1],final_subtriangle[...,2],c='green') #the final subtriangle (which I think produces an unreasonably small angle??)
#colour vertices progressively so I'm certain where they start / end:
color_list = ['black','brown','blue','yellow']
index = 0
for vertex in problematic_polygon_array:
    color = color_list[index]
    ax.scatter(vertex[...,0],vertex[...,1],vertex[...,2],color=color,s=80)
    index += 1
ax.elev = 0.
#testing lat / lon conventions for comparison with JPL document
cartesian_coord = numpy.array([[0,0,1]])
ax.scatter(cartesian_coord[...,0],cartesian_coord[...,1],cartesian_coord[...,2],c='green',s=100)
spherical_test_coord = voronoi_utility.convert_cartesian_array_to_spherical_array(cartesian_coord)
spherical_theta_modified = numpy.zeros((1,3))
spherical_theta_modified[0,0] = spherical_test_coord[0,0]
spherical_theta_modified[0,1] = spherical_test_coord[0,1] - (math.pi / 2.) 
spherical_theta_modified[0,2] = spherical_test_coord[0,2] 
theta_adjusted_cartesian_test = voronoi_utility.convert_spherical_array_to_cartesian_array(spherical_theta_modified)
ax.scatter(theta_adjusted_cartesian_test[...,0],theta_adjusted_cartesian_test[...,1],theta_adjusted_cartesian_test[...,2],s=300,c='purple')
spherical_theta_modified


Out[120]:
array([[ 1.        , -1.57079633,  0.        ]])