15,123,509 members
Articles / Game Development
Article
Posted 12 Nov 2012

90.1K views
62 bookmarked

# Delaunay Triangulation For Fast Mesh Generation

Rate me:
Practical on the fly fast mesh generation from arbitrary points.

## Introduction

Ten years ago, computing meshes for surfaces in real time for surfaces wasn't realistic, and having a customizable source code module wasn't available either. There are faster versions, but they are large implementations and they are hard to read and modify.  With modern computers and modern languages, not only can simple meshes be generated in real time, the code can be written in ways easy to understand. It is the authors opinion that computational geometry and computer vision are entering a new age where real time processing is realistic for a growing set of problems.  The code is designed to be reused such that a vertex can be used to both generate the mesh, but also be apart of other data structures and track other aspects at a given point.

To form a surface that can be shaded by a graphics card, the goal is to create a list of vertexes, and a list of triplet indexes of the vertexes that form counter clockwise triangles. Whether it is to visualize a differential equation solved by finite element modeling (FEM) or to view a shape of an object in 3D, the first step is to form the triangle mesh.  Usually selecting a set of points is easy, For simple objects like a cylinder or sphere, the point generation leads directly to triangle generation. Generating a mesh from an arbitrary set of points is where Delaunay's Triangulation proves valuable.  The goal of the code wasn't to compute the ideal triangulation, but instead produce a “good enough” solution for most practical 3D or FEM problems and do so quickly. Thus the code recursively improves the triangulation after adding each point to either a recursion count is reached or it is an ideal Delaunay triangulation. This limits how long a mesh can take to compute and avoids infinite recursion should there be a software bug.

## Background

Delaunay's Triangulation is named for Boris Delaunay. The general idea is to form a mesh where each triangle's three points lie on the edge of a circle that doesn't contain any other point. This means given any two adjacent triangle's (quadrilateral) the sum of the angles opposite the dividing line are less than 180 degrees.

Triangle mesh a,b,d and b, c, d don't meet the Delaunay condition. But b,c,a and c,d,a do. This is because angle bcd + dab > 180, but abc+cda < 180. This forces the mesh to have triangles that tend to be as close to evenly spaced as possible.

The Bowyer Watson algorithm iteratively adds a point at a time to a mesh to transform it from one Delaunay mesh to another. The general algorithm works by deleting all triangles whose circumscribed circle contains the new point and reform the mesh. The code provided takes a simpler approach.

From top to bottom, the triangles are flipped to meet the Delaunay condition with 3 points:

This does the same thing but with seven points. Note how at first the triangles are mostly slivers and other triangles are huge.

After one level of flipping per insertion the meshing is more regular (below).

After a second round of flipping, the effects of the first flip are then flipped as needed (below).

## Design

### Mesh Generation

The Bower-Watson algorithm is a bit complex because:

• Adding one point may cause a set of triangles to be deleted, and reformed.
• Reforming the set at once explores forming many possible triangulations.
• One must search for the affected triangle's.
• There isn't a pattern for adjusting the mesh in an ordered fashion.

The “hull” based techniques that are also commonly used spiral out from a point and reform the triangulation from the inside out. This approach is nice because it takes an ordered approach to computing and recomputing the mesh. The disadvantage as the software is recomputing the circle for each triangle and sorting through the triangles in real time.

The code provided takes a hybrid approach that forms a linked list of triangles. Each triangle consists of three vertexes, and each side, there is a pointer to the adjacent triangle (or null). The code has the following advantages:

• Each triangle has a link to its neighbors.
• Recursion explores outward if a triangle is altered, its neighbors may be affected (no circles).
• Triangles aren't not generally deleted, they are added or modified.
1. Compute the bounding rectangle:
2. Note that a->b->c is counter clockwise, so if it had a face rendered by say direct3D, it would face out from the page. Triangle c->b->a is clockwise so it's face is toward the page. As long as the direction is maintained, the start point isn't relevant, so bdc is the same as dcb or cbd. This concept in abstract algebra is often called commutation.

In memory at this point:

Triangle 1 = a,b,c, and side b,c = triangle 2.

Triangle 2 = bdc, and c,b = triangle 1.

3. Find the triangle containing the new point (for example, abc), and replace the old triangle 1 with three new triangles.

4. Repair the pointers:

5. abe with side be should point to ebc.

Abe with side ea should point to aec.

Etc.

Then repair those old adjacent triangles that used to point to abc:  bdc side bc used to point to abc, make it point to ebc.

6. Examine each triangle and it's outward adjacent:. For example bec shares a side with bdc.

7. Examine the angles across from the adjacent sides. If greater than 180 degrees, flip the adjacent line.

9. Reform the old triangles into the new triangles, fix the triangle edges.

10. For each triangle and it's adjacent triangle repeat steps 5-8 recursively.

### Triangle Contains a Point

The key to the algorithm is finding the correct triangle to modify.

Classically, a triangle contains a point if it is on the interior side. If a line passes through two points, the line is often expressed as y = mx + b, where m is the slope, and b is the x=0 intercept. The area above the line can be expressed as y > mx+b. So, I a triangle has edge s, r and the other point is t. Then if test point p, is on the same side of the inequality as the point t, then the test point p is on the same side of the line as the other vertex.

This means to compute if a point p is inside, the code compares to see if it's on the same side as the vertex for an edge. That is a lot of computations. Caching the information speeds up the comparison but not enough.  Most triangles don't contain the point, and the point is not anywhere near the triangle. To accelerate the computations, only points that are near the triangle need to be tested using the classical method. The software computes the distance between each vertex and the triangle's center. If the point is not closer than this distance it does not need to be tested.

### Acceleration

The triangle objects cache information calculated on demand to form the mesh. When a triangle is "flipped" with it's neighbor, two of it's side's change but two of the points remain. When a vertex moves, it invalidates the existing information about those edges.

For example, the center of a triangle changes if any vertex changes, setting points A, B, C sets the calculated to false, asking for the center calculates it ans sets the flag.

C#
```Vertex m_Center;

public Vertex Center
{
get
{
if (m_centerComputed) return m_center;
m_center = new Vertex(
(A.X + B.X + C.X) / 3f,
(A.Y + B.Y + C.Y) / 3f,
(A.Z + B.Z + C.Z) / 3f);

float delta = m_center.DeltaSquared(A);
float tmp = m_center.DeltaSquared(B);
delta = delta > tmp ? delta : tmp;
tmp = m_center.DeltaSquared(C);
delta = delta > tmp ? delta : tmp;
FarthestFromCenter = delta;
m_centerComputed = true;

return m_center;
}
} ```

Once the center is computed the software can reuse this information and it is only then recalculated on demand.

## Using the code

The form provides a demonstration that generates a semi-clustered seeded set of points. The mesh class has the following API:

• `Points` – a list of Vertex objects used in the mesh.
• `Facets` – a list of triangle facets of the mesh.
• `Bounds` – Rectangular bounds used to form the initial square region.
• `Recursion` – How many times to iterate out from the initial changed triangle. For most meshes under 100 the odds of needing to recursively re-triangulate more than 5 or 6 times is very low. This allows a user to balance computation time against the quality of the mesh.
• `Compute()` - Compute from a list of vertexes the mesh.
• `Append()` - Add in a new vertex to the triangulation.
• `Setup()` - Set the bounds.
• `Draw()` - Draw to a `System.Drawing.Graphics` object the mesh for debug.
• `GetVertexIndicies()` - Get the set of indexes (3 per triangle) that form the counter clockwise facets.

Example:

C#
```List<Vertex>
set = new List<Vertex>();
Random.Random r = new Random.Random(seed);
for(int i = 0; i < count; i++)
{
}
Mesh m = new Mesh();
m.Recursion = 6;
m.Compute(set, new RectangleF(0, 0, 640, 480));

m.Draw(g, 0, 0, 640, 680);
List<Vertex> p = m.Points;               // Points for directX vertex buffer.
int
[] indicies = m.GetVertexIndicies(); // Indexes for directX surface. ```

## Performance

On my laptop (a decent I7, w/Windows 7) the algorithm runs fast enough over small meshes to be done at video speeds. Each point is a search for the triangle that needs to be split into three, so the time grows with that search. My goal was to use this in another project where I needed to be able to reuse the code, maintain it, and have it run reasonably fast. While this is relatively fast, there are others that are faster, if you need something super fast OpenCV provides a version of the algorithm.

```50 Vertexes 5 msec
100 Vertexes 11 msec
500 Vertexes 129 msec
1k Vertexes 380 msec
5k Vertexes 7.8 seconds
10k Vertexes 31 seconds
60 Hz 16 msec 150 vertexes
30 Hz 33 msec 230 vertexes

20 Hz 50 msec 280 vertexes```

## Share

Phil is a Principal Software developer focusing on weird yet practical algorithms that run the gamut of embedded and desktop (PID loops, Kalman filters, FFTs, client-server SOAP bindings, ASIC design, communication protocols, game engines, robotics).

In his personal life he is a part time mad scientist, full time dad, and studies small circle jujitsu, plays guitar and piano.

 First Prev Next
 SDK version , How to Run this program? Member 1537596429-Sep-21 8:48 Member 15375964 29-Sep-21 8:48
 My vote of 5 tugrulGtx4-May-21 9:47 tugrulGtx 4-May-21 9:47
 3D Triangulation. Can this be used? MikeChatz26-Feb-19 22:28 MikeChatz 26-Feb-19 22:28
 Can use it for 3d triangulation? Владимир Добржанский20-Jul-18 5:56 Владимир Добржанский 20-Jul-18 5:56
 not really fast ilovesc26-Jun-15 21:12 ilovesc 26-Jun-15 21:12
 if calulate points count more than 10k, it will be very slow.
 Re: not really fast Member 1260270224-Jan-17 21:07 Member 12602702 24-Jan-17 21:07
 little bug Member 1137558724-Jan-15 5:49 Member 11375587 24-Jan-15 5:49
 Re: little bug HoshiKata18-Oct-15 8:28 HoshiKata 18-Oct-15 8:28
 delaunay triangulation Member 107112591-Apr-14 11:42 Member 10711259 1-Apr-14 11:42
 delaunay triangulation Member 1071125930-Mar-14 16:54 Member 10711259 30-Mar-14 16:54
 Re: delaunay triangulation HoshiKata31-Mar-14 5:29 HoshiKata 31-Mar-14 5:29
 Re: delaunay triangulation Member 1071125931-Mar-14 10:19 Member 10711259 31-Mar-14 10:19
 Thank you! taonv8914-Nov-13 19:26 taonv89 14-Nov-13 19:26