Triangulation

  

 

Input description: A set of points or a polyhedon.

Problem description: Partition the interior of the point set or polyhedron into triangles.

Discussion: Triangulation is a fundamental problem in computational geometry, because the first step in working with complicated geometric objects is to break them into simple geometric objects. The simplest geometric objects are triangles in two dimensions, and tetrahedra in three. Classical applications of triangulation include finite element analysis and computer graphics.    

   A particularly interesting application of triangulation is surface or function interpolation. Suppose that we have sampled the height of a mountain at a certain number of points. How can we estimate the height at any point q in the plane? If we project the points on the plane, and then triangulate them, the triangulation completely partitions the plane into regions. We can estimate the height of q by interpolating among the three points of the triangle that contains it. Further, this triangulation and the associated height values define a surface of the mountain suitable for graphics rendering.

In the plane, a triangulation is constructed by adding nonintersecting chords between the vertices until no more such chords can be added. Specific issues arising in triangulation include:

• Does the shape of the triangles in your triangulation matter? - There are usually many different ways to partition your input into triangles. Consider a set of n points in convex position in the plane. The simplest way to triangulate them would be to add to the convex hull diagonals from the first point to all of the others. However, this has the tendency to create skinny triangles.  

  If the shape of the triangles matters for your application, you will usually want to avoid skinny triangles, or equivalently, small angles in the triangulation. The Delaunay triangulation of a point set minimizes the maximum angle over all possible triangulations. This isn't exactly what we are looking for, but it is pretty close, and the Delaunay triangulation has enough other interesting properties (including that it is dual to the Voronoi diagram) to make it the quality triangulation of choice. Further, it can be constructed in time, with implementations described below.   

• What dimension are we working in? - As always, three-dimensional problems are harder than two-dimensional problems. The three-dimensional generalization of triangulation involves partitioning the space into four-vertex tetrahedra by adding nonintersecting faces. One important difficulty is that for certain polyhedra there is no way to tetrahedralize the interior without adding extra vertices. Further, it is NP-complete to decide whether such a tetrahedralization exists, so we should not feel afraid to add extra vertices to simplify our problem.  
• What constraints does the input have? - When we are triangulating a polygon or polyhedra, we are restricted to adding chords that do not intersect any of the boundary facets. In general, we may have a set of obstacles or constraints that cannot be intersected by inserted chords. The best such triangulation is likely to be the so-called constrained Delaunay triangulation. Implementations are described below.  
• Are you allowed to add extra points? - When the shape of the triangles does matter, it might pay to strategically add a small number of extra ``Steiner'' points to the data set to facilitate the construction of a triangulation (say) with no small angles. As discussed above, there may be no triangulation possible for certain polyhedra without adding Steiner points.  

To construct a triangulation of a convex polygon in linear time, just pick an arbitrary starting vertex v and insert chords from v to each other vertex in the polygon.   Because the polygon is convex, we can be confident that none of the boundary edges of the polygon will be intersected by these chords and that all of them lie within the polygon. The simplest algorithm for constructing general polygon triangulations tries each of the possible chords and inserts them if they do not intersect a boundary edge or previously inserted chord. There are practical algorithms that run in time and theoretically interesting algorithms that run in linear time. See the implementations and notes below for details.

Implementations: Triangle, by Jonathan Shewchuk of Carnegie-Mellon University, is a C language code that generates Delaunay triangulations, constrained Delaunay triangulations (forced to have certain edges), and quality-conforming Delaunay triangulations (which avoid small angles by inserting extra points).     It has been widely used for finite element analysis and other applications and is fast and robust. Triangle is the first thing I would try if I needed a two-dimensional triangulation code. Although Triangle is available at http://www.cs.cmu.edu/ quake/triangle.html, it is copyrighted by the author and may not be sold or included in commercial products without a license.

GEOMPACK is a suite of Fortran 77 codes by Barry Joe of the University of Alberta, for 2- and 3-dimensional triangulation and convex decomposition problems.    In particular, it does both Delaunay triangulation and convex decompositions of polygonal and polyhedral regions, as well as arbitrary-dimensional Delaunay triangulations. They can be obtained from ftp://ftp.cs.ualberta.ca/pub/geompack.

Steve Fortune is the author of a widely used 2D code for Voronoi diagrams and Delaunay triangulations, written in C. This code is smaller and probably simpler to work with than either of the above, if all you need is the Delaunay triangulation of points in the plane. It is based on Fortune's own sweepline algorithm [For87] for Voronoi diagrams and is available from Netlib (see Section ) at http://netlib.bell-labs.com/netlib/voronoi/index.html.  

O'Rourke [O'R94] provides asymptotically slow implementations in C of polygon triangulation (in ) and Delaunay triangulation (in ). These will be unusable for more than modest numbers of points, but see Section if interested. See Section .

Algorithm 624 [Ren84] of the Collected Algorithms of the ACM is a Fortran implementation of triangulation for surface interpolation. See Section . A linear-time implementation for triangulating a planar map is included with LEDA (see Section ).   

Higher-dimensional Delaunay triangulations are a special case of higher-dimensional convex hulls, and Qhull [BDH97] appears to be the convex hull code of choice for general dimensions (i.e. three dimensions and beyond). It is written in C, and it can also construct Voronoi vertices, furthest-site Voronoi vertices, and half-space intersections. Qhull has been widely used in scientific applications and has a well-maintained home page at http://www.geom.umn.edu/software/qhull/.   

Notes: After a long search, Chazelle [Cha91] discovered a linear-time algorithm for triangulating a simple polygon. This algorithm is sufficiently hopeless to implement that it qualifies more as an existence proof. The first algorithm for polygon triangulation was given by [GJPT78]. An algorithm by Tarjan and Van Wyk [TW88] followed before Chazelle's result. Expositions on polygon and point set triangulation include [O'R94, PS85].

Linear-time algorithms for triangulating monotone polygons have been long known [GJPT78] and are the basis of algorithms for triangulating simple polygons. A polygon is monotone when there exists a direction d such that any line with slope d intersects the polygon in at most two points.  

A heavily studied class of optimal triangulations seeks to minimize the total length of the chords used. The computational complexity of constructing this minimum weight triangulation is a long-standing open problem in computational geometry, so the interest has shifted to provably good approximation algorithms. The minimum weight triangulation of a convex polygon can be found in time using dynamic programming, as discused in Section .   

Related Problems: Voronoi diagrams (see page ), polygon partitioning (see page ).