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


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

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]:
array([[ 0.81325255, -0.10395777, -0.57254962],
       [ 0.45141247, -0.58900654,  0.67029701],
       [-0.629159  , -0.77472842, -0.06288739],
       [-0.44992442,  0.10452284,  0.88692897],
       [ 0.        ,  0.        ,  0.        ]])

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]:
(1996, 3)

In [118]:
array_circumcenter_coords


Out[118]:
array([[-0.26579721,  0.3158127 , -0.29206454],
       [-0.38607866,  0.26836967,  0.17417512],
       [-0.25290039,  0.35967617,  0.24814029],
       ..., 
       [-0.05847558, -0.34502138,  0.35752067],
       [-0.01743171, -0.36129696,  0.34650686],
       [-0.01695518, -0.354884  ,  0.35305303]])

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]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f29c2f50750>

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]:
(1996, 3)

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)



In [290]:
#attempt to use the Delaunay tetrahedralization neighbour information to connect the Voronoi vertices around the generators (i.e., to produce the Voronoi regions)
#still a work in progress

def filter_tetrahedron_to_triangle(current_tetrahedron_coord_array):
    current_triangle_coord_array = [] #initialize as a list
    for row in current_tetrahedron_coord_array: #ugly to use for loop for this, but ok for now!
        if row[0] == 0 and row[1] == 0 and row[2] == 0: #filter out origin row
            continue
        else:
            current_triangle_coord_array.append(row)
    current_triangle_coord_array = np.array(current_triangle_coord_array)
    return current_triangle_coord_array

array_tetrahedra = tri.points[tri.simplices]
num_generators = len(tri.points) - 1 #don't want to work with the dummy generator at the origin
gen_counter = 0
list_of_sublists_of_tetrahedron_indices_surrounding_each_generator = []

for generator in tri.points[:-1]: #don't want to work with the dummy generator at the origin (which is the last point in the generator set)
    #print 'generator:', generator
    #identify the simplices (tetrahedra minus origin = triangles) that contain a given generator as a vertex
    indices_of_triangles_surrounding_generator = np.unique(np.where(tri.points[tri.simplices] == generator)[0])
    #print 'number of triangles surrounding generator:', len(indices_of_triangles_surrounding_generator)
    #print 'indices_of_triangles_surrounding_generator:', indices_of_triangles_surrounding_generator
    #pick any one of the triangles surrounding the generator and pick a non-generator vertex
    first_tetrahedron_index = indices_of_triangles_surrounding_generator[0]
    first_tetrahedron = array_tetrahedra[first_tetrahedron_index]
    #print 'first_tetrahedron:', first_tetrahedron
    first_triangle = filter_tetrahedron_to_triangle(first_tetrahedron)
    #print 'first_triangle:', first_triangle #first triangle has index 0 in the indices_of_triangles_surrounding_generator
    #pick one of the two non-generator vertices in the first triangle
    indices_non_generator_vertices_first_triangle = np.unique(np.where(first_triangle != generator)[0])
    #print 'indices_non_generator_vertices_first_triangle:', indices_non_generator_vertices_first_triangle
    
    ordered_list_tetrahedron_indices_surrounding_current_generator = [first_tetrahedron_index] 
    
    #determine the appropriate ordering of Voronoi vertices to close the Voronoi region (polygon) by traversing the Delaunay neighbour data structure from scipy
    vertices_remaining = len(indices_of_triangles_surrounding_generator) - 1
    #choose the neighbour opposite the first non-generator vertex of the first triangle
    neighbour_tetrahedral_index = tri.neighbors[first_tetrahedron_index][indices_non_generator_vertices_first_triangle[0]]
    #print 'neighbour_tetrahedral_index:', neighbour_tetrahedral_index
    ordered_list_tetrahedron_indices_surrounding_current_generator.append(neighbour_tetrahedral_index)
    #print 'ordered_list_tetrahedron_indices_surrounding_current_generator:', ordered_list_tetrahedron_indices_surrounding_current_generator
    vertices_remaining -= 1
    
    #for all subsequent triangles it is the common non-generator vertex with the previous neighbour that should be used to propagate the connection chain to the following neighbour
    #the common vertex with the previous neighbour is the the vertex of the previous neighbour that was NOT used to locate the current neighbour
    #since there are only two candidate vertices on the previous neighbour and I've chosen to use the vertex with index 0, the remaining vertex on the previous neighbour is the non-generator vertex with index 1
    common_vertex_coordinate = first_triangle[indices_non_generator_vertices_first_triangle[1]]
    #print 'first common_vertex_coordinate:', common_vertex_coordinate
    while vertices_remaining > 0:
        #print 'common_vertex_coordinate:', common_vertex_coordinate
        current_tetrahedron_index = ordered_list_tetrahedron_indices_surrounding_current_generator[-1]
        current_tetrahedron_coord_array = array_tetrahedra[current_tetrahedron_index]
        #it seems that it is not quite guaranteed that the last row of the tetrahedron is the origin point (though this is usually the case) -- will need a more robust filter to remove the origin and produce the relevant triangular simplex
        current_triangle_coord_array = filter_tetrahedron_to_triangle(current_tetrahedron_coord_array)
        #print 'current_triangle_coord_array:', current_triangle_coord_array
        #print 'current_tetrahedron_coord_array', current_tetrahedron_coord_array
        indices_candidate_vertices_current_triangle_excluding_generator = np.unique(np.where(current_triangle_coord_array != generator)[0])
        #print 'indices_candidate_vertices_current_triangle_excluding_generator:', indices_candidate_vertices_current_triangle_excluding_generator
        array_candidate_vertices = current_triangle_coord_array[indices_candidate_vertices_current_triangle_excluding_generator]
        #if len(indices_candidate_vertices_current_triangle_excluding_generator) > 2:
        #    print 'violating array:', array_candidate_vertices
        #    print 'generator', generator
        #    print 'current_triangle_coord_array', current_triangle_coord_array
        #    print 'current_tetrahedron_coord_array', current_tetrahedron_coord_array
        #print 'array_candidate_vertices:', array_candidate_vertices
        #print 'current_tetrahedron_coord_array:', current_tetrahedron_coord_array
        current_tetrahedron_index_for_neighbour_propagation = np.unique(np.where(current_tetrahedron_coord_array == common_vertex_coordinate)[0])
        #print 'current_tetrahedron_index_for_neighbour_propagation:', current_tetrahedron_index_for_neighbour_propagation
        next_tetrahedron_index_surrounding_generator = tri.neighbors[current_tetrahedron_index][current_tetrahedron_index_for_neighbour_propagation][0]
        if next_tetrahedron_index_surrounding_generator == -1: #debug -- trying to deal with case of 'boundary' simplex -- whatever that means on the sphere!!
            #print 'problematic array_candidate_vertices:', array_candidate_vertices
            #print 'related common_vertex_coordinate:', common_vertex_coordinate
            full_list_neighbour_indices_current_tetrahedron = tri.neighbors[current_tetrahedron_index]
            for neighbour_index in full_list_neighbour_indices_current_tetrahedron:
                if neighbour_index != -1 and neighbour_index not in ordered_list_tetrahedron_indices_surrounding_current_generator:
                    next_tetrahedron_index_surrounding_generator = neighbour_index
        #print 'next_tetrahedron_index_surrounding_generator:', next_tetrahedron_index_surrounding_generator
        common_vertex_coordinate = array_candidate_vertices[array_candidate_vertices != common_vertex_coordinate] #for the next iteration
        ordered_list_tetrahedron_indices_surrounding_current_generator.append(next_tetrahedron_index_surrounding_generator)
        #print 'ordered_list_tetrahedron_indices_surrounding_current_generator:', ordered_list_tetrahedron_indices_surrounding_current_generator
        vertices_remaining -= 1
    list_of_sublists_of_tetrahedron_indices_surrounding_each_generator.append(ordered_list_tetrahedron_indices_surrounding_current_generator)
    gen_counter += 1
    print 'generator', gen_counter, 'of', num_generators


generator 1 of 1000
generator 2 of 1000
generator 3 of 1000
generator 4 of 1000
generator 5 of 1000
generator 6 of 1000
generator 7 of 1000
generator 8 of 1000
generator 9 of 1000
generator 10 of 1000
generator 11 of 1000
generator 12 of 1000
generator 13 of 1000
generator 14 of 1000
generator 15 of 1000
generator 16 of 1000
generator 17 of 1000
generator 18 of 1000
generator 19 of 1000
generator 20 of 1000
generator 21 of 1000
generator 22 of 1000
generator 23 of 1000
generator 24 of 1000
generator 25 of 1000
generator 26 of 1000
generator 27 of 1000
generator 28 of 1000
generator 29 of 1000
generator 30 of 1000
generator 31 of 1000
generator 32 of 1000
generator 33 of 1000
generator 34 of 1000
generator 35 of 1000
generator 36 of 1000
generator 37 of 1000
generator 38 of 1000
generator 39 of 1000
generator 40 of 1000
generator 41 of 1000
generator 42 of 1000
generator 43 of 1000
generator 44 of 1000
generator 45 of 1000
generator 46 of 1000
generator 47 of 1000
generator 48 of 1000
generator 49 of 1000
generator 50 of 1000
generator 51 of 1000
generator 52 of 1000
generator 53 of 1000
generator 54 of 1000
generator 55 of 1000
generator 56 of 1000
generator 57 of 1000
generator 58 of 1000
generator 59 of 1000
generator 60 of 1000
generator 61 of 1000
generator 62 of 1000
generator 63 of 1000
generator 64 of 1000
generator 65 of 1000
generator 66 of 1000
generator 67 of 1000
generator 68 of 1000
generator 69 of 1000
generator 70 of 1000
generator 71 of 1000
generator 72 of 1000
generator 73 of 1000
generator 74 of 1000
generator 75 of 1000
generator 76 of 1000
generator 77 of 1000
generator 78 of 1000
generator 79 of 1000
generator 80 of 1000
generator 81 of 1000
generator 82 of 1000
generator 83 of 1000
generator 84 of 1000
generator 85 of 1000
generator 86 of 1000
generator 87 of 1000
generator 88 of 1000
generator 89 of 1000
generator 90 of 1000
generator 91 of 1000
generator 92 of 1000
generator 93 of 1000
generator 94 of 1000
generator 95 of 1000
generator 96 of 1000
generator 97 of 1000
generator 98 of 1000
generator 99 of 1000
generator 100 of 1000
generator 101 of 1000
generator 102 of 1000
generator 103 of 1000
generator 104 of 1000
generator 105 of 1000
generator 106 of 1000
generator 107 of 1000
generator 108 of 1000
generator 109 of 1000
generator 110 of 1000
generator 111 of 1000
generator 112 of 1000
generator 113 of 1000
generator 114 of 1000
generator 115 of 1000
generator 116 of 1000
generator 117 of 1000
generator 118 of 1000
generator 119 of 1000
generator 120 of 1000
generator 121 of 1000
generator 122 of 1000
generator 123 of 1000
generator 124 of 1000
generator 125 of 1000
generator 126 of 1000
generator 127 of 1000
generator 128 of 1000
generator 129 of 1000
generator 130 of 1000
generator 131 of 1000
generator 132 of 1000
generator 133 of 1000
generator 134 of 1000
generator 135 of 1000
generator 136 of 1000
generator 137 of 1000
generator 138 of 1000
generator 139 of 1000
generator 140 of 1000
generator 141 of 1000
generator 142 of 1000
generator 143 of 1000
generator 144 of 1000
generator 145 of 1000
generator 146 of 1000
generator 147 of 1000
generator 148 of 1000
generator 149 of 1000
generator 150 of 1000
generator 151 of 1000
generator 152 of 1000
generator 153 of 1000
generator 154 of 1000
generator 155 of 1000
generator 156 of 1000
generator 157 of 1000
generator 158 of 1000
generator 159 of 1000
generator 160 of 1000
generator 161 of 1000
generator 162 of 1000
generator 163 of 1000
generator 164 of 1000
generator 165 of 1000
generator 166 of 1000
generator 167 of 1000
generator 168 of 1000
generator 169 of 1000
generator 170 of 1000
generator 171 of 1000
generator 172 of 1000
generator 173 of 1000
generator 174 of 1000
generator 175 of 1000
generator 176 of 1000
generator 177 of 1000
generator 178 of 1000
generator 179 of 1000
generator 180 of 1000
generator 181 of 1000
generator 182 of 1000
generator 183 of 1000
generator 184 of 1000
generator 185 of 1000
generator 186 of 1000
generator 187 of 1000
generator 188 of 1000
generator 189 of 1000
generator 190 of 1000
generator 191 of 1000
generator 192 of 1000
generator 193 of 1000
generator 194 of 1000
generator 195 of 1000
generator 196 of 1000
generator 197 of 1000
generator 198 of 1000
generator 199 of 1000
generator 200 of 1000
generator 201 of 1000
generator 202 of 1000
generator 203 of 1000
generator 204 of 1000
generator 205 of 1000
generator 206 of 1000
generator 207 of 1000
generator 208 of 1000
generator 209 of 1000
generator 210 of 1000
generator 211 of 1000
generator 212 of 1000
generator 213 of 1000
generator 214 of 1000
generator 215 of 1000
generator 216 of 1000
generator 217 of 1000
generator 218 of 1000
generator 219 of 1000
generator 220 of 1000
generator 221 of 1000
generator 222 of 1000
generator 223 of 1000
generator 224 of 1000
generator 225 of 1000
generator 226 of 1000
generator 227 of 1000
generator 228 of 1000
generator 229 of 1000
generator 230 of 1000
generator 231 of 1000
generator 232 of 1000
generator 233 of 1000
generator 234 of 1000
generator 235 of 1000
generator 236 of 1000
generator 237 of 1000
generator 238 of 1000
generator 239 of 1000
generator 240 of 1000
generator 241 of 1000
generator 242 of 1000
generator 243 of 1000
generator 244 of 1000
generator 245 of 1000
generator 246 of 1000
generator 247 of 1000
generator 248 of 1000
generator 249 of 1000
generator 250 of 1000
generator 251 of 1000
generator 252 of 1000
generator 253 of 1000
generator 254 of 1000
generator 255 of 1000
generator 256 of 1000
generator 257 of 1000
generator 258 of 1000
generator 259 of 1000
generator 260 of 1000
generator 261 of 1000
generator 262 of 1000
generator 263 of 1000
generator 264 of 1000
generator 265 of 1000
generator 266 of 1000
generator 267 of 1000
generator 268 of 1000
generator 269 of 1000
generator 270 of 1000
generator 271 of 1000
generator 272 of 1000
generator 273 of 1000
generator 274 of 1000
generator 275 of 1000
generator 276 of 1000
generator 277 of 1000
generator 278 of 1000
generator 279 of 1000
generator 280 of 1000
generator 281 of 1000
generator 282 of 1000
generator 283 of 1000
generator 284 of 1000
generator 285 of 1000
generator 286 of 1000
generator 287 of 1000
generator 288 of 1000
generator 289 of 1000
generator 290 of 1000
generator 291 of 1000
generator 292 of 1000
generator 293 of 1000
generator 294 of 1000
generator 295 of 1000
generator 296 of 1000
generator 297 of 1000
generator 298 of 1000
generator 299 of 1000
generator 300 of 1000
generator 301 of 1000
generator 302 of 1000
generator 303 of 1000
generator 304 of 1000
generator 305 of 1000
generator 306 of 1000
generator 307 of 1000
generator 308 of 1000
generator 309 of 1000
generator 310 of 1000
generator 311 of 1000
generator 312 of 1000
generator 313 of 1000
generator 314 of 1000
generator 315 of 1000
generator 316 of 1000
generator 317 of 1000
generator 318 of 1000
generator 319 of 1000
generator 320 of 1000
generator 321 of 1000
generator 322 of 1000
generator 323 of 1000
generator 324 of 1000
generator 325 of 1000
generator 326 of 1000
generator 327 of 1000
generator 328 of 1000
generator 329 of 1000
generator 330 of 1000
generator 331 of 1000
generator 332 of 1000
generator 333 of 1000
generator 334 of 1000
generator 335 of 1000
generator 336 of 1000
generator 337 of 1000
generator 338 of 1000
generator 339 of 1000
generator 340 of 1000
generator 341 of 1000
generator 342 of 1000
generator 343 of 1000
generator 344 of 1000
generator 345 of 1000
generator 346 of 1000
generator 347 of 1000
generator 348 of 1000
generator 349 of 1000
generator 350 of 1000
generator 351 of 1000
generator 352 of 1000
generator 353 of 1000
generator 354 of 1000
generator 355 of 1000
generator 356 of 1000
generator 357 of 1000
generator 358 of 1000
generator 359 of 1000
generator 360 of 1000
generator 361 of 1000
generator 362 of 1000
generator 363 of 1000
generator 364 of 1000
generator 365 of 1000
generator 366 of 1000
generator 367 of 1000
generator 368 of 1000
generator 369 of 1000
generator 370 of 1000
generator 371 of 1000
generator 372 of 1000
generator 373 of 1000
generator 374 of 1000
generator 375 of 1000
generator 376 of 1000
generator 377 of 1000
generator 378 of 1000
generator 379 of 1000
generator 380 of 1000
generator 381 of 1000
generator 382 of 1000
generator 383 of 1000
generator 384 of 1000
generator 385 of 1000
generator 386 of 1000
generator 387 of 1000
generator 388 of 1000
generator 389 of 1000
generator 390 of 1000
generator 391 of 1000
generator 392 of 1000
generator 393 of 1000
generator 394 of 1000
generator 395 of 1000
generator 396 of 1000
generator 397 of 1000
generator 398 of 1000
generator 399 of 1000
generator 400 of 1000
generator 401 of 1000
generator 402 of 1000
generator 403 of 1000
generator 404 of 1000
generator 405 of 1000
generator 406 of 1000
generator 407 of 1000
generator 408 of 1000
generator 409 of 1000
generator 410 of 1000
generator 411 of 1000
generator 412 of 1000
generator 413 of 1000
generator 414 of 1000
generator 415 of 1000
generator 416 of 1000
generator 417 of 1000
generator 418 of 1000
generator 419 of 1000
generator 420 of 1000
generator 421 of 1000
generator 422 of 1000
generator 423 of 1000
generator 424 of 1000
generator 425 of 1000
generator 426 of 1000
generator 427 of 1000
generator 428 of 1000
generator 429 of 1000
generator 430 of 1000
generator 431 of 1000
generator 432 of 1000
generator 433 of 1000
generator 434 of 1000
generator 435 of 1000
generator 436 of 1000
generator 437 of 1000
generator 438 of 1000
generator 439 of 1000
generator 440 of 1000
generator 441 of 1000
generator 442 of 1000
generator 443 of 1000
generator 444 of 1000
generator 445 of 1000
generator 446 of 1000
generator 447 of 1000
generator 448 of 1000
generator 449 of 1000
generator 450 of 1000
generator 451 of 1000
generator 452 of 1000
generator 453 of 1000
generator 454 of 1000
generator 455 of 1000
generator 456 of 1000
generator 457 of 1000
generator 458 of 1000
generator 459 of 1000
generator 460 of 1000
generator 461 of 1000
generator 462 of 1000
generator 463 of 1000
generator 464 of 1000
generator 465 of 1000
generator 466 of 1000
generator 467 of 1000
generator 468 of 1000
generator 469 of 1000
generator 470 of 1000
generator 471 of 1000
generator 472 of 1000
generator 473 of 1000
generator 474 of 1000
generator 475 of 1000
generator 476 of 1000
generator 477 of 1000
generator 478 of 1000
generator 479 of 1000
generator 480 of 1000
generator 481 of 1000
generator 482 of 1000
generator 483 of 1000
generator 484 of 1000
generator 485 of 1000
generator 486 of 1000
generator 487 of 1000
generator 488 of 1000
generator 489 of 1000
generator 490 of 1000
generator 491 of 1000
generator 492 of 1000
generator 493 of 1000
generator 494 of 1000
generator 495 of 1000
generator 496 of 1000
generator 497 of 1000
generator 498 of 1000
generator 499 of 1000
generator 500 of 1000
generator 501 of 1000
generator 502 of 1000
generator 503 of 1000
generator 504 of 1000
generator 505 of 1000
generator 506 of 1000
generator 507 of 1000
generator 508 of 1000
generator 509 of 1000
generator 510 of 1000
generator 511 of 1000
generator 512 of 1000
generator 513 of 1000
generator 514 of 1000
generator 515 of 1000
generator 516 of 1000
generator 517 of 1000
generator 518 of 1000
generator 519 of 1000
generator 520 of 1000
generator 521 of 1000
generator 522 of 1000
generator 523 of 1000
generator 524 of 1000
generator 525 of 1000
generator 526 of 1000
generator 527 of 1000
generator 528 of 1000
generator 529 of 1000
generator 530 of 1000
generator 531 of 1000
generator 532 of 1000
generator 533 of 1000
generator 534 of 1000
generator 535 of 1000
generator 536 of 1000
generator 537 of 1000
generator 538 of 1000
generator 539 of 1000
generator 540 of 1000
generator 541 of 1000
generator 542 of 1000
generator 543 of 1000
generator 544 of 1000
generator 545 of 1000
generator 546 of 1000
generator 547 of 1000
generator 548 of 1000
generator 549 of 1000
generator 550 of 1000
generator 551 of 1000
generator 552 of 1000
generator 553 of 1000
generator 554 of 1000
generator 555 of 1000
generator 556 of 1000
generator 557 of 1000
generator 558 of 1000
generator 559 of 1000
generator 560 of 1000
generator 561 of 1000
generator 562 of 1000
generator 563 of 1000
generator 564 of 1000
generator 565 of 1000
generator 566 of 1000
generator 567 of 1000
generator 568 of 1000
generator 569 of 1000
generator 570 of 1000
generator 571 of 1000
generator 572 of 1000
generator 573 of 1000
generator 574 of 1000
generator 575 of 1000
generator 576 of 1000
generator 577 of 1000
generator 578 of 1000
generator 579 of 1000
generator 580 of 1000
generator 581 of 1000
generator 582 of 1000
generator 583 of 1000
generator 584 of 1000
generator 585 of 1000
generator 586 of 1000
generator 587 of 1000
generator 588 of 1000
generator 589 of 1000
generator 590 of 1000
generator 591 of 1000
generator 592 of 1000
generator 593 of 1000
generator 594 of 1000
generator 595 of 1000
generator 596 of 1000
generator 597 of 1000
generator 598 of 1000
generator 599 of 1000
generator 600 of 1000
generator 601 of 1000
generator 602 of 1000
generator 603 of 1000
generator 604 of 1000
generator 605 of 1000
generator 606 of 1000
generator 607 of 1000
generator 608 of 1000
generator 609 of 1000
generator 610 of 1000
generator 611 of 1000
generator 612 of 1000
generator 613 of 1000
generator 614 of 1000
generator 615 of 1000
generator 616 of 1000
generator 617 of 1000
generator 618 of 1000
generator 619 of 1000
generator 620 of 1000
generator 621 of 1000
generator 622 of 1000
generator 623 of 1000
generator 624 of 1000
generator 625 of 1000
generator 626 of 1000
generator 627 of 1000
generator 628 of 1000
generator 629 of 1000
generator 630 of 1000
generator 631 of 1000
generator 632 of 1000
generator 633 of 1000
generator 634 of 1000
generator 635 of 1000
generator 636 of 1000
generator 637 of 1000
generator 638 of 1000
generator 639 of 1000
generator 640 of 1000
generator 641 of 1000
generator 642 of 1000
generator 643 of 1000
generator 644 of 1000
generator 645 of 1000
generator 646 of 1000
generator 647 of 1000
generator 648 of 1000
generator 649 of 1000
generator 650 of 1000
generator 651 of 1000
generator 652 of 1000
generator 653 of 1000
generator 654 of 1000
generator 655 of 1000
generator 656 of 1000
generator 657 of 1000
generator 658 of 1000
generator 659 of 1000
generator 660 of 1000
generator 661 of 1000
generator 662 of 1000
generator 663 of 1000
generator 664 of 1000
generator 665 of 1000
generator 666 of 1000
generator 667 of 1000
generator 668 of 1000
generator 669 of 1000
generator 670 of 1000
generator 671 of 1000
generator 672 of 1000
generator 673 of 1000
generator 674 of 1000
generator 675 of 1000
generator 676 of 1000
generator 677 of 1000
generator 678 of 1000
generator 679 of 1000
generator 680 of 1000
generator 681 of 1000
generator 682 of 1000
generator 683 of 1000
generator 684 of 1000
generator 685 of 1000
generator 686 of 1000
generator 687 of 1000
generator 688 of 1000
generator 689 of 1000
generator 690 of 1000
generator 691 of 1000
generator 692 of 1000
generator 693 of 1000
generator 694 of 1000
generator 695 of 1000
generator 696 of 1000
generator 697 of 1000
generator 698 of 1000
generator 699 of 1000
generator 700 of 1000
generator 701 of 1000
generator 702 of 1000
generator 703 of 1000
generator 704 of 1000
generator 705 of 1000
generator 706 of 1000
generator 707 of 1000
generator 708 of 1000
generator 709 of 1000
generator 710 of 1000
generator 711 of 1000
generator 712 of 1000
generator 713 of 1000
generator 714 of 1000
generator 715 of 1000
generator 716 of 1000
generator 717 of 1000
generator 718 of 1000
generator 719 of 1000
generator 720 of 1000
generator 721 of 1000
generator 722 of 1000
generator 723 of 1000
generator 724 of 1000
generator 725 of 1000
generator 726 of 1000
generator 727 of 1000
generator 728 of 1000
generator 729 of 1000
generator 730 of 1000
generator 731 of 1000
generator 732 of 1000
generator 733 of 1000
generator 734 of 1000
generator 735 of 1000
generator 736 of 1000
generator 737 of 1000
generator 738 of 1000
generator 739 of 1000
generator 740 of 1000
generator 741 of 1000
generator 742 of 1000
generator 743 of 1000
generator 744 of 1000
generator 745 of 1000
generator 746 of 1000
generator 747 of 1000
generator 748 of 1000
generator 749 of 1000
generator 750 of 1000
generator 751 of 1000
generator 752 of 1000
generator 753 of 1000
generator 754 of 1000
generator 755 of 1000
generator 756 of 1000
generator 757 of 1000
generator 758 of 1000
generator 759 of 1000
generator 760 of 1000
generator 761 of 1000
generator 762 of 1000
generator 763 of 1000
generator 764 of 1000
generator 765 of 1000
generator 766 of 1000
generator 767 of 1000
generator 768 of 1000
generator 769 of 1000
generator 770 of 1000
generator 771 of 1000
generator 772 of 1000
generator 773 of 1000
generator 774 of 1000
generator 775 of 1000
generator 776 of 1000
generator 777 of 1000
generator 778 of 1000
generator 779 of 1000
generator 780 of 1000
generator 781 of 1000
generator 782 of 1000
generator 783 of 1000
generator 784 of 1000
generator 785 of 1000
generator 786 of 1000
generator 787 of 1000
generator 788 of 1000
generator 789 of 1000
generator 790 of 1000
generator 791 of 1000
generator 792 of 1000
generator 793 of 1000
generator 794 of 1000
generator 795 of 1000
generator 796 of 1000
generator 797 of 1000
generator 798 of 1000
generator 799 of 1000
generator 800 of 1000
generator 801 of 1000
generator 802 of 1000
generator 803 of 1000
generator 804 of 1000
generator 805 of 1000
generator 806 of 1000
generator 807 of 1000
generator 808 of 1000
generator 809 of 1000
generator 810 of 1000
generator 811 of 1000
generator 812 of 1000
generator 813 of 1000
generator 814 of 1000
generator 815 of 1000
generator 816 of 1000
generator 817 of 1000
generator 818 of 1000
generator 819 of 1000
generator 820 of 1000
generator 821 of 1000
generator 822 of 1000
generator 823 of 1000
generator 824 of 1000
generator 825 of 1000
generator 826 of 1000
generator 827 of 1000
generator 828 of 1000
generator 829 of 1000
generator 830 of 1000
generator 831 of 1000
generator 832 of 1000
generator 833 of 1000
generator 834 of 1000
generator 835 of 1000
generator 836 of 1000
generator 837 of 1000
generator 838 of 1000
generator 839 of 1000
generator 840 of 1000
generator 841 of 1000
generator 842 of 1000
generator 843 of 1000
generator 844 of 1000
generator 845 of 1000
generator 846 of 1000
generator 847 of 1000
generator 848 of 1000
generator 849 of 1000
generator 850 of 1000
generator 851 of 1000
generator 852 of 1000
generator 853 of 1000
generator 854 of 1000
generator 855 of 1000
generator 856 of 1000
generator 857 of 1000
generator 858 of 1000
generator 859 of 1000
generator 860 of 1000
generator 861 of 1000
generator 862 of 1000
generator 863 of 1000
generator 864 of 1000
generator 865 of 1000
generator 866 of 1000
generator 867 of 1000
generator 868 of 1000
generator 869 of 1000
generator 870 of 1000
generator 871 of 1000
generator 872 of 1000
generator 873 of 1000
generator 874 of 1000
generator 875 of 1000
generator 876 of 1000
generator 877 of 1000
generator 878 of 1000
generator 879 of 1000
generator 880 of 1000
generator 881 of 1000
generator 882 of 1000
generator 883 of 1000
generator 884 of 1000
generator 885 of 1000
generator 886 of 1000
generator 887 of 1000
generator 888 of 1000
generator 889 of 1000
generator 890 of 1000
generator 891 of 1000
generator 892 of 1000
generator 893 of 1000
generator 894 of 1000
generator 895 of 1000
generator 896 of 1000
generator 897 of 1000
generator 898 of 1000
generator 899 of 1000
generator 900 of 1000
generator 901 of 1000
generator 902 of 1000
generator 903 of 1000
generator 904 of 1000
generator 905 of 1000
generator 906 of 1000
generator 907 of 1000
generator 908 of 1000
generator 909 of 1000
generator 910 of 1000
generator 911 of 1000
generator 912 of 1000
generator 913 of 1000
generator 914 of 1000
generator 915 of 1000
generator 916 of 1000
generator 917 of 1000
generator 918 of 1000
generator 919 of 1000
generator 920 of 1000
generator 921 of 1000
generator 922 of 1000
generator 923 of 1000
generator 924 of 1000
generator 925 of 1000
generator 926 of 1000
generator 927 of 1000
generator 928 of 1000
generator 929 of 1000
generator 930 of 1000
generator 931 of 1000
generator 932 of 1000
generator 933 of 1000
generator 934 of 1000
generator 935 of 1000
generator 936 of 1000
generator 937 of 1000
generator 938 of 1000
generator 939 of 1000
generator 940 of 1000
generator 941 of 1000
generator 942 of 1000
generator 943 of 1000
generator 944 of 1000
generator 945 of 1000
generator 946 of 1000
generator 947 of 1000
generator 948 of 1000
generator 949 of 1000
generator 950 of 1000
generator 951 of 1000
generator 952 of 1000
generator 953 of 1000
generator 954 of 1000
generator 955 of 1000
generator 956 of 1000
generator 957 of 1000
generator 958 of 1000
generator 959 of 1000
generator 960 of 1000
generator 961 of 1000
generator 962 of 1000
generator 963 of 1000
generator 964 of 1000
generator 965 of 1000
generator 966 of 1000
generator 967 of 1000
generator 968 of 1000
generator 969 of 1000
generator 970 of 1000
generator 971 of 1000
generator 972 of 1000
generator 973 of 1000
generator 974 of 1000
generator 975 of 1000
generator 976 of 1000
generator 977 of 1000
generator 978 of 1000
generator 979 of 1000
generator 980 of 1000
generator 981 of 1000
generator 982 of 1000
generator 983 of 1000
generator 984 of 1000
generator 985 of 1000
generator 986 of 1000
generator 987 of 1000
generator 988 of 1000
generator 989 of 1000
generator 990 of 1000
generator 991 of 1000
generator 992 of 1000
generator 993 of 1000
generator 994 of 1000
generator 995 of 1000
generator 996 of 1000
generator 997 of 1000
generator 998 of 1000
generator 999 of 1000
generator 1000 of 1000

In [291]:
#plot to confirm that the above algorithm is working for a single generator and its surrounding Voronoi vertices:
fig_test_ordering_single_voronoi_polygon = plt.figure()
ax = fig_test_ordering_single_voronoi_polygon.add_subplot('111', projection = '3d')
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
voronoi_region_vertices_ordered = array_Voronoi_vertices[ordered_list_tetrahedron_indices_surrounding_current_generator]
print voronoi_region_vertices_ordered
polygon = Poly3DCollection([voronoi_region_vertices_ordered], alpha = 0.5)
polygon.set_color('blue')
ax.add_collection3d(polygon)
ax.scatter(voronoi_region_vertices_ordered[...,0], voronoi_region_vertices_ordered[...,1], voronoi_region_vertices_ordered[...,2], c = 'k', s = 10)
ax.scatter(generator[...,0], generator[...,1], generator[...,2], c= 'red', edgecolor = 'none')
ax.azim = 0
ax.elev = -180
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_zlim((0.4,0.8))
ax.set_ylim((0.0,0.4))


[[-0.35963455  0.16525321  0.91834327]
 [-0.37405408  0.11836398  0.91982255]
 [-0.46530628  0.06476555  0.88277714]
 [-0.50211722  0.08889887  0.86021816]
 [-0.46863676  0.15046071  0.87048329]
 [-0.39790958  0.17871525  0.89984933]]
Out[291]:
(0.0, 0.4)

So, I appear to have successfully constructed the Voronoi cell for the first generator using the algorithm Ross proposed, with proper ordering of Voronoi vertices. Now need to produce the overall spherical Voronoi diagram and probably also select a subset of the Voronoi regions for visual inspection in this manner.


In [295]:
#plot the spherical Voronoi diagram
fig_spherical_Voronoi = plt.figure()
ax = fig_spherical_Voronoi.add_subplot('111', projection = '3d')
import matplotlib.colors as colors
for list_tetrahedron_indices_around_current_generator in list_of_sublists_of_tetrahedron_indices_surrounding_each_generator:
    random_color = colors.rgb2hex(sp.rand(3))
    current_voronoi_region_vertices_ordered = array_Voronoi_vertices[list_tetrahedron_indices_around_current_generator]
    polygon = Poly3DCollection([current_voronoi_region_vertices_ordered], alpha = 1.0)
    polygon.set_color(random_color)
    ax.add_collection3d(polygon)
    
ax.set_xlim((-1,1))
ax.set_ylim((-1,1))
ax.set_zlim((-1,1))
fig_spherical_Voronoi.set_size_inches(8,8)



In [297]:
#check that the surface area reconstituted from the full set of Voronoi polygons is close to the theoretical surface area of the a unit sphere:
total_surface_area = 0
for list_tetrahedron_indices_around_current_generator in list_of_sublists_of_tetrahedron_indices_surrounding_each_generator:
    current_voronoi_region_vertices_ordered = array_Voronoi_vertices[list_tetrahedron_indices_around_current_generator]
    total_surface_area += voronoi_utility.calculate_surface_area_of_a_spherical_Voronoi_polygon(current_voronoi_region_vertices_ordered, 1.0)

import math
theoretical_surface_area = 4.0 * math.pi 
print 'total_surface_area', total_surface_area
print 'theoretical_surface_area', theoretical_surface_area
print 'percent surface area reconstitution:', (total_surface_area / theoretical_surface_area) * 100.


total_surface_area 12.5663706144
theoretical_surface_area 12.5663706144
percent surface area reconstitution: 100.0

Can't get any better than that!


In [298]:
#just to be absolutely sure, up next I'll want to plot a subset of the Voronoi regions and their associated generators, but I can't imagine there are any self-intersecting Voronoi polygons now that we have full reconstitution

In [311]:
fig_panel_zoom_voronoi_regions = plt.figure()
subplot_number = 1
array_random_indices = np.random.randint(random_coordinate_array.shape[0], size = 30)
for generator_index in array_random_indices:
    ax = fig_panel_zoom_voronoi_regions.add_subplot(6,5,subplot_number, projection = '3d')
    list_tetrahedron_indices_around_current_generator = list_of_sublists_of_tetrahedron_indices_surrounding_each_generator[generator_index]
    current_voronoi_region_vertices_ordered = array_Voronoi_vertices[list_tetrahedron_indices_around_current_generator]
    generator = random_coordinate_array[generator_index]
    ax.scatter(generator[...,0], generator[...,1], generator[...,2], c= 'k', edgecolor = 'none')
    polygon = Poly3DCollection([current_voronoi_region_vertices_ordered], alpha = 0.4)
    polygon.set_color('blue')
    ax.add_collection3d(polygon)
    ax.set_xlim((generator[...,0].min() - 0.1, generator[...,0].max() + 0.1))
    ax.set_ylim((generator[...,1].min() - 0.1, generator[...,1].max() + 0.1))
    ax.set_zlim((generator[...,2].min() - 0.1, generator[...,2].max() + 0.1))
    subplot_number += 1
    
    
fig_panel_zoom_voronoi_regions.set_size_inches(20,10)


The results look excellent -- I don't see any indication of issues with the randomly-selected Voronoi regions, after accounting for 3D perspectives, etc.


In [338]:
#now test implementation of above improvements into the main module
voronoi_instance = voronoi_utility.Voronoi_Sphere_Surface(random_coordinate_array[:-1,...], 1.0) #allow the main module to add the zero point back in

In [339]:
dictionary_voronoi_polygon_vertices = voronoi_instance.voronoi_region_vertices_spherical_surface()

In [342]:
#plot to confirm that the overall spherical Voronoi diagram generated by the main module looks reasonable
fig_main_module_spherical_Voronoi = plt.figure()
ax = fig_main_module_spherical_Voronoi.add_subplot('111', projection = '3d')
for generator_index, voronoi_region in dictionary_voronoi_polygon_vertices.iteritems():
    #print voronoi_region
    random_color = colors.rgb2hex(sp.rand(3))
    polygon = Poly3DCollection([voronoi_region], alpha = 1.0)
    polygon.set_color(random_color)
    ax.add_collection3d(polygon)
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
ax.set_zlim(-1,1)
fig_main_module_spherical_Voronoi.set_size_inches(8,8)



In [343]:
#overall spherical Voronoi diagram produced by the main module look excellent using the new algorithm
#now check the reconstituted surface area
dictionary_surface_areas_per_generator_Voronoi_region = voronoi_instance.voronoi_region_surface_areas_spherical_surface()

In [344]:
reconstituted_surface_area = sum(dictionary_surface_areas_per_generator_Voronoi_region.itervalues())
percent_surface_area_reconstituted = (reconstituted_surface_area / (4 * math.pi)) * 100.0
print 'percent_surface_area_reconstituted:', percent_surface_area_reconstituted


percent_surface_area_reconstituted: 100.0

In [345]:
#again, excellent, now check a random subset of the polygons up close:
fig_main_module_polygon_check = plt.figure()
subplot_number = 1
for generator_index in array_random_indices:
    ax = fig_main_module_polygon_check.add_subplot(6,5,subplot_number, projection = '3d')
    voronoi_region_vertices_ordered = dictionary_voronoi_polygon_vertices[generator_index]
    generator = random_coordinate_array[generator_index]
    ax.scatter(generator[...,0], generator[...,1], generator[...,2], c= 'k', edgecolor = 'none')
    polygon = Poly3DCollection([voronoi_region_vertices_ordered], alpha = 0.4)
    polygon.set_color('blue')
    ax.add_collection3d(polygon)
    ax.set_xlim((generator[...,0].min() - 0.1, generator[...,0].max() + 0.1))
    ax.set_ylim((generator[...,1].min() - 0.1, generator[...,1].max() + 0.1))
    ax.set_zlim((generator[...,2].min() - 0.1, generator[...,2].max() + 0.1))
    subplot_number += 1
    
    
fig_main_module_polygon_check.set_size_inches(20,10)



In [347]:
#looks excellent, now I think the main module (with its improved implementation) is ready for testing with real (experimental / computational) data

In [ ]:


In [ ]:


In [ ]: