The ComputationalGeometry package has a DelaunayTriangulation[] function. It returns a list of points connected to each point, ordered counterclockwise.


showTriangulation[tri_, opt : OptionsPattern[]] := 
 Graphics[GraphicsComplex[points, Line /@ Thread /@ tri], opt]

points = Tuples[Range[0, 5, 1], 2];

tri = DelaunayTriangulation[points];

Mathematica graphics

Question: How can I obtain a list of all triangles instead?

My first naive first try doesn't work correctly:

makeTriangles[points_] := 
 Flatten[Function[{p, list}, 
    Prepend[#, p] & /@ Partition[list, 2, 1, {1, 1}]] @@@ 
   DelaunayTriangulation[points], 1] (* doesn't work *)

GraphicsComplex[points, Line@makeTriangles[points]] // Graphics

Mathematica graphics

Why does it give an incorrect result?

This function simply takes all points $\{A, B, C, D, \ldots\}$ connected to a point $P$, and constructs the triangles $PAB, PBC, PCD, \ldots$. Since $A,B,C,...$ are in counterclockwise order, I assumed all these would be valid triangles. But take the following case:

Mathematica graphics

$PAB$ will not be a valid triangle, even though $ABC$ are in counterclockwise order.

  • 1
    $\begingroup$ I already realized that I don't actually need this for my work, but I though some might find the question interesting. $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 14:29

4 Answers 4


It might be easier to use TriangularSurfacePlot3D to find the Delaunay triangulation of the points. For example,

triangles[points_] := Module[{pl},
  pl = TriangularSurfacePlot[ArrayPad[points, {{0, 0}, {0, 1}}]];
  Cases[pl, Polygon[a_] :> Flatten[(Position[points, #[[{1, 2}]]] & /@ a)], 

Graphics[GraphicsComplex[points, {EdgeForm[Black], FaceForm[],    

produces this:

Mathematica graphics

If we use the approach in the original question, each valid triangle should appear exactly trice in the list, so a way to extract the triangles from the result of DelaunayTriangulation would be

triangles2[points_] := Module[{tr, triples},
  tr = DelaunayTriangulation[points];
  triples = Flatten[Function[{v, list},
        (* account for nodes with connectivity 2 or less *)
        1, {},
        2, {Flatten[{v, list}]}, 
        _, {v, ##} & @@@ Partition[list, 2, 1, {1, 1}]
    ] @@@ tr, 1];
  Cases[GatherBy[triples, Sort], a_ /; Length[a] == 3 :> a[[1]]]]


I've moved Needs["ComputationalGeometry`"] outside of the function definition of triangles. As Szabolcs correctly remarked in his comments, putting Needs["ComputationalGeometry`"] inside the definition will cause shadowing problems because of the creation of the symbol Global`DelaunayTriangulation

Edit 2

Apparently ListDensityPlot uses a Delaunay triangulation as well, and is much faster than TriangularSurfacePlot, so the first part of this answer could be made much more efficient by rewriting it according to

triangles1[points_] := Module[{pl},
  pl = ListDensityPlot[ArrayPad[points, {{0, 0}, {0, 1}}]];
  Cases[pl, Polygon[a_] :> a, Infinity][[1]]]

Note that since ListDensityPlot returns a GraphicsComplex and it keeps the order of the points the same, the index lists for the polygons can extracted directly from the plot without having to lookup the indices of the vertices in points.

  • $\begingroup$ A comment on putting the Needs inside the function definition: this will only work if you explicitly include the context names in the symbols, i.e. use ComputationalGeometry`DelaunayTriangulation. This is because the context of symbols is determined at parse time. When the function definition is evaluated in a fresh kernel, the package is not loaded yet, so there's no DelaunayTriangulation symbol yet. But as soon as Mathematica sees this symbol in the definition, it creates it---in the Global` context. When the function is called for the first time, ... $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 19:29
  • $\begingroup$ ... the package will be loaded, but ComputationalGeometry`DelaunayTriangulation will be shadowed by Global`DelaunayTriangulation. $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 19:30
  • $\begingroup$ Please see this question (finally found it). $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 19:37
  • $\begingroup$ @Szabolcs Thanks, I didn't realise about Needs[...] inside function definitions. I guess I missed it because I had already loaded the package in the same session. $\endgroup$
    – Heike
    Commented Jan 19, 2012 at 23:05
  • $\begingroup$ Can you edit the answer please? $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 23:31

In Mathematica 10, you can use DelaunayMesh on a set of points. This returns a MeshRegion. You can use MeshCoordinates to return a list of coordinates of the points (should be the same as the initial set of points) and then MeshCells to return the triangles.

See Interactive Computational Geometry for more details.



It seems that the simplest way to find Delaunay triangles is to read them off from Voronoi diagram. Here is only two lines of code:


listSow[{pt_, list_}] := Scan[Sow[pt, #] &, list];

DelaunayTriangles[pts_] := Select[
Last@Reap@Scan[listSow, Last@VoronoiDiagram[pts]], (Length[#] == 3) &]

Application is as follows:

npts = 100;
pts = Transpose@{RandomReal[{0, 2}, npts], RandomReal[{0, 1}, npts]};

{{FaceForm@None, EdgeForm@LightGray,Polygon/@ DelaunayTriangles[pts]}, {Red, Point@Range@npts}}]

enter image description here

Obsolete part

I just would like to notice that the part of Heike's answer, there DelaunayTriangulation is used, is not completely correct. Let me demonstrate it with a simple example:

pts = {{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}, {1/2, 1/(2 Sqrt[3])}};  
delaunayGraph = DelaunayTriangulation[pts];

Let's now find the graph edges

edges = Join@@Table[Thread[{v, Select[delaunayGraph[[v, 2]], # > v &]}], {v, 1, Length[pts] - 1}]; 

and depict them

Graphics[{Line@pts[[#]] & /@ edges, {Red, PointSize -> Medium, Point[pts]},  
FontSize -> 14, MapIndexed[Inset[First@#2, #1, {Right, Bottom}] &, pts]}]

enter image description here

It is easy to see that the triangle {1,2,3} appears also thrice in Heike's list.

(* {{1, 2, 4}, {1, 4, 3}, {1, 3, 2}, {2, 3, 4}} *)

Probably, the triangular convex hull is the only such example. In any case here is my "brute force" solution:

area[v1_, v2_] := Det[Subtract @@@ {v1, v2}]  

findSmallest[{a_, b_}, list_] := First@SortBy[list, Abs@area[pts[[{#, b}]], pts[[{a, b}]]] &]

findPartners[{a_, b_}, list_] := findSmallest[{a, b}, #] & /@  
GatherBy[list, Sign@area[pts[[{#, b}]], pts[[{a, b}]]] &] 

tri[a_, b_] := Sequence @@ (Sort[{a, #, b}] & /@ 
 findPartners[{a, b}, Intersection[delaunayGraph[[a, 2]], delaunayGraph[[b, 2]]]]);

triples = Union[tri @@@ edges]  
(* {{1, 2, 4}, {1, 3, 4}, {2, 3, 4}} *)

The idea is as follows: firstly, for an every edge {a,b} we find additional vertices attached to it (they form triangle with the edge), see Intersection in tri. Then we gather these vertices according to parts of the plane divided by the edge linear continuation, see GatherBy in findPartners. Finally, we find the triangle of the smallest area in each part of the plane, see SortBy in findSmallest.


How about

triangles[vertices_] := 
    DeleteDuplicates[Sort /@ Flatten[Partition[#[[2]], 2, 1] /.
                     {x_Integer, y_Integer} :> {#[[1]], x, y} & /@ vertices, 1]]

(using Partition[] to make lists of triangles with three vertices instead of four)

In[95]:= triangles[tri]
Out[95]= {{1,7,8},{1,2,8},{2,8,9},{2,3,9},{3,9,10},{3,4,10},{4,10,11},{4,5,11},{5,11,12},{5,6,12},{7,13,14},{7,8,14},{8,14,15},{8,9,15},{9,15,16},{9,10,16},{10,16,17},{10,11,17},{11,17,18},{11,12,18},{13,19,20},{13,14,20},{14,20,21},{14,15,21},{15,21,22},{15,16,22},{16,22,23},{16,17,23},{17,23,24},{17,18,24},{19,25,26},{19,20,26},{20,26,27},{20,21,27},{21,27,28},{21,22,28},{22,28,29},{22,23,29},{23,29,30},{23,24,30},{25,31,32},{25,26,32},{26,32,33},{26,27,33},{27,33,34},{27,28,34},{28,34,35},{28,29,35},{29,35,36},{29,30,36}}

Graphics[{FaceForm[None], EdgeForm[Thickness[Medium]],
   (* the triangles *) GraphicsComplex[points, Polygon[triangles[tri]]],
   (* and their labels *) MapIndexed[Text[#2[[1]], Mean[#1 /. Thread[Range@Length@tri -> points]]] &, triangles[tri]]}]

Mathematica graphics

  • $\begingroup$ This is the same what I did (with minor variations) and has the same problem. Not all of the points that are connected to a given point need to be connected together, not even if they are adjacent in the counterclockwise order. $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 17:15
  • $\begingroup$ Just test on points = Tuples[{0, 1}, 2]; tri = {{1, {2, 3, 4}}, {2, {1, 4}}, {3, {4, 1}}, {4, {2, 1, 3}}}; DelaunayTriangulationQ[points, tri]. This is a manually constructed example, but it is a valid triangulation, and I think there's no guarantee something like this won't be returned by DelaunayTriangulation[]. $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 17:15
  • $\begingroup$ I should have explained a bit better what was wrong with my original (flawed) implementation. I have now updated the question to avoid confusion. $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 17:32
  • $\begingroup$ I tried the DelaunayTriangulation directly on your points, and I get a different triangulation than your example: same points in each list, but the last list {{4,{2,1,3}} is a cyclic permutation in the numbers of the first list {{1,{3,4,2}}). And this one works with my code, while yours does not. I wonder if it possible to prove whether DT will always return lists whose repeated elements are cyclic permutations only. (If that makes any sense!) $\endgroup$
    – JxB
    Commented Jan 19, 2012 at 21:27
  • $\begingroup$ Yes, I constructed the example by hand to show that your function is not robust. This data is also a valid triangulation, and the points are still in counterclockwise order. I cannot trust that DelaunayTriangulation won't ever return such a result because according to the documentation any cyclic permutation of the points is a possible output (i.e. it mentions no restrictions). $\endgroup$
    – Szabolcs
    Commented Jan 19, 2012 at 21:32

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.