Thanks CJ. I've made your specific suggested changes. I added a new task for improving the concision of the variable names to the to-do list above.

Cool. :) I think I can fix the problem with the docstrings (and maybe some PEP8 issues along the way).

@tylerjereddy a few comments on open items on the list:

add SphericalVoronoi to release notes / reference guide (?!) as indicated in development FAQ

Addition to the reference guide is already achieved by your edit to the docstring of spatial/__init__.py. And I wouldn't worry about the release notes, we can do that afterwards. So this box can be checked.

the organization of functions defined outside the SphericalVoronoi class may have to be adjusted ...

This isn't necessary, your code will take very little time to import. Everything looks good as is.

there are almost certainly good prospects for additional vectorization of some components of the code and / or usage of Cython for critical loops; I'm less familiar with the latter but open to giving it a go if encouraged

Do you think that your code is unacceptably slow (for real-world size problems or compared to another implementation)? I see some for-loops that might be quite slow, but only optimize at the end if you think it's needed.

The plotting example works for me even though it's not pretty. Mayavi is indeed a too heavy dependency.

@rgommers It seems that some of the (essential and much appreciated) refactoring Nikolai has done has caused a pretty substantial loss in performance based on my informal benchmarks:

The Jupyter notebook used to produce the crude benchmarks has been made public (https://github.com/tylerjereddy/bench_spherical_voronoi/blob/master/assess_spherical_Voronoi_algorithm.ipynb).

So, I think it would be reasonable to expect the performance of the code in the scipy fork to roughly match the performance of the code from my repo (https://github.com/tylerjereddy/py_sphere_Voronoi). At the moment, I don't think it is acceptable, even for my personal use at ~20,000 input points maximum.

Why is the scipy fork performing slower? Speculation probably isn't productive without line profiling (as usual). However, the new SphericalVoronoi class automatically generates unsorted Voronoi regions when it is initialized (mostly for consistency with things like scipy.spatial.ConvexHull or Voronoi where the primary use data structure is automatically generated when the class is instantiated). However, my suspicion is that splitting the general Voronoi vertex generation and Voronoi region vertex sorting operations in this fashion has had some cost (they are more 'integrated' in old code). [the challenges of sorting the Voronoi region vertices are discussed in the initial presentation of the algorithm to the Python community: https://www.youtube.com/watch?v=gxNa9BD5CnQ]

I'm going to have to read the refactored code carefully. One question--if I can manage to restore the performance is there an accepted way to incorporate a "performance test" to prevent future "performance regressions" of this nature? I imagine unit tests aren't quite the right place for various reasons--execution time could change with computers / other parts of code improving, etc., though I suppose one could decorate the "unit test" with the appropriate flags such that it only runs when a comprehensive suite is desired.

Hi @tylerjereddy, a quick reply about benchmarking (without having time to look into the details here now): we'd love to have proper benchmarks for any new feature that's added. Those live in https://github.com/scipy/scipy/tree/master/benchmarks and are easy to add.

@rgommers It also seems that removing import scipy from the Examples section of the docstring is causing an issue when calling scipy.rand() in one of the Travis tests, even though the docs build for me locally. The control flow does import from scipy before reaching that point, so that's a bit strange if scipy is supposed to be imported by default.

scipy.rand isn't a scipy but a numpy function: try scipy.rand is np.random.rand

There is nothing in the scipy namespace that should be used (except test()), it's only there for historical reasons. See http://docs.scipy.org/doc/scipy/reference/api.html#guidelines-for-importing-functions-from-scipy

That is a massive difference in performance indeed, looks like it needs fixing.

Ok, I've restored (and actually improved upon) the performance:

I have a few more comments and concerns that I'll try to summarize later.

• Performance

  • The Competition: There seems to be a paucity of open source alternatives for calculating spherical Voronoi diagrams. There's an interactive implementation written by Jason Davies (https://www.jasondavies.com/maps/voronoi/) using the javascript D3 library. As this implementation appears to allow progressive (interactive) insertion of data points (and also proceeds via the Convex Hull), it will probably have the same quadratic time complexity as our algorithm (see Table 1 below).

    Benchmarking would take at least some effort as only interactive interaction is allowed (can't feed in array data to get output--which is part of the motivation for developing something more flexible like this in scipy). His example with all major airports in the world (https://www.jasondavies.com/maps/voronoi/airports/) cites about 3000 data points, which is easily within the well-performing region of our implementation.

    There's a closed source implementation in Wolfram (http://demonstrations.wolfram.com/VoronoiDiagramOnASphere/), but this also only appears to allow for interaction rather than feeding in data and returning vertices, regions, etc. It also briefly describes an algorithm that exploits the same properties as ours, and would likely have the same time complexity.

    Apparently, the bleeding edge version of CGAL should be capable of calculating spherical Voronoi diagrams, and indeed it was a former CGAL core developer that provided key insight into the method for sorting the Voronoi regions in my early attempts. I think setting up CGAL properly with C++ libraries can be tricky though. Last time I checked, Qhull made a brief comment about how the spherical Voronoi diagram could be computed, but didn't offer this functionality. If I had to guess, I would suspect that the CGAL implementation of a similar algorithm would run at least a bit faster than ours with all the effort they put into C++ optimization, etc. I can't find a reference to it in their documentation though.

  • Theory: There are faster (loglinear) algorithms recently published in the mathematics literature for spherical Voronoi diagrams (Table 1). I made several attempts to incorporate the stereographic projection paradigm, without much success. In the future, a faster implementation would be nice, but these can be really hard to implement. At least if we have a working slower algorithm it will serve as a useful unit test for correct results when attempting to improve with a fancier algorithm.

  • Generator limit: It appears that the upper limit on the number of testable generators (input data points) in my benchmarks relates to the amount of physical memory used (the space complexity). Although I had made some initial attempts to 'benchmark' the space complexity, this seems to be a bit harder to do well. Suffice it to say that it certainly seems clear that at around 70,000 generators a machine with 16 GB of RAM starts to choke. I routinely deal with ~20,000 generators when working with viruses. Hard to say what other application domains will require--certainly we can't produce a spherical Voronoi diagram with the coordinates of all ~3 trillion trees on our planet. Nonetheless, there are < 200 countries in the world for looking at epidemiology-related things on the planet surface, and surely many practical applications involve < 70,000 generators.

Table 1: Spherical Voronoi Algorithms

• Unit testing:

  • The unit tests have changed quite a lot from my original implementation. I think Nikolai felt there were too many 'functional' tests rather than proper unit tests. That was probably true. It is noteworthy that one of the removed unit tests verified an unproved (but empirically observed) property of spherical Voronoi diagrams--the upper limit on the number of Voronoi vertices is 2n - 4, where n is the number of generators.

    A similar relationship (2n - 5) has been formally proven in the plane [Theorem 4.12 in 'Discrete and Computational Geometry' (2011) Devadoss and O'Rourke], but not for a sphere (yet). Anyway, good to have a record of this here if that is proven in the future. Certainly looks true.

• Example in Docstring: I'm a bit confused by the new example that was coded in (see my annotations below). Actually, it may very well still be correct, but perhaps @niknow can confirm that these issues are just the result of the perspective / plotting limitations?!

Beyond that, I guess there's still a few check points in the original list above, although they may be less important now.

Regarding the Example in Docstring: I can confirm that there are visibility issues with the generators and the Voronoi vertices. If you rotate the plot in the interactive view, then sometimes generators or Voronoi vertices become visible/invisible depending on the current perspective on the plot. All in all plotting spherical Voronoi diagrams in this way is certainly not ideal. Using matplotlibs built in Poly3DCollection and scatter is only meant as a temporary solution.

Regarding the 2n-4 unittest: If one plays around with the interactive version from https://www.jasondavies.com/maps/voronoi/ , one can see that the claim that a spherical Voronoi diagram always has exactly v=2n-4 vertices is wrong in general:

This picture shows two diagrams with the same number n of generator points, but the right one has one Voronoi vertex more than the left one. However, the left one is a pretty special case, because the degree of that vertex in the graph is four (would our algorithm actually detect that?).
Since the Voronoi vertices are the projections of the circumcenters of the simplices of the Delaunay triangulation, the number of Voronoi vertices is bounded by the number of simplices of a Delaunay triangulation. This is bounded by the McMullen`s upper bound theorem (see for instance https://en.wikipedia.org/wiki/Upper_bound_theorem), which in our case gives v <= 2n-2. So I guess one could always check that.

The fact that the better test v=2n-4 seems to be correct in all your test cases is probably because the problem described above is unlikely to occur in practice and the following should be true: Let v be the number of vertices of a spherical Voronoi diagram. Assume that the union over all faces of the diagram covers the whole sphere and that every vertex has degree 3. Then v=2n-4.
(This can be seen as follows: Remove one face from the Voronoi diagram and project the remaining diagram onto the plane via stereographic projection centered at the generator of the face you removed. The result will be a planar graph with f=n-1 bounded faces. Since all vertices have degree 3 and summing over the degrees of all vertices gives two times the number e of edges, we obtain 2e = 3v. Substituting this into Eulers formula v-e+f=1 for planar graphs gives the result.)

Consequently: If we test the condition v=2n-4 on diagrams, which do not fall into this very special category, I think this is actually okay. Alternatively, we could test for v <= 2n-2. What do you think?

@niknow Just to clarify, the upper limit should be 2n-4 Voronoi vertices, not the exact number of vertices. Does the exception exceed the upper limit [I can't tell if the diagram with one less vertex is at the 2n-4 limit already, or somewhere below it--so adding one more vertex doesn't seem to rule the upper limit out clearly]? If it exceeds the upper limit, that would be an issue, but if it is <= 2n-4, the unit test still holds as far as I can tell. The textbook also defines the planar version (2n-5) as an upper limit only. If we have a concrete example where 2n-4 is exceeded in a spherical Voronoi diagram then yeah maybe v <= 2n-2 could be used, although maybe not crucial.

If one plays around with the interactive version from https://www.jasondavies.com/maps/voronoi/ , one can see that the claim that a spherical Voronoi diagram always has exactly v=2n-4 vertices is wrong in general

Isn't this what is meant by cooking up a pathological input which is not in (the spherical equivalent of) general position?

We don't really have any checks for the spherical equivalent of general position in the input data. It is not difficult to cook up pathological input that would cause issues.

@argriffing I ended up rebasing on the command line (didn't see any 'button' for this based on the already-rebased branch created above?). Had to add one more commit to clean up the testing module after the rebase for some reason (I had to resolve a few conflicts along the way). I suspect the latter commit shouldn't have been necessary with a perfect conflict resolution workflow using git rebase --continue as I did, but it seems fine now.

After recompiling / installing scipy the full suite of scipy.spatial tests passed. Checking the "Files Changed" tab seems to indicate that the major source of conflicts Thanks.txt was fairly cleanly resolved (i.e., no obliterated contributions there).

didn't see any 'button'

Now I see that this PR does not have its own branch -- the changes are directly to master. Maybe that has something to do with it.

afaik, github ui doesn't have helpers for rebasing things, you have to do it manually (as done above)

Thanks, merged based on above comments.

Hello,

thank you for you scipy.spatial-SphericalVoronoi function. I also see your videos about this. It's really useful for my research project. But i meet some problems and I'd like to ask you.

I see your functions are visualized in matplotlib. I change it to mayavi. There's no big difference but the visualization.

The problem is followed.

1. I tried your example, firstly to get the spherical voronoi vertices, then to produce the polygon and faces. I treat the sv.vertices as verticesList and sv.regions as facesList to get a new mesh, and then uses the half-edge method to see the mesh faces. The sv.regions have conflict with it, they seem not always in clockwise or anticlockwise, so i have to recognize and redefine the faces again. Here is what i get.
   
   I think maybe it's better to produce the sv.regions compatible with half-edge.

2. As you said in the video, there is a problem for the vertices which are located only on one half of the sphere. There exists some wrong spherical voronoi vertices on the sphere. Is there a way to avoid it?
   The simple example is projecting a planar mesh on hyperbolic paraboloid surface, getting the gauss images vertices of each faces' normal vectors, using these spherical vertices to get a spherical voronoi. But the result is weird.
   
   I tried it in rhino, the result seems better
   
   The same problem also appear in similar example that generator vertices in hemisphere. Is there a good way to handle this?

3. Also above example, if the generator vertices are multiple, the algorithm can not work, right?
   The gauss image has multiple vertices, so i have already delete one of them to get the spherical voronoi.

Thank you for anyone who cares about this problem and answer it!

Glad to hear people are finding a use for this code! I agree there's some room for improvement in some areas--I have a few WIP PRs spattered around & I even tried to write a small grant recently for 3D plotting things, but hasn't taken off yet.

I'll try to add a few short answers to queries below, for now. Certainly open to collaborating on improvements.

1. You're suggesting that we should guarantee that the vertices of regions are always in one orientation (i.e., clockwise only)? The docs indicate that sorting can produce either orientation at the moment--maybe consistency would be an improvement if the code isn't too much more complicated

2. Handling pathological generators can be tricky -- the function tries to detect a few things like duplicates & mismatches between generators and the radius of the sphere. I've not found a suitable literature reference for the case of all generators in a single hemisphere -- Question: Spherical Voronoi Diagram for Generators in a single hemisphere #8859 remains unresolved -- I think this is the same issue.

3. Do you mean that you have duplicate generators, or generators that are extremely close to each other? This is definitely something that SphericalVoronoi aims to detect and raise an error for -- have unique points is a pretty basic requirement ("general position") for many algorithms in computational geometry, as you likely know.

If you have useful literature references or you'd like to open a PR that you think can improve behavior that is likely to be welcomed.

Glad to hear people are finding a use for this code! I agree there's some room for improvement in some areas--I have a few WIP PRs spattered around & I even tried to write a small grant recently for 3D plotting things, but hasn't taken off yet.

I'll try to add a few short answers to queries below, for now. Certainly open to collaborating on improvements.

1. You're suggesting that we should guarantee that the vertices of regions are always in one orientation (i.e., clockwise only)? The docs indicate that sorting can produce either orientation at the moment--maybe consistency would be an improvement if the code isn't too much more complicated
2. Handling pathological generators can be tricky -- the function tries to detect a few things like duplicates & mismatches between generators and the radius of the sphere. I've not found a suitable literature reference for the case of all generators in a single hemisphere -- Question: Spherical Voronoi Diagram for Generators in a single hemisphere #8859 remains unresolved -- I think this is the same issue.
3. Do you mean that you have duplicate generators, or generators that are extremely close to each other? This is definitely something that SphericalVoronoi aims to detect and raise an error for -- have unique points is a pretty basic requirement ("general position") for many algorithms in computational geometry, as you likely know.

If you have useful literature references or you'd like to open a PR that you think can improve behavior that is likely to be welcomed.

Thank you very much for your reply. These answers are helpful for me. In my case, it's not hard to orient the directions of each regions. You are definitely right that deduplicating generators is necessary for it. I will think about to solve the hemisphere problem and if i get a good result or read some good references, i will share it here.

Thanks again for your help and your contributions to the spherical voronoi problem!

Hello,

i think maybe projecting the generators(in a hemisphere) to a plane from the south(or north) polar, then use the scipy.spatial.Voronoi to compute planar voronoi, after getting the voronoi vertices and regions(faces), projecting them back to the hemisphere in a same way, then to get a spherical voronoi. In my case, it seems work.

But be aware that the generators should not be multiple and choose the opposite polar of the hemisphere.

The process is just like the following figure in the paper
https://ac.els-cdn.com/S0925772102000779/1-s2.0-S0925772102000779-main.pdf?_tid=ea2fe42a-4150-4cd8-8fb6-a1e0e0447b70&acdnat=1539864627_ff7c6e697b49f202b26ff22e51750220

But i don't quite understand if there is relation with this paper. I didn't find if there are some materials about it and if someone has already proved the feasibility. Hope someone may add such findings here.

Thank you!

Ah yeah, I saw that paper -- couldn't quite get the algorithm working, but I suppose it may indeed circumvent the hemisphere vulnerability. The algorithms are all loglinear time complexity, but I've not seen anyone successfully implement that one in the wild.

Hello,

i'd like to know where can i get the Source code of scipy.spatial.Voronoi? The source code in scipy.org webpage can't be found.

I am working with this function and want to get only one vertex when two vor.vertices are near to each other, like they have same value in decimal 2.

Hoping someone can help me. A little hurry.

Thanks a lot.

@WWmore the source code is in this repo, please use the search function. Voronoi is here: https://github.com/scipy/scipy/blob/master/scipy/spatial/qhull.pyx#L2399

Also, please don't add questions on unrelated PRs - best to use Stackoverflow or the scipy-user mailing list.