17
$\begingroup$

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

Example:

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

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

tri = DelaunayTriangulation[points];
showTriangulation[tri]

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.

|improve this question
$\endgroup$
  • 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 Jan 19 '12 at 14:29
15
$\begingroup$

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

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

Graphics[GraphicsComplex[points, {EdgeForm[Black], FaceForm[],    
  Polygon[triangles[points]]}]]

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

Needs["ComputationalGeometry`"];
triangles2[points_] := Module[{tr, triples},
  tr = DelaunayTriangulation[points];
  triples = Flatten[Function[{v, list},
      Switch[Length[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]]]]

Edit

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.

|improve this answer
$\endgroup$
  • $\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 Jan 19 '12 at 19:29
  • $\begingroup$ ... the package will be loaded, but ComputationalGeometry`DelaunayTriangulation will be shadowed by Global`DelaunayTriangulation. $\endgroup$ – Szabolcs Jan 19 '12 at 19:30
  • $\begingroup$ Please see this question (finally found it). $\endgroup$ – Szabolcs Jan 19 '12 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 Jan 19 '12 at 23:05
  • $\begingroup$ Can you edit the answer please? $\endgroup$ – Szabolcs Jan 19 '12 at 23:31
6
$\begingroup$

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.

|improve this answer
$\endgroup$
4
$\begingroup$

Update

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:

Needs["ComputationalGeometry`"]  

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

Graphics@GraphicsComplex[pts,
{{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:

Needs["ComputationalGeometry`"]  
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.

triangles2[pts]  
(* {{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.

|improve this answer
$\endgroup$
3
$\begingroup$

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

|improve this answer
$\endgroup$
  • $\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 Jan 19 '12 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 Jan 19 '12 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 Jan 19 '12 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 Jan 19 '12 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 Jan 19 '12 at 21:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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