Delaunay Triangulation of Points Scattered Over Triangular Domain

[zekeriya.sari]

I want to triangulate a set of points that are scattered over a two-dimensional triangular domain. For this purpose, I have been looking into a couple of packages including VoronoiDelaunay, GeometryBasics, GeometricalPredicates, Meshing. As the plotting tool, I want to use Makie .

I started with GeometryBasics and saw that it is possible to mesh a rectangle.

using GeometryBasics
using Makie

rect = Rect(0., 0., 1., 1.)
msh = GeometryBasics.mesh(rect)
scene = Makie.mesh(msh, color=1:length(msh.position))
wireframe!(msh, linewidth=2)
scatter!(msh.position, markersize=10)

But the problem here is that I only have two triangles. To get more fine grained grid of points, it is possible to use Tesselation.

using GeometryBasics

using Makie

rect = Rect(0., 0., 1., 1.)

msh = GeometryBasics.mesh(Tesselation(rect, (5,5)))

scene = Makie.mesh(msh, color=1:length(msh.position))

wireframe!(msh, linewidth=2)

scatter!(msh.position, markersize=10)

The grid looks fine, but the problem is that the plane is rectangular, but not triangular. I tried to construct a triangle and use Tesselation.

using GeometryBasics
using Makie

outerpnts = [Point(0., 0), Point(1., 0.), Point(0.5, 1.)]
tri = Triangle(outerpnts...)
Tesselation(tri)

But I got a method error, since tri is not of type GeometryPrimitive

julia> Tesselation(tri)
ERROR: MethodError: no method matching Tesselation(::GeometryBasics.Ngon{2, Float64, 3, Point2{Float64}})
Closest candidates are:
Tesselation(::Any, ::Integer) at /home/sari/.julia/packages/GeometryBasics/inb5O/src/interfaces.jl:66
Tesselation(::GeometryPrimitive{Dim, T}, ::Tuple{Vararg{var"#s10", N}} where var"#s10"<:Integer) where {Dim, T, N} at /home/sari/.julia/packages/GeometryBasics/inb5O/src/interfaces.jl:61
Stacktrace:
[1] top-level scope
@ REPL[16]:1
[2] eval
@ ./boot.jl:360 [inlined]

Then, I thought that triangulation can be constructed from a set of points that are inside of a triangle

using GeometryBasics
using Makie

# Construct a triangle
outerpnts = [Point(0., 0), Point(1., 0.), Point(0.5, 1.)]
tri = Triangle(outerpnts...)

# Generate points that are inside of the triangle
innerpnts = Point[]
while length(innerpnts) ≤ 10
p = Point(rand(2)...)
p ∈ tri && push!(innerpnts, p)
end

# Mesh the inner and outer points
totalmsh = GeometryBasics.mesh(append!(outerpnts, innerpnts))

# Plot mesh
mesh(totalmsh, color=1 : length(totalmsh.position))
wireframe!(totalmsh, linewidth=2)
scatter!(totalmsh.position, markersize=10)

Here is what I got. The triangulation is not complete. Besides, the triangulation is not correct since some of the triangles overlap.

So I thought that VoronoiDelaunay should be the package I should use. I could construct a Delaunay tesselation but I had to use Gadfly as the plotting package. Also, again the region is not triangular.

using VoronoiDelaunay
using Gadfly

tess = DelaunayTessellation()
width = max_coord - min_coord
pts = Point2D[Point(min_coord + rand() * width, min_coord + rand() * width) for i in 1:100]
push!(tess, pts)

x, y = getplotxy(delaunayedges(tess))
set_default_plot_size(15cm, 15cm)
plot(x=x, y=y, Geom.path)

Any suggestions for the triangulation of two-dimensional points scattered over triangular meshes?

[jkbest2]

You might try TriangleMesh.jl. We’ve been using it to develop SPDEPrecisionMatrices.jl. The readme includes an example of creating a mesh.

[leethargo]

There’s another Julia wrapper of TRIANGLE in Triangulate.jl, which I found easier to use and is actively maintained.

[zekeriya.sari]

Thanks for the advice @jkbest2 . Based on the readme you referred to, here is the complete script for the solution.

using Makie 
using GeometryBasics 
using Clustering
import TriangleMesh

# Construct a triangle 
outerpnts = [Point(0., 0), Point(1., 0.), Point(0.5, 1.)]
tri = Triangle(outerpnts...)

# Generate points that are inside of the triangle
innerpnts = Point[]
while length(innerpnts) ≤ 2000
    p = Point(rand(2)...)
    p ∈ tri && push!(innerpnts, p)
end
innerpntsmat = [getindex.(innerpnts, 1) getindex.(innerpnts, 2)]
innernodes = collect(kmeans(innerpntsmat', 100).centers')
outernodes = [getindex.(outerpnts, 1) getindex.(outerpnts, 2)]
msh = TriangleMesh.create_mesh([outernodes; innernodes])

# Convert TriangleMesh.TriMesh to GeometryBasics.Mesh
pts = [Point(val[1], val[2]) for val in eachcol(msh.point)]
fcs = [TriangleFace(val[1], val[2], val[3]) for val in eachcol(msh.cell)]
msh = GeometryBasics.Mesh(pts, fcs)

# Plot mesh
scn = Makie.mesh(msh, color = 1 : length(msh.position))
Makie.wireframe!(msh) 

[zekeriya.sari]

Thank you for the advice @leethargo . Here is the result I got when I used Triangulate.

using Makie 
using GeometryBasics 
using Clustering
using Triangulate

# Construct a triangle 
outerpnts = [Point(0., 0), Point(1., 0.), Point(0.5, 1.)]
tri = Triangle(outerpnts...)

# Generate points that are inside of the triangle
innerpnts = Point[]
while length(innerpnts) ≤ 2000
    p = Point(rand(2)...)
    p ∈ tri && push!(innerpnts, p)
end
innerpntsmat = [getindex.(innerpnts, 1) getindex.(innerpnts, 2)]
innernodes = collect(kmeans(innerpntsmat', 100).centers')
outernodes = [getindex.(outerpnts, 1) getindex.(outerpnts, 2)]

# Triangulate the points.
triin = TriangulateIO() 
triin.pointlist = [outernodes; innernodes]'
msh, _ = triangulate("vcDQ", triin)

# Convert TriangleMesh.TriMesh to GeometryBasics.Mesh
pts = [Point(val[1], val[2]) for val in eachcol(msh.pointlist)]
fcs = [TriangleFace(val[1], val[2], val[3]) for val in eachcol(msh.trianglelist)]
msh = GeometryBasics.Mesh(pts, fcs)

# Plot mesh
scn = Makie.mesh(msh, color = 1 : length(msh.position))
Makie.wireframe!(msh)