How to speed up the function DelaunayTriangulation?

Question

First define a function meshGrid to generate some points:

meshGrid[{x1_, x2_, y1_, y2_}, h0_] := 
  With[{yh0 = h0*Sqrt[3.]/2}, 
    Array[{(#1 - 1)*h0 + x1 + (1 + (-1)^#2) h0/4, (#2 - 1)*yh0 + 
        y1} &, Ceiling@{(x2 - x1)/h0, (y2 - y1)/yh0}]]~Flatten~1;

p = meshGrid[{-1, 1, -1, 1}, 0.05];

The computing time of DelaunayTriangulation:

Needs["ComputationalGeometry`"];
DelaunayTriangulation[p] // Timing // First

In my computer it gives 18.533s

Matlab saves much time if does the same thing:

====================================Update=====================================

@halirutan really made a great attempt to point the way, but I failed to compile and didn't get the right answer, maybe I should learn something before. Here I find another way in this blogpost, which also relates to Qhull but easier to implement. You can get more infomation from here. Before changing anything, two files need to be downloaded, one is mPower, from which we need is mPower.m, another one is Qhull.

You can get the rest steps from that blog, only step two is worthy of note:

step 2: download qhull for windows, you may need to change the name, and put it into the folder C:\qhull. Then Copy all the *.exe files in bin folder and paste them in folder qhull, errors will occur without this step.

Answer 1

Short answer

Yes, it is possible to speed up the Delaunay-triangulation and make it as fast as it is in Matlab.

If you are not afraid of some setup-work, then one possibility is to use a package which calls a c-implementation of the Delaunay-triangulation. One package I know is qh-math which is available in the Wolfram-library:

This package includes source code and support files needed to create a MathLink-based interface to the Qhull library (http://www.qhull.org) algorithm for Delaunay Triangulation. The sources are based on work done originally by Alban Tsui at the Imperial College of Science, Technology and Medicine.

And btw, this is exactly what Matlab is using: http://www.qhull.org/html/qh-faq.htm#math

Usage

I assume the program qh-math.exe is located in my download-folder. For your system you have to change this in the Install call. The usage is pretty easy. First you Install the MathLink program and after this you can call qDelaunayTriangulation[..] like a normal Mathemtatica function:

lnk = Install["/home/patrick/Downloads/qh-math/qh-math.exe"];

And then you can triangulate your points

meshGrid[{x1_, x2_, y1_, y2_}, h0_] := 
  With[{yh0 = h0*Sqrt[3.]/2}, 
    Array[{(#1 - 1)*h0 + x1 + (1 + (-1)^#2) h0/4, (#2 - 1)*yh0 + 
        y1} &, Ceiling@{(x2 - x1)/h0, (y2 - y1)/yh0}]]~Flatten~1;
p = meshGrid[{-1, 1, -1, 1}, 0.05];

{t, del} = AbsoluteTiming[qDelaunayTriangulation[p]];

On my machine this took only t=0.032471 seconds. The the result looks nice

Graphics[MapIndexed[{ColorData[29, First[#2]], 
    Polygon[#1]} &, (Part[p, #] & /@ del)]]

Please note that the output is different from DelaunayTriangulation. This version really gives a triangle index list like {{5, 6, 2}, {10, 7, 4}, {1, 5, 6},....

Freshly compiled qh-math.exe for Windows

Due to the great efforts of @Oleksandr R. we have now compiled versions of qh-math.exe and all the commandline tools from qhull. Please download a zip with all files for your system:

• qhull.zip for 64bit Windows
• qhull.zip for 32bit Windows

Compiling your own qh-math

I'm on Linux here and since there is no executable program included I had to compile it by myself. Since it can happen, that your program does not work (it's kind of old) you may have to compile it for your machine too. Therefore, I explain it step by step

Compiling: First you download the archive with the sources and unpack it. The following steps all takes place in the terminal. On Windoze you may want to do this in Visual Studio or with Cygwin.

First I store the path-name to my dev-directory of MathLink in a variable

MROOT="/usr/local/Wolfram/Mathematica/8.0/SystemFiles/Links/MathLink/\
 DeveloperKit/Linux-x86-64/CompilerAdditions/"

Then I had to install the qhull development files. Here, I could use my package manager, while on other systems you may need to download and install it from the home page of qhull.

sudo apt-get install libqhull-dev

Then you go into the unpacked folder of qh-math and use mprep of Mathematica to process the template file

$MROOT/mprep -o qh-math.tm.c qh-math.tm

Now you can compile the sources into a MathLink program

gcc -I${MROOT} -L${MROOT} -I/usr/include/qhull -lqhull -lML64i3 -lm \
 -lpthread -lrt -lstdc++ qh-math.c qh-math.tm.c 

If you use a recent version of qhull you have to rename the variable in qh-math.c

char qh_version[] = "qh-math.c 2000/7/6";

into maybe qh_version_blub. Otherwise it clashes with a definition in the qhull lib.

The final MathLink program qh-math.exe is now ready to use in this directory.

Answer 2

Try:

Graphics`Mesh`MeshInit[]
mo = Delaunay[mypts];
faces = mo@"FaceCoordinates"

then Polygon[faces] is the triangulation.

Answer 3

As J.M. suggested in a comment, the Delaunay function in Graphics`Mesh can be used for this, though my method for getting the data out is somewhat unsatisfying.

The Delaunay function will take a list of points in 2D and return a MeshObject. Looking at the InputForm of the MeshObject one can see a list of the original points and a list of integer triplets corresponding to the points making up each triangle. Perhaps one of the functions in the package will extract this data from the MeshObject but I couldn't find it. Nor could I extract the relevant parts using normal expression manipulations. I resorted to converting the InputForm to a string in order to replace MeshObject with List.

The code is not as fast as the other solutions, but it is much quicker than the ComputationalGeometry triangulation. It runs in about 0.3 seconds on my PC.

p = meshGrid[{-1, 1, -1, 1}, 0.05];

Graphics`Mesh`MeshInit[];

extractTriangulation[mo_] := 
 ToExpression[StringReplace[ToString[InputForm[mo]], "MeshObject" -> "List"]][[3, 2]]

del = extractTriangulation@Delaunay@p;

Graphics[{EdgeForm[Black], FaceForm[], GraphicsComplex @@ MapAt[Polygon, del, 2]}]

Answer 4

Fast Delaunay triangulations are sorely missing from Mathematica. I really hope the next version will have them.

To complement the other excellent answers, here's how to access MATLAB's CGAL based Delaunay triangulator (delaunay) from Mathematica:

First, install the MATLink package. Now load it and set up the delaunay function:

Needs["MATLink`"]
OpenMATLAB[]

delaunay = Composition[Round, MFunction["delaunay"]];

Now you can simply evaluate delaunay[points] to get the triangles, represented as point indices.

Note: in MATLAB, delaunay, which is based on CGAL, is faster than delaunayn which is based on Qhull.

_{Disclaimer: I'm one of the authors of MATLink. The reason why I posted this answer is that I have referenced this question more than once and I wanted to give one more practically usable and fast method.}

Answer 5

Version 10 (which is at the moment publicly accessible through the Programming Cloud) includes the DelaunayMesh and VoronoiMesh functions.

They work for 1D, 2D or 3D point sets. Example:

points = RandomReal[1, {100,2}];

mesh = DelaunayMesh[points]

Now various properties of the triangulation can be queried using mesh["some-property"]. Available properties can be listed using mesh["Properties"].

Answer 6

I found the answer by Szabolcs really useful in finding the answer to the same question, so I thought I'd just leave the exact way to do it here in case someone stumbles across it like I did. The speed boost is very good - I'd been relying on this function to get the output I needed and for a list of about 3000 points it was taking about 20 seconds on my machine using DelaunayTriangulation and about 0.1 for the below equivalent!

points=RandomReal[1,{100,2}];
vertexadjacency = DelaunayMesh[points]["VertexVertexConnectivity"]