Handling Spherical Polygons

[tylerjereddy]

It would be really nice if matplotlib could handle the plotting of spherical polygons. I (and several other volunteers) have been working to integrate spherical Voronoi diagram code into scipy (pull request: scipy/scipy#5232), but actually plotting the Voronoi regions on the surface of a sphere is problematic:

Possible applications include: visualization for search and rescue operations on the surface of the earth, global network coverage, epidemiology, airspace territory assignment (i.e., https://www.jasondavies.com/maps/voronoi/airports/), computational virology (my research), and many other sphere-based analyses.

I am not sure how difficult this would be and if it can be done in a 'lightweight' fashion, or if this really falls into the domain of mayavi and similar 3D-specialized packages.

[mdboom]

Astronomy could really use such a tool as well.

Just to confirm the problem. The problem is that the edges of the polygons are being plotted as straight line segments and not spherical arcs. Is that correct?

I wonder if there is a standard transformation/approximation from spherical arcs to bezier curves in the 2D projected plane. That would obviate the need for doing the interpolation on our side -- we could just let the renderer/vector file format handle it.

[WeatherGod]

It would seem that you have already implemented it in some rudimentary form
or another given your attached plot. I am always willing to incorporate new
plotting methods that makes such plots easier. For plotting them more
accurately with proper curves, that will be trickier because mplot3d does
not currently utilize the transform stack. I have made a few attempts at
converting it, but my python-fu is not strong...

On Wed, Oct 21, 2015 at 11:36 AM, Tyler Reddy notifications@github.com
wrote:

It would be really nice if matplotlib could handle the plotting of
spherical polygons. I (and several other volunteers) have been working to
integrate spherical Voronoi diagram code into scipy (pull request:
scipy/scipy#5232 scipy/scipy#5232), but
actually plotting the Voronoi regions on the surface of a sphere is
problematic:

[image: sample_voronoi]
https://cloud.githubusercontent.com/assets/7903078/10641479/55470a0c-7811-11e5-824f-4f7ea54c8469.png

Possible applications include: visualization for search and rescue
operations on the surface of the earth, global network coverage,
epidemiology, airspace territory assignment (i.e.,
https://www.jasondavies.com/maps/voronoi/airports/), computational
virology (my research), and many other sphere-based analyses.

I am not sure how difficult this would be and if it can be done in a
'lightweight' fashion, or if this really falls into the domain of mayavi
and similar 3D-specialized packages.

—
Reply to this email directly or view it on GitHub
#5294.

[tylerjereddy]

@mdboom Yes, so something like this:

polygon = Poly3DCollection([voronoi_region],alpha=1.0)

will plot a planar polygon rather than one that follows the spherical arc. One would want a way for the surface curvature to be more reasonably approximated in the plot given the ordered vertices of the spherical polygon in an input array for example.

[WeatherGod]

The problem is going to be a bit trickier than projecting the path of the
polygon boundary. Imagine a polygon that represents a hemisphere, and
rotate it such that half of it is out of view. Projecting the boundaries of
the polygon isn't enough because that would never account for the "buldge"
of the face of the polygon. The projection has to be all the way down to
the "view" space, I think. I am completely unaware of the math needed to
accomplish that.

On Wed, Oct 21, 2015 at 11:53 AM, Tyler Reddy notifications@github.com
wrote:

@mdboom https://github.com/mdboom Yes, so something like this:

polygon = Poly3DCollection([voronoi_region],alpha=1.0)

will plot a planar polygon rather than one that follows the spherical arc.
One would want a way for the surface curvature to be more reasonably
approximated in the plot given the ordered vertices of the spherical
polygon in an input array for example.

—
Reply to this email directly or view it on GitHub
#5294 (comment)
.

[WeatherGod]

Note that just because I am unaware of the math needed to accomplish it
doesn't mean that it doesn't exist or that it isn't trivial. It just means
that 3D projections are not my domain.

On Wed, Oct 21, 2015 at 11:58 AM, Benjamin Root ben.v.root@gmail.com
wrote:

The problem is going to be a bit trickier than projecting the path of the
polygon boundary. Imagine a polygon that represents a hemisphere, and
rotate it such that half of it is out of view. Projecting the boundaries of
the polygon isn't enough because that would never account for the "buldge"
of the face of the polygon. The projection has to be all the way down to
the "view" space, I think. I am completely unaware of the math needed to
accomplish that.

On Wed, Oct 21, 2015 at 11:53 AM, Tyler Reddy notifications@github.com
wrote:

@mdboom https://github.com/mdboom Yes, so something like this:

polygon = Poly3DCollection([voronoi_region],alpha=1.0)

will plot a planar polygon rather than one that follows the spherical
arc. One would want a way for the surface curvature to be more reasonably
approximated in the plot given the ordered vertices of the spherical
polygon in an input array for example.

—
Reply to this email directly or view it on GitHub
#5294 (comment)
.

[tylerjereddy]

Yeah, I just thought I'd share the idea and see what people thought. I know the differential geometer I've been working with (Nikolai) might have a few ideas, so maybe he will chime in at some point.

[tylerjereddy]

There does seem to be a very mathematical / detailed discussion about doing this with Mathematica here: http://mathematica.stackexchange.com/questions/78705/plot-a-partition-of-the-sphere-given-vertices-of-polygons

I'll try to take a closer look at that discussion soon.

[niknow]

Hi,
I'm one of the colleagues @tylerjereddy mentioned. The reason why our plot of spherical Voronoi regions is not really good indeed is that Poly3DCollection plots a Euclidean 2D polygon in Euclidean 3D. If your Voronoi diagram has a lot of generators (i.e. if the polygons are relatively small compared to the sphere) this looks already fairly nice (see Tylers picture above). But if the polygons get bigger, the result looks less and less spherical (see the ''Example in Docstring'' picture from scipy/scipy#5232))

I already tried out various solutions to fix this problem, which I just uploaded to

https://github.com/niknow/plot-spherical-polygon/blob/master/plot_revision.py

but I'm not happy with any of these, which is why I did'nt start a pull request yet.

My general idea was: If the polygon is so big such that you can no longer ignore the difference between a Euclidean 2D plot in 3D and a true spherical plot, just split up the polygon into smaller ones. This can be done by calculating convex combinations of the vertices (see plotdata-function) and projecting them onto the sphere. But this is a very crude attempt and the resulting plot is very slow. Judging from

http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html

shouldn't it be possible to do this using plot_trisurf? As you can see from my examples in plot_revision.py I tried that out, but there seems to be a problem with the ordering of the triangles in the triangulation (probably resulting from the fact that trisurf does'nt really understand that we want a spherical triangulation). Also, the mesh grid might look more ''spherical'' if one first projects the polygon to 2D using polar coordinates, generate the convex combinations in 2D and then project back. (One has to be careful though, since the polygon might be cut in half by the polar coordinates.)
I was busy with some other papers recently, but I'll think about it again in more detail as soon as I have more time. If anybody has a good idea how to continue this, I'm of course very happy to discuss it.

[mdboom]

@tylerjereddy: Do you have any code you could point me to to produce the above plot? I don't care so much about the triangulation itself as much as the output data right before you pass it to matplotlib. I have some ideas I'd like to play with, and it would be helpful to hit the ground running with some real data.

[tylerjereddy]

@mdboom I've produced a 'minimal' Jupyter notebook (https://github.com/tylerjereddy/matplotlib_spherical_polygon_example/blob/master/spherical_polygon_matplotlib.ipynb) that should remove any of the confusing Voronoi stuff and just reduce it down to spherical polygons and matplotlib.

As you can see, I did that by pickling the spherical polygon coordinate arrays into a small file that can be cloned from that repo along with the notebook. Hope that helps?

Let me know if there's something I can do to help.

[mdboom]

Thanks. That's exactly what I had in mind.

[tylerjereddy]

@mdboom It seems that this may not be too hard to do in Python based on my discussion with @rougier (core developer of glumpy): https://twitter.com/GlumpyProject/status/699247950896226304

He leveraged the Python tool written by @prideout here: http://prideout.net/blog/?p=44

[rougier]

The solution is a simple subdivision of the patch. The code is:

from euclid import *

def subdivide(verts, faces):
    """Subdivide each triangle into four triangles, pushing verts to the unit sphere"""
    triangles = len(faces)
    for faceIndex in xrange(triangles):

        # Create three new verts at the midpoints of each edge:
        face = faces[faceIndex]
        a,b,c = (Vector3(*verts[vertIndex]) for vertIndex in face)
        verts.append((a + b).normalized()[:])
        verts.append((b + c).normalized()[:])
        verts.append((a + c).normalized()[:])

        # Split the current triangle into four smaller triangles:
        i = len(verts) - 3
        j, k = i+1, i+2
        faces.append((i, j, k))
        faces.append((face[0], i, k))
        faces.append((i, face[1], j))
        faces[faceIndex] = (k, j, face[2])

    return verts, faces

The "trick" is the normalized call that push new vertices onto the unit sphere (verts are supposed to be on the unit sphere). Several calls may be needed depending on the original number of faces.

[pwolfram]

👍

[niknow]

Thank you all very much for your input. I've tried to wrap up this discussion in pull request #6248. I'm happy for any feedback in the discussion thread of the PR.

[tylerjereddy]

I quite like the idea of using spherical linear interpolation (Slerp), which as the article notes is quite common in 3D graphics. I wrote a Cython version of slerp in a WIP scipy PR, and for one pathological input scenario you can get a much nicer result by interpolating between vertices in this manner:

Before:

After (with a bit of azimuthal rotation as well -- see note below):

One advantage here might be that we could just intercept the vertices, interpolate between them, and then just allow the regular matplotlib machinery to do its thing with more vertices. Would perhaps just need an appropriate keyword for slerp_count when initializing i.e., Poly3DCollection -- that also leaves the performance and quality of the interpolation up to the end user.

I started working on a PR (ok, just for a few minutes), when I noticed there didn't seem to be Cython code in matplotlib (!), and doing the interpolation iteration in Cython is quite a useful speedup, especially if you wanted thousands of additional points between many vertices.

One issue that perhaps all approaches will have though is the one raised by @WeatherGod -- take a look at the sensitivity of "perceived quality" of the polygon to the azimuthal angle (much uglier here):

Here's the crude code I used to intercept the vertices of the original spherical polygon and slerp them into many interpolated vertices:

100 # try applying spherical linear interpolation to improve plot
102 N = dual_poly.shape[0]
103 n_int = 900
104 interpolated_polygon = np.zeros((N * n_int, 3),dtype=np.float64)
105 
106 counter = 0
107 for i in range(N):
108     if i == (N-1):
109         next_index = 0
110     else:
111         next_index = i + 1
112 
113     interpolated_polygon[counter:(counter + n_int), ...] = slerp(dual_poly[i],
114                                                                dual_poly[next_index],
115                                                                n_int)
116     counter += n_int
117     
118 
119 
120 # plot the suspected Dual Pole containing
121 # polygon for visual inspection
122 polygon = Poly3DCollection([interpolated_polygon], alpha=1.0)
123 ax.add_collection3d(polygon)
124 
125 ax.scatter(dual_poly[...,0],
126            dual_poly[...,1],
127            dual_poly[...,2],                      
128            c='k',                                 
129            s=20)

[tylerjereddy]

Of course, there would be other things to consider too -- a polygon with consecutive antipodal vertices has an infinite number of possible great circle arcs so that would have to be rejected, and there's also I suppose a need to specify the radius of the sphere.

[WeatherGod]

This is really neat. Again, there are limits to what mplot3d can do since it isn't a real 3d renderer, the example above is an interesting form of the Escher effect that I have seen hints of before in old bug reports, but I have never been able to figure out.before. But with your example, it totally makes sense now. When the 3D vertices get projected into 2D space, a knot forms with the 2d line string, forming a bow-tie appearance. For the moment, we won't be including cython code into the codebase. Now, if there is anything we can do to make it easier for an external package to provide nice-looking curved surfaces, let us know. Perhaps it might make sense to have something better within mplot3d one day for general polygon surfaces.

…

On Thu, Feb 22, 2018 at 8:18 PM, Tyler Reddy ***@***.***> wrote: Of course, there would be other things to consider too -- a polygon with consecutive antipodal vertices has an infinite number of possible great circle arcs so that would have to be rejected, and there's also I suppose a need to specify the radius of the sphere. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#5294 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AARy-G0f_hOuP8i4BOV7-jeAwy9YK9-Rks5tXhHJgaJpZM4GTE2P> .

[pwolfram]

I'm really glad to see the discussion and progress here-- this will be cool if it can become a standard part of matplotlib!

[github-actions]

This issue has been marked "inactive" because it has been 365 days since the last comment. If this issue is still present in recent Matplotlib releases, or the feature request is still wanted, please leave a comment and this label will be removed. If there are no updates in another 30 days, this issue will be automatically closed, but you are free to re-open or create a new issue if needed. We value issue reports, and this procedure is meant to help us resurface and prioritize issues that have not been addressed yet, not make them disappear. Thanks for your help!