Use DelaunayTriangulation.jl for tricontourf

Project.toml

@@ -11,6 +11,7 @@ ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
+DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -33,20 +34,19 @@ MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Match = "7eb4fadd-790c-5f42-8a69-bfa0b872bfbf"
 MathTeXEngine = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
-MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca"
 Observables = "510215fc-4207-5dde-b226-833fc4488ee2"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Packing = "19eb6ba3-879d-56ad-ad62-d5c202156566"
 PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043"
 PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RelocatableFolders = "05181044-ff0b-4ac5-8273-598c1e38db00"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 Showoff = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
 SignedDistanceFields = "73760f76-fbc4-59ce-8f25-708e95d2df96"
-PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StableHashTraits = "c5dd0088-6c3f-4803-b00e-f31a60c170fa"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -63,6 +63,7 @@ ColorSchemes = "3.5"
 ColorTypes = "0.8, 0.9, 0.10, 0.11"
 Colors = "0.9, 0.10, 0.11, 0.12"
 Contour = "0.5, 0.6"
+DelaunayTriangulation = "0.6.2"
 Distributions = "0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
 DocStringExtensions = "0.8, 0.9"
 FFMPEG = "0.2, 0.3, 0.4"
@@ -81,17 +82,16 @@ LaTeXStrings = "1.2"
 MakieCore = "=0.6.3"
 Match = "1.1"
 MathTeXEngine = "0.5"
-MiniQhull = "0.4"
 Observables = "0.5.3"
 OffsetArrays = "1"
 Packing = "0.5"
 PlotUtils = "1"
 PolygonOps = "0.1.1"
+PrecompileTools = "1.0"
 RelocatableFolders = "0.1, 0.2, 0.3, 1.0"
 Setfield = "1"
 Showoff = "0.3, 1.0.2"
 SignedDistanceFields = "0.4"
-PrecompileTools = "1.0"
 StableHashTraits = "0.3"
 StatsBase = "0.31, 0.32, 0.33"
 StatsFuns = "0.9, 1.0"

ReferenceTests/Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0"
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"

ReferenceTests/src/ReferenceTests.jl

@@ -27,6 +27,7 @@ using Colors
 using LaTeXStrings
 using GeometryBasics
 using DelimitedFiles
+using DelaunayTriangulation
 
 basedir(files...) = normpath(joinpath(@__DIR__, "..", files...))
 loadasset(files...) = FileIO.load(assetpath(files...))

ReferenceTests/src/tests/examples2d.jl

@@ -676,6 +676,73 @@ end
     f
 end
 
+@reference_test "tricontourf with boundary nodes" begin
+    n = 20
+    angles = range(0, 2pi, length = n+1)[1:end-1]
+    x = [cos.(angles); 2 .* cos.(angles .+ pi/n)]
+    y = [sin.(angles); 2 .* sin.(angles .+ pi/n)]
+    z = (x .- 0.5).^2 + (y .- 0.5).^2 .+ 0.5.* RNG.randn.()
+
+    inner = [n:-1:1; n] # clockwise inner 
+    outer = [(n+1):(2n); n+1] # counter-clockwise outer
+    boundary_nodes = [[outer], [inner]]
+    f, ax, _ = tricontourf(x, y, z; boundary_nodes)
+    scatter!(x, y, color = z, strokewidth = 1, strokecolor = :black)
+    f
+end
+
+@reference_test "tricontourf with boundary nodes and edges" begin
+    curve_1 = [
+    [(0.0, 0.0), (5.0, 0.0), (10.0, 0.0), (15.0, 0.0), (20.0, 0.0), (25.0, 0.0)],
+    [(25.0, 0.0), (25.0, 5.0), (25.0, 10.0), (25.0, 15.0), (25.0, 20.0), (25.0, 25.0)],
+    [(25.0, 25.0), (20.0, 25.0), (15.0, 25.0), (10.0, 25.0), (5.0, 25.0), (0.0, 25.0)],
+    [(0.0, 25.0), (0.0, 20.0), (0.0, 15.0), (0.0, 10.0), (0.0, 5.0), (0.0, 0.0)]
+    ]  
+    curve_2 = [
+        [(4.0, 6.0), (4.0, 14.0), (4.0, 20.0), (18.0, 20.0), (20.0, 20.0)],
+        [(20.0, 20.0), (20.0, 16.0), (20.0, 12.0), (20.0, 8.0), (20.0, 4.0)],
+        [(20.0, 4.0), (16.0, 4.0), (12.0, 4.0), (8.0, 4.0), (4.0, 4.0), (4.0, 6.0)]
+    ] 
+    curve_3 = [
+        [(12.906, 10.912), (16.0, 12.0), (16.16, 14.46), (16.29, 17.06),
+        (13.13, 16.86), (8.92, 16.4), (8.8, 10.9), (12.906, 10.912)]
+    ] 
+    curves = [curve_1, curve_2, curve_3]
+    points = [
+        (3.0, 23.0), (9.0, 24.0), (9.2, 22.0), (14.8, 22.8), (16.0, 22.0),
+        (23.0, 23.0), (22.6, 19.0), (23.8, 17.8), (22.0, 14.0), (22.0, 11.0),
+        (24.0, 6.0), (23.0, 2.0), (19.0, 1.0), (16.0, 3.0), (10.0, 1.0), (11.0, 3.0),
+        (6.0, 2.0), (6.2, 3.0), (2.0, 3.0), (2.6, 6.2), (2.0, 8.0), (2.0, 11.0),
+        (5.0, 12.0), (2.0, 17.0), (3.0, 19.0), (6.0, 18.0), (6.5, 14.5),
+        (13.0, 19.0), (13.0, 12.0), (16.0, 8.0), (9.8, 8.0), (7.5, 6.0),
+        (12.0, 13.0), (19.0, 15.0)
+    ]
+    boundary_nodes, points = convert_boundary_points_to_indices(curves; existing_points=points)
+    edges = Set(((1, 19), (19, 12), (46, 4), (45, 12)))
+
+    x = getx.(points)
+    y = gety.(points)
+    z = (x .- 1) .* (y .+ 1) 
+    f, ax, _ = tricontourf(x, y, z, boundary_nodes = boundary_nodes, edges = edges, levels = 30)
+    f
+end
+
+@reference_test "tricontourf with provided triangulation" begin
+    θ = [LinRange(0, 2π * (1 - 1/19), 20); 0]
+    xy = Vector{Vector{Vector{NTuple{2,Float64}}}}()
+    cx = [0.0, 3.0]
+    for i in 1:2
+        push!(xy, [[(cx[i] + cos(θ), sin(θ)) for θ in θ]])
+        push!(xy, [[(cx[i] + 0.5cos(θ), 0.5sin(θ)) for θ in reverse(θ)]])
+    end
+    boundary_nodes, points = convert_boundary_points_to_indices(xy)
+    tri = triangulate(points; boundary_nodes=boundary_nodes, check_arguments=false)
+    z = [(x - 3/2)^2 + y^2 for (x, y) in each_point(tri)]
+
+    f, ax, tr = tricontourf(tri, z, colormap = :matter)
+    f
+end
+
 @reference_test "contour labels 2D" begin
     paraboloid = (x, y) -> 10(x^2 + y^2)
 

ReferenceTests/src/tests/refimages.jl

@@ -10,6 +10,7 @@ using ReferenceTests.LaTeXStrings
 using ReferenceTests.DelimitedFiles
 using ReferenceTests.Test
 using ReferenceTests.Colors: RGB, N0f8
+using ReferenceTests.DelaunayTriangulation
 using Makie: Record, volume
 
 @testset "primitives" begin

docs/Project.toml

@@ -6,6 +6,7 @@ Chain = "8be319e6-bccf-4806-a6f7-6fae938471bc"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 ColorVectorSpace = "c3611d14-8923-5661-9e6a-0046d554d3a4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

docs/examples/plotting_functions/tricontourf.md

@@ -74,6 +74,101 @@ f
 ```
 \end{examplefigure}
 
+Constrained triangulations can also be used by passing boundary nodes or edges as keyword arguments into `tricontourf`, rather than having to identify the triangles manually as above. The form of the boundary nodes should conform with the specification in DelaunayTriangulation.jl. For example, the above annulus could have been plotted as follows.
+
+\begin{examplefigure}{svg = true}
+```julia
+using CairoMakie
+
+using Random
+Random.seed!(123)
+
+n = 20
+angles = range(0, 2pi, length = n+1)[1:end-1]
+x = [cos.(angles); 2 .* cos.(angles .+ pi/n)]
+y = [sin.(angles); 2 .* sin.(angles .+ pi/n)]
+z = (x .- 0.5).^2 + (y .- 0.5).^2 .+ 0.5.*randn.()
+
+inner = [n:-1:1; n] # clockwise inner 
+outer = [(n+1):(2n); n+1] # counter-clockwise outer
+boundary_nodes = [[outer], [inner]]
+f, ax, _ = tricontourf(x, y, z; boundary_nodes)
+scatter!(x, y, color = z, strokewidth = 1, strokecolor = :black)
+f
+```
+\end{examplefigure}
+
+Boundary nodes make it possible to support more complicated regions, possibly with holes, than is possible by only providing points themselves.
+
+\begin{examplefigure}{svg = true}
+```julia
+using CairoMakie, DelaunayTriangulation
+CairoMakie.activate!() # hide 
+
+## Start by defining the boundaries, and then convert to the appropriate interface 
+curve_1 = [
+    [(0.0, 0.0), (5.0, 0.0), (10.0, 0.0), (15.0, 0.0), (20.0, 0.0), (25.0, 0.0)],
+    [(25.0, 0.0), (25.0, 5.0), (25.0, 10.0), (25.0, 15.0), (25.0, 20.0), (25.0, 25.0)],
+    [(25.0, 25.0), (20.0, 25.0), (15.0, 25.0), (10.0, 25.0), (5.0, 25.0), (0.0, 25.0)],
+    [(0.0, 25.0), (0.0, 20.0), (0.0, 15.0), (0.0, 10.0), (0.0, 5.0), (0.0, 0.0)]
+] # outer-most boundary: counter-clockwise  
+curve_2 = [
+    [(4.0, 6.0), (4.0, 14.0), (4.0, 20.0), (18.0, 20.0), (20.0, 20.0)],
+    [(20.0, 20.0), (20.0, 16.0), (20.0, 12.0), (20.0, 8.0), (20.0, 4.0)],
+    [(20.0, 4.0), (16.0, 4.0), (12.0, 4.0), (8.0, 4.0), (4.0, 4.0), (4.0, 6.0)]
+] # inner boundary: clockwise 
+curve_3 = [
+    [(12.906, 10.912), (16.0, 12.0), (16.16, 14.46), (16.29, 17.06),
+    (13.13, 16.86), (8.92, 16.4), (8.8, 10.9), (12.906, 10.912)]
+] # this is inside curve_2, so it's counter-clockwise 
+curves = [curve_1, curve_2, curve_3]
+points = [
+    (3.0, 23.0), (9.0, 24.0), (9.2, 22.0), (14.8, 22.8), (16.0, 22.0),
+    (23.0, 23.0), (22.6, 19.0), (23.8, 17.8), (22.0, 14.0), (22.0, 11.0),
+    (24.0, 6.0), (23.0, 2.0), (19.0, 1.0), (16.0, 3.0), (10.0, 1.0), (11.0, 3.0),
+    (6.0, 2.0), (6.2, 3.0), (2.0, 3.0), (2.6, 6.2), (2.0, 8.0), (2.0, 11.0),
+    (5.0, 12.0), (2.0, 17.0), (3.0, 19.0), (6.0, 18.0), (6.5, 14.5),
+    (13.0, 19.0), (13.0, 12.0), (16.0, 8.0), (9.8, 8.0), (7.5, 6.0),
+    (12.0, 13.0), (19.0, 15.0)
+]
+boundary_nodes, points = convert_boundary_points_to_indices(curves; existing_points=points)
+edges = Set(((1, 19), (19, 12), (46, 4), (45, 12)))
+
+## Extract the x, y 
+x = getx.(points)
+y = gety.(points)
+z = (x .- 1) .* (y .+ 1) 
+f, ax, _ = tricontourf(x, y, z, boundary_nodes = boundary_nodes, edges = edges, levels = 30)
+f
+```
+\end{examplefigure}
+
+You can also provide a `Triangulation` object itself from DelaunayTriangulation.jl into `tricontourf`.
+
+\begin{examplefigure}{svg = true}
+```julia
+using CairoMakie, DelaunayTriangulation
+CairoMakie.activate!() # hide 
+
+using Random
+Random.seed!(1234)
+
+θ = [LinRange(0, 2π * (1 - 1/19), 20); 0]
+xy = Vector{Vector{Vector{NTuple{2,Float64}}}}()
+cx = [0.0, 3.0]
+for i in 1:2
+    push!(xy, [[(cx[i] + cos(θ), sin(θ)) for θ in θ]])
+    push!(xy, [[(cx[i] + 0.5cos(θ), 0.5sin(θ)) for θ in reverse(θ)]])
+end
+boundary_nodes, points = convert_boundary_points_to_indices(xy)
+tri = triangulate(points; boundary_nodes=boundary_nodes, check_arguments=false)
+z = [(x - 3/2)^2 + y^2 for (x, y) in each_point(tri)]
+
+f, ax, tr = tricontourf(tri, z, colormap = :matter)
+f
+```
+\end{examplefigure}
+
 #### Relative mode
 
 Sometimes it's beneficial to drop one part of the range of values, usually towards the outer boundary.

src/Makie.jl

@@ -49,7 +49,7 @@ import ImageIO
 import FileIO
 import SparseArrays
 import TriplotBase
-import MiniQhull
+import DelaunayTriangulation as DelTri 
 import Setfield
 import REPL
 

src/basic_recipes/tricontourf.jl

@@ -1,10 +1,12 @@
 struct DelaunayTriangulation end
 
 """
+    tricontourf(triangles::Triangulation, zs; kwargs...)
     tricontourf(xs, ys, zs; kwargs...)
 
-Plots a filled tricontour of the height information in `zs` at horizontal positions `xs`
-and vertical positions `ys`.
+Plots a filled tricontour of the height information in `zs` at over the points in the 
+provided triangulation, or alternatively at the horizontal positions `xs` and 
+vertical positions `ys`.
 
 ## Attributes
 
@@ -14,7 +16,8 @@ and vertical positions `ys`.
 - `mode = :normal` sets the way in which a vector of levels is interpreted, if it's set to `:relative`, each number is interpreted as a fraction between the minimum and maximum values of `zs`. For example, `levels = 0.1:0.1:1.0` would exclude the lower 10% of data.
 - `extendlow = nothing`. This sets the color of an optional additional band from `minimum(zs)` to the lowest value in `levels`. If it's `:auto`, the lower end of the colormap is picked and the remaining colors are shifted accordingly. If it's any color representation, this color is used. If it's `nothing`, no band is added.
 - `extendhigh = nothing`. This sets the color of an optional additional band from the highest value of `levels` to `maximum(zs)`. If it's `:auto`, the high end of the colormap is picked and the remaining colors are shifted accordingly. If it's any color representation, this color is used. If it's `nothing`, no band is added.
-- `triangulation = DelaunayTriangulation()`. The mode with which the points in `xs` and `ys` are triangulated. Passing `DelaunayTriangulation()` performs a delaunay triangulation. You can also pass a preexisting triangulation as an `AbstractMatrix{<:Int}` with size (3, n), where each column specifies the vertex indices of one triangle.
+- `boundary_nodes = nothing`: Boundary constraints for the triangulation, in case the triangulation is constructed from DelaunayTriangulation()` above. Boundary nodes should match the specification in DelaunayTriangulation.jl. Note that these are not used if a triangulation is manually provided.
+- `edges = nothing`: Constrained edges for the triangulation, in case the triangulation is constructed from DelaunayTriangulation()` above. The edges should not intersect each other edge at an interior. Note that these are not used if a triangulation is manually provided.
 
 ### Generic
 
@@ -41,12 +44,36 @@ $(ATTRIBUTES)
         nan_color = :transparent,
         inspectable = theme(scene, :inspectable),
         transparency = false,
-        triangulation = DelaunayTriangulation()
+        triangulation = DelaunayTriangulation(),
+        boundary_nodes = nothing, 
+        edges = nothing,
     )
 end
 
-function Makie.convert_arguments(::Type{<:Tricontourf}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, z::AbstractVector{<:Real})
-    map(x -> elconvert(Float32, x), (x, y, z))
+function Makie.used_attributes(::Type{<:Tricontourf}, ::AbstractVector{<:Real}, ::AbstractVector{<:Real}, ::AbstractVector{<:Real})
+    return (:boundary_nodes, :edges, :triangulation)
+end
+
+function Makie.convert_arguments(::Type{<:Tricontourf}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, z::AbstractVector{<:Real}; 
+    boundary_nodes=nothing, 
+    edges=nothing, 
+    triangulation=DelaunayTriangulation())
+    z = elconvert(Float32, z)
+    points = [x'; y']
+    if triangulation isa DelaunayTriangulation
+        tri = DelTri.triangulate(points; boundary_nodes=boundary_nodes, edges=edges, check_arguments=false)
+    else
+        # Wrap user's provided triangulation into a Triangulation. Their triangulation must be such that DelTri.add_triangle! is defined. 
+        if typeof(triangulation) <: AbstractMatrix{<:Int} && size(triangulation, 1) != 3 
+            triangulation = triangulation' 
+        end
+        tri = DelTri.Triangulation(points)
+        triangles = DelTri.get_triangles(tri)
+        for τ in DelTri.each_solid_triangle(triangulation)
+            DelTri.add_triangle!(triangles, τ)
+        end
+    end
+    return (tri, z)
 end
 
 function compute_contourf_colormap(levels, cmap, elow, ehigh)
@@ -90,8 +117,14 @@ function compute_highcolor(eh, cmap)
     end
 end
 
-function Makie.plot!(c::Tricontourf{<:Tuple{<:AbstractVector{<:Real},<:AbstractVector{<:Real},<:AbstractVector{<:Real}}})
-    xs, ys, zs = c[1:3]
+function Makie.plot!(c::Tricontourf{<:Tuple{<:DelTri.Triangulation, <:AbstractVector{<:Real}}})
+    tri, zs = c[1:2]
+    xs = lift(tri) do tri 
+        return [DelTri.getx(p) for p in DelTri.each_point(tri)]
+    end 
+    ys = lift(tri) do tri 
+        return [DelTri.gety(p) for p in DelTri.each_point(tri)]
+    end
 
     c.attributes[:_computed_levels] = lift(c, zs, c.levels, c.mode) do zs, levels, mode
         return _get_isoband_levels(Val(mode), levels, vec(zs))
@@ -131,7 +164,7 @@ function Makie.plot!(c::Tricontourf{<:Tuple{<:AbstractVector{<:Real},<:AbstractV
         lows = levels[1:end-1]
         highs = levels[2:end]
 
-        trianglelist = compute_triangulation(triangulation, xs, ys)
+        trianglelist = compute_triangulation(triangulation)
         filledcontours = filled_tricontours(xs, ys, zs, trianglelist, levels)
 
         levelcenters = (highs .+ lows) ./ 2
@@ -154,10 +187,10 @@ function Makie.plot!(c::Tricontourf{<:Tuple{<:AbstractVector{<:Real},<:AbstractV
         return
     end
 
-    onany(calculate_polys, c, xs, ys, zs, c._computed_levels, is_extended_low, is_extended_high, c.triangulation)
+    onany(calculate_polys, c, tri, xs, ys, zs, c._computed_levels, is_extended_low, is_extended_high)
     # onany doesn't get called without a push, so we call
     # it on a first run!
-    calculate_polys(xs[], ys[], zs[], c._computed_levels[], is_extended_low[], is_extended_high[], c.triangulation[])
+    calculate_polys(xs[], ys[], zs[], c._computed_levels[], is_extended_low[], is_extended_high[], tri[])
 
     poly!(c,
         polys,
@@ -175,16 +208,8 @@ function Makie.plot!(c::Tricontourf{<:Tuple{<:AbstractVector{<:Real},<:AbstractV
     )
 end
 
-function compute_triangulation(::DelaunayTriangulation, xs, ys)
-    vertices = [xs'; ys']
-    return MiniQhull.delaunay(vertices)
-end
-
-function compute_triangulation(triangulation::AbstractMatrix{<:Int}, xs, ys)
-    if size(triangulation, 1) != 3
-        throw(ArgumentError("Triangulation matrix must be of size (3, n) but is of size $(size(triangulation))."))
-    end
-    triangulation
+function compute_triangulation(tri)
+    return [T[j] for T in DelTri.each_solid_triangle(tri), j in 1:3]'
 end
 
 # FIXME: TriplotBase augments levels so here the implementation is just repeated without that step