# Create triangular mesh from random list of points

I have a list of points. I would like to take these points and create a mesh of triangles from them, making sure triangles don't overlap. So here's a list of points:

p0 = Transpose[{RandomReal[{0, 10}, {100}], RandomReal[{0, 12}, {100}]}];
ListPlot@p0 Now I've managed to take the first point in the list, find its two nearest neighbours and construct a triangle from this:

Graphics[{Thick, Red, Line[Append[Nearest[p0, p0[[1, All]], 3], p0[[1, All]]]]}] What I'd not like to do is from this starting point keep connecting points to make a mesh where all the triangle corners are at one of the points and no triangles overlap. Any suggestions on how to do this?

--Below may be irrelevant--

I imagine DelaunayTriangulation might come into this, however I'm not sure how. Also when I run it I don't understand what it returns:

DelaunayTriangulation[p0] // MatrixForm • p.s. It doesn't really seem productive to forbid new users from submitting images and posting more than 2 links. It just makes it harder to explain the question! So I apologise for the links to images and the third link at the bottom which requires removing [take this out] for the link to work :)
– Tom
Mar 27, 2013 at 18:29
• Don't worry about the images, people here are usually helpful on solve this issue. I've shortened your Q as it cointained some irrelevant overhead, please feel free to roll back the edit if you find it incorrect. Mar 27, 2013 at 19:11
• Delaunay is the way to go. You will need to use those index lists to create the list of triangle segments. Mar 27, 2013 at 19:15
• You now have two upvotes and with them, the vaunted power to post images! Please post responsibly. :) Mar 27, 2013 at 19:17
• @MarkMcClure haha fantastic. Just wrapping my brain round your answer. Think I have enough information to figure it out. Thanks :)
– Tom
Mar 27, 2013 at 19:54

First, you can generate your random points like so:

SeedRandom;
pts = RandomReal[{0, 12}, {100, 2}];


The DelaunayTriangulation command returns an adjacency list representation of the triangulation.

Needs["ComputationalGeometry"];
dt = DelaunayTriangulation[pts];
dt // Column This says that the first point should be connected to the 2nd, the 24th, etc. Given {u, {v1,v2,v3,___}}, we need a toPairs function to form {{u,v1},{u,v2},{u,v3},___}. We then need to map toPairs onto the triangulation and Flatten that result one level. This is all accomplished as follows.

toPairs[{m_, ns_List}] := Map[{m, #} &, ns];
edges = Flatten[Map[toPairs, dt], 1];


Finally, we visualize using a GraphicsComplex.

Graphics[GraphicsComplex[pts, {Line[edges],
Red, PointSize[Large], Point[pts]}]] • Amazing. Any chance you could briefly describe each step and I'll try and understand what's going on? Thank you. (Each step from "toPairs" onwards)
– Tom
Mar 27, 2013 at 19:20
• @ThomasJebbSturges Does the edit help? Mar 27, 2013 at 19:35
• This may be useful too: mathematica.stackexchange.com/questions/277/… Mar 27, 2013 at 19:47
• @MarkMcClure It was tough but I've managed to comprehend it. Go me. What a concise way of doing it! Thank you kindly :)
– Tom
Mar 27, 2013 at 20:19
• @ThomasJebbSturges No problem! Mar 27, 2013 at 20:21

There are some new functions in Mathematica 10 that make this very easy:

r = {{-6, 6}, {-6, 6}};
pts = RandomSample[Permutations[Range[-5, 5], {2}], 10];
Grid[{
{"The sites", "Delaunay trianguation", "Voronoi diagram"},
{
Graphics[{Red, Point[pts]}, PlotRange -> r],
Show[dm = DelaunayMesh[ pts], Graphics[{Red, Point[pts]}],
PlotRange -> r],
Show[VoronoiMesh[ pts], Graphics[{Red, Point[pts]}], PlotRange -> r]
}
}, Frame -> All]

MeshCoordinates[ dm ]
MeshCells[ dm , 2]
MeshCells[ dm , 2][[ All, 1]] MeshRegion

{{-5., -4.}, {3., -4.}, {5., -2.}, {4., 0.}, {-4., -1.},
{-3., 2.}, {0., -1.}, {3., -5.}, {2., 4.}, {5., -5.}}

{Polygon[{5, 1, 7}], Polygon[{6, 7, 9}], Polygon[{7, 6, 5}],
Polygon[{9, 7, 4}], Polygon[{1, 8, 7}], Polygon[{8, 10, 2}],
Polygon[{2, 3, 4}], Polygon[{3, 2, 10}], Polygon[{2, 4, 7}],
Polygon[{8, 2, 7}]}

{{5, 1, 7}, {6, 7, 9}, {7, 6, 5}, {9, 7, 4}, {1, 8, 7}, {8, 10,
2}, {2, 3, 4}, {3, 2, 10}, {2, 4, 7}, {8, 2, 7}}


So you use DelaunayMesh to create a MeshRegion from the point set, and then you can use MeshCells as shown to get the triangles. MeshCells gives you triples of indexes into the MeshCoordinates.

I took the above code from Interactive Computational Geometry (disclaimer - I am the author).

It seems you are asking for the Delaunay triangulation.

There's a function for this in the Computational Geometry package, which Mark described.

Another, usually much faster option is using ListDensityPlot:

ldp = ListDensityPlot[ArrayPad[p0, {0, {0, 1}}], Mesh -> All,
ColorFunction -> (White &)] You can extract the polygons from this graphic if needed.

Cases[ldp, Polygon[idx_] :> idx, Infinity]


This will return the triangles as point index triplets.

You can also use the undocumented function ListDensityPlot relies on, if you wish.

• This looks like a neat way of doing it. Thank you for the response. It took me long enough to wrap my small brain around Marks answer so I'll stick with what I (now) know!
– Tom
Mar 27, 2013 at 20:21
• Hmm... I can't decide whether I like this or not. I guess I better upvote, just in case. Mar 27, 2013 at 20:22
• @Mark I don't like it, but it's definitely much faster than using the package function. I used it for a while with this, because DelaunayTriangulation was too slow. Mar 27, 2013 at 20:25
• @Szabolcs It's just that ListDensityPlot is about the last thing I would've thought of. Back in the day, this returned a DensityGraphics object. I guess it now returns a GraphicsComplex with polygons. I use the GraphicsMesh stuff all the time, though. Mar 27, 2013 at 20:32
• @MarkMcClure It's probably better to use GraphicsMesh instead. After all the structure returned by ListDensityPlot` is also undocumented and it may change at any time. So my effort to avoid mentioning undocumented functions was a bit misguided I guess. Mar 27, 2013 at 20:55