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]:
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()
In [8]:
array_facet_normals
In [9]:
generator_Voronoi_region_dictionary,dictionary_sorted_Voronoi_point_coordinates_for_each_generator = voronoi_instance.Voronoi_polygons_spherical_surface()
In [10]:
dictionary_sorted_Voronoi_point_coordinates_for_each_generator
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
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]:
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()
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]:
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]:
In [60]:
numpy.diff(derivative_vector_3D,axis=0)
Out[60]:
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]: