Fast, robust construction of 2D Delaunay and Voronoi tessellations on generic point types. Implementation follows algorithms described in the Arepo paper and used (for e.g.) in the Illustris Simulation. License: MIT. Bug reports welcome!
Incrementally insert points to a valid Delaunay tessallation, while restoring Delaunayhood by flipping triangles. Point location (i.e. which triangle should it divide into three) is accelerated by spatial sorting. Spatial sorting allows to add points which are close in space thus walking the tesselation is fast. Initial tessalletion includes two triangles built by 4 points which are outside of the allowed region for users. These "external" triangles are skipped when iterating over Delaunay/Voronoy edges. Fast and robust predicates are provided by the GeometricalPredicates package. Benchmarks suggest this package is a bit faster than CGAL, see here benchmark of an older version which is also a bit slower than current.
min_coord <= x <= max_coordwhere
max_coord=2.0-2eps(Float64). Note this is a bit different than what is required by the
Define and push individual points like this:
using VoronoiDelaunay tess = DelaunayTessellation() push!(tess, Point(1.5, 1.5))
creation of points is explained in the GeometricalPredicates package documentation.
Pushing arrays of points is more efficient:
width = max_coord - min_coord a= Point2D[Point(min_coord+rand()*width, min_coord+rand()*width) for i in 1:100] push!(tess, a)
notice care taken for correct range of coordinates.
max_coord are defined in the package. We can further optimize by giving a
sizehint at time of construction:
tess = DelaunayTessellation(100)
or at any later point:
Delaunay tesselations need at least 3 points to be well defined. Voronoi need 4. Remember this when iterating or plotting. Iterating over Delaunay edges is done like this:
i = 0 for edge in delaunayedges(tess) i += 1 # or, do something more useful :) end
DelaunayEdge contains two points a and b, they can be accesse with
Iterating over Voronoi edges is similar:
i = 0 for edge in voronoiedges(tess) i += 1 # or, do something more useful :) end
VoronoiEdge is a bit different than a
Point2D and not the generators, as they have different coordinates. To get the generators use
getgenb(edge) these give the relevant
AbstractPoint2D which were used to create the edge.
If the generators are not needed when iterating over the Voronoi edges (e.g. when plotting) then a more efficient way to iterate is:
i = 0 e=Nothing for edge in voronoiedgeswithoutgenerators(tess) i += 1 # do something more useful here :) end
edge is a
VoronoiEdgeWithoutGenerators, the points
b can be accessed as usual.
Iterating over Delaunay triangles:
i = 0 for delaunaytriangle in tess i += 1 # or, do something more useful :) end
delaunaytriangle here is of type
DelaunayTriagle which is s subtype of
AbstractNegativelyOrientedTriangle. To get the generators of this triangle use the
getc methods. You can do all other operations and predicate tests on this triangle as explained in GeometricalPredicates
Locating a point, i.e. finding the triangle it is inside:
t = locate(tess, Point(1.2, 1.3))
if the point is outside of the tessellation then
isexternal(t) == true holds. This is good for type stability, at least better than returning a
Nothing. It is assumed that the point we want to locate is actually in the allowed points region. Performance is best when locating points close to each other (this is also why spatial sorting is used). Future versions may implement a hierarchal approach for fast random locations.
Navigating from a triangle to its neighbours is done like this:
t = movea(tess, t) # move to the direction infront of generator a t = moveb(tess, t) # move to the direction infront of generator b t = movec(tess, t) # move to the direction infront of generator c
The following retrieves a couple of vectors ready to plot Voronoi edges:
x, y = getplotxy(voronoiedges(tess))
and for Delaunay edges:
x, y = getplotxy(delaunayedges(tess))
Now plotting can be done with your favorite plotting package, for e.g.:
using Gadfly plot(x=x, y=y, Geom.path)
To make a nice looking plot remember to limit the axes and aspect ratio. For e.g.:
set_default_plot_size(15cm, 15cm) plot(x=x, y=y, Geom.path, Scale.x_continuous(minvalue=1.0, maxvalue=2.0), Scale.y_continuous(minvalue=1.0, maxvalue=2.0))
You can create a tesselation from an image, just like the tesselation of the
julia logo at the top of this README. This was created from a png with
import Images: imread img = imread("julia.png") tess = from_image(img, 25000)
23 days ago