# How to speed up the function DelaunayTriangulation?

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.

-
Now if only there were a way to leverage GraphicsMeshDelaunay[]... – J. M. Oct 22 '12 at 11:28
In this specific case it could probably be done a lot faster given the regular nature of your points. Is that of any value to you? – Mr.Wizard Oct 22 '12 at 11:32
@J.M. don't leave us guessing! – Mr.Wizard Oct 22 '12 at 11:33
@J.M. I can't find any infomation about in help, how to use it – withparadox2 Oct 22 '12 at 11:34
@Mr.Wizard I will use DelaunayTriangulation in many different situations,can't guarantee regular pattern, so I want to speed it up – withparadox2 Oct 22 '12 at 11:43

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:

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/\


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 \


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.

-
There is an exe file for windows,how can I implement it in my program? I have never used MathLink before, is there any tutorials for me to follow – withparadox2 Oct 22 '12 at 13:42
very pretty triangles. – s0rce Oct 22 '12 at 18:21
@OleksandrR. I have version 2009.1-1 here. Except for the adjustment of the one variable name I did not change anything. Although, there were warnings about types during compilation which maybe fatal on other machines. Chat?. – halirutan Oct 23 '12 at 0:48
Why isn't this the default and integrated with MMA? – s0rce Oct 24 '12 at 1:51
A correction: MATLAB uses Qhull for delaunayn and related functions: N-dimensional Delaunay triangulations. For 2D (and maybe for 3D) it uses CGAL, which is MUCH faster than Qhull: benchmark delaunay vs delanayn in MATLAB. There are libraries which are again considerably faster than CGAL, e.g. Triangle. – Szabolcs Apr 4 '13 at 17:40

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}];

-

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"].

-

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.

-
Unfortunately I only have matlab2009 to try, but it seems I can only use the function once before the link crashes, and only with less than 1500 points. Is that a known current limitation ? (Mac OSX Mountain Lion, Mathematica 8) – lalmei Oct 6 '13 at 17:16
@lalmei I'm not sure what you mean by "the link crashes", but why don't you send me an email (matlink.m@gmail.com) with a complete description of the problem and complete version info (Mma, MATLAB, MatLink, your OS, 32/64 bit, student version or not)? We don't really support that MATLAB version but it may be possible to fix the problem. – Szabolcs Oct 7 '13 at 1:39
@lalmei Please just send an email to the address above, with all details included, so we can debug and fix this issue. – Szabolcs Oct 7 '13 at 16:35
@lalmei I haven't heard from you since your comment. Did the problem go away for you? – Szabolcs Oct 15 '13 at 19:51

Try:

GraphicsMeshMeshInit[]
mo = Delaunay[mypts];
faces = mo@"FaceCoordinates"


then Polygon[faces] is the triangulation.

-
Have you read the answer of @Simon? He uses kind of the same approach. Can you explain in detail, how your solution is different? – halirutan Oct 24 '12 at 4:27
@hal, OTOH, Ulises provides a slightly more convenient way to extract the coordinates of the triangulation. After executing dt = Delaunay[meshGrid[{-1, 1, -1, 1}, 0.05]];, one can now do either Graphics[GraphicsComplex[dt["Coordinates"], {EdgeForm[Black], FaceForm[], Polygon[dt["Faces"]]}]] to get a picture like the one in Simon's answer, or Graphics[MapIndexed[{ColorData[29, First[#2]], Polygon[#1]} &, dt["FaceCoordinates"]]] to get a picture like the one in your answer. – J. M. Oct 24 '12 at 5:12
Ulises, if you could be kind enough to share and elaborate on the properties of a MeshObject[] that can be exposed, your answer would be certainly more useful... – J. M. Oct 24 '12 at 5:13
@J.M., I've discovered a few: dt /@ {"MeshElements", "Faces", "Edges", "Coordinates", "VertexCoordinates", "EdgeCoordinates", "FaceCoordinates", "Dimension"}. – Simon Woods Oct 24 '12 at 8:29
@Simon, neat! How'd you manage to spelunk these? (I only found "Faces" and "Coordinates" through trial and error.) – J. M. Oct 24 '12 at 9:48

As J.M. suggested in a comment, the Delaunay function in GraphicsMesh 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];

GraphicsMeshMeshInit[];

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

del = extractTriangulation@Delaunay@p;

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


-
I can't get your method to work with version 7. FWIW typically you can use patterns to extract data from pseudo-atomic objects, e.g. SparseArray[Range@9] /. _[data__] :> {data}. – Mr.Wizard Oct 23 '12 at 15:18
@Mr.Wizard, GraphicsMeshMeshInit[] loads the package, so Delaunay can be used without a context. what happens in v7, does Delaunay[p] give you anything at all? – Simon Woods Oct 23 '12 at 15:40
I didn't realize that. I get a MeshObject but InputForm only returns MeshObject[None, None]. Interestingly I see something if I use GeometryPlot[Delaunay@p] so there is something there, but neither InputForm nor /. _[data___]` appear to access it in v7. – Mr.Wizard Oct 23 '12 at 15:43