Delaunay Triangles

For one of my projects, I needed the so-called Delaunay triangulation of a set of points.

Loosely put, the Delaunay triangulation is the most efficient way to draw triangles between pairs of points. Each point is connected by lines to its closest neighbours, in such a way that all line parts form triangles, and do not intersect otherwise. No triangles overlap; in fact, the surface is completely covered with one nice layer of different triangular tiles.

The Delaunay triangulation was invented in 1934 by, and named after, the Russian mathematician Boris Nikolaevich Delaunay (1890-1980). It has a lot of applications in science and computer graphics. It is often used in the graphic representation of geometrically irregularly distributed data—think weather maps or altitude maps. Its 3D-variant is important in creating virtual worlds for video games, among many other things.

Although at first glance, obtaining the Delaunay triangulation seems to be almost trivial, in fact it's a quite complicated task, the more so if you want to do it efficiently for greater numbers of points. Quite a lot of ingenious algorithms have been devised, and the field is the subject of ongoing research.

I searched the Internet for Delaunay triangulation (and the closely related Voronoi or Direchlet diagrams), and found loads of algorithms. However, most of them were unsatisfactory. Some consisted of hundreds or even thousands of lines of incomprehensive spaghetti code. Some were much too heavy, sporting lots of options to create triangulations with special constraints. Others were optimized in a way that obfuscated the underlying method completely.

So, I ended up in creating my own implementation. I present it here because, perhaps, someday someone else will be in need of a Delaunay triangulation. But also because I think it's a nice and intriguing problem in itself. Moreover, it shows some nice tricks with the Standard Template Library (STL).

The Algorithm

The overall algorithm leans heavily on an important property of Delaunay triangulations, that is: Apart from the vertices, there are no other points on or inside the circumscribed circle of any triangle. In other words, all circumscribed circles are empty.

Knowing that, and thinking hard, you can devise the following way to add one vertex to an already existing triangulation (in pseudo-code).

add_vertex(vertex)
{
   for (each triangle)
   {
      if (vertex is inside triangle's circumscribed circle)
      {
         store triangle's edges in edgebuffer
         remove triangle
      }
   }
   remove all double edges from edgebuffer,
      keeping only unique ones
   for (each edge in edgebuffer)
   {
      form a new triangle between edge and vertex
   }
}

In pictures, it looks like this:

To insert a new vertex... ... first, find which triangles have a circumscribed circle encompassing the vertex.
Remove those triangles, but remember their edges. Remove double edges, keeping only the unique ones.
Form new triangles between the remaining edges and the new vertex... ... and, finally, put those new triangles back.

Now, you just have to find a way to start the whole process. Somehow, you should create a valid Delaunay triangulation to begin with. You then successively add all your vertices.

This valid initial triangulation is easy to find. You just make a big "super triangle" that encompasses all your vertices. Of course, it means that superfluous triangles will be formed, but they can be removed afterwards easily.

The complete algorithm thus amounts to the following in pseudo-code.

triangulate()
{
   create supertriangle and add it to the triangulation
   for (each vertex)
   {
      add_vertex(vertex)
   }
   for (each triangle)
   {
      if (one or more vertices stem from supertriangle)
      {
         remove triangle
      }
   }
}

Optimizations

The algorithm can be optimized by pre-sorting the vertices along the horizontal axis. It is then possible to successively render triangles "completed," as soon as their circumscribed circles are completely to the left of the current vertex. I implemented this optimization. Its effect is hardly noticeable for small numbers of vertices, but for larger numbers it is dramatic. I do the sorting simply by using the std::set container, which keeps its object sorted automatically.

I also pre-calculated the center and the radius of the circumscribed circles. These optimizations are less important, because they don't enhance the fundamental algorithm. However, they do speed it up considerably.

Not Perfect

The algorithm is not perfect. In particular, it doesn't handle situations with more than three vertices lying on one circle very well. I found it out the hard way, after investigating the comment of Roger Labbe, below. Admittedly, I made a mistake in the first version by incorrectly assuming that no two triangles could have the same circumcircle. I corrected that, but there is still a problem in the unlikely situation that a lot of vertices happen to lie on the same circle.

Implementation

I implemented the whole thing in a C++ class Delaunay, that can be used in MFC and non-MFC environments. To keep my code clean and well structured, I used several features of the STL, such as function objects and seldom-used algorithms such as std::remove_copy_if().

I won't guarantee that it will fit seamlessly in any conceivable project, but it certainly will be easier to use than some of the other implementations around. Also, please notice that this isn't the best implementation possible. In particular, it doesn't deal with contours and holes in the set of vertices. As such, it is only a starting point.

Demo

The demo is nothing special. It's an MFC application, displaying a Delaunay triangulation of random points. Just for fun, I used some of those nifty STL function objects in the demo code as well.

References

As stated, there is much information about Delaunay triangulation on the Internet. Much of it, though, is of an academic nature, exploring the many interesting mathematical properties. Some of the more accessible sites are:



Downloads

Comments

  • Improvement for convex result

    Posted by Wolfgang Ortmann on 06/13/2013 06:02am

    The solution mentioned by airproject increases the probability to get the convex hull as border of the result. For a real solution of that problem one has to go further and use the "infinite point" Pi instead of the supertriangle. Than you can start with two points P1 and P2 and the two "infinite" Triangles P1P2Pi and P2P1Pi, consisting of two finite and the infinite point. The algorithm remains nearly the same. The circumcirle of the infinite triangles is the straight line determined by the two finite points and inside means "to the right of that line". An implementation of the Delaunay triangulation based on the original idea and including my improvement is part of the Version 7.0Beta of the image analysis library ICE, available on http://www.inf-cv.uni-jena.de/Lehrstuhl/Software/ICE/Download.html (LGPL). Regards, Wolfgang (wolfgang.ortmann@uni-jena.de)

    Reply
  • Improvement for convex result

    Posted by airproject on 04/09/2013 05:52am

    As some already noticed, there is an issue with this algorithm. In some particular cases of vertices positions, the resulting set of triangles might create concave structure. The cause of such behaviour is the initial supertriangle. In this implementation, it is reduced to the size that is just enough to enclose all vertices. This makes it highly probable to have a significant triangles formed with the supertriangle vertices, and have them deleted in the finalizing step. If the supertriangle was much bigger than the spread of the input vertces (dx, dy), the algorithm would create the triangles in a convex fashion. Simple solution found in Gilles Dumoulin's implementation is to initialize the supertriangle as follows: REAL dMax = (dx dy) ? dx : dy; REAL xMid = (xMax + xMin) / 2.0; REAL yMid = (yMax + yMin) / 2.0; vSuper[0] = Vertex(xMid - 20 * dMax, yMid - dMax); vSuper[1] = Vertex(xMid, yMid + 20 * dMax); vSuper[2] = Vertex(xMid + 20 * dMax,yMid - dMax); This mod should sort out the convexity issue.

    Reply
  • artist

    Posted by Denyse Le Blanc on 01/05/2013 10:23am

    One day while just rinsing my eyes in admiration with Robert Delaunay, I came upon Delaunay triangulation. Could you send me a download to triangulate pictures? thank you Denyse

    Reply
  • Worked great

    Posted by Jan Ekholm on 08/22/2012 12:12am

    I implemented the algorithm just to test how the Delaunay algorithm works and it seems to work perfectly. The only small issue is the need for the "super triangle(s)". You do not always have those bounding triangles when you start. Anyway, the visualization is really great and made it absolutely trivial to implement!

    Reply
  • performance time

    Posted by posthumus on 06/29/2012 03:09am

    hey the time displayed at the bottom-left corner of the screen...... is it the time of the -- total process - ( triangulation + display) or only triangulation

    Reply
  • Posted by Francois Bilodeau on 04/08/2012 12:10pm

    Real great and so fast. I compile it with Visual Studio and it worked But I want to integrate Delaunay.cpp with a project in QT I am using the QT4.8 sdk by nokia and compiling with minGW Everything goes well but for 2 lines using the same function remove_if in you code around line 273: tIterator itEnd = remove_if(workset.begin(), workset.end(), triangleIsCompleted(itVertex, output, vSuper)); I get this message from minGW: c:\qtsdk\mingw\bin\..\lib\gcc\mingw32\4.4.0\include\c++\bits\stl_algo.h:-1: In function '_FIter std::remove_if(_FIter, _FIter, _Predicate) [with _FIter = std::_Rb_tree_const_iterator, _Predicate = triangleIsCompleted]': c:\qtsdk\mingw\bin\..\lib\gcc\mingw32\4.4.0\include\c++\bits\stl_algo.h:1161: erreur : passing 'const triangle' as 'this' argument of 'triangle Reply

  • Problem with VIsual Studio 2010

    Posted by PabloLiberman on 08/18/2011 02:34am

    Hi, I try to use this demo in VIsual Studio 2010 and the remove_if makes bugs, somebody try to run it in vs2010? (it makes compilation bug)

    • I solved a problem.

      Posted by Becky on 01/28/2013 11:10pm

      I solved a problem. This problem occurred because that VC2010's and VC2008 are differ. so, we should add some code that operator '=' overriding in class triangle. const triangle& operator=(const triangle &tri;) const { if (this != &tri;) { memcpy ((void*)this, &tri;, sizeof(tri)); } return *this; } have a good time :)

      • Thank you so much!

        Posted by Michal on 02/02/2013 07:59am

        Great solution! have a good time :)

        Reply
      Reply
    • pdg

      Posted by Francois Bilodeau on 04/08/2012 01:14pm

      Hello, same problem with QT and minGW microsoft does not use original algorithm.h As a way around, modify the code like this: /* FB */ tIterator itEnd; /* FB */ for ( itEnd = workset.begin(); itEnd != workset.end(); itEnd++) /* FB */ triangleIsCompleted(itVertex, output, vSuper); /* FB */ -- itEnd; //FB tIterator itEnd = remove_if(workset.begin(), workset.end(), triangleIsCompleted(itVertex, output, vSuper)); edgeSet edges; // A triangle is 'hot' if the current vertex v is inside the circumcircle. // Remove all hot triangles, but keep their edges. // FB itEnd = remove_if(workset.begin(), itEnd, vertexIsInCircumCircle(itVertex, edges)); /* FB */ tIterator itAt; /* FB */ for ( itAt = workset.begin(); itAt != itEnd; itAt++) /* FB */ vertexIsInCircumCircle(itVertex, edges); /* FB */ itEnd = --itAt; workset.erase(itEnd, workset.end()); // remove_if doesn't actually remove; we have to do this explicitly.

      Reply
    Reply
  • polygon with points inside

    Posted by Oliver Neuse on 09/17/2010 04:13pm

    It might be, I don't understand anything (then forget this comment), but...
    try this:
    input:
    -3, -3
     3, -3
     3,  3
    -3, 3
    -1, 1
    -1, -1
     0,  0
    then you reliably get 7 triangles, instead of 8
    regards from germany

    • polygon with points inside

      Posted by Oliver Neuse on 09/18/2010 11:28am

      sorry, points were wrong, too late: -3, -2 3, -2 3, 2 -3, 2 -1, 1 -1, -1 0, 0

      Reply
    Reply
  • How to make the surface have a convex boundary?

    Posted by xiaoxiao914 on 08/03/2006 12:10am

    I noticed that your result surface may not have a convex boundary. Could you tell me how to modify the algorithm to make the surface have a convex boundary? I wish you can understand me.:-) love from China

    Reply
  • Steve Sloan -software engineering article ca 1978

    Posted by profray on 02/17/2006 11:03am

    Steve wrote a nice article on Delauney tesselation and included a source fortran code with it. I still use it as a first step in a three step program. Second step is smoothing the triangle via Clough-Tocher finite element plates (Salkauskus in his book about 1980 used it), then a brute force contouring. I still have the source and executables, produces autocad dfx output files. Rich Ray (ray@engr.sc.edu)

    Reply
  • Loading, Please Wait ...

Leave a Comment
  • Your email address will not be published. All fields are required.

Top White Papers and Webcasts

  • IBM Worklight is a mobile application development platform that lets you extend your business to mobile devices. It is designed to provide an open, comprehensive platform to build, run and manage HTML5, hybrid and native mobile apps.

  • On-demand Event Event Date: October 23, 2014 Despite the current "virtualize everything" mentality, there are advantages to utilizing physical hardware for certain tasks. This is especially true for backups. In many cases, it is clearly in an organization's best interest to make use of physical, purpose-built backup appliances rather than relying on virtual backup software (VBA - Virtual Backup Appliances). Join us for this webcast to learn why physical appliances are preferable to virtual backup appliances, …

Most Popular Programming Stories

More for Developers

Latest Developer Headlines

RSS Feeds