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:

DelaunayTriangulation[p] // Timing // First

In my computer it gives 18.533s

Matlab saves much time if does the same thing:

enter image description here


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

  • $\begingroup$ Now if only there were a way to leverage Graphics`Mesh`Delaunay[]... $\endgroup$ Commented Oct 22, 2012 at 11:28
  • $\begingroup$ 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? $\endgroup$
    – Mr.Wizard
    Commented Oct 22, 2012 at 11:32
  • $\begingroup$ @J.M. don't leave us guessing! $\endgroup$
    – Mr.Wizard
    Commented Oct 22, 2012 at 11:33
  • $\begingroup$ @J.M. I can't find any infomation about in help, how to use it $\endgroup$ Commented Oct 22, 2012 at 11:34
  • 1
    $\begingroup$ @Mr.Wizard I will use DelaunayTriangulation in many different situations,can't guarantee regular pattern, so I want to speed it up $\endgroup$ Commented Oct 22, 2012 at 11:43

6 Answers 6


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


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)]]

Mathematica graphics

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:

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


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.

  • 2
    $\begingroup$ 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 $\endgroup$ Commented Oct 22, 2012 at 13:42
  • $\begingroup$ very pretty triangles. $\endgroup$
    – s0rce
    Commented Oct 22, 2012 at 18:21
  • $\begingroup$ @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?. $\endgroup$
    – halirutan
    Commented Oct 23, 2012 at 0:48
  • 4
    $\begingroup$ Why isn't this the default and integrated with MMA? $\endgroup$
    – s0rce
    Commented Oct 24, 2012 at 1:51
  • 4
    $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Commented Apr 4, 2013 at 17:40


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

then Polygon[faces] is the triangulation.

  • $\begingroup$ Have you read the answer of @Simon? He uses kind of the same approach. Can you explain in detail, how your solution is different? $\endgroup$
    – halirutan
    Commented Oct 24, 2012 at 4:27
  • $\begingroup$ @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. $\endgroup$ Commented Oct 24, 2012 at 5:12
  • $\begingroup$ 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... $\endgroup$ Commented Oct 24, 2012 at 5:13
  • 6
    $\begingroup$ @J.M., I've discovered a few: dt /@ {"MeshElements", "Faces", "Edges", "Coordinates", "VertexCoordinates", "EdgeCoordinates", "FaceCoordinates", "Dimension"}. $\endgroup$ Commented Oct 24, 2012 at 8:29
  • 1
    $\begingroup$ @Simon, neat! How'd you manage to spelunk these? (I only found "Faces" and "Coordinates" through trial and error.) $\endgroup$ Commented Oct 24, 2012 at 9:48

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


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

del = extractTriangulation@Delaunay@p;

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

enter image description here

  • 1
    $\begingroup$ 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}. $\endgroup$
    – Mr.Wizard
    Commented Oct 23, 2012 at 15:18
  • $\begingroup$ @Mr.Wizard, Graphics`Mesh`MeshInit[] loads the package, so Delaunay can be used without a context. what happens in v7, does Delaunay[p] give you anything at all? $\endgroup$ Commented Oct 23, 2012 at 15:40
  • 1
    $\begingroup$ 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. $\endgroup$
    – Mr.Wizard
    Commented Oct 23, 2012 at 15:43

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:


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.

  • $\begingroup$ 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) $\endgroup$
    – lalmei
    Commented Oct 6, 2013 at 17:16
  • $\begingroup$ @lalmei I'm not sure what you mean by "the link crashes", but why don't you send me an email ([email protected]) 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. $\endgroup$
    – Szabolcs
    Commented Oct 7, 2013 at 1:39
  • $\begingroup$ @lalmei Please just send an email to the address above, with all details included, so we can debug and fix this issue. $\endgroup$
    – Szabolcs
    Commented Oct 7, 2013 at 16:35
  • $\begingroup$ @lalmei I haven't heard from you since your comment. Did the problem go away for you? $\endgroup$
    – Szabolcs
    Commented Oct 15, 2013 at 19:51

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


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!

vertexadjacency = DelaunayMesh[points]["VertexVertexConnectivity"]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.