Constrained Construction of Planar Delaunay Triangulations without Flipping

Full Text

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 (18051859) 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 (1868-1908) when he extended the Dirichlet tessellations to higher dimensions and Boris Delaunay (18901980) 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: Voronoi diagram in a plane 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. Figure 2: Delaunay triangulation for the Voronoi sites in Figure 1 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. 24. 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: f : D ® V (1) f bijective mapping for duality (set of ordered pairs) D elements of the Delaunay triangulation (domain of the mapping) V elements of the Voronoi diagram (target of the mapping) The bijective mapping function is defined as follows: f (primal node n ) = dual region r f -1(dual region r ) = primal node n (2) i i i i f (primal edge e ) = dual arc a f -1(dual arc a ) = primal edge e k k k k f (primal cell c ) = dual vertex v f -1(dual vertex v ) = primal cell c m m m m Figure 3: Duality between a Voronoi diagram and a Delaunay triangulation Figure 3 shows mapping (2) for the Voronoi diagram in figure 1 and the Delaunay triangulation in figure 2. The reversal of the inclusion relationship is illustrated for cell c2 in the primal subdivision, which contains nodes n1, n4 and n5. Mapping f yields the following images: f(c2 ) = v2 f(n1) = r1 f(n4 ) = r4 f(n5 ) = r5 (3) Regions r1, r4 and r5 , which are the duals of the nodes n1, n4 and n5 contained in cell c2 , all contain vertex v2 , which is the dual of cell c2 . The reversal of the inclusion is illustrated by edge e8 in the primal subdivision, which contains nodes n4 and n5. Mapping f yields the following images: f(e8 ) = a8 f(n4 ) = r4 f(n5 ) = r5 (4) Regions r4 and r5 , which are the duals of nodes n4 and n5 contained in edge e8 in the primal subdivision, all contain arc a8 , which is the dual of edge e8. 25. 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 figure 4. Figure 4: Graphical construction of a Delaunay triangulation Let the coordinates of a set of sites be given. Following the empty circle concept of Delaunay, an empty circle passing through three of the sites is chosen by inspection. The midpoint v1 of this circle is a Voronoi vertex. The edges of the triangle are Delaunay edges and the sites at the corners of the triangle are Delaunay nodes. The three bisectors of the edges pass through v1 and contain Voronoi arcs with v1 as end point. The left diagram in figure 5 shows the empty circle and the right diagram the subdivision of the plane into three point sets, each of which is nearest to one of the sites of the initial Delaunay triangle. The second vertex of the Delaunay arcs radiating from v1 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 n1n4 is a1a2 . Let the bisector leaving v1 in the direction of point a2 be followed. This bisector is the boundary that separates points that are nearer to site n1 from points that are nearer to site n4 . After this arc has reached the point of intersection with bisector f1 f 2 , it enters a set of points which are closer to site n2 than to sites n1 and n4 . The point of intersection is therefore the second vertex v2 of the arc. If the arc leaving v1 in the direction of point e2 is followed, none of the bisectors in figure 4 is intersected. This arc is infinite. Similarly, the arc leaving v1 in the direction of point b2 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. Figure 5: Initial step of the graphical construction The vertices, which are constructed on the bisectors leaving the first vertex, are collected in a list. After all bisectors at node v1 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 n1n4 is an edge because it is bisected by arc v1v2 . Segment n2n3 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 ak overlays the bisector of edge ek . Nodes ns and nt of edge ek are the duals of regions rs and rt whose common boundary is arc ak . Vertices vg and vh of arc ak are the duals of cells cg and ch whose common boundary is edge ek . 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. 26. 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 x1 of the sites in the first group are less than the coordinates x1 off all other sites, the coordinates x1 of the sites in the second group are less than the coordinates x1 of the remaining sites, and so on. Figure 6 shows a given set of 6 sites ni 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 x2 of the three nodes differ significantly. It is assumed that the sites are not collinear. Special methods are required if the sites are collinear. Figure 6: Construction of a Delaunay triangulation by partitioning and merging 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 n1 and n5 . 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 b1|3 and b1|5 , where bk|m denotes the bisector of the line segment connecting sites nk and nm. The bisectors intersect in a Voronoi vertex (not shown in D3) which is the midpoint of a circle through sites n1,n3 and n5. The triangle with corners n1, n3 and n5 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 310 = 59049 sites, 10 - 1 = 9 quired. The complexity of the method is O(n log(n)). levels of merging are re- Figure 7: Flipping of edges in the merge operation 27. Construction of Delaunay Triangulations by Flipping Figure 7 shows that part of the merging process can be achieved by selecting a quadrangle and replacing its diagonal by the other diagonal. In diagram D2 of figure 7 diagonal e1 | 4 can be replaced by diagonal e2|5 in the quadrangle with nodes n1, n5 ,n4 and n2 , where ek|m is the edge connecting nodes nk and nm. In the literature, this operation is referred to as edge flipping in the quadrangle. Diagram D4 in the same figure shows that flipping can be accompanied by the construction of an additional edge. 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. Angle vector: Different triangulations are compared by means of an angle vector. Let the number of triangles in the triangulation be m. The 3m angles of the triangles are arranged in ascending order in in a vector t which is called the angle vector of the triangulation. The angle vector t of a triangulation T is compared row-wise with the angle vector t of another triangulation T of the same node set. The angles in vectors t and t are denoted by ai and ai . The comparison starts in with the entries in row 0. The first pair of angles (ak ,ak ) which are not equal determines the lexicographic order of the two vectors. If ak > ak , then vector t is larger than vector t. a1 a2 ... ak ... a3m t = t = t > t Þ " (a = a ) Ù a > a (5) a1 a2 ... ak ... a3m i<k i i k k 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 t . c c b b a d a d edge ab edge cd Figure 8: Flipping edge ab to edge cd in a quadrilateral abcd t ³ t Þ edge ab is regular (6) If abcd is a rectangle, both edges are regular. Otherwise, one edge is regular and the other edge is irregular. Proof: Figure 9 shows two triangulations for a quadrilateral with nodes nk . The triangles formed by diagonal n1n3 are shown in diagram (a), the triangles formed by diagonal n2n4 in diagram (b). Diagram (a) contains the circumscribed circle of triangle n1n2 n3 . Diagram (b) contains the circumscribed circle for triangle n1n2 n4 . The extended diagonals n1 n3 and n2 n4 intersect the circumscribed circles in point u. It will be shown that one of the edges n1n3 and n2 n4 is regular and the other irregular. Edge n2n4 of triangulation (b) is regular if the six angles a + e,b,h,c,d + f ,g of the two triangles in (b) are larger than the smallest angle j in triangulation (a). The geometry of the triangles and of the circumscribed circle of triangle n1n2 n3 yields the following relations between the angles in triangulation (a): e + x f + y = c x + d = h (7) = b y + a = g Similarly, the geometry of the triangles and the circumscribed circle of triangle n1n2 n4 yield the following relations between the angles in triangulation (b): c - s g - t = e s + d = a t + f = h (8) = b u d a x c e h g a e h g a a t b y b f b h u b e f c d c s d Figure 9: Regularity test for the diagonals of a convex quadrilateral It follows from (8) that angles b + c and g + h cannot be the smallest angle in triangulation (a): b + c g + h = e + f + x + y = a + d + x + y > e (9) > a Each of the remaining four angles a,d,e,f can be the smallest angle in triangulation (a). It will now be shown that in all four cases, each of the six angles a + e,b,h,c,d + f ,g in triangulation (b) is larger than the smallest angle in triangulation (a). Case 1: j = a < d,e,f : (10) Case 2 : j = d < a,e,f : (11) a + e > b = a f + y > a + y > a a + e > b = d + e > f + y > d d + y > d h c d + f g = s + d > = e + x > > d > = t + a > d > a a + x > a a a h c d + f g = s + d > = e + x > > d = t + a > d d + x > d t + d > d Case 3 : j = e < a,d,f : (12) Case 4 : j = f < a,d,e : (13) a + e > b = d + e > f + y > e e + y > e a + e > b = f + e > f f + y > f h c d + f g = s + d > = e + x > > e + f > = t + a > s + e > e e e t + e > e h c d + f g = s + d > = e + x > > f = t + a > s + f > f f + x > f t + f > f The proof shows that edge n2 n4 is regular. Generally, the regular edge in a quadrilateral can be recognized by two attributes: Empty circle: Edge n1 n3 in triangulation (a) is irregular. Node n4 lies inside the circumscribed circle and the circle is therefore not an empty circle. Edge n2 n4 in triangulation (b) is regular. Node n3 lies outside the circumscribed circle, which does not contain any nodes and is therefore an empty circle. Shortest diagonal: The regular diagonal n2 n4 is shorter than the irregular diagonal n1n3 . This follows from triangulation (a) where triangles n1n2 n3 and n1n2u are congruent, such that chords n1 n3 andn2u are of equal length. Node n4 is an inner point of chord n2u such that the length of n2n4 is less than the length of n2 u, and therefore less than the length of n1n3 . 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. Figure 10: Flipping edges in a the total flip algorithm 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 a1 is inserted and edges a1n1,a1n2 and a1n3 are constructed. Because the circles passing through (n1,n2,a1),(n2,n3,a1) and (n3,n1,a1) are empty, the three edges are regular and there is no flipping. Node a2 is inserted and edges a1a2,n2 a2 and n3 a2 are constructed. The circles passing through (a1,a2,n3 ) and (n2,n3,a2 ) are empty such that edges a1n3 and n2 n3 are regular. The circle passing through (a1,a2 ,n2 ) contains node n1 such that edge a1n2 must be flipped. Node n3 is inserted and edges (n1a3 ),(n2 a3 ) and (a2 a3 ) are constructed. The circles passing through (n1,a3,n2 ) and (n2,a2,a3 ) are empty such that edges n1 n2 and n2 a2 are regu- lar. The circle passing through (a2,a3,n1) contains nodes a1 and n3 such that edge n1a2 must be flipped. The procedure can be continued by inserting additional nodes. Figure 11: Insertion of nodes a1, a2 and a3 in the envelope triangle (n1 n2 n3 ) 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. 28. 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 O(n2 ). 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: o The triangulation is restricted to a subdomain of the convex hull, which is called its face. o Edges are prescribed between pairs of nodes in the given set. Figure 12: Description of a face as the intersection of the interiors of polygons A face is defined by means of a set of polygons as shown in figure 12. The oriented polygons are not self-intersecting and do not intersect or touch each other. The outer polygon contains all of the inner polygons. An inner polygon does not contain any other polygon. The area of the face is defined as the intersection of the inner point sets of the polygons. The edges which are prescribed in the face may not intersect each other or a polygon of the face. Figure 13: Arrow data structure for a triangulation 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 D1, whose circumscribed circle is empty, and A is an arbitrary node outside the circumscribed circle, then a triangle D2 can be constructed with node A and two of the nodes of set { R,S,T } as corners, such that the circumscribed circle of D2 does not contain the third node in the set. Figure 14 shows an example of triangles D1 and D2 and their circumscribed circles. Figure 14: Given triangle D1 and constructed triangle for the triangulation theorem Proof: The midpoint M1 of the circumscribed circle of the given triangle D1 is chosen as origin of a Cartesian coordinate system (x,y). The corners of triangle D1 are named R,S,T in the anti-clockwise direction. Choose an edge ST of triangle D1 such that the third corner R and point A are located on opposite sides of the line containing edge ST, as shown in figure 15. Axis x points from origin M1 to the midpoint M2 of edge ST. Figure 15: Geometric construction for the proof of the triangulation theorem The midpoint of the circle which passes through points A,T and S is denoted by M2. The radii of the circles with midpoints M1 and M3 are denoted by r1and r2 , the distance between points R and M3 by z. Point R is located outside the circumscribed circle of triangle ATS if the condition z > r2 is satisfied. The rule of Pythagoras is used in the following derivations. z2 = 3 1 1 (x - x )2 + y2 = x2 - 2x x + (x2 + y2 ) 3 1 3 1 1 = x2 - 2x x § r2 3 1 3 1 r2 = (x - x )2 + y2 2 3 2 2 = x2 - 2x x + (x2 + y2 ) 3 2 3 2 2 = x2 - 2x x § r2 3 2 3 1 z2 - r2 = 2 x (x - x ) (14) 2 3 2 1 The counter-clockwise orientation of boundary RST and the choice of the coordinate system assure that x2 > x1 . Because edge ST was chosen such that points R and A lie on opposite sides of the line passing through points S and T, and because point A is located outside the circumscribed circle of triangle RST, coordinate x3 of the midpoint of the circumscribed circle of triangle AST is positive. It follows from expression (14) that z that point R lies outside the circumscribed circle AST. This proves the theorem. > r2 , such Basic construction step: Let a domain, which is to be triangulated, be described by a set of nodes N = {n0, ...,nm-1) and a set of arrows E = {e0,...,eq-1) on its boundary. Assume without loss of generality that arrow e0 points from node n0 to node n1 . A node nv Î N - {n0,n1} is to be determined such that the circumscribed circle of triangle n0 n1nv is empty. A Cartesian local coordinate system with axes s and t is constructed with the midpoint U of arrow n0 n1 as origin. Unit basis vector g1 points in the direction from n0 to n1, and unit basis vector g2 is constructed such that g1,g2 and the normal vector n of the plane form a right hand system. Figure 16: Geometric construction for the basic triangulation step 1 g = x1 - x0 = c -d := g = (15) x1 - x0 d 2 c 1 x1 - x0 (x - x )2 + (y - y )2 1 0 1 0 y1 - y0 A node nk Î N - {n0,n1} is chosen and the circle through nodes n0,n1 and nk with midpoint Mk is constructed as shown in figure 16. Because nk must be an interior point of the face, only points to the left of axis s are considered. Node nk is therefore rejected if the area of triangle n0 n1nk 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 e0 , is an interior point. Let the vector from the global origin G to midpoint Mk be denoted by mk . The distance from midpoint Mk to nodes n0 and nk is equal. This leads to the following equation for the circle: ( m - x )T ( m - x ) = ( m · x )T ( m · x ) k 0 k 0 k k k k 2(x · x )T m = xT x · xT x (16) k 0 k k k 0 0 Vector mk from origin G to midpoint Mk equals the sum of vector xU from G to the midpoint U of arrow e0, and vector tkg2, where tk is the t-coordinate of midpoint Mk : m = x o t g = 1 (x o x ) + t g (17) k u k 2 2 0 1 k 2 Substitution of mk from equation (17) into equation (16) yields an expression for tk : (x - x )T (x · x ) t = k 0 k 1 (18) k 2(x k 0 2 · x )T g Algorithm: The triangulation is initiated by choosing a boundary arrow e0 of the outer polygon with nodes n0 and n1, and a node n2 Î N - {n0,n1} is selected. Coordinate t2 is computed with expression (18). Triangle n0 n1n2 is called the trial triangle. The next test point n3 is chosen in set N - {n0,n1,n2 } and its coordinate t3 is computed. The values of t2 and t3 determine whether triangle n0 n1n3 replaces n0 n1n2 as trial triangle. The test is illustrated in figure 17: · If t3 > t2 the circumscribed circle of triangle n0 n1n3 is not empty because it contains node n2 . Triangle n0 n1n2 remains the trial triangle. · If t2 > t3 the circumscribed circle of triangle n0 n1n2 is not empty because it contains node n3. Triangle n0 n1n3 is empty and becomes the new trial triangle. · If t2 = t3 the four nodes lie on a common circle. Either of the triangles can be chosen as trial triangle. Figure 17: Empty circle test for a trial triangle The procedure is repeated for the remaining nodes in set N. After the node with the smallest coordinate tk 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 e0 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. Figure 18: Automatic creation of triangle t2 after the construction of triangle t1 Lens: The number of nodes, for which the coordinate tk is computed, is reduced significantly by introducing a lens. Only the nodes which are located inside the lens are considered in the search. Point n2 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 n0 or n1 as first trial node n2 . The complexity of the construction of a quadtree with n entries is O(n log(n)). Figure 19: Lens for trial triangle n0 n1nk Example: 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. Figure 20: Triangulation of measuring stations at a lake 1. 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.

About the authors

Vera V Galishnikova

Peoples’ Friendship University of Russia (RUDN University)

Email: galishni@gmail.com
Associate Professor, Director of the Department of Architecture and civil engineering, Engineering Academy, RUDN University. Research Interests: Computational Civil Engineering, Building information modeling, Topological computer models of buildings, Computational geometry, Computational mechanics of complex steel structural systems - latticed plates and shells, thin-walled plate and plate-rod structures. Nonlinear finite element analysis of space frames. Nonlinear stability of structures 6 Miklukho-Maklaya St., Moscow, 117198, Russian Federation

Peter Jan Pahl

Technische Universität Berlin

Email: pahl@ifb.bv.tuberlin.de
Prof. Dr. Dr. h. c. mult., Department of Civil Engineering, Technical University Berlin (TUB). Research Interests: Mathematical modeling and optimization of comple[ structural systems, Computational Civil Engineering, Building information modeling, Topological computer models of buildings, Computational geometry, Computational mechanics of complex steel structural systems - latticed plates and shells, thin-walled plate and plate-rod structures. Nonlinear finite element analysis of space frames. Nonlinear stability of structures 17 Juni Str., 135, 10623, Berlin, Germany

References