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

**URL:** <https://discourse.julialang.org/t/ann-delaunaytriangulation-jl-v0-5-constrained-triangulation-and-mesh-refinement/97508>\
**Category:** Package Announcements\
**Tags:** package, announcement\
**Created:** [April 15, 2023, 8:26am UTC](https://discourse.julialang.org/t/ann-delaunaytriangulation-jl-v0-5-constrained-triangulation-and-mesh-refinement/97508 "2023-04-15T08:26:57Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [April 15, 2023, 8:26am UTC](https://discourse.julialang.org/t/ann-delaunaytriangulation-jl-v0-5-constrained-triangulation-and-mesh-refinement/97508/1 "2023-04-15T08:26:57Z")

</div>

Hi all,

I have just released v0.5 of my package [DelaunayTriangulation.jl](https://github.com/DanielVandH/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](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/predicates/) are implemented with [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl), and many predicates have been extended from ExactPredicates.jl.
- [Unconstrained](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/triangulations/unconstrained/) and [constrained](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/triangulations/constrained/) triangulations. Support is provided for many types of domains, as given in the docs.
- [Mesh refinement](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/triangulations/refinement/), with support for custom angle and area constraints.
- Dynamic point insertion, point deletion, and segment insertion, amongst many other [operations](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/operations/).
- Computation of convex hulls, either [from the triangulation itself](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/data_structures/convex_hull/) or using a [Graham scan](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/other_features/convex_hull/).
- [Triangulation of convex polygons](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/tri_algs/convex/).
- [Efficient point location](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/other_features/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](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/other_features/pole_of_inaccessibility/), i.e. the point in a polygon that is furthest from a boundary (see e.g. [this blogpost](https://blog.mapbox.com/a-new-algorithm-for-finding-a-visual-center-of-a-polygon-7c77e6492fbc)).
- [Fully customisable interface](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/interface/interface/) for defining geometric primitives.
- [Simple iteration over mesh elements, including points, edges, or triangles](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/data_structures/triangulation/).
- Computation of [statistics](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/data_structures/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](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/triangulations/refinement/) 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.

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

```

 ![julia_constrained_refinement](https://global.discourse-cdn.com/julialang/original/3X/3/f/3fb7f4cfd3044beacb68f1b86c720ecee55bb7b6.png)

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.

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

```

 ![tassy_constrained_refinement](https://global.discourse-cdn.com/julialang/original/3X/e/6/e6202b17c4cfcdc387a236a02a09595732048143.png)

---

<div class="post-metadata">

**Author:** ![cormullion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cormullion/32/49131_2.png) [@cormullion](https://discourse.julialang.org/u/cormullion)\
**Post date:** [April 15, 2023, 11:07am UTC](https://discourse.julialang.org/t/ann-delaunaytriangulation-jl-v0-5-constrained-triangulation-and-mesh-refinement/97508/2 "2023-04-15T11:07:27Z")

</div>

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

 ![Screenshot 2023-04-15 at 12.05.16](https://global.discourse-cdn.com/julialang/original/3X/e/d/ed4344a15724751459bc5cf73604b4fe432abf98.png)

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [April 15, 2023, 11:08am UTC](https://discourse.julialang.org/t/ann-delaunaytriangulation-jl-v0-5-constrained-triangulation-and-mesh-refinement/97508/3 "2023-04-15T11:08:37Z")

</div>

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

[https://user-images.githubusercontent.com/95613936/232210266-615e08bd-d4f2-4a40-849e-e832778a4b71.mp4](https://user-images.githubusercontent.com/95613936/232210266-615e08bd-d4f2-4a40-849e-e832778a4b71.mp4)

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [April 24, 2023, 6:02pm UTC](https://discourse.julialang.org/t/ann-delaunaytriangulation-jl-v0-5-constrained-triangulation-and-mesh-refinement/97508/4 "2023-04-24T18:02:54Z")

</div>

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](https://github.com/DanielVandH/DelaunayTriangulation.jl) 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.

```julia
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)
end

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/a/8a200a2a728d6df1e3984d1879f51f021a5728f0.jpeg)

And a nice animation to finish off.

[https://user-images.githubusercontent.com/95613936/234076094-41489505-8fe2-431e-a1a2-e8b25ecfdb8c.mp4](https://user-images.githubusercontent.com/95613936/234076094-41489505-8fe2-431e-a1a2-e8b25ecfdb8c.mp4)