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:
- Geometric predicates are implemented with ExactPredicates.jl, and many predicates have been extended from ExactPredicates.jl.
- Unconstrained and constrained triangulations. Support is provided for many types of domains, as given in the docs.
- Mesh refinement, with support for custom angle and area constraints.
- Dynamic point insertion, point deletion, and segment insertion, amongst many other operations.
- Computation of convex hulls, either from the triangulation itself or using a Graham scan.
- Triangulation of convex polygons.
- Efficient point location on convex geometries, even with interior holes. Partial support exists for non-convex geometries, although it is much slower and not perfect.
- Computation of the pole of inaccessibility, i.e. the point in a polygon that is furthest from a boundary (see e.g. this blogpost).
- Fully customisable interface for defining geometric primitives.
- Simple iteration over mesh elements, including points, edges, or triangles.
- Computation of statistics over individual triangular elements and over a complete triangulation.
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)
resize_to_layout!(fig)
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)]
reverse!(tassy)
unique!(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)
