17

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

enter image description here

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

enter image description here

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

image

VividD
  • 3,660
  • 4
  • 26
  • 42
Tom
  • 3,416
  • 1
  • 18
  • 33
  • 1
    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 '13 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. – István Zachar Mar 27 '13 at 19:11
  • Delaunay is the way to go. You will need to use those index lists to create the list of triangle segments. – Daniel Lichtblau Mar 27 '13 at 19:15
  • 1
    You now have two upvotes and with them, the vaunted power to post images! Please post responsibly. :) – Mark McClure Mar 27 '13 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 '13 at 19:54

3 Answers3

19

First, you can generate your random points like so:

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

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

Needs["ComputationalGeometry`"];
dt = DelaunayTriangulation[pts];
dt // Column

enter image description here

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

enter image description here

Mark McClure
  • 32,469
  • 3
  • 103
  • 161
  • 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 '13 at 19:20
  • @ThomasJebbSturges Does the edit help? – Mark McClure Mar 27 '13 at 19:35
  • This may be useful too: http://mathematica.stackexchange.com/questions/277/how-to-get-actual-triangles-from-delaunaytriangulation – Szabolcs Mar 27 '13 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 '13 at 20:19
  • @ThomasJebbSturges No problem! – Mark McClure Mar 27 '13 at 20:21
  • Slight mod to assigning edges indicies edges = DeleteDuplicates[Sort /@ Flatten[Map[toPairs, dt], 1]]; halves the number of links, more importantly it ensures a set of unique links (that you visualize and) that OP will be needing in the iterations of his algorithm. – BoLe Apr 11 '13 at 12:40
14

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]

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

DelaunayTraingulation

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

Jim Arlow
  • 301
  • 3
  • 2
11

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

Mathematica graphics

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.

Community
  • 1
Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
  • 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 '13 at 20:21
  • 1
    Hmm... I can't decide whether I like this or not. I guess I better upvote, just in case. – Mark McClure Mar 27 '13 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. – Szabolcs Mar 27 '13 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 Graphics``Mesh` stuff all the time, though. – Mark McClure Mar 27 '13 at 20:32
  • @MarkMcClure It's probably better to use Graphics`Mesh 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. – Szabolcs Mar 27 '13 at 20:55
  • @Szabolcs Hey, I upvoted! I was joking when I said I wasn't sure if I liked it. – Mark McClure Mar 27 '13 at 21:50
  • @Mark But I do agree that while it's practical, it's not a nice way ... I was always annoyed by having to extract data from the output of plotting functions because there's no public API to access the underlying methods ... For example I'd love to have access to the adaptive sampling algorithm used by two-variable plotting functions. I'd also love to have a documented and fast Delaunay triangulator. There would be so many applications ... I don't want to rely too much on undocumented things for code that I may need to use a few years later. – Szabolcs Mar 27 '13 at 21:54
  • @Mark Here's yet another example where I used plotting functions for numerical calculations. It would be so valuable to have direct access to the underlying numerical methods and to be able to tune them and use them for calculations, not only for plotting. – Szabolcs Mar 27 '13 at 21:57
  • Great! It's a speed-up of 150 and more (164 for 5k points on my comp). This would make the algorithm OP is coding (recalculating triangulation over and over) useful actually. Little Q though: is anything in Graphics`Mesh even faster (for I am not familiar with the package; for example, how can I get a list of functions there?)? – BoLe Apr 11 '13 at 14:22
  • ?Graphics(backtick)Mesh(backtick)* // I think I get the function list part. – BoLe Apr 11 '13 at 14:37
  • @BoLe According to my benchmarking there's a problem with the complexity of the algorithm used for triangulation by Mma's plotting functions: it doesn't look like $O(n \log n)$ which it should be ... Around 10000-30000 points you'll hit a wall and you won't be able to go much higher. If you need to triangulate a lot, I suggest you look for ways to connect to external libraries. Try e.g. this (though Qhull isn't too speedy either) – Szabolcs Apr 11 '13 at 14:41
  • Wov. The Qhull qDelaunayTriangulation is faster for another order. My timing (5k points; same output, ie. triangle indicies): DelaunayTriangulation - 133 s; via ListDensityPlot - 0.784 s; Qhull - 0.049 s. That are speed-up factors 170 from first route to the second and another 16 to the third (or 2700 from Mma's built-in to the linked Qhull). Now you're mentioning there're even faster triangulators, like Triangle; they could be linked into Mma in the same way (compilation, MathLink), right? – BoLe Apr 11 '13 at 16:05