10/17/2014

1 day ago

78 commits

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.

- Due to numerical restrictions the point coordinates must be within
`min_coord <= x <= max_coord`

where`min_coord=1.0+eps(Float64)`

and`max_coord=2.0-2eps(Float64)`

. Note this is a bit different than what is required by the`GeometricalPredicates`

package. - The following features are not implemented, but are in the TODO list; In order of priority: centroid tessellations (Lloy's method), Weighted generators (both power and sum), bounding, maybe restricting. Hierarchal tessellations for fast random locatings; Distributed tessellation construction. 3D. Order of priority may change of course :)

```
Pkg.add("VoronoiDelaunay")
```

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. `min_coord`

and `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:

```
sizehint(tess, 100)
```

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
```

a `DelaunayEdge`

contains two points a and b, they can be accesse with `geta(edge)`

and `getb(edge)`

.
Iterating over Voronoi edges is similar:

```
i = 0
for edge in voronoiedges(tess)
i += 1
# or, do something more useful :)
end
```

a `VoronoiEdge`

is a bit different than a `DelaunayEdge`

: here `a`

and `b`

are `Point2D`

and not the generators, as they have different coordinates. To get the generators use `getgena(edge)`

and `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
```

here `edge`

is a `VoronoiEdgeWithoutGenerators`

, the points `a`

and `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 `geta`

, `getb`

, and `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 `from_file`

(see `examples/img_to_vorono.jl`

):

```
import Images: imread
img = imread("julia.png")
tess = from_image(img, 25000)
```