Delaunay Triangularization

Hello!

I’d like to be able to obtain the following output from MATLAB but in Julia. Here is a toy example in MATLAB:

>> x = rand(5,1)

x =

    0.8147
    0.9058
    0.1270
    0.9134
    0.6324

>> y = rand(5,1)

y =

    0.0975
    0.2785
    0.5469
    0.9575
    0.9649

>> tri = delaunay(x,y)

tri =

     2     3     1
     5     3     2
     5     2     4

Is there any way to obtain the output (i.e. the tri variable) in Julia?

Any help would be greatly appreciated.

Check out https://github.com/JuliaPDE/SurveyofPDEPackages#grids

2 Likes

Thank you for the link. I tried the VoronoiDelaunay package but that doesn’t seem to do what I want (probably I don’t know how) and the MiniQhull package requires Linux unless I’m reading this wrong.

My ultimate goal is to complete a spatial econometrics package and I need to be able to build a contiguity based weight matrix.

Yes, MiniQhull is designed for linux at this moment. Perhaps it also works for mac OS but we have not tested it in this scenario…

Have you looked at GeoStats.jl? Depending on what exactly you’re trying to do, it might be there already or maybe an opportunity to combine efforts?

I’ve looked at GeoStats and that currently doesn’t have what I need. I think what I need is very simple and the MiniQHull would work (theoretically) but I’m on Windows.

Like I mentioned earlier, my ultimate goal is to produce a spatial econometrics library for Julia. Most of the code is written (think alpha version) and I’d like to include other options for uses to create spatial weight matrices.

You can do it easily with GMT.jl, but it’s a non-negligible ovehead to have as a dependency

using GMT

triangulate(rand(5,2))
Array{GMT.GMTdataset,1} with 1 segments
First segment DATA
4×3 Array{Float64,2}:
 2.0  0.0  1.0
 4.0  3.0  1.0
 1.0  3.0  2.0
 4.0  1.0  0.0

Thank you for the reply. I’m really looking for a native Julia “Dealunay” function if possible so that users can avoid installing extra software.

I may take a crack at coding one myself and given I can’t leave home, it may be sooner than later!

FYI: I found the following implementation in Python:

Delaunay in Python

I do not know Python but this might be a good place to start.

The link “Delaunay in Python” doesn’t work. However Delaunay triangulation is implemented in the SciPy package (using the Qhull library). For me, the simplest thing would be calling the SciPy from Julia.

Thank you for taking a look at this code. I do know about SciPy, but I’d like to avoid having to call Python code. Nothing against Python, but I want to avoid the “you have to have additional stuff for Julia to work blah, blah”.

All I really need is a very simple 2D Delaunay routine that will take latitude and longitude coordinates and output the vertices like in my MATLAB example above.

If someone can point me to existing code that works, I can try and adapt to Julia. It may not be exhaustive but I think it’s valuable nonetheless.

MiniQhull is a really nice package, thanks for the great work. But what is missing for it to work on Windows?

A bit late to this thread but I’m wondering why the VoronoiDelaunay doesn’t work for you, assuming you don’t need to do any 3D tetrahedralizations, just 2d Delaunay triangulations. At work we’ve been using VoronoiDelaunay to do triangulations of unstructured 2d point clouds. It’s worked great for us. Here’s a little code snippet:

using GeometricalPredicates, VoronoiDelaunay
mutable struct Point2DI <: AbstractPoint2D
  _x::Float64
  _y::Float64
  _i::Int64
end

Point2DI(x::Float64, y::Float64) = Point2DI(x, y, 0)

import GeometricalPredicates.getx
import GeometricalPredicates.gety
import VoronoiDelaunay.isexternal

getx(p::Point2DI) = p._x
gety(p::Point2DI) = p._y
geti(p::Point2DI) = p._i

function isexternal(p::Point2DI)
  getx(p) < VoronoiDelaunay.min_coord || getx(p) > VoronoiDelaunay.max_coord
end

x = randn(10); y = randn(10)

n = length(x)
d0 = VoronoiDelaunay.min_coord
dd = VoronoiDelaunay.max_coord - VoronoiDelaunay.min_coord
dn = VoronoiDelaunay.max_coord
x0,xn = extrema(x)
y0,yn = extrema(y)
dx = xn-x0
dy = yn-y0
u = d0 .+ ((x .- x0) ./ dx) .* dd
v = d0 .+ ((y .- y0) ./ dy) .* dd
P = Point2DI[ Point2DI(ui, vi, i) for (ui, vi, i) in zip(u, v, 1:n) ]
push!(DT, P)

That should give you a delaunay triangulation object DT that you can use for subsequent calculations such as interpolation etc.

2 Likes

Honestly, I don’t know… I have never used windows for technical computing. In any case, we are using cmake.jl to compile the C functions in the project. Thus, this compilation step should be portable. My doubt is about the mechanism we use to find the compiled C objects from Julia, I don’t know if it is portable…

We would be happy if someone is ready to help with the windows extension!

2 Likes

I have open an issue in MiniQhul about the windows and Mac OS extension

2 Likes

Probably not really important, but for what it’s worth, here is the correct link.

I have been investigating this issue and it appears that the open source program Octave uses Qhull. I have performed my toy example above and it works as expected. The Octave code calls Qhull and in the MinGW64 bin folder it appears that the developers of Octave have compiled the Qhull code for use on Windows.

I’m not sure if this helps but I thought I’d bring it up. I can provide more information if anyone is interrested.

Thank you.

What is DT? I tried DT = DelaunayTessellation() but got

julia> push!(DT, P)
+ 
ERROR: MethodError: no method matching push!(::DelaunayTessellation2D{Point2D}, ::Array{Point2DI,1})
Closest candidates are:
  push!(::Any, ::Any, !Matched::Any) at abstractarray.jl:2252
  push!(::Any, ::Any, !Matched::Any, !Matched::Any...) at abstractarray.jl:2253
  push!(!Matched::Array{Any,1}, ::Any) at array.jl:940

Hey,
I used this one in a paper: [just in case you want to see something more: ]

using GR
n, tri = GR.delaunay(x, y)

A complete MWE is as follows:

using Random
using GR: delaunay
Random.seed!(123)
x, y = rand(10), rand(10)
n, tri = delaunay(x,y)
tri
13×3 Array{Int64,2}:
  2  3   6
  4  6   3
  2  6   1
  9  3   8
  9  4   3
  4  9   5
  5  9   8
  7  1   6
  7  5   1
  6  4   7
  4  5   7
 10  5   8
  1  5  10
1 Like

Sorry about that. Looks like I left a line out of my example. The key line I left out was:

DT = DelaunayTessellation2D{Point2DI}(n)

My example is using a custom point type Point2DI. When using a custom point type you need to specify that point type when initialising the DelaunayTessellation. When you do DT = DelaunayTessellation() it returns a DelaunayTessellation2D{Point2D}, rather than a DelaunayTessellation2D{Point2DI}

Discourse won’t let me edit my original post so here’s the updated self-contained example that actually works now:

using GeometricalPredicates, VoronoiDelaunay
mutable struct Point2DI <: AbstractPoint2D
  _x::Float64
  _y::Float64
  _i::Int64
end

Point2DI(x::Float64, y::Float64) = Point2DI(x, y, 0)

import GeometricalPredicates.getx
import GeometricalPredicates.gety
import VoronoiDelaunay.isexternal

getx(p::Point2DI) = p._x
gety(p::Point2DI) = p._y
geti(p::Point2DI) = p._i

function isexternal(p::Point2DI)
  getx(p) < VoronoiDelaunay.min_coord || getx(p) > VoronoiDelaunay.max_coord
end

x = randn(10); y = randn(10)

n = length(x)
DT = DelaunayTessellation2D{Point2DI}(n)
d0 = VoronoiDelaunay.min_coord
dd = VoronoiDelaunay.max_coord - VoronoiDelaunay.min_coord
dn = VoronoiDelaunay.max_coord
x0,xn = extrema(x)
y0,yn = extrema(y)
dx = xn-x0
dy = yn-y0
u = d0 .+ ((x .- x0) ./ dx) .* dd
v = d0 .+ ((y .- y0) ./ dy) .* dd
P = Point2DI[ Point2DI(ui, vi, i) for (ui, vi, i) in zip(u, v, 1:n) ]
push!(DT, P)

Note that you’re not required to define your own point type, you can use the built in Point2D type. I just decided to illustrate using a custom point type in my example.

2 Likes