Up until now I have been using the excellent Triangle.net library to carry out my trianglulation for models I generate myself. Models I buy from turbosquid are already triangulated so dont require any attention.

I’ve been running into problems trying to constrain the trinangulation algorithm to force certain triangle edges to be used; for instance in the case of roads, I definitely want the road edge line to be used irrespective if the triangulation algorithm could make a better (more equalateral) triangle.

The Triangle.net library does support a level of constraints, but its quite poorly documented and I wanted to understand much better the whole concept of Delaunay Triangulation and write a post on the practical implementation. There are quite a few C and C++ libraries out there, but little explanation of how they work or why they work the way they do – most of them concentrate on performance optimisation tricks, but I wanted just to be able to understand the basic process before getting into performance.

**Step 1 – Definitions**

A Delaunay Triangulation defines a set of triangles which join a set of points such that all the points belong to a triangle, and no triangle edges cross. A Delaunay Triangulation also fulfills the criteria that the triangles are as equalateral as possible. This is important in graphics programming because “long thin” triangles create visual artefacts and suffer from mathematical inaccuracies when carrying out hit tests etc. So the more equalateral the triangle is, the better it looks, the more evenly it tessellates and the easier it is to decimate (make more or less dense accordingly)

Heres a Delaunay triangulation of a mesh of 8 points.

**Step 2 – The Method**

Starting with an arbitrary set of 2D Points, collect them into a list (or other array type) structure.

Calculate an additional 3 points (we’ll name “Pf”) which, when made into a triangle, are further away from all the original points than they are to each other. This is called the Super Triangle – you could theoretically use infinity points to describe this, but you’ll get floating point inaccuracies from this, so just calculate the bounding box of the original point set and double it. This forms the first Triangle.

For the first Point (arbitrarily point 1) insert it into the super Triangle and split it into three resultant triangles. You now have three triangles and have used up one of your Points.

Each of your Triangles must be stored in some kind of array, along with its Circumcircle; the center and radius of a circle on whose circumference all the points of the triangle fall. Its clear that them more equalateral the Triangle the tighter the Circumcircle will be to the Triangle area – if the triangle is very un-equalateral, the Circumcircle will need to be much larger.

For every subsequent Point

- Gather all the Triangles where the Point falls within their Circumcircle – i.e. the distance of point to the Circumcircle center is less than the Circumcircle radius. (this is where a lot of optimisation is concentrated on – for my case, I just iterate all the Triangles testing each one in turn). We’ll call this set of Triangles T0.
- For all the triangles in T0 work out which of their edges are internal – i.e. shared between other members of the set T0. These edges we’ll call E0.
- Delete all T0 triangles from the set of all Triangles.
- Join all the points of T0 that do not describe any of the E0 edges to the Point we are inserting, to form a set of new Triangles T1.
- Add all the new Triangles T1 to the set of all remaining triangles.

This process is illustrated having selected Point 0 and found two of the existing triangle circumcircles overlap Point 0; the green lines show the internal shared edges of those triangles (these get deleted). Consequently we create new triangles 9-0-10, 9-1-0, 8-1-0 and 8-0-10

The key to the algorithm is that you just keep repeating the same basic move until you run out of points. You will end up with this;

Once you’ve processed every Point.

- Delete all Triangles that have as one of their points a point from the set “Pf” (8,9 and 10 in this example)
- What you have left is the Delaunay triangulation of your points. Almost….

Whereas the mesh is nice and is quite likely to be a Delaunay mesh, there are circumstances where the algorithm will join two triangles along an edge where a a shorter alternative edge existed (in this case the quad 0,1,7,6 has a rather acute angle on the corner 0,1,7, and would be better expressed as triangles 0,6,7 and 0,7,1. For a high quality “true” Delaunay mesh we should

- Look at each pair of triangles that share a side.
- If that resultant polygon is convex (you can test this by checking all the internal angles – they should all be less than 180 degrees – if they are not, a concavity exists) then calculate the length of both possible configurations of shared edge;
- If the alternative shared edge is shorter than the existing one, flip the triangle definitions around so they share the shorter edge, delete the existing triangles, and insert the new pair.

**Step 3 – Constraints**

The above method will generate a convex polygon – but you have no control over the triangles generated. In many cases (not least of which is defining a concave hull) there are certain edges you want to pre-define which must be present in the eventual mesh, even though they are a sub-optimal Delaunay mesh. For example there may be some reason why we need triangle 3,4,6 to be represented in the eventual mesh, even though it would never get created using Delaunay

- Once triangulation has finished, process a list of Constrained Edges that you want to be in your eventual mesh.
- For each Constrained Edge (CE1) generate a circle and identify all the Triangles whose circumcircle intersects it. (this is an optimisation step – you could simply iterate all triangles). This is set T2.
- For each Triangle in T2 test each edge and if it crosses CE1 its add it to set T3 and delete the Triangle from the main mesh. In this case triangles 3,7,4 and 7,4,6 both intersect the edge 3,6 and get deleted from the mesh.
- Use a divide-and-conquor algorithm to generate a mesh joining all the triangles in T3 with the Constrained Edge CE1.

Ok – that last step was a bit of a mouthful. Whats a divide-and-conquor algorithm ? All the points of triangles T3 must lie either on one side or other of the CE1 line (line 3,6 in this example). So divide the points into their two constituent piles (Pile1 and Pile2) and add the two points of the constraint line CE1. Since both sets came from a Delaunay triangulation they will form two seperate convex polygons with edge CE1 forming one of the outside edges – just Delaunay both polygons and add the resultant triangles back into the main mesh.

This is hard to picture with just one triangle, so heres a better illustration; the blue line was a Contraint on a Delaunay mesh – all the red points belonged to triangles that intersected the blue line. The points are sorted into Pile1 (red) and Pile2 (yellow). Each set is triangulated using standard Delaunay to produce two convex polygons, but seperately. The subsequent mesh of both halves is no longer Delaunay, but honours the Constraint line.

**Step 3 – Data, Optimizations.**

You will need to work with double-precision data for all the mathematical elements of this algorithm – floats are just too imprecise. Many examples are quite prescriptive about the kinds of list structures to be used to efficiently parse the many triangles and points you will generate – you can experiment with these, but from the point of view of learning the method, I just stuck with a simple List<> and a couple of others; a Triangle class and a Vector2Double struct. I didn’t bother storing the triangle adjacencies, but did store the Triangle circumcircles in the Triangle class.

In terms of optimisations, a very good method of accellerating the initial Delaunay triangulation is to process all the points in spatial order – either top to bottom or left to right. When each triangle is generated and stored you can check the list of triangles to see if their circumcircles can need be considered any longer for future points – all future points will be further away in either the X or Y coordinate (depending on your choice of spatial order) so triangles whose circumcircles didn’t intersect for the current point and are left/below it, cannot now match any future point, so you no longer need to worry about them. This means the set of triangle/circle matches to do gets smaller and smaller for each subsequent point.

If that resultant polygon in convex == “If that resultant polygon is convex”?

Maybe. Great explanation. I love your blog

Thanks for the correction. Glad you enjoy the writing.