Working with Delaunay Triangulations

Definition of Delaunay Triangulation

Delaunay triangulations are widely used in scientific computing in many diverse applications. While there are numerous algorithms for computing triangulations, it is the favorable geometric properties of the Delaunay triangulation that make it so useful.

The fundamental property is the Delaunay criterion. In the case of 2-D triangulations, this is often called the empty circumcircle criterion. For a set of points in 2-D, a Delaunay triangulation of these points ensures the circumcircle associated with each triangle contains no other point in its interior. This property is important. In the illustration below, the circumcircle associated with T1 is empty. It does not contain a point in its interior. The circumcircle associated with T2 is empty. It does not contain a point in its interior. This triangulation is a Delaunay triangulation.

The triangles below are different. The circumcircle associated with T1 is not empty. It contains V3 in its interior. The circumcircle associated with T2 is not empty. It contains V1 in its interior. This triangulation is not a Delaunay triangulation.

Delaunay triangles are said to be “well shaped” because in fulfilling the empty circumcircle property, triangles with large internal angles are selected over ones with small internal angles. The triangles in the non-Delaunay triangulation have sharp angles at vertices V2 and V4. If the edge {V2, V4} were replaced by an edge joining V1 and V3, the minimum angle would be maximized and the triangulation would become a Delaunay triangulation. Also, the Delaunay triangulation connects points in a nearest-neighbor manner. These two characteristics, well-shaped triangles and the nearest-neighbor relation, have important implications in practice and motivate the use of Delaunay triangulations in scattered data interpolation.

While the Delaunay property is well defined, the topology of the triangulation is not unique in the presence of degenerate point sets. In two dimensions, degeneracies arise when four or more unique points lie on the same circle. The vertices of a square, for example, have a nonunique Delaunay triangulation.

The properties of Delaunay triangulations extend to higher dimensions. The triangulation of a 3-D set of points is composed of tetrahedra. The next illustration shows a simple 3-D Delaunay triangulation made up of two tetrahedra. The circumsphere of one tetrahedron is shown to highlight the empty circumsphere criterion.

A 3-D Delaunay triangulation produces tetrahedra that satisfy the empty circumsphere criterion.

Creating Delaunay Triangulations

MATLAB^® provides two ways to create Delaunay triangulations:

• The functions delaunay and delaunayn

• The delaunayTriangulation class

The delaunay function supports the creation of 2-D and 3-D Delaunay triangulations. The delaunayn function supports creating Delaunay triangulations in 4-D and higher.

Tip

Creating Delaunay triangulations in dimensions higher than 6-D is generally not practical for moderate to large point sets due to the exponential growth in required memory.

The delaunayTriangulation class supports creating Delaunay triangulations in 2-D and 3-D. It provides many methods that are useful for developing triangulation-based algorithms. These class methods are like functions, but they are restricted to work with triangulations created using delaunayTriangulation. The delaunayTriangulation class also supports the creation of related constructs such as the convex hull and Voronoi diagram. It also supports the creation of constrained Delaunay triangulations.

In summary:

• The delaunay function is useful when you only require the basic triangulation data, and that data is sufficiently complete for your application.

• The delaunayTriangulation class offers more functionality for developing triangulation-based applications. It is useful when you require the triangulation and you want to perform any of these operations:

  • Search the triangulation for triangles or tetrahedra enclosing a query point.

  • Use the triangulation to perform a nearest-neighbor point search.

  • Query the triangulation's topological adjacency or geometric properties.

  • Modify the triangulation to insert or remove points.

  • Constrain edges in the triangulation—this is called a constrained Delaunay triangulation.

  • Triangulate a polygon and optionally remove the triangles that are outside of the domain.

  • Use the Delaunay triangulation to compute the convex hull or Voronoi diagram.

Using the delaunay and delaunayn functions

The delaunay and delaunayn functions take a set of points and produce a triangulation in matrix format. Refer to Triangulation Matrix Format for more information on this data structure. In 2-D, the delaunay function is often used to produce a triangulation that can be used to plot a surface defined in terms of a set of scattered data points. In this application, it’s important to note that this approach can only be used if the surface is single-valued. For example, it could not be used to plot a spherical surface because there are two z values corresponding to a single (x, y) coordinate. A simple example demonstrates how the delaunay function can be used to plot a surface representing a sampled data set.

This example shows how to use the delaunay function to create a 2-D Delaunay triangulation from the seamount data set. A seamount is an underwater mountain. The data set consists of a set of longitude (x) and latitude (y) locations, and corresponding seamount elevations (z) measured at those coordinates.

Load the seamount data set and view the (x, y) data as a scatter plot.

load seamount
plot(x,y,'.','markersize',12)
xlabel('Longitude'), ylabel('Latitude')
grid on

Construct a Delaunay triangulation from this point set and use triplot to plot the triangulation in the existing figure.

tri = delaunay(x,y);
hold on, triplot(tri,x,y), hold off

Add the depth data (z) from seamount to lift the vertices and create the surface. Create a new figure and use trimesh to plot the surface in wireframe mode.

figure
hidden on
trimesh(tri,x,y,z)
xlabel('Longitude'),ylabel('Latitude'),zlabel('Depth in Feet');

If you want to plot the surface in shaded mode, use trisurf instead of trimesh.

A 3-D Delaunay triangulation also can be created using the delaunay function. This triangulation is composed of tetrahedra.

This example shows how to create a 3-D Delaunay triangulation of a random data set. The triangulation is plotted using tetramesh, and the FaceAlpha option adds transparency to the plot.

rng('default')
X = rand([30 3]);
tet = delaunay(X);
faceColor  = [0.6875 0.8750 0.8984];
tetramesh(tet,X,'FaceColor', faceColor,'FaceAlpha',0.3);

MATLAB provides the delaunayn function to support the creation of Delaunay triangulations in dimension 4-D and higher. Two complementary functions tsearchn and dsearchn are also provided to support spatial searching for N-D triangulations. See Spatial Searching for more information on triangulation-based search.

Using the delaunayTriangulation Class

The delaunayTriangulation class provides another way to create Delaunay triangulations in MATLAB. While delaunay and delaunayTriangulation use the same underlying algorithm and produce the same triangulation, delaunayTriangulation provides complementary methods that are useful for developing Delaunay-based algorithms. These methods are like functions that are packaged together with the triangulation data into a container called a class. Keeping everything together in a class provides a more organized setup that improves ease of use. It also improves the performance of triangulation-based searches such as point-location and nearest-neighbor. delaunayTriangulation supports incremental editing of the Delaunay triangulation. You also can impose edge constraints in 2-D.

Triangulation Representations introduces the triangulation class, which supports topological and geometric queries for 2-D and 3-D triangulations. A delaunayTriangulation is a special kind of triangulation. This means you can perform any triangulation query on a delaunayTriangulation in addition to the Delaunay-specific queries. In more formal MATLAB language terms, delaunayTriangulation is a subclass of triangulation.

This example shows how to create, query, and edit a Delaunay triangulation from the seamount data using delaunayTriangulation. The seamount data set contains (x, y) locations and corresponding elevations (z) that define the surface of the seamount.

Load and triangulate the (x, y) data.

load seamount
DT = delaunayTriangulation(x,y)

DT = 
  delaunayTriangulation with properties:

              Points: [294x2 double]
    ConnectivityList: [566x3 double]
         Constraints: []

The Constraints property is empty because there aren't any imposed edge constraints. The Points property represents the coordinates of the vertices, and the ConnectivityList property represents the triangles. Together, these two properties define the matrix data for the triangulation.

The delaunayTriangulation class is a wrapper around the matrix data, and it offers a set of complementary methods. You access the properties in a delaunayTriangulation in the same way you access the fields of a struct.

Access the vertex data.

DT.Points;

Access the connectivity data.

DT.ConnectivityList;

Access the first triangle in the ConnectivityList property.

DT.ConnectivityList(1,:)

ans = 1×3

   230   205   152

delaunayTriangulation provides an easy way to index into the ConnectivityList property matrix.

Access the first triangle.

DT(1,:)

ans = 1×3

   230   205   152

Examine the first vertex of the first triangle.

DT(1,1)

ans = 
230

Examine all the triangles in the triangulation.

DT(:,:);

Indexing into the delaunayTriangulation output, DT, works like indexing into the triangulation array output from delaunay. The difference between the two are the extra methods that you can call on DT (for example, nearestNeighbor and pointLocation).

Use triplot to plot the delaunayTriangulation. The triplot function is not a delaunayTriangulation method, but it accepts and can plot a delaunayTriangulation.

triplot(DT); 
axis equal
xlabel('Longitude'), ylabel('Latitude')
grid on

Alternatively, you could use triplot(DT(:,:), DT.Points(:,1), DT.Points(:,2)); to get the same plot.

Use the delaunayTriangulation method, convexHull, to compute the convex hull and add it to the plot. Since you already have a Delaunay triangulation, this method allows you to derive the convex hull more efficiently than a full computation using convhull.

hold on
k = convexHull(DT);
xHull = DT.Points(k,1);
yHull = DT.Points(k,2);
plot(xHull,yHull,'r','LineWidth',2); 
hold off

You can incrementally edit the delaunayTriangulation to add or remove points. If you need to add points to an existing triangulation, then an incremental addition is faster than a complete retriangulation of the augmented point set. Incremental removal of points is more efficient when the number of points to be removed is small relative to the existing number of points.

Edit the triangulation to remove the points on the convex hull from the previous computation.

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on

% The convex hull topology duplicates the start and end vertex.
% Remove the duplicate entry.
k(end) = [];

% Now remove the points on the convex hull.
DT.Points(k,:) = []

DT = 
  delaunayTriangulation with properties:

              Points: [274x2 double]
    ConnectivityList: [528x3 double]
         Constraints: []

% Plot the new triangulation.
hold on
triplot(DT); 
hold off

There is one vertex that is just inside the boundary of the convex hull that was not removed. The fact that it is interior to the hull can be seen using the Zoom-In tool in the figure. You could plot the vertex labels to determine the index of this vertex and remove it from the triangulation. Alternatively, you can use the nearestNeighbor method to identify the index more readily.

The point is close to location (211.6, -48.15). Use the nearestNeighbor method to find the nearest vertex.

vertexId = nearestNeighbor(DT, 211.6, -48.15)

vertexId = 
50

Now remove that vertex from the triangulation.

DT.Points(vertexId,:) = []

DT = 
  delaunayTriangulation with properties:

              Points: [273x2 double]
    ConnectivityList: [525x3 double]
         Constraints: []

Plot the new triangulation.

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on
hold on
triplot(DT); 
hold off

Add points to the existing triangulation. Add 4 points to form a rectangle around the triangulation.

Padditional = [210.9 -48.5; 211.6 -48.5; ...
     211.6 -47.9; 210.9 -47.9];
DT.Points(end+(1:4),:) = Padditional

DT = 
  delaunayTriangulation with properties:

              Points: [277x2 double]
    ConnectivityList: [548x3 double]
         Constraints: []

Close all existing figures.

close all

Plot the new triangulation.

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on
hold on
triplot(DT); 
hold off

You can edit the points in the triangulation to move them to a new location. Edit the first of the additional point set (the vertex ID 274).

DT.Points(274,:) = [211 -48.4];

Close all existing figures.

close all

Plot the new triangulation

figure
plot(xHull,yHull,'r','LineWidth',2); 
axis equal
xlabel('Longitude'),ylabel('Latitude')
grid on
hold on
triplot(DT); 
hold off

Use the a method of the triangulation class, vertexAttachments, to find the attached triangles. Since the number of triangles attached to a vertex is variable, the method returns the attached triangle IDs in a cell array. You need braces to extract the contents.

attTris = vertexAttachments(DT,274);
hold on
triplot(DT(attTris{:},:),DT.Points(:,1),DT.Points(:,2),'g')
hold off

delaunayTriangulation also can be used to triangulate points in 3-D space. The resulting triangulation is composed of tetrahedra.

This example shows how to use a delaunayTriangulation to create and plot the triangulation of 3-D points.

rng('default')
P = rand(30,3);
DT = delaunayTriangulation(P)

DT = 
  delaunayTriangulation with properties:

              Points: [30x3 double]
    ConnectivityList: [102x4 double]
         Constraints: []

faceColor  = [0.6875 0.8750 0.8984];
tetramesh(DT,'FaceColor', faceColor,'FaceAlpha',0.3);

The tetramesh function plots both the internal and external faces of the triangulation. For large 3-D triangulations, plotting the internal faces might be an unnecessary use of resources. A plot of the boundary might be more appropriate. You can use the freeBoundary method to get the boundary triangulation in matrix format. Then pass the result to trimesh or trisurf.

Constrained Delaunay Triangulation

The delaunayTriangulation class allows you to constrain edges in a 2-D triangulation. This means you can choose a pair of points in the triangulation and constrain an edge to join those points. You can picture this as “forcing” an edge between one or more pairs of points. The following example shows how edge constraints can affect the triangulation.

The triangulation below is a Delaunay triangulation because it respects the empty circumcircle criterion.

Triangulate a set of points with an edge constraint specified between vertex V1 and V3.

Define the point set.

P = [2 4; 6 1; 9 4; 6 7];

Define a constraint, C, between V1 and V3.

C = [1 3];
DT = delaunayTriangulation(P,C);

Plot the triangulation and add annotations.

triplot(DT)

% Label the vertices.
hold on
numvx = size(P,1);
vxlabels = arrayfun(@(n) {sprintf('V%d', n)}, (1:numvx)');
Hpl = text(P(:,1)+0.2, P(:,2)+0.2, vxlabels, 'FontWeight', ...
  'bold', 'HorizontalAlignment','center', 'BackgroundColor', ...
  'none');
hold off


% Use the incenters to find the positions for placing triangle labels on the plot.
hold on
IC = incenter(DT);
numtri = size(DT,1);
trilabels = arrayfun(@(P) {sprintf('T%d', P)}, (1:numtri)');
Htl = text(IC(:,1),IC(:,2),trilabels,'FontWeight','bold', ...
'HorizontalAlignment','center','Color','blue');
hold off

% Plot the circumcircle associated with the triangle, T1.
hold on
[CC,r] = circumcenter(DT);
theta = 0:pi/50:2*pi;
xunit = r(1)*cos(theta) + CC(1,1);
yunit = r(1)*sin(theta) + CC(1,2);
plot(xunit,yunit,'g');
axis equal
hold off

The constraint between vertices (V1, V3) was honored, however, the Delaunay criterion was invalidated. This also invalidates the nearest-neighbor relation that is inherent in a Delaunay triangulation. This means the nearestNeighbor search method provided by delaunayTriangulation cannot be supported if the triangulation has constraints.

In typical applications, the triangulation might be composed of many points, and a relatively small number of edges in the triangulation might be constrained. Such a triangulation is said to be locally non-Delaunay, because many triangles in the triangulation might respect the Delaunay criterion, but locally there might be some triangles that do not. In many applications, local relaxation of the empty circumcircle property is not a concern.

Constrained triangulations are generally used to triangulate a nonconvex polygon. The constraints give us a correspondence between the polygon edges and the triangulation edges. This relationship enables you to extract a triangulation that represents the region. The following example shows how to use a constrained delaunayTriangulation to triangulate a nonconvex polygon.

Define and plot a polygon.

figure()
axis([-1 17 -1 6]);
axis equal
P = [0 0; 16 0; 16 2; 2 2; 2 3; 8 3; 8 5; 0 5];
patch(P(:,1),P(:,2),'-r','LineWidth',2,'FaceColor',...
     'none','EdgeColor','r');
 
% Label the points.
hold on
numvx = size(P,1);
vxlabels = arrayfun(@(n) {sprintf('P%d', n)}, (1:numvx)');
Hpl = text(P(:,1)+0.2, P(:,2)+0.2, vxlabels, 'FontWeight', ...
  'bold', 'HorizontalAlignment','center', 'BackgroundColor', ...
  'none');
hold off

Create and plot the triangulation together with the polygon boundary.

figure()
subplot(2,1,1);
axis([-1 17 -1 6]);
axis equal

P = [0 0; 16 0; 16 2; 2 2; 2 3; 8 3; 8 5; 0 5];
DT = delaunayTriangulation(P);
triplot(DT)
hold on; 
patch(P(:,1),P(:,2),'-r','LineWidth',2,'FaceColor',...
     'none','EdgeColor','r');
hold off

% Plot the standalone triangulation in a subplot.
subplot(2,1,2);
axis([-1 17 -1 6]);
axis equal
triplot(DT)

This triangulation cannot be used to represent the domain of the polygon because some triangles cut across the boundary. You need to impose a constraint on the edges that are cut by triangulation edges. Since all edges have to be respected, you need to constrain all edges. The steps below show how to constrain all the edges.

Enter the constrained edge definition. Observe from the annotated figure where you need constraints (between (V1, V2), (V2, V3), and so on).

C = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];

In general, if you have N points in a sequence that define a polygonal boundary, the constraints can be expressed as C = [(1:(N-1))' (2:N)'; N 1];.

Specify the constraints when you create the delaunayTriangulation.

DT = delaunayTriangulation(P,C);

Alternatively, you can impose constraints on an existing triangulation by setting the Constraints property: DT.Constraints = C;.

Plot the triangulation and polygon.

figure('Color','white')
subplot(2,1,1);
axis([-1 17 -1 6]);
axis equal
triplot(DT)
hold on; 
patch(P(:,1),P(:,2),'-r','LineWidth',2, ...
     'FaceColor','none','EdgeColor','r'); 
hold off

% Plot the standalone triangulation in a subplot.
subplot(2,1,2);
axis([-1 17 -1 6]);
axis equal
triplot(DT)

The plot shows that the edges of the triangulation respect the boundary of the polygon. However, the triangulation fills the concavities. What is needed is a triangulation that represents the polygonal domain. You can extract the triangles within the polygon using the delaunayTriangulation method, isInterior. This method returns a logical array whose true and false values that indicate whether the triangles are inside a bounded geometric domain. The analysis is based on the Jordan Curve theorem, and the boundaries are defined by the edge constraints. The ith triangle in the triangulation is considered to be inside the domain if the ith logical flag is true, otherwise it is outside.

Now use the isInterior method to compute and plot the set of domain triangles.

% Plot the constrained edges in red.
figure('Color','white')
subplot(2,1,1);
plot(P(C'),P(C'+size(P,1)),'-r','LineWidth', 2);
axis([-1 17 -1 6]);

% Compute the in/out status.
IO = isInterior(DT);
subplot(2,1,2);
hold on;
axis([-1 17 -1 6]);

% Use triplot to plot the triangles that are inside.
% Uses logical indexing and dt(i,j) shorthand
% format to access the triangulation.
triplot(DT(IO, :),DT.Points(:,1), DT.Points(:,2),'LineWidth', 2)
hold off;

Triangulation of Point Sets Containing Duplicate Locations

The Delaunay algorithms in MATLAB construct a triangulation from a unique set of points. If the points passed to the triangulation function, or class, are not unique, the duplicate locations are detected and the duplicate point is ignored. This produces a triangulation that does not reference some points in the original input, namely the duplicate points. When you work with the delaunay and delaunayn functions, the presence of duplicates may be of little consequence. However, since many of the queries provided by the delaunayTriangulation class are index based, it is important to understand that delaunayTriangulation triangulates and works with the unique data set. Therefore, indexing based on the unique point set is the convention. This data is maintained by the Points property of delaunayTriangulation.

The following example illustrates the importance of referencing the unique data set stored within the Points property when working with delaunayTriangulation:

rng('default')
P = rand([25 2]);
P(18,:) = P(8,:)
P(16,:) = P(6,:)
P(12,:) = P(2,:)

DT = delaunayTriangulation(P)

When the triangulation is created, MATLAB issues a warning. The Points property shows that the duplicate points have been removed from the data.

DT = 

  delaunayTriangulation with properties:

              Points: [22x2 double]
    ConnectivityList: [31x3 double]
         Constraints: []

If for example, the Delaunay triangulation is used to compute the convex hull, the indices of the points on the hull are indices with respect to the unique point set, DT.Points. Therefore, use the following code to compute and plot the convex hull:

K = DT.convexHull();
plot(DT.Points(:,1),DT.Points(:,2),'.');
hold on
plot(DT.Points(K,1),DT.Points(K,2),'-r');

If the original data set containing the duplicates were used in conjunction with the indices provided by delaunayTriangulation, then the result would be incorrect. The delaunayTriangulation works with indices that are based on the unique data set DT.Points. For example, the following would produce an incorrect plot, because K is indexed with respect to DT.Points and not P:

K = DT.convexHull();
plot(P(:,1),P(:,2),'.');
hold on
plot(P(K,1),P(K,2),'-r');

It’s often more convenient to create a unique data set by removing duplicates prior to creating the delaunayTriangulation. Doing this eliminates the potential for confusion. This can be accomplished using the unique function as follows:

rng('default')
P = rand([25 2]);
P(18,:) = P(8,:)
P(16,:) = P(6,:)
P(12,:) = P(2,:)

[~, I, ~] = unique(P,'first','rows');
I = sort(I);
P = P(I,:);
DT = delaunayTriangulation(P)  % The point set is unique

Related Topics

• Spatial Searching