Use DelaunayTriangulation.jl for tricontourf (#2896)

* Use DelaunayTriangulation.jl * Example wasn't passing edges * Update to v0.6.1 with removed constprop * Change DelaunayTriangulation to 0.6.2 and fix doc code * Add DelaunayTriangulation to doc depdencies * Support passing Triangulation objects into tricontourf * Fix scoping issue * Make the default form of tricontourf tricontourf(tri::DelTri.Triangulation, zs) * Add reference tests * Fix missing DelaunayTriangulation dep * Remove MiniQhull dependency * Accidentally left in deved packages * PrecompileTools is needed * Remove boundary_nodes and edges kwargs * Update the DelaunayTriangulation compat * small docs fixes --------- Co-authored-by: Simon <sdanisch@protonmail.com> Co-authored-by: Julius Krumbiegel <julius.krumbiegel@gmail.com> Co-authored-by: Julius Krumbiegel <22495855+jkrumbiegel@users.noreply.github.com>
MakieOrg · Jun 29, 2023 · c8039e2 · c8039e2
1 parent d93343d
commit c8039e2
Show file tree

Hide file tree

Showing 9 changed files with 217 additions and 25 deletions.
diff --git a/Project.toml b/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"
@@ -34,7 +35,6 @@ 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"
@@ -64,6 +64,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, 0.7"
 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"
@@ -83,7 +84,6 @@ MacroTools = "0.5"
 MakieCore = "=0.6.3"
 Match = "1.1"
 MathTeXEngine = "0.5"
-MiniQhull = "0.4"
 Observables = "0.5.3"
 OffsetArrays = "1"
 Packing = "0.5"

diff --git a/ReferenceTests/Project.toml b/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"

diff --git a/ReferenceTests/src/ReferenceTests.jl b/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...))

diff --git a/ReferenceTests/src/tests/examples2d.jl b/ReferenceTests/src/tests/examples2d.jl
@@ -685,6 +685,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]]
+    tri = DelaunayTriangulation.triangulate([x'; y'], boundary_nodes = boundary_nodes)
+    f, ax, _ = tricontourf(tri, z)
+    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)))
+
+    tri = triangulate(points; boundary_nodes = boundary_nodes, edges = edges, check_arguments = false)
+    z = [(x - 1) * (y + 1) for (x, y) in each_point(tri)]
+    f, ax, _ = tricontourf(tri, z, 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)
 

diff --git a/ReferenceTests/src/tests/refimages.jl b/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

diff --git a/docs/Project.toml b/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"

diff --git a/docs/examples/plotting_functions/tricontourf.md b/docs/examples/plotting_functions/tricontourf.md
@@ -44,6 +44,8 @@ f
 
 #### Triangulation modes
 
+Manual triangulations can be passed as a 3xN matrix of integers, where each column of three integers specifies the indices of the corners of one triangle in the vector of points.
+
 \begin{examplefigure}{svg = true}
 ```julia
 using CairoMakie
@@ -74,6 +76,107 @@ f
 ```
 \end{examplefigure}
 
+By default, `tricontourf` performs unconstrained triangulations.
+Greater control over the triangulation, such as allowing for enforced boundaries, can be achieved by using [DelaunayTriangulation.jl](https://github.com/DanielVandH/DelaunayTriangulation.jl) and passing the resulting triangulation as the first argument of `tricontourf`.
+For example, the above annulus can also be plotted as follows:
+
+\begin{examplefigure}{svg = true}
+```julia
+using CairoMakie
+using DelaunayTriangulation
+CairoMakie.activate!() # hide
+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]]
+points = [x'; y']
+tri = triangulate(points; boundary_nodes = boundary_nodes)
+f, ax, _ = tricontourf(tri, z;
+    axis = (; aspect = 1, title = "Constrained triangulation\nvia DelaunayTriangulation.jl"))
+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
+using 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 
+tri = triangulate(points; boundary_nodes = boundary_nodes, edges = edges, check_arguments = false)
+z = [(x - 1) * (y + 1) for (x, y) in each_point(tri)] # note that each_point preserves the index order
+f, ax, _ = tricontourf(tri, z, levels = 30; axis = (; aspect = 1))
+f
+```
+\end{examplefigure}
+
+\begin{examplefigure}{svg = true}
+```julia
+using CairoMakie
+using 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)] # note that each_point preserves the index order
+
+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.

diff --git a/src/Makie.jl b/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
 import MacroTools

diff --git a/src/basic_recipes/tricontourf.jl b/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 the horizontal positions `xs` and 
+vertical positions `ys`. A `Triangulation` from DelaunayTriangulation.jl can also be provided instead of `xs` and `ys`
+for specifying the triangles, otherwise an unconstrained triangulation of `xs` and `ys` is computed.
 
 ## Attributes
 
@@ -14,7 +16,7 @@ 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.
+- `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, or as a `Triangulation` from DelaunayTriangulation.jl.
 
 ### Generic
 
@@ -41,12 +43,33 @@ $(ATTRIBUTES)
         nan_color = :transparent,
         inspectable = theme(scene, :inspectable),
         transparency = false,
-        triangulation = DelaunayTriangulation()
+        triangulation = DelaunayTriangulation(),
+        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 (:triangulation,)
+end
+
+function Makie.convert_arguments(::Type{<:Tricontourf}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, z::AbstractVector{<:Real}; 
+    triangulation=DelaunayTriangulation())
+    z = elconvert(Float32, z)
+    points = [x'; y']
+    if triangulation isa DelaunayTriangulation
+        tri = DelTri.triangulate(points)
+    elseif !(triangulation isa DelTri.Triangulation)
+        # 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 +113,8 @@ 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]
 
     c.attributes[:_computed_levels] = lift(c, zs, c.levels, c.mode) do zs, levels, mode
         return _get_isoband_levels(Val(mode), levels, vec(zs))
@@ -117,7 +140,7 @@ function Makie.plot!(c::Tricontourf{<:Tuple{<:AbstractVector{<:Real},<:AbstractV
     polys = Observable(PolyType[])
     colors = Observable(Float64[])
 
-    function calculate_polys(xs, ys, zs, levels::Vector{Float32}, is_extended_low, is_extended_high, triangulation)
+    function calculate_polys(triangulation, zs, levels::Vector{Float32}, is_extended_low, is_extended_high)
         empty!(polys[])
         empty!(colors[])
 
@@ -131,7 +154,10 @@ function Makie.plot!(c::Tricontourf{<:Tuple{<:AbstractVector{<:Real},<:AbstractV
         lows = levels[1:end-1]
         highs = levels[2:end]
 
-        trianglelist = compute_triangulation(triangulation, xs, ys)
+        xs = [DelTri.getx(p) for p in DelTri.each_point(triangulation)] # each_point preserves indices
+        ys = [DelTri.gety(p) for p in DelTri.each_point(triangulation)]
+
+        trianglelist = compute_triangulation(triangulation)
         filledcontours = filled_tricontours(xs, ys, zs, trianglelist, levels)
 
         levelcenters = (highs .+ lows) ./ 2
@@ -154,10 +180,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, 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(tri[], zs[], c._computed_levels[], is_extended_low[], is_extended_high[])
 
     poly!(c,
         polys,
@@ -175,16 +201,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