# Delaunay Triangulation of Points Scattered Over Triangular Domain

**URL:** <https://discourse.julialang.org/t/delaunay-triangulation-of-points-scattered-over-triangular-domain/54653>\
**Category:** General Usage\
**Tags:** geometry\
**Created:** [February 4, 2021, 9:30pm UTC](https://discourse.julialang.org/t/delaunay-triangulation-of-points-scattered-over-triangular-domain/54653 "2021-02-04T21:30:24Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![zekeriya.sari](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zekeriya.sari/32/13695_2.png) [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)\
**Post date:** [February 4, 2021, 9:30pm UTC](https://discourse.julialang.org/t/delaunay-triangulation-of-points-scattered-over-triangular-domain/54653/1 "2021-02-04T21:30:24Z")

</div>

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.

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

```

 ![rectangle_with_two_triangles](https://global.discourse-cdn.com/julialang/original/3X/1/d/1dd4c0ca6b78b839406e13aef828b728e1ec729d.jpeg)

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`.

```julia

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)

```

 ![rectangle](https://global.discourse-cdn.com/julialang/original/3X/6/9/6971ed0e340851aa4d1608852573828f69be555c.png)  
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`.

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

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

 ![rectangular_mesh](https://global.discourse-cdn.com/julialang/original/3X/6/4/644b590a0f1f3eaaacf1e32649a00255213f38ca.png)

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.

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

```

 ![voronoi](https://global.discourse-cdn.com/julialang/original/3X/b/a/bac88e371a58a3904d10d46d6852c1e51ca2934f.png)

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

---

<div class="post-metadata">

**Author:** ![jkbest2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jkbest2/32/7350_2.png) [@jkbest2](https://discourse.julialang.org/u/jkbest2)\
**Post date:** [February 4, 2021, 11:36pm UTC](https://discourse.julialang.org/t/delaunay-triangulation-of-points-scattered-over-triangular-domain/54653/2 "2021-02-04T23:36:10Z")

</div>

You might try [TriangleMesh.jl](https://github.com/konsim83/TriangleMesh.jl). We’ve been using it to develop [SPDEPrecisionMatrices.jl](https://github.com/ElOceanografo/SPDEPrecisionMatrices.jl). The readme includes an example of creating a mesh.

---

<div class="post-metadata">

**Author:** ![leethargo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leethargo/32/6004_2.png) [@leethargo](https://discourse.julialang.org/u/leethargo)\
**Post date:** [February 5, 2021, 10:04am UTC](https://discourse.julialang.org/t/delaunay-triangulation-of-points-scattered-over-triangular-domain/54653/3 "2021-02-05T10:04:21Z")

</div>

There’s another Julia wrapper of `TRIANGLE` in [Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl), which I found easier to use and is actively maintained.

---

<div class="post-metadata">

**Author:** ![zekeriya.sari](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zekeriya.sari/32/13695_2.png) [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)\
**Post date:** [February 5, 2021, 9:28pm UTC](https://discourse.julialang.org/t/delaunay-triangulation-of-points-scattered-over-triangular-domain/54653/4 "2021-02-05T21:28:07Z")

</div>

Thanks for the advice @jkbest2 . Based on the [readme](https://github.com/ElOceanografo/SPDEPrecisionMatrices.jl) you referred to, here is the complete script for the solution.

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

```

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

---

<div class="post-metadata">

**Author:** ![zekeriya.sari](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zekeriya.sari/32/13695_2.png) [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)\
**Post date:** [February 5, 2021, 10:01pm UTC](https://discourse.julialang.org/t/delaunay-triangulation-of-points-scattered-over-triangular-domain/54653/5 "2021-02-05T22:01:48Z")

</div>

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

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

```

 ![triangle2](https://global.discourse-cdn.com/julialang/original/3X/1/c/1cdad5351d40b772a9f97a1f20873ca3143df52e.jpeg)