[ANN] DelaunayTriangulation.jl v0.5: Constrained triangulation and mesh refinement

Hi all,

I have just released v0.5 of my package DelaunayTriangulation.jl, a package for computing Delaunay triangulations of planar point sets. This is probably the last version before v1.0, I don’t expect much of the interface to change at all.

This is a pretty big change, as with it comes the possibility of generating constrained Delaunay triangulations and a method for mesh refinement, compared to v0.1 which only had the ability to compute unconstrained Delaunay triangulations. It’s especially nice as it avoids the need to have rely so heavily on Gmsh just to have constrained triangulations. On the way to v0.5, many features have become available as I have needed them for the main triangulation algorithms, some of them being:

As a demonstration of just how much the package can do, I show below an example of triangulating and refining a multipolygon defining the Julia logo. (The definition of the points has been truncated, see the mesh refinement section of the docs for the complete definition of the points.) This definition just shows boundary nodes, but constrained segments could also be provided with the edges keyword argument of triangulate (or afterwards with add_edge!; once again, see the docs). Much simpler examples are given in the docs - I hope they are sufficiently exhaustive.

using DelaunayTriangulation, CairoMakie 

## Defining the points
C = (15.7109521325776, 33.244486807457)
D = (14.2705719699703, 32.8530791545746)
E = (14.3, 27.2)
--- truncated input ---
S5 = (26.5782996794556, 29.7503874690082)
T5 = (26.7603175886393, 30.3384453294476)
U5 = (27.3203726938197, 30.7024811478149)

## Defining the point set and the boundary segments
J_curve = [[C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, C]]
U_curve = [[T, U, V, W, Z, A1, B1, C1, D1, E1, F1, G1, H1, I1, J1, K1, L1, M1, N1, O1, T]]
L_curve = [[P1, Q1, R1, S1, P1]]
I_curve = [[T1, U1, V1, W1, T1]]
A_curve_outline = [[
    K5, W3, Z3, A4, B4, C4, D4, E4, F4, G4, H4, I4, J4, K4, L4, M4, N4,
    O4, P4, Q4, R4, S4, T4, U4, V4, W4, Z4, A5, B5, C5, D5, E5, F5, G5,
    H5, I5, J5, K5]]
A_curve_hole = [[L5, M5, N5, O5, P5, Q5, R5, S5, T5, U5, L5]]
dot_1 = [[Z1, A2, B2, C2, D2, E2, F2, G2, H2, I2, J2, Z1]]
dot_2 = [[Z2, A3, B3, C3, D3, E3, F3, G3, H3, I3, J3, Z2]]
dot_3 = [[K2, L2, M2, N2, O2, P2, Q2, R2, S2, T2, U2, V2, W2, K2]]
dot_4 = [[K3, L3, M3, N3, O3, P3, Q3, R3, S3, T3, U3, V3, K3]]
curves = [J_curve, U_curve, L_curve, I_curve, A_curve_outline, A_curve_hole, dot_1, dot_2, dot_3, dot_4]
nodes, points = convert_boundary_points_to_indices(curves)

## Now triangulate
tri = triangulate(points; boundary_nodes=nodes, check_arguments=false) # about 0.002s
orig_tri = deepcopy(tri) # for plotting later

## Refine
A = get_total_area(tri)
stats = refine!(tri; min_angle=26.45, max_area=0.005A / 9) # about 0.02

## Show some statistics for a demonstration 
fig = Figure(fontsize=33)
ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", title=L"(a):$ $ Original", width=400, height=400, titlealign=:left)
triplot!(ax, orig_tri, show_convex_hull=false) # orig_tri was deepcopy(triangulate(pts)) from before 
ax = Axis(fig[1, 2], xlabel=L"x", ylabel=L"y", title=L"(b):$ $ Refined", width=400, height=400, titlealign=:left)
triplot!(ax, tri, show_convex_hull=false)
areas = get_all_stat(stats, :area) ./ (0.005A)
angles = getindex.(get_all_stat(stats, :angles), 1) # 1 is the smallest
ax = Axis(fig[2, 1], xlabel=L"A/A(\Omega)", ylabel=L"$ $Count", title=L"(c):$ $ Area histogram", width=400, height=400, titlealign=:left)
hist!(ax, areas, bins=0:0.001:0.1)
ax = Axis(fig[2, 2], xlabel=L"\theta_{\min}", ylabel=L"$ $Count", title=L"(d):$ $ Angle histogram", width=400, height=400, titlealign=:left)
hist!(ax, rad2deg.(angles), bins=0:0.2:40)
vlines!(ax, [26.45], color=:red)

Also, as a much more practical example, here’s an example of triangulating a boundary of Tasmania; the tassy.txt file comes from a download image of Tasmania that I then traced using ImageJ.

using DelimitedFiles
tassy = readdlm("tassy.txt")
ymax = @views maximum(tassy[:, 2])
tassy = [(x, ymax - y) for (x, y) in eachrow(tassy)]
push!(tassy, tassy[begin])
boundary_nodes, points = convert_boundary_points_to_indices(tassy)
tri = triangulate(points; boundary_nodes=boundary_nodes)
A = get_total_area(tri)
stats = refine!(tri; max_area=1e-3A)


Useful package, and very well documented, thanks! I will find some non-serious uses for it… :slight_smile:


In addition to the main post, here’s a nice demonstration of the incremental nature of the algorithm.


I decided that before v1.0 I would like to have an initial version of a Voronoi tessellation. So, v0.6 that I’ve just released now also includes Voronoi tessellations, with support for clipping to a convex hull and smoothing via Lloyd’s algorithm (i.e. computation of a centroidal Voronoi tessellation). I was hoping to be able to clip to any geometry, but I couldn’t seem to figure out some of the complicated details - my implementation is quite complicated, likely more than it needs to be. Anyway!

Here are some examples of the tessellation in action. The algorithms all start from an initial Delaunay triangulation; no implementation like Fortune’s algorithm is given.

using DelaunayTriangulation, CairoMakie, StableRNGs
rng = StableRNG(999888771)
pts = 25randn(rng, 2, 500)
tri = triangulate(pts; rng)
vorn = voronoi(tri)
clip_vorn = voronoi(tri, true)
smooth_vorn = centroidal_smooth(clip_vorn; rng)

cmap = Makie.cgrad(:jet)
colors = get_polygon_colors(vorn, cmap)
fig = Figure(fontsize=38, resolution=(1540, 485))
for (j, vor) in enumerate((vorn, clip_vorn, smooth_vorn))
    ax = Axis(fig[1, j], width=400, height=400)
    voronoiplot!(ax, vor, strokecolor=:red, strokewidth=0.2, polygon_color=colors, markersize=4)
    xlims!(ax, -100, 100)
    ylims!(ax, -100, 100)

And a nice animation to finish off.