Finding points in a Delaunay triangulation

Hello everyone,

Is there a specific package/function that computes the id’s of the simplices that contain a given (probably large) set of points in a n-d Delaunay triangulation? Something equivalent to Matlab’s tsearchn function.

I am using the Delaunay.jl package but found no useful function.

Thank you all!

We have some search algorithms in Meshes.jl that you can take a look: Neighbor search · Meshes.jl

Currently refactoring some of these, but the stable version may already have what you need?

Check the test suite for examples of usage.

Thanks for your answer. I went through the package but didn’t find what I was looking for. Thank you anyway!!

This thread [ANN] Announcing CGAL.jl might provide some more leads.

Check locate in GitHub - JuliaGeometry/VoronoiDelaunay.jl: Fast and robust Voronoi & Delaunay tessellation creation with Julia

1 Like

Howdy!

I’m afraid CGAL.jl can’t provide a solution for OP’s general issue. The mappings I’ve exposed only pertain to 2D Delaunay Triangulations. More work could be done to make 3D ones available with some effort, but even then it seems too specific regarding what they’re looking for.

AFAIK CGAL does not have generic n-D Delaunay Triangulation algorithms, but it might just be a matter of exploring the plethora of packages deeper. I’d love to be proven wrong!

1 Like

Hello there!
After searching, it seems that there are no available packages for this specific topic yet. I shall wait for upcoming updates!
Thank you all.

@jmaffi please provide a MWE with what you expect as the input and output. The definition of the problem is unclear still.

Hello all! I’m sorry for the delayed answer… quite busy months.
I apologise if the original question was not clear enough. Hope this example will help.

Say that I have a [1 3]x[1 3] grid in 2D space:

x = [1 2 3;1 2 3;1 2 3]
3×3 Matrix{Int64}:
 1  2  3
 1  2  3
 1  2  3

y = [1 1 1;2 2 2; 3 3 3]
3×3 Matrix{Int64}:
 1  1  1
 2  2  2
 3  3  3

The complete list of points is:

9×2 Matrix{Int64}:
 1  1
 1  2
 1  3
 2  1
 2  2
 2  3
 3  1
 3  2
 3  3

A Delaunay triangulation would have the following connectivity list:

8×3 Matrix{Int64}:
 1  4  2
 4  5  2
 2  5  3
 5  8  6
 3  5  6
 4  7  5
 6  8  9
 5  7  8

So, a total of 8 triangles connecting all 9 points of the grid. I plotted it in Matlab:
delaunay

If I want to know which triangle holds the point say [2.4, 1.2], I expect a function that would take this point as an input and provide the number 6 as a result, since the 6th entry of the connectivity list is the triangle defined by points [4 7 5], which are the [2 1], [3 1] and [2 2] respectively.

I am looking for such a function, not only for 2D cases, but for a general n-D triangulation. The Matlab equivalent is tsearchn.

I hope it’s clearer this time.

Thank you all once again.

Yes, it is much clearer now @jmaffi.

This is a solution with Meshes.jl:

using Meshes, MeshViz

import GLMakie as Mke

grid = CartesianGrid(Point(1.,1.), Point(3.,3.), dims=(2,2))
mesh = triangulate(grid) # triangles of your example

viz(mesh, showfacets = true)

# K-nearest search with K=1
searcher = KNearestSearch(mesh, 1)

search(Point(2.4,1.2), searcher)
1-element view(::Vector{Int64}, 1:1) with eltype Int64:
 3

This means that the triangle 3 is the closest to the point of interest:

mesh[3]
Triangle{2,Float64}
  └─Point(2.0, 1.0)
  └─Point(3.0, 1.0)
  └─Point(3.0, 2.0)
2 Likes

@juliohm thanks a lot for your answer.

Is it possible that triangulate works only for 2D grids? That is what I understand from the documentation. I tried creating the following 3D grid:

using Meshes, MeshViz

grid = CartesianGrid(Point(1,0,0), Point(10,5,5), dims = (10,6,6))
10×6×6 CartesianGrid{3,Float64}
  minimum: Point(1.0, 0.0, 0.0)
  maximum: Point(10.0, 5.0, 5.0)
  spacing: (0.9, 0.8333333333333334, 0.8333333333333334)

import GLMakie as Mke
viz(grid, showfacets = true)

This is what I get, as expected:

cube

However, if I try to triangulate it:

mesh = triangulate(grid)

ERROR: MethodError: no method matching triangulate(::Hexahedron{3, Float64, Vector{Point3}})
Closest candidates are:
  triangulate(::Box{2}) at C:\Users\JM\.julia\packages\Meshes\mqr1d\src\discretization.jl:104
  triangulate(::Triangle{Dim, T} where {Dim, T}) at C:\Users\JM\.julia\packages\Meshes\mqr1d\src\discretization.jl:106
  triangulate(::Quadrangle{Dim, T} where {Dim, T}) at C:\Users\JM\.julia\packages\Meshes\mqr1d\src\discretization.jl:108
  ...
Stacktrace:
 [1] triangulate(mesh::CartesianGrid{3, Float64})
   @ Meshes C:\Users\JM\.julia\packages\Meshes\mqr1d\src\discretization.jl:128
 [2] top-level scope
   @ none:1

I don’t know if the package does not compute Delaunay triangulations for higher dimensions, or if I’m using the function in a wrong way. By the error message, it seems as if the code interprets the CartesianGrid as an Hexahedron (which is not the case).

So I would need to create a 3-D or even 4-D Delaunay triangulations and use the KNearestSearch function you suggested (or other).

I appreciate your help.

Regards,

@jmaffi you cannot split an hexahedron into triangles. Hexahedron are 3D geometries and triangles are 2D geometries. This is what the error message is saying.

You could split the hexahedrons into tetrahedrons but that is not implemented yet. Contributions are welcome.

Yes, but the grid is not a hexahedron at all, it’s a rectangular prism, as you can see in the previous image.
This 3D grid can be successfully triangulated into tetrahedra.
The result with Matlab is the following:

prism

And who is saying that the grid is a hexahedron? The error message is clear: it is saying that in order to triangulate a grid you need to triangulate its elements, which are hexahedra. Hence, the error.

That is what I said, except that it is not called a triangulation anymore, it is a “tetrahedrization”. More generally, you are after simplex decompositions of grids.

Ah, I did not get you the first time! Thank you for clarifying.
Guess I’ll be waiting for updates.
Regards,