Voronoi diagrams

[mo271]

a feature request: Voronoi diagrams, see https://groups.google.com/d/topic/sage-devel/JgqyInSu8S8/discussion

A list of things that could be included in an implementation of Voronoi diagrams could be:

• the graph structure of the Voronoi diagram and it's dual, the Delaunay triangulation
• generalizations
  • Voronoi diagrams with weights (possibly empty regions!)
  • Voronoi diagrams with respect to other metrics (regions not necessary convex Polyhedra anymore)
  • farthest points Voronoi diagrams
• a closest pair method for the list of points/ nearest neighbor
• a plotting routine in 2d and 3d

Perhaps in the file sage/geometry/triangulation/point_configuration.py

CC: @vbraun @jplab @sagetrac-tmonteil @pelegm

Component: geometry

Keywords: voronoi, delaunay, days84

Author: Moritz Firsching

Branch/Commit: 5dc6bd5

Reviewer: Jean-Philippe Labbé

Issue created by migration from https://trac.sagemath.org/ticket/13517

[mo271]

a first try

[mo271]

comment:1

Attachment: trac_#13517_voronoi_first_try.patch.gz

I added a patch just to see if I understood correctly how patching and the mercurial export etc works..

[vbraun]

comment:2

The patch is correct ;-)

You should format the examples in the docstring as in http://www.sagemath.org/doc/developer/conventions.html#docstring-markup-with-rest-and-sphinx, then they become tests that you can run with sage -t sage/geometry/voroni_diagram.py.

If you want to lay the path for more functionality you should create a class that holds the initial data (the points) and caches the polyhedral cells. Then other methods, for example a plot() method, can be easily added without having to pass lists of polyhedra around.

Also you don't have to hard-code floating-point doubles as the ring, you could allow exact rationals at the same time. This requires you to figure out what a good ring for the input points is, of course. You could just put them into a PointConfiguration and use

sage: pc = PointConfiguration([(1,2), (1,2/3)])
sage: pc.base_ring()
Rational Field
sage: pc = PointConfiguration([(1,2), (1,0.3)])
sage: pc.base_ring()
Real Field with 53 bits of precision

[mo271]

comment:3

Thanks for the hints, I will get cracking on this.

Replying to @vbraun:

The patch is correct ;-)

You should format the examples in the docstring as in http://www.sagemath.org/doc/developer/conventions.html#docstring-markup-with-rest-and-sphinx, then they become tests that you can run with sage -t sage/geometry/voroni_diagram.py.

If you want to lay the path for more functionality you should create a class that holds the initial data (the points) and caches the polyhedral cells. Then other methods, for example a plot() method, can be easily added without having to pass lists of polyhedra around.

Also you don't have to hard-code floating-point doubles as the ring, you could allow exact rationals at the same time. This requires you to figure out what a good ring for the input points is, of course. You could just put them into a PointConfiguration and use

sage: pc = PointConfiguration([(1,2), (1,2/3)])
sage: pc.base_ring()
Rational Field
sage: pc = PointConfiguration([(1,2), (1,0.3)])
sage: pc.base_ring()
Real Field with 53 bits of precision

[mo271]

Attachment: trac_#13517_new_class_initial_version.patch.gz

first working version

[mo271]

comment:5

I added a first version, which provides the basic functionality

[miguelmarco]

comment:6

Thanks for the work, i think it is a nice addition.

A couple of comments though.

-I don't understand why your patch adds the file voronoi_diagram.py.out. Is that made on purpose or is it a mistake?

-Unless i am missing something, i think that allowing exact fields (like rationals, or maybe even algebraic reals) should work just the same way. There is really no need to hard code RDF, it would be better to use the same field as the given entry. This should be a really trivial change: just define a self._field as self._points.base_ring(), and then replace RDF by self._field in the rest of your code.

[mo271]

reviewed

[mo271]

comment:7

Attachment: trac_#13517_new_class_reviewed.patch.gz

thanks for looking at the file.

• the file voronoi_diagram.py.out was indeed a mistake, the new patch does not include that file.

• I introduced a self._base_ring variable to handle QQ and RDF depending on input, as you suggested. (Now if the base_ring of the PointConfiguration is a subring of the rationals, it uses QQ, otherwise RDF)

[miguelmarco]

comment:9

Would it be a lot of trouble to introduce also the posibility of using algebraic extensions of QQ as base ring? Arithmetic in those cases may be slow, but i can imagine some situations where it would be important to know the vertices of the diagram in an exact form.

[mo271]

comment:10

Replying to @miguelmarco:

Would it be a lot of trouble to introduce also the possibility of using algebraic extensions of QQ as base ring? Arithmetic in those cases may be slow, but i can imagine some situations where it would be important to know the vertices of the diagram in an exact form.

I don't see how one could easily adapt the algorithm for algebraic extensions of QQ for the following reason:
The Polyhedron-method only works with QQ and RDF.

If you really want support for extensions of QQ, either the algorithm of finding the Voronoi-diagram has to be changed, or one has adapt the Polyhedron method. Both are nontrivial tasks, which might make a good upgrade later on.

[kcrisman]

comment:11

See also

• https://sage.math.clemson.edu:34567/home/pub/388/ which seems to have an implementation of Voronoi diagrams
• http://pypi.python.org/pypi/voropy/0.0.1
• http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html
• We include qhull somewhere, which is also part of an optional package in R
• Some information at various posts at this blog, though I don't know how detailed they are
  And so forth. And Delaunay triangulations exist in matplotlib too, I think... would it be faster to use that, especially since they fixed the problem we had in Clarify some comments concerning the matplotlib Delaunay code in list_plot3d.py #13167?

[nathanncohen]

comment:12

Hellooooooooooooo !!

Just thought I would add a link toward d3.js here. It's a javascript library that can be used to plot a lot of things, and there's an example of a Voronoi Diagram on its website :

http://bl.ocks.org/mbostock/4060366

I wrote a patch (#14953) that uses this library to draw graphs :-)

Nathann

[jplab]

comment:67

There is going to be a conflict with #22496. I guess rebasing will just do the thing. Let's wait some time that the "rc" versions pass...

[jplab]

comment:68

Dear Moritz,

I guess you can now rebase and once the bot gives the green light, I'll give positive review again.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

ef78175	Voronoi diagram merge
7a7905d	added to new reference index

[sagetrac-git]

Changed commit from 90477d4 to 7a7905d

[mo271]

comment:71

Ok, just rebasing didn't qute work, I did another merge (and squashed). In any case I hope the bot doesn't complain now.

I added it to the misc section in the discrete geometry reference.

Let's see what the bot thinks.

[jplab]

comment:73

Hi Moritz,

There is a doctest failure:

File "voronoi_diagram.py", line 197, in sage.geometry.voronoi_diagram.VoronoiDiagram.regions
Failed example:
    V = VoronoiDiagram([[1, 3, .3], [2, -2, 1], [-1, 2, -.1]]); V.regions()
Expected:
    {P(1.00000000000000, 3.00000000000000, 0.300000000000000): A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 1 vertex, 2 rays, 1 line,
    P(2.00000000000000, -2.00000000000000, 1.00000000000000): A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 1 vertex, 2 rays, 1 line,
    P(-1.00000000000000, 2.00000000000000, -0.100000000000000): A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 1 vertex, 2 rays, 1 line}
Got:
    {P(1.00000000000000, 3.00000000000000, 0.300000000000000): A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 1 vertex, 2 rays, 1 line,
     P(-1.00000000000000, 2.00000000000000, -0.100000000000000): A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 1 vertex, 2 rays, 1 line,
     P(2.00000000000000, -2.00000000000000, 1.00000000000000): A 3-dimensional polyhedron in RDF^3 defined as the convex hull of 1 vertex, 2 rays, 1 line}
**********************************************************************

I guess the output should fix the order. Or change the doctest.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

fc3c49a

make doctest independent of order in dict

[sagetrac-git]

Changed commit from 7a7905d to fc3c49a

[tscrim]

comment:75

It is better if you use:

            sage: V.regions() == {P[0]: Polyhedron(base_ring=RDF, lines=[(-RDF(0.375), RDF(0.13888888890000001), RDF(1.5277777779999999))],
            ....:                                                 rays=[(RDF(9), -RDF(1), -RDF(20)), (RDF(4.5), RDF(1), -RDF(25))],
            ....:                                                 vertices=[(-RDF(1.1074999999999999), RDF(1.149444444), RDF(9.0138888890000004))]),
            ....:                 P[1]: Polyhedron(base_ring=RDF, lines=[(-RDF(0.375), RDF(0.13888888890000001), RDF(1.5277777779999999))],
            ....:                                                 rays=[(RDF(9), -RDF(1), -RDF(20)), (-RDF(2.25), -RDF(1), RDF(2.5))],
            ....:                                                 vertices=[(-RDF(1.1074999999999999), RDF(1.149444444), RDF(9.0138888890000004))]),
            ....:                 P[2]: Polyhedron(base_ring=RDF, lines=[(-RDF(0.375), RDF(0.13888888890000001), RDF(1.5277777779999999))],
            ....:                                                 rays=[(RDF(4.5), RDF(1), -RDF(25)), (-RDF(2.25), -RDF(1), RDF(2.5))],
            ....:                                                 vertices=[(-RDF(1.1074999999999999), RDF(1.149444444), RDF(9.0138888890000004))])}

[sagetrac-git]

Changed commit from fc3c49a to 5dc6bd5

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

5dc6bd5

nicer linebreaks

[mo271]

comment:77

Thanks for the tip, Travis! I changed it accordingly.

I hope all is well now.

[mo271]

comment:79

Dear JP, I will set it on positive review on your behalf, as soon as the bot does not complain about the white space changes

[vbraun]

Changed branch from public/13517 to 5dc6bd5