CONSTRAINED CONSTRUCTION OF PLANAR DELAUNAY TRIANGULATIONS WITHOUT FLIPPING

The construction of Voronoi diagrams and Delaunay triangulations finds wide application in many branches of science. Delaunay triangulations have properties which make them more desirable than other triangulations for the same node set. Delaunay has characterized his triangulations by the empty circle property. The partitioning and flipping methods which have been developed for digital construction of Voronoi diagrams and Delaunay triangulations only make indirect use of this property. A novel method of construction is proposed, which is based directly on the empty circle property of Delaunay. The geometry of the steps of the algorithm is simple and can be grasped intuitively. The method can be applied to constrained triangulations, in which a triangulation domain and some of the edges are prescribed. A data structure for triangulations of concave and multiply-connected domains is presented which permits convenient specification of the constraints and the triangulation. The method is readily implemented, efficient and robust.


Voronoi Diagrams and Delaunay Triangulations
Over centuries, outstanding scientists have considered the following problem: Given a set of sites in a plane, partition the plane into regions in such a way that all points of a region are at least as near to a particular site as to any other site of the set. This problem has arisen in many applications [1] [2] [3]. Johannes Keppler (1571-1639) encountered Voronoi diagrams when he considered the densest packing of spheres, René Descartes (1596-1650) when he investigated the distribution of matter relative to fixed stars, Johann Dirichlet (1805-1859) when he constructed integer lattices, John Snow (1813-1853) when he related the outbreak of cholera in London to the location of water pumps, Georges Voronoi (1868Voronoi ( -1908 when he extended the Dirichlet tessellations to higher dimensions and Boris Delaunay (1890Delaunay ( -1980 when he generalized Voronoi diagrams and their duals to irregularly spaced sites in ddimensional space. Many additional examples of Voronoi diagrams have been reported, for example in gold mining, crystallography, metallurgy and meteorology. Today diagrams which relate sets of nearest points to sites in a plane are known as Voronoi diagrams, sometimes also as Dirichlet tessellations. Extensions of the nearness concept to weighted tessellations by Dirichlet [4] have been studied by Aurenhammer [5] and are today known as power diagrams. Extensions to power diagrams are not treated in this paper.
The concept of triangulation has been related to the concept of sets of nearest points by Voronoi [2]. The relationship was deepened when Delaunay [3] distinguished the triangulation, which is derived from a Voronoi diagram, from other triangulations for the same set of sites by defining the empty circle property as follows: A triangulation is Delaunay if, and only if, the interior of the circumscribed circles none of its triangles contains sites. Triangulations with the empty circle property are today called Delaunay triangulations. They have applications in many branches of science, for example in computational geometry and in the networks for discrete methods of physical analysis such as the finite element method.  Figure 1 shows an example of a set of sites in a plane and the corresponding Voronoi diagram. The straight line segment that is a boundary between two Voronoi regions is called a Voronoi arc. An intersection of Voronoi arcs is called a Voronoi vertex. A Voronoi circle can be drawn with a Voronoi vertex as midpoint such that it passes through the sites of the Voronoi regions that meet at the Voronoi vertex. Voronoi diagrams have distinctive properties [6]. Each Voronoi region is convex. The diagram contains bounded and unbounded regions. At least three arcs and at least three regions meet at a Voronoi vertex, and the distances of a Voronoi vertex from the sites of the regions which meet at the vertex are equal. Figure 2 shows the Delaunay triangulation for the set of sites in figure 1. If a Voronoi circle passes through exactly three sites, these sites are the corners of a Delaunay triangle and are called Delaunay nodes. If the circle passes through more than 3 sites, the sites define a convex polygon that can be decomposed into Delaunay triangles. The decomposition is not unique. The straight line segment connecting two nodes of a Delaunay triangle is called a Delaunay edge. The Delaunay triangulation is the convex hull of the set of sites in the Voronoi diagram. While the Delaunay triangulation is not the only possible triangulation of the sites of the Voronoi diagram, it is the triangulation for which the magnitude of the smallest angle of all triangles is at a maximum. This property is advantageous for discrete analysis, for example by the finite element method, where the use of nearly equilateral triangles increases the accuracy of the solution. GEOMETRICAL   A triangulation for a set of sites is called unconstrained, if neither the boundary of the area covered by the triangles is prescribed, nor any of its interior edges. The boundary of an unconstrained triangulation is the convex hull of the given set of sites. The two main classes of algorithms, which have been developed for unconstrained triangulation, make use of partitioning of the set of sites and of edge flipping respectively.
If a subdomain of the area covered by the convex hull is prescribed for the triangulation, or if any of its internal edges is prescribed, the triangulation is called constrained. The goal of the research, which is reported in this paper, is to determine whether it is promising to extend one of the two conventional methods to include constraints, or whether the problem can be solved more advantageously with a novel approach. For this purpose, the relationships between the existing methods of unconstrained triangulation are analyzed. It will be shown that they are all based on the same basic property, which was discovered by Delaunay. A novel method of constrained analysis, which is based directly on this property, will be presented.

Duality of Voronoi Diagrams and Delaunay Triangulations
Voronoi diagrams and Delaunay triangulations are related by the concept of duality. Duality is a concept which is encountered in mathematics in many different forms. In the context of Voronoi diagrams and Delaunay triangles, duality is a bijective mapping which reverses the inclusion relationships. This mapping will now be defined.
A Delaunay triangulation is a planar subdivision whose element set D consists of nodes, edges and triangles. It is arbitrarily named the primal subdivision of the plane, and its elements are called primal elements. The Voronoi diagram is a planar subdivision whose element set V consists of vertices, arcs and regions. It is called the dual subdivision of the plane and its elements are called dual elements. Consider a mapping f of the set of elements D of the Delaunay triangulation to the set of elements V of the Voronoi diagram:

Graphical Construction of Delaunay Triangulations
The historic development from Voronoi diagrams to Delaunay triangulations suggests an intuitive approach to the construction of Delaunay triangulations. Using paper, pencil and the capability of human vision to analyze shape, the graphic construction proceeds as shown in  The second vertex of the Delaunay arcs radiating from 1 v is determined next. For this purpose the bisectors of the line segments connecting the given sites are constructed as shown in figure 4. For example, the bisector of line segment 1 4 n n is 1 2 a a . Let the bisector leaving 1 v in the direction of point 2 a be followed. This bisector is the boundary that separates points that are nearer to site 1 n from points that are nearer to site 4 n . After this arc has reached the point of intersection with bisector 1 2 f f , it enters a set of points which are closer to site 2 n than to sites 1 4 n and n . The point of intersection is therefore the second vertex 2 v of the arc. If the arc leaving 1 v in the direction of point 2 e is followed, none of the bisectors in figure 4 is intersected. This arc is infinite. Similarly, the arc leaving 1 v in the direction of point 2 b is infinite. The bisectors which are overlaid by the constructed arcs are removed, because the Voronoi regions are convex, such that each bisector can overlay only one arc of the Voronoi diagram. The vertices, which are constructed on the bisectors leaving the first vertex, are collected in a list. After all bisectors at node 1 v have been treated, one of the vertices in the list is chosen and removed from the list. The remaining bisectors leaving this vertex are traversed in the direction of the line segment which they bisect. If they intersect another bisector, a Voronoi vertex is constructed at the point of intersection. It is the second vertex of the arc and is added to the list. If none of the remaining bisectors is intersected, the arc is infinite. The bisector is removed from the diagram in figure 4. The procedure ends if the vertex list is empty after the bisectors at a vertex have been treated.
The line segment connecting two nodes is an edge of the triangulation if it is bisected by one of the Voronoi arcs that are constructed. For example, segment 1 4 n n is an edge because it is bisected by arc 1 2 v v . Segment 2 3 n n is not an edge of the Delaunay diagram because it is not bisected by any of the arcs that have been constructed.
The graphical construction shows how the duality mapping in figure 3 can be derived. Arc k a overlays the bisector of edge k e . Nodes s t n and n of edge k e are the duals of regions s t r and r whose common boundary is arc k a . Vertices g h v and v of arc k a are the duals of cells g c and h c whose common boundary is edge k e . The end of an infinite arc is the dual of the infinite area consisting of the difference of the area of the infinite plane and the areas of the Delaunay triangles.

Construction of Delaunay Triangulations by Partitioning
The graphical method of construction in section 3 suffers from the disadvantage that the inter-section of the bisector leaving a Voronoi vertex must be tested with all other bisectors of the diagram. While human vision readily solves this problem, digital testing of the bisections is expensive. Algorithms which reduce the problem by partitioning the set of sites into subsets, whose Voronoi diagram can be determined by inspection, and merge these diagrams pairwise until only the desired Voronoi diagram remains, are called "divide and conquer" algorithms. An algorithm of this type by Guibas [7] has been used by many researchers as basis for diverse modified algorithms of similar type, which differ primarily in the merge operation. Many contri-butions to the topic have been made, for example by Shamos [8] and Skvortsov [9] [10]. A modified version of the Guibas algorithm is presented in figure 6.
The given sites are decomposed into sets of two or three sites, such that coordinates 1 x of the sites in the first group are less than the coordinates 1 x off all other sites, the coordinates 1 x of the sites in the second group are less than the coordinates 1 x of the remaining sites, and so on. Figure 6 shows a given set of 6 sites i n in diagram D1 and their decomposition into two groups of 3 sites each. The Delaunay triangle for each group of 3 sites is constructed by inspection. The triangles can be highly distorted if the coordinates 2 x of the three nodes differ significantly. It is assumed that the sites are not collinear. Special methods are required if the sites are collinear. In order to merge the two groups, a base line segment connecting a node of the left group with a node of the right group is determined such that all nodes lie on or above the base line. In diagram D2, the base line connects nodes 1 5 n and n . The gap between the two triangles in diagram D2 is triangulated such that s single triangulation results, which is convex.

Diagram
D3 shows that this is achieved by constructing bisectors 1|3 n , n and n is a Delaunay triangle because the circle is empty. Diagrams D4 to D6 show the analogous construction of the other Delaunay triangles for the merge Operation. If there are more than 6 nodes, additional groups are formed in diagram D2 and there will be more than a single level of merge operations. For a set of 10 3 59049  sites, 10 1 9   levels of merging are required. The complexity of the method is O(n log(n)).  The concept of edge flipping was first presented by Lawson [11] and developed further by Joe [12], de Loera [13] and many others. The basic concept is to construct an arbitrary triangulation for the given set of nodes and then to flip edges until the circumscribed circles of all triangles are empty as required by Delaunay. The convergence of this procedure is proved by studying the effect of edge flipping on a triangulation which is not yet Delaunay.
Triangulation T is considered to be better than triangulation T if vector t is larger than vector . t The optimal triangulation maximizes the smallest angle that occurs in the triangles.

Regular edge:
In order to prove the existence of optimal triangulations, the effect of flipping edges in an existing triangulation is studied. Figure 8 shows the effect of flipping an edge in a convex quadrilateral of a triangulation. Let t be the angle vector for triangles abc and adb in the triangulation with edge ab. Let t be the angle vector for the triangulation with edge cd. Edge ab is called regular if vector t is lexicographically larger than or equal to vector . Proof: Figure 9 shows two triangulations for a quadrilateral with nodes k n . The triangles formed by diagonal 1 3 n n are shown in diagram (a), the triangles formed by diagonal 2 4 n n in diagram (b). Diagram (a) contains the circumscribed circle of triangle 1 2 3 n n n . Diagram (b) contains the circumscribed circle for triangle 1 2 4 n n n . The extended diagonals 1 3 2 4 n n and n n intersect the circumscribed circles in point u. It will be shown that one of the edges 1 3 n n and 2 4 n n is regular and the other irregular.
The proof shows that edge 2 4 n n is regular. Generally, the regular edge in a quadrilateral can be recognized by two attributes: Empty circle: Edge 1 3 n n in triangulation (a) is irregular. Node 4 n lies inside the circumscribed circle and the circle is therefore not an empty circle. Edge 2 4 n n in triangulation (b) is regular. Node 3 n lies outside the circumscribed circle, which does not contain any nodes and is therefore an empty circle.
Shortest diagonal: The regular diagonal 2 4 n n is shorter than the irregular diagonal 1 3 n n . This follows from triangulation (a) where triangles 1 2 3 n n n and 1 2 n n u are congruent, such that chords 1 3 n n and 2 n u are of equal length. Node 4 n is an inner point of chord 2 n u such that the length of 2 4 n n is less than the length of 2 n u, and therefore less than the length of 1 3 n n .

Total flip algorithm:
The total flip algorithm is based on the property derived above, that a diagonal of a quadrangle formed by two neighboring triangles is a regular edge if it is shorter than the other diagonal of the quadrangle. The algorithm is initiated by constructing an arbitrary triangulation for the given node set. The quadrilaterals in the triangulation are traversed. If the diagonal of a convex quadrilateral is not a regular edge, it is flipped such that it becomes regular. Flipping of a diagonal in a quadrilateral can destroy the regularity in adjoining quadrilaterals as shown in figure 10. In triangulation (a) the diagonal be of quadrilateral abde is regular, whereas diagonal bd of quadrilateral bcde is irregular. Diagonal bd is flipped to create triangulation (b). In triangulation (b) the diagonal be in quadrilateral abcd is irregular, whereas diagonal ce in quadrilateral bcde is regular. Diagonal be is flipped to create triangulation (c). Diagonals ac and ce are regular in triangulation (c), which is the Delaunay triangulation.
The traversal of the quadrilaterals of the triangulation must be continued until all diagonals are regular (the circumscribed circles of all triangles are empty). Unfortunately, the algorithm converges slowly for large node sets and is therefore not suitable for application.

Incremental flip algorithm:
An envelope triangle is constructed whose interior contains all nodes of the given set. The given nodes are added to the triangle one at a time, as shown in figure 11. After each node addition, edge flipping is applied until the circumscribed circles of all triangles in the triangulation are empty. Joe [12] and Edelsbrunner [14] have shown that the incremental flip algorithm always converges. Figure 11 shows an example of the edge flipping in the envelope triangle after nodes have been added. Node 1 a is inserted and edges 1 1 1 2 1 3 a n ,a n and a n are constructed. Because the circles passing through 1  (a ,a ,n ) and (n ,n ,a ) are empty such that edges 1 3 2 3 a n and n n are regular. The circle passing through 1 2 2 (a ,a ,n ) contains node 1 n such that edge 1 2 a n must be flipped. Node 3 n is inserted and edges 1 3 n n and n a are regu-lar. The circle passing through 2 3 1 (a ,a ,n ) contains nodes 1 3 a and n such that edge 1 2 n a must be flipped. The procedure can be continued by inserting additional nodes. Figure 11: Insertion of nodes 1 2 3 a , a and a in the envelope triangle 1 2 3 (n n n ) The insertion of a node and construction of edges inside a triangle can make edges of neighboring triangles irregular. These edges must be identified and tested. A data structure is required in which the edges, whose regularity has not yet been tested, are administered. Baudson [15] has shown that this data structure is complex and expensive. In addition, not every edge, which should be flipped, is flippable. If the union of the triangles adjacent to the edge is not a convex quadrangle, the edge is not flippable. It can, however, be shown [16] that the incremental flip algorithm converges because the search for locally irregular edges is finite and does not cycle, and because at least one of the locally irregular edges is flippable.

Constrained Delaunay Triangulation
Evaluation and novel concept: The analysis of the existing methods of unconstrained Delaunay triangulation has shown that the methods are all based on the empty circle property which Delaunay discovered. The methods differ in the application of the property. The fundamental problem is evident in the graphical construction in section 3. Because all n sites must be traversed, and the number of bisectors to be tested per site is proportional to n, the complexity of the graphical method is 2 O(n ). Partitioning of the set of the sites reduces the complexity of the triangulation to O(nlog(n)), but the initial set of triangles is highly deformed and special cases must be considered to make the algorithms robust. The merging operation can be based directly on the empty circle property as shown, though Lawson [11] proceeded differently by constructing and merging Voronoi diagrams. The flipping, which is required in the merge operations of the partitioning method, is the primary operation of the second class, the flip triangulation algorithms. Their theory is well developed, but special effort is required to make the method robust, and complex data structures and procedures must be developed when then method is implemented. As shown above, the flip triangulation algorithm also uses the empty circle property.
On the basis of this evaluation, it was decided to base the constrained Delaunay triangulation on two foundations. The first is a data structure which permits convenient and explicit specification of the constraints and is adaptable as the triangulation develops. The second is direct application of the Delaunay empty circle property based on the triangulation theorem.
Data structure: Unconstrained Delaunay triangulation of a node set is confined to a domain of the plane, which is the convex hull of the node set. The triangulation is called constrained if it is restricted by one or both of the following conditions: -The triangulation is restricted to a subdomain of the convex hull, which is called its face.
-Edges are prescribed between pairs of nodes in the given set.  A data structure is required to relate the triangulation algorithm to the prescribed constraints. It is created by splitting each edge of the triangulation into two half-edges [7], which are called arrows. The arrow objects of the closed polygon on the boundary of a triangle are connected by a pointer next, as shown in figure 13. The arrows leaving a node of the triangulation are ordered by their dihedral angle about the normal vector of the plane. The arrows of a dihedral cycle are connected by a pointer rot, as shown in figure 13. This data structure can also be used to describe the dual Voronoi diagram. Quadedge data structures [17] have been developed which combine the description of the Delaunay and Voronoi subdivisions for a node set.
Triangulation theorem: If nodes R, S, T are the corners of a triangle 1 D , whose circumscribed circle is empty, and A is an arbitrary node outside the circumscribed circle, then a triangle 2 D can be constructed with node A and two of the nodes of set {R,S,T } as corners, such that the circumscribed circle of 2 D does not contain the third node in the set.     is to be determined such that the circumscribed circle of triangle 0 1 v n n n is empty.
A Cartesian local coordinate system with axes s and t is constructed with the midpoint U of arrow 0 1 n n as origin. Unit basis vector 1 g points in the direction from 0 1 n to n , and unit basis vector 2 g is constructed such that 1 2 , g g and the normal vector n of the plane form a right hand system. : A node k 0 1 n N {n ,n }   is chosen and the circle through nodes 0 1 k n ,n and n with midpoint k M is constructed as shown in figure 16. Because k n must be an interior point of the face, only points to the left of axis s are considered. Node k n is therefore rejected if the area of triangle 0 1 k n n n is negative. Depending on the shape of the face, not all points to the left of axis s are interior points. However, the point left to axis s, which is nearest to edge 0 , e is an interior point.
Let the vector from the global origin G to midpoint k M be denoted by k . m The distance from midpoint k M to nodes 0 k n and n is equal. This leads to the following equation for the circle: Vector k m from origin G to midpoint k M equals the sum of vector U x from G to the midpoint U of arrow 0 , e and vector k 2 t , g where k t is the t-coordinate of midpoint k M : Substitution of k m from equation (17) into equation (16) yields an expression for k t : -If 3 2 t t  the circumscribed circle of triangle 0 1 3 n n n is not empty because it contains node 2 n .Triangle 0 1 2 n n n remains the trial triangle.
-If 2 3 t t  the circumscribed circle of triangle 0 1 2 n n n is not empty because it contains node 3 n .Triangle 0 1 3 n n n is empty and becomes the new trial triangle.
-If 2 3 t t  the four nodes lie on a common circle. Either of the triangles can be chosen as trial triangle. The procedure is repeated for the remaining nodes in set N. After the node with the smallest coordinate k t has been determined, the corresponding triangle is added to the triangulation. This operation modifies the boundary of the face. Nodes and edges, which no longer are part of the area of the face that has not yet been triangulated, are removed. The first edge without a neighboring triangle which succeeds 0 e in the new boundary is treated in a similar manner. The process is repeated and terminates after the outer polygon has been traversed. Each of the inner polygons is then treated by analogy.
After the first traversal of the outer polygon and all inner polygons, the number of inner polygons may differ from the initial specification. If the modified face is not completely triangulated, a second traversal is initiated with the outer polygon. The traversals are repeated until the face is completely triangulated. The algorithm is illustrated in figure 20.
The addition of one triangle by the described method described can lead to the automatic creation of another triangle. It is not guaranteed that the automatically created triangle has an empty circumscribed circle. This deviation from the empty circle property is due to the constraints that have been specified, and cannot be avoided. The simplest case of automatic generation of triangles is shown in figure 18 with only 4 nodes.

Lens:
The number of nodes, for which the coordinate k t is computed, is reduced significantly by introducing a lens. Only the nodes which are located inside the lens are considered in the search. Point 2 n is chosen freely. After a trial triangle has been determined, the circumferential circle of the triangle determines the lens as shown in figure 19. With each improvement of the trial triangle, the size of the lens decreases. The points within the lens can be determined at low cost, for example by storing the node set in a quadtree [18]. The quadtree can be used to choose one of the neighbors of nodes 0 1 n or n as first trial node 2 n . The complexity of the construction of a quadtree with n entries is O(n log(n)).  Figure 20 shows a set of measuring stations at a lake with an island. The stations are treated as nodes of a constrained Delaunay triangulation, which is to be restricted to the water surface. Figure 20 illustrates the growth of the triangulation in the traversals of the polygons of its boundary. The time required for the triangulation and its graphical presentation on a conventional laptop computer is hardly noticeable.

Results and Conclusions
It has been shown that the complexity of the partitioning and the flip algorithms, and some of the problems encountered with their robustness, can be avoided by returning to the empty circle property of Delaunay and using the bisector concept illustrated in the graphical example to construct constrained triangulations incrementally. The complexity of the novel method depends on the complexity of the process which determines the subset of the given nodes that lies within a lens. Quadtrees are one of the possible tools for this task.
The extension of the method presented in this paper to polyhedral subdivisions of threedimensional space containing the points nearest to given sites is a promising area of research. This will require the replacement of the empty circle property by an empty sphere property, which has also been discovered by Delaunay. The half-edge structure for the description of the triangulation and its constraints will be replaced by a half-face structure, where the dihedral cycles at the nodes are replaced by dihedral cycles at the edges. The three-dimensional cells will be described by the intersection of an oriented outer polyhedron with a set of inner polyhedra. The circular lens will be replaced by a spherical lens, and quadtrees will be replaced by octrees to collect the nodes located inside the lens.
The data structure required for the description of three-dimensional Voronoi diagrams and Delaunay triangulations can also be useful for the description of polyhedral subdivisions of space when original buildings are mapped to digital models in building information modelling (BIM). This is an area of great practical importance and intensive current research.