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

In [73]:
#What exactly is wrong with the arrangements of vertices in the above polygon on the unit sphere that prevents a proper surface area calculation?
#One of them does look a bit 'out of plane--' take a closer look
problematic_polygon_spherical_polar_coords = voronoi_utility.convert_cartesian_array_to_spherical_array(problematic_polygon_array)
problematic_polygon_spherical_polar_coords


Out[73]:
array([[ 0.99146433,  0.4232636 ,  1.87706301],
       [ 0.99146433,  0.34894733,  2.14924571],
       [ 0.99146433,  0.31446395,  1.33345602],
       [ 0.99146433,  0.36111677,  1.18564273]])

In [86]:
#the radial distances are the same, so the polygon appears to be completely on the surface
voronoi_utility.calculate_surface_area_of_a_spherical_Voronoi_polygon(problematic_polygon_array,1.0)


subtriangle vertex coords: [ 0.1316095   0.32464041  0.92751769] [-0.12278101  0.38828208  0.90397089] [-0.18533492  0.28384049  0.9317119 ]
a,b,c side lengths on subtriangle: 0.223145325401 0.371087531608 0.322828413589
Euclidean edge lengths (debug): 0.124862328241 0.319587218516 0.263285482943
current vertex inner angle (degrees): 84.191989221
subtriangle vertex coords: [-0.12278101  0.38828208  0.90397089] [-0.18533492  0.28384049  0.9317119 ] [ 0.07210294  0.29806975  0.94284522]
a,b,c side lengths on subtriangle: 0.318516021867 0.286683715775 0.223145325401
Euclidean edge lengths (debug): 0.258071064662 0.218241224301 0.124862328241
current vertex inner angle (degrees): 61.5331597412
subtriangle vertex coords: [-0.18533492  0.28384049  0.9317119 ] [ 0.07210294  0.29806975  0.94284522] [ 0.1316095   0.32464041  0.92751769]
a,b,c side lengths on subtriangle: 0.19647632473 0.371087531608 0.318516021867
Euclidean edge lengths (debug): 0.0669474706899 0.319587218516 0.258071064662
current vertex inner angle (degrees): 89.5275516286
subtriangle vertex coords: [ 0.07210294  0.29806975  0.94284522] [ 0.1316095   0.32464041  0.92751769] [-0.12278101  0.38828208  0.90397089]
a,b,c side lengths on subtriangle: 0.322828413589 0.286683715775 0.19647632473
Euclidean edge lengths (debug): 0.263285482943 0.218241224301 0.0669474706899
current vertex inner angle (degrees): 61.9839286961
theta: 5.18775783859
n: 4
lower theta boundary: 6.28318530718
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-86-4407d9456f13> in <module>()
      1 #the radial distances are the same, so the polygon appears to be completely on the surface
----> 2 voronoi_utility.calculate_surface_area_of_a_spherical_Voronoi_polygon(problematic_polygon_array,1.0)

/Users/treddy/github_projects/py_sphere_Voronoi/voronoi_utility.py in calculate_surface_area_of_a_spherical_Voronoi_polygon(array_ordered_Voronoi_polygon_vertices, sphere_radius)
     18     print 'lower theta boundary:', (n - 2) * math.pi
     19     surface_area_Voronoi_polygon_on_sphere_surface = (theta - ((n - 2) * math.pi)) * (sphere_radius ** 2)
---> 20     assert surface_area_Voronoi_polygon_on_sphere_surface > 0, "Surface areas of spherical polygons should be > 0 but got: {SA}".format(SA=surface_area_Voronoi_polygon_on_sphere_surface)
     21     #print 'surface_area_Voronoi_polygon_on_sphere_surface:', surface_area_Voronoi_polygon_on_sphere_surface
     22     return surface_area_Voronoi_polygon_on_sphere_surface

AssertionError: Surface areas of spherical polygons should be > 0 but got: -1.09542746859

In [81]:
#work with the final subtriangle (i.e., plot above), which seems to produce an unreasonably small angle by visual inspection
final_subtriangle = numpy.array([[ 0.07210294 ,0.29806975 ,0.94284522],[ 0.1316095  ,0.32464041 ,0.92751769],[-0.12278101 ,0.38828208 ,0.90397089]])

In [85]:
#I'm pretty sure b should be longer than a in the final subtriangle!!

In [91]:
math.acos(numpy.dot(numpy.array([0,0,1]),numpy.array([1,0,0])))


Out[91]:
1.5707963267948966

In [ ]:


In [30]:
#rough work to assess possible issues with surface area reconstitution in one of my unit tests
#%load -r 238-247 test_voronoi_utility.py

In [75]:
unit_sphere_surface_area = 4 * math.pi
circumdiameter = 2.0
u, v = numpy.mgrid[0.01:2*numpy.pi:15j, 0.6:numpy.pi-0.6: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))
input_sphere_coordinate_array = numpy.zeros((150,3))
input_sphere_coordinate_array[...,0] = x.ravel()
input_sphere_coordinate_array[...,1] = y.ravel()
input_sphere_coordinate_array[...,2] = z.ravel()

In [76]:
fig7 = plt.figure()
fig7.set_size_inches(12,12)
ax = fig7.add_subplot(111, projection='3d')
ax.scatter(input_sphere_coordinate_array[...,0],input_sphere_coordinate_array[...,1],input_sphere_coordinate_array[...,2],c='red',s=60,linewidth=1.0,alpha=0.9)
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
ax.set_zlim(-1,1)


Out[76]:
(-1, 1)

In [84]:
#working to improve the selection of random points on the unit sphere for my unit testing, etc., to less pathological case (i.e., remove pole bias!!)
test_random_coords_on_sphere = voronoi_utility.generate_random_array_spherical_generators(4000,1.0)

In [85]:
test_random_coords_on_sphere[0:10]


Out[85]:
array([[-0.46232822, -0.86759182,  0.18313126],
       [-0.25819432,  0.48899412, -0.83319892],
       [-0.28662818,  0.90375812, -0.31790808],
       [-0.2363103 ,  0.78704402,  0.56984133],
       [ 0.39804765,  0.83641857, -0.3767785 ],
       [ 0.02374909,  0.31747978,  0.9479676 ],
       [-0.40173819,  0.8753041 ,  0.26916382],
       [ 0.77293988,  0.44538366,  0.451882  ],
       [ 0.61905806,  0.78492307,  0.0257466 ],
       [ 0.46852946, -0.50506803, -0.72483545]])

In [86]:
fig8 = plt.figure()
fig8.set_size_inches(19,19)
ax = fig8.add_subplot(111, projection='3d')
ax.scatter(test_random_coords_on_sphere[...,0],test_random_coords_on_sphere[...,1],test_random_coords_on_sphere[...,2],c='red')


Out[86]:
<mpl_toolkits.mplot3d.art3d.Patch3DCollection at 0x118b51710>

In [91]:
#this looks more reasonable by visual inspection anyway
import voronoi_utility

In [46]:
# trying to assess the issue preventing proper Voronoi diagram data structure generation with large sphere radii or sphere radii with many decimal places
prng_2 = numpy.random.RandomState(556)
large_sphere_radius = 1.0
cartesian_coord_array_large_radius = voronoi_utility.generate_random_array_spherical_generators(50,large_sphere_radius,prng_2)
random_dist_voronoi_instance = voronoi_utility.Voronoi_Sphere_Surface(cartesian_coord_array_large_radius,large_sphere_radius)
#dictionary_voronoi_region_vertices = random_dist_voronoi_instance.voronoi_region_vertices_spherical_surface()
facet_coordinate_array_Delaunay_triangulation = voronoi_utility.produce_triangle_vertex_coordinate_array_Delaunay_sphere(random_dist_voronoi_instance.hull_instance)
#print 'facet_coordinate_array_Delaunay_triangulation:', voronoi_utility.convert_cartesian_array_to_spherical_array(facet_coordinate_array_Delaunay_triangulation)
array_Voronoi_vertices = voronoi_utility.produce_array_Voronoi_vertices_on_sphere_surface(facet_coordinate_array_Delaunay_triangulation,large_sphere_radius,random_dist_voronoi_instance.sphere_centroid)
def spherical_distance(u,v):
    #print 'u:', voronoi_utility.convert_cartesian_array_to_spherical_array(u)
    #print 'v:', voronoi_utility.convert_cartesian_array_to_spherical_array(v)
    #dot_product = numpy.dot(u,v)
    #print 'dot_product:',dot_product
    spherical_u = voronoi_utility.convert_cartesian_array_to_spherical_array(u)
    spherical_v = voronoi_utility.convert_cartesian_array_to_spherical_array(v)
    spherical_u[0] = 1.0 #normalizing
    spherical_v[0] = 1.0
    normalized_u = voronoi_utility.convert_spherical_array_to_cartesian_array(spherical_u)
    normalized_v = voronoi_utility.convert_spherical_array_to_cartesian_array(spherical_v)
    normalized_dot_product = numpy.dot(normalized_u,normalized_v)
    #print 'normalized_dot_product:', normalized_dot_product
    spherical_distance = math.acos(normalized_dot_product)
    #print 'spherical_distance:', spherical_distance
    return spherical_distance


fig9 = plt.figure()
fig9.set_size_inches(19,19)
ax = fig9.add_subplot(111, projection='3d')
#ax.scatter(cartesian_coord_array_large_radius[...,0],cartesian_coord_array_large_radius[...,1],cartesian_coord_array_large_radius[...,2],c='red')
ax.scatter(array_Voronoi_vertices[...,0],array_Voronoi_vertices[...,1],array_Voronoi_vertices[...,2],c='red')

distance_matrix_Voronoi_vertices_to_generators = scipy.spatial.distance.cdist(array_Voronoi_vertices,cartesian_coord_array_large_radius,spherical_distance)



In [ ]:


In [45]:
distance_matrix_Voronoi_vertices_to_generators.sort(axis=1)
for Voronoi_vertex_generator_distances in distance_matrix_Voronoi_vertices_to_generators:
    print Voronoi_vertex_generator_distances[0:5]


[ 0.39300873  0.39300873  0.39300873  0.40371148  0.42374955]
[ 0.35944917  0.35944917  0.35944917  0.39266285  0.56853184]
[ 0.31767752  0.31767752  0.31767752  0.51769914  0.82013499]
[ 0.46523859  0.46523859  0.46523859  0.59924043  0.66456415]
[ 0.55263489  0.55263489  0.55263489  0.61610348  0.63059547]
[ 0.37101463  0.37101463  0.37101463  0.3865098   0.65877482]
[ 0.40436468  0.40436468  0.40436468  0.46711492  0.56578136]
[ 0.25173723  0.25173723  0.25173723  0.42434607  0.54613719]
[ 0.39011454  0.39011454  0.39011454  0.42813428  0.43636739]
[ 0.10812033  0.10812033  0.10812033  0.46231734  0.52155696]
[ 0.31116015  0.31116015  0.31116015  0.41364399  0.43317231]
[ 0.55875526  0.55875526  0.55875526  0.63746426  0.64467896]
[ 0.33764806  0.33764806  0.33764806  0.43086644  0.46177352]
[ 0.28618995  0.28618995  0.28618995  0.40451276  0.45789433]
[ 0.58249837  0.58249837  0.58249837  0.59768207  0.59880257]
[ 0.5813974   0.5813974   0.5813974   0.59156372  0.59816389]
[ 0.42580529  0.42580529  0.42580529  0.73107598  0.75166369]
[ 0.51683797  0.51683797  0.51683797  0.5687968   0.65576761]
[ 0.48545962  0.48545962  0.48545962  0.67248896  0.7264545 ]
[ 0.51451136  0.51451136  0.51451136  0.56697043  0.57430318]
[ 0.30285062  0.30285062  0.30285062  0.32818801  0.41386589]
[ 0.33340065  0.33340065  0.33340065  0.33473419  0.34203413]
[ 0.50463553  0.50463553  0.50463553  0.51224513  0.54385749]
[ 0.30958462  0.30958462  0.30958462  0.32141347  0.41904017]
[ 0.28186335  0.28186335  0.28186335  0.38541421  0.47220081]
[ 0.38022985  0.38022985  0.38022985  0.41112275  0.41263727]
[ 0.4034992   0.4034992   0.4034992   0.70571299  0.7357063 ]
[ 0.38400445  0.38400445  0.38400445  0.46292929  0.47650834]
[ 0.18724437  0.18724437  0.18724437  0.33498789  0.61904958]
[ 0.50110262  0.50110262  0.50110262  0.51388241  0.53724399]
[ 0.32350278  0.32350278  0.32350278  0.53640255  0.62991573]
[ 0.33565514  0.33565514  0.33565514  0.52704816  0.70208716]
[ 0.25368923  0.25368923  0.25368923  0.33284299  0.49466969]
[ 0.40187743  0.40187743  0.40187743  0.42979115  0.64951106]
[ 0.29817926  0.29817926  0.29817926  0.54758084  0.66508063]
[ 0.47819976  0.47819976  0.47819976  0.50738201  0.59557462]
[ 0.46410535  0.46410535  0.46410535  0.52365517  0.63439219]
[ 0.42692749  0.42692749  0.42692749  0.50124783  0.55862237]
[ 0.39523333  0.39523333  0.39523333  0.50659283  0.56104333]
[ 0.41280904  0.41280904  0.41280904  0.47056611  0.49283385]
[ 0.26989633  0.26989633  0.26989633  0.29202492  0.61986761]
[ 0.28814159  0.28814159  0.28814159  0.31154062  0.62265917]
[ 0.33380508  0.33380508  0.33380508  0.33447017  0.34304404]
[ 0.5126155   0.5126155   0.5126155   0.53399     0.65103485]
[ 0.46008682  0.46008682  0.46008682  0.55768814  0.61500266]
[ 0.29849746  0.29849746  0.29849746  0.4411112   0.53121361]
[ 0.4740187   0.4740187   0.4740187   0.55634291  0.58003934]
[ 0.5156411   0.5156411   0.5156411   0.64545089  0.6504358 ]
[ 0.53380547  0.53380547  0.53380547  0.61734111  0.62721063]
[ 0.29879678  0.29879678  0.29879678  0.32724508  0.38634793]
[ 0.33355045  0.33355045  0.33355045  0.34354787  0.35434993]
[ 0.32544819  0.32544819  0.32544819  0.33734302  0.35270237]
[ 0.20317536  0.20317536  0.20317536  0.64245867  0.78006517]
[ 0.32572079  0.32572079  0.32572079  0.41085416  0.52009216]
[ 0.28384049  0.28384049  0.28384049  0.46126881  0.49243966]
[ 0.51943436  0.51943436  0.51943436  0.67674066  0.72565607]
[ 0.47568598  0.47568598  0.47568598  0.52151516  0.73124843]
[ 0.49377283  0.49377283  0.49377283  0.51228752  0.72944067]
[ 0.41614012  0.41614012  0.41614012  0.66246431  0.82601986]
[ 0.42639417  0.42639417  0.42639417  0.42924579  0.46824383]
[ 0.4087893   0.4087893   0.4087893   0.42079796  0.44154011]
[ 0.40608922  0.40608922  0.40608922  0.42206025  0.43859067]
[ 0.31557499  0.31557499  0.31557499  0.36318019  0.55993745]
[ 0.41450034  0.41450034  0.41450034  0.55650709  0.59497751]
[ 0.47275964  0.47275964  0.47275964  0.49886204  0.70431489]
[ 0.35797193  0.35797193  0.35797193  0.50854009  0.75007451]
[ 0.51207439  0.51207439  0.51207439  0.6143982   0.63466821]
[ 0.42571231  0.42571231  0.42571231  0.43015315  0.46363779]
[ 0.39168177  0.39168177  0.39168177  0.45127987  0.48211646]
[ 0.2468822   0.2468822   0.2468822   0.31981897  0.507467  ]
[ 0.36198615  0.36198615  0.36198615  0.36996294  0.51155619]
[ 0.36579486  0.36579486  0.36579486  0.37166383  0.50179494]
[ 0.35140977  0.35140977  0.35140977  0.36093531  0.47062622]
[ 0.2173299   0.2173299   0.2173299   0.32610821  0.35807475]
[ 0.28241814  0.28241814  0.28241814  0.29450632  0.56483082]
[ 0.30353639  0.30353639  0.30353639  0.3750688   0.45749794]
[ 0.2628028   0.2628028   0.2628028   0.45576495  0.47756619]
[ 0.46958451  0.46958451  0.46958451  0.51251852  0.6066728 ]
[ 0.29132597  0.29132597  0.29132597  0.34760527  0.36134803]
[ 0.40687486  0.40687486  0.40687486  0.47829277  0.50655546]
[ 0.34183169  0.34183169  0.34183169  0.39849896  0.49178277]
[ 0.34973014  0.34973014  0.34973014  0.35344054  0.39524212]
[ 0.34163311  0.34163311  0.34163311  0.36346905  0.39765666]
[ 0.45277325  0.45277325  0.45277325  0.6019323   0.61203126]
[ 0.41890022  0.41890022  0.41890022  0.60720038  0.67604181]
[ 0.46422943  0.46422943  0.46422943  0.6231506   0.75935796]
[ 0.2660881   0.2660881   0.2660881   0.50208911  0.64142315]
[ 0.2943888   0.2943888   0.2943888   0.38715973  0.43950416]
[ 0.33329247  0.33329247  0.33329247  0.42570454  0.52282429]
[ 0.1399132   0.1399132   0.1399132   0.46766012  0.60102234]
[ 0.43413988  0.43413988  0.43413988  0.60260308  0.63659857]
[ 0.34452659  0.34452659  0.34452659  0.53249492  0.60796803]
[ 0.34207856  0.34207856  0.34207856  0.37601783  0.56380808]
[ 0.31878035  0.31878035  0.31878035  0.39295736  0.52631828]
[ 0.44813881  0.44813881  0.44813881  0.46031964  0.55065575]
[ 0.4490223   0.4490223   0.4490223   0.45638449  0.5465128 ]

Some of the above distance arrays are reasonable (3 equidistant generators) while others have a single closest generator, which shouldn't happen for a proper Voronoi diagram. Could this be sensitive to distances being calculated using a Euclidean metric rather than the true spherical surface metric?!


In [31]:
-2.95 + 3.14


Out[31]:
0.18999999999999995

In [10]:
#more debugging
prng_2 = numpy.random.RandomState(556)
large_sphere_radius = 1.0
cartesian_coord_array_large_radius = voronoi_utility.generate_random_array_spherical_generators(500,large_sphere_radius,prng_2)
random_dist_voronoi_instance = voronoi_utility.Voronoi_Sphere_Surface(cartesian_coord_array_large_radius,large_sphere_radius)
dictionary_surface_areas = random_dist_voronoi_instance.voronoi_region_surface_areas_spherical_surface()
sum(dictionary_surface_areas.itervalues()) / (4 * math.pi * (large_sphere_radius**2))


Out[10]:
0.96429087482093701

In [14]:
sphere_radius = 87.0
coordinates_on_sphere_1 = numpy.array([[0,0,1],[1,0,0]]) * sphere_radius
voronoi_utility.calculate_haversine_distance_between_spherical_points(coordinates_on_sphere_1[0],coordinates_on_sphere_1[1],sphere_radius)


lambda_1: 0.0
lambda_2: 1.57079632679
phi_1: 0.0
phi_2: 0.0
spherical_distance: 136.659280431
Out[14]:
136.659280431156