Download Delaunay.zip
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.
- Compute the bounding rectangle:
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.
Find the triangle containing the new point (for example, abc), and replace the old triangle 1 with three new triangles.
Repair the pointers:
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.
-
Examine each triangle and it's outward adjacent:. For example bec shares a side with bdc.
Examine the angles across from the adjacent sides. If greater than 180 degrees, flip the adjacent line.
If the triangle is not flipped skip to step 10.
Reform the old triangles into the new triangles, fix the triangle edges.
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.
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:
List<Vertex>
set = new List<Vertex>();
Random.Random r = new Random.Random(seed);
for(int i = 0; i < count; i++)
{
set.Add(new Vertex(r.Float(0,width), r.Float(0, height), 0));
}
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;
int
[] indicies = m.GetVertexIndicies();
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
