Vid Domiter

vid.domiter@uni-mb.si

University of Maribor

Faculty of Electrical and Computer Science

Laboratory for
Geometric Modelling and Multimedia Algorithms

Smetanova 17, SI-2000 Maribor, Slovenia

1. Introduction

2. The algorithm

2.1. Insertion of points

2.2. Insertion of edges

2.2.1 Removal of intersected triangles

2.2.2 Triangulation of pseudo-polygons

3. Results

4. Conclusion

This paper presents an algorithm for obtaining a constrained Delaunay triangulation from a given planar graph.
The main advantage towards other algorithms is that I use an efficient Žalik's algorithm, using a plane
subdivison for obtaining a Delaunay triangulation. It is used for insertion of points into existing
triangulation. The other part of algorithm presents a method for inserting edges, already proposed by Anglada.
The algorithm is fast and efficient and therefore appropriate for GIS applications.

**Keywords: **Constrained Delaunay triangulation, two-level subdivision, computational geometry, GIS

Triangulation of set of points on the plane presents a partition of an area, surrounded by a convex hull. Triangles
are composed by vertices, given as an input. But if the triangles fulfill additional condition of an empty circumcircle,
they present a Delaunay triangulation (DT) [9]. I will go a step further, as many
have taken before, by constructing a constrained Delaunay triangulation (CDT). The input is now a planar graph
*G* which consists of vertices *V*
and non-crossing edges *E* (*G = {V, E}*). The edges *E* are actually edges of future triangles,
obtained by triangulation.

**2. The algorithm**

The presented algorithm belongs to a group of incremental insertion algorithms. Its structure consists of two parts:
**2.1. Insertion of points**

The bottleneck of incremental DT algorithm is the search for a triangle, into which the currently integrating
point falls. Žalik and Kolingerova [9] have transformed the problem into
finding the closest point, which takes less processing time. Three main steps of the algorithm are:
*initialized*. The proposed structure is
two-level uniform plane sudivision (2LUPS). The plane, which contains all points, is divided into cells
(the first level).
In case of non-uniformly distributed input points some cells can become overpopulated. Therefore an additional
variable *Treshold*) is defined. When the number of points in cell exceeds *Treshold*, cells
are divided again (the second level). 2LUPS structure is a combination of uniform subdivision (the number of
cells equals in both, x and y direction),
used for the second level, and adaptive subdivision at the first level, which uses
heuristics [9].
The most appropriate cell division in second level forms 16 new cells.
**2.2. Insertion of edges**

The currently inserting edge is defined by vertices that are already a part of triangulation.
Figure 4
shows an example of the whole process of edge insertion and in (Figure 7)
we can
see the pseudocode.
The algorithm follows two main steps:
**2.2.1. Removal of intersected triangles**

At the beginning, the first triangle must be located. This is done in constant time, because we know
which vertices define our edge, thus we can immediately locate the starting point.
Next, we must determine which triangle is the first cut by an edge. We know that one vertex
can belong to more than one triangle. So during the process of insertion of points, triangles are built,
and every point is equiped with a list of surrounding triangles.
Figure 5
shows the searching method. For each overpassed triangle, intersection
between inserted edge and edge laying opposite to starting point, is checked. If they do intersect,
our search is completed.
**Figure 7: **Algorithm for edge insertion
**2.2.2 Triangulation of pseudo-polygons**

Triangulation of pseudo-polygons (see pseudocode in Figure 10) is based
on strategy "divide and conquer" (Figure 8).
The main idea is to triangulate the both upper and lower polygons and link them when finished. As
pseudocode shows, we go through all vertices of polygon and check the empty circumcircle criterium for
each possible triangle, formed by visited vertex and edge. The checking ends when criterium is satisfied. If
currently selected point is *p*, then *p* divides the polygon into left and right
sub-polygons, which are again recursively triangulated. Recursion stops when a polygon, containing only
three vertices, is reached. Trivially, a triangle is constructed and added to triangulation (it must be linked
with its neighbouring triangles). These steps are repeated until
the whole polygon is triangulated.
**Figure 10: **Algorithm for triangulation of pseudo-polygon.

**3. Results**

Several tests have been made on AMD Celeron 1.7 GHz processor with 256 Mbytes of memory on disposal.
The operating system was Windows XP.

** Table 1: ** Meassured CPU time for tested polygons
Examples in Table 1b,c,e,f present real input data and examples in
Table 1a,d are artificially formed
and can be seen in Figure 11.

**4. Conclusion**

The algorithm copes well with real data and presents a good groundwork for advanced research. The advantages towards other
CDT algorithms are especially noticeable in cases of non-uniformally distributed input data. In future, I intend to
overcome some disadvantages of the algorithm. When edges are long, relatively to whole triangulation, triangles surrounding
them turn out very thin. These triangles are unwanted. One of the ideas is using a polyline that uses triangle edges and tries
to fit the inserted edge. Such modification provides a next challenge in my future work.

**References**

CDT algorithms can be divided in the same way as DT algorithms, into three main groups:

**Algorithms based on "Divide and conquer" strategy**

As usually, we divide the input set into smaller sets. They can be divided using different approaches, like ribbons or areas. When the sets are small enough, they are trivial to triangulate. The last and the most difficult step is merging two sets together. Those algorithms are quite demanding to implement, however, they are convenient for parallel processing (Hardwick [5] A good example was proposed by Chew [2]. A succesive refinement of the Delaunay triangulation was described by Ruppert [6], where additional Steiner points are used.**Sweep-Line algorithms**

They use an imaginary sweep-line which divides a working area into two sub-areas. If we watch the area behind the sweep-line, we can see that it has already been triangulated. Yet, the area before it still waits to be processed. Each step adds a new triangle to existing triangulation. This is achieved by connecting a new point or an edge to the boundary of the current triangulation considering the Delaunay criterium. Very popular algorithm was introduced by Fortune [3]. His algorithm was developed to construct a Voronoi diagram [1] - a dual graph of Delaunay triangulation. Shewchuk [7] presented a successful algorithm for constructing higher-dimensional CDT.**Incremental algorithms**

Are the most popular algorithms today and probably the simplest to implement. They build triangulation gradually, by inserting new vertices or edges. Every step preserves and ensures the rule of empty circumcircle (Guibas, Knuth and Sharir [4], Žalik and Kolingerova [9], Vigo [8]). There are two groups of incremental algorithms:- incremental search algorithms,
- incremental insertion algorithms.

This article presents an algorithm, which is a combination of two known algorithms (DT algorithm using a nearest-point paradigm, presented by Žalik and Kolingerova [9] and an incremental CDT algorithm explained by Anglada [8]). Sections 2.1 and 2.2 describe each idea in details. At the end, a conclusion, based on results will be given.

- Initialization,
- Triangulation,
- Finalisation.

Due to the insertional nature of this algorithm, determination of artificial triangle is to be done first. The triangle has to be constructed large enough, so all points fall inside of it.

Insertion of the first point is handled under initialization. The result is a big triangle, split in three smaller ones. This step is trivial, because triangulation is definitely a DT (Figure 1a).

**Figure 1:** Inserted point *p _{i}* splits one triangle and forms three new triangles (a) or splits two and
forms four new triangles (b), when it falls on the common edge.

*Triangulation* is the main process. It consists of two steps.
The first step is finding a triangle containing inserted point. This step uses a 2LUPS searching structure and a
searching mechanism, usually performed by so-called spiral algorithm, which is used here[9].
First, the point is inserted into a cell of 2LUPS structure. Then, a search for the closest vertex is started
(Figure 2).
Here, we can see the idea presented by \v Zalik and Kolingerova [9].

**Figure 2:** Finding the closest point.

In Figure 2, inserted point is *p _{i}*. A distance calculate a distance

After the closest vertex is found, it is possible, that it belongs to one of the triangles containing the inserted point. If not, the next closest vertex has to be found. This is easy, because the candidates are previously examined vertices.

The second step is splitting the corresponding triangle into three or four sub-triangles (Figure 1a and Figure 1b), which are then checked for empty circumcircle criterium (Figure 3). If a triangle fails the checking, it must be legalized. This means that a common edge of the tested triangle and its edge-neighbour must be flipped to obtain two new triangles. Again, they must be tested. The procedure is implemented recursively.

**Figure 3:** Triangle fails Delaunay criterium and edge flip must be done.

The final step is *finalization*. Its basic task is a removal of triangles containing the vertices of
artificial triangle, formed in *initialization*.
The algorithm works with time complexity of *O(n ^{1,1})* for most cases. The worst time
complexity is

- Removal of triangles cut by edge
*ab*to gain an empty area around the edge (Figure 4b,c) ,
- Triangulation of the previously created area, so called pseudo-polygons (Figure 4d) .

**Figure 4:** Situation before insertion of edge *ab*(a), removal of intersected triangles (b) forms two pseudo-polygons
P_{U} and P_{L} (c) and triangulation of P_{U} and P_{L} (d).

**Figure 5: **Cycling through surrounding triangles until the first triangle is found, in our example triangle *t _{1}*.

Now, ground for triangle removal has been arranged. Moving through triangles is simple
(Figure 6). To perform this
task, the proposed functions *OpposedTriangle* and *OpposedVertex* are implemented.

**Figure 6: **Moving through triangles is stressed by arrows.

Procedure InsertEdgeCDT(T:CDT, ab:Edge){Precondition: a,b in(T) and ab not in(T)}Find the triangle t in(T) that contains a and is cut by ab P_{U}:=EmptyList P_{L}:=EmptyList v:=a While v not in t do t_{seq}:=OpposedTriangle(t,v) v_{seq}:=OpposesdVertex(t_{seq},t) If v_{seq}above the edge ab then AddList(P_{U},v_{seq}) v:=Vertex shared by t and t_{seq}above ab Else If v_{seq}below the edge ab then AddList(P_{L}, v_{seq}) v:=Vertex shared by t and t_{seq}below ab Else {v_{seq}on the edge ab} InsertEdgeCDT(T, v_{seq}b) a:=v_{seq}break EndIf Remove t from T t:=tseq EndWhile TriangulatePseudoPolygon(P_{U},ab,T) TriangulatePseudoPolygon(P_{L},ab,T) Reconstitute the triangle adjacencies of T Add edge ab to T Mark the edge ab from T as fixed EndProc

Let us look at the triangles *T _{papbpc}* and

When triangle is deleted, all links must be updated. During this process, all points above and below edge
are stored separately in two lists,
*P _{U}* and

**Figure 8: **Triangulation of pseudo-polygon: Triangle *T _{abc}* fulfills Delaunay criterium, so
point

When inputing data, it can sometimes occur, that a point falls directly on the edge. So pseudocode is
slightly improved, compared to Anglada's [8], because it can handle such cases. Point *{p}*
appears to split the edge *ab*, so the algorithm presumes the same. Recursively, it calls a
procedure for edge insertion again, yet with a new edge *pb* (edge *ab* is now treated as
*ap+pb*). Comparisson is shown in Figure 9.

**Figure 9: **Example shows a situation, where point *p* falls on the edge *ab*.
In picture a), the situation is not handled, therefore the CDT is incorrect. Picture b) shows our refinement,
where edge *ab* is splitted to new edges *ap* and *pb*), both treated separately. CDT is correct.

Procedure TriangulatePseudoPolygon(P:VertexList, ab:Edge, T:CDT) If P has more than one element then c:=First vertex of P For each vertex v $\epsilon$ P do If v E CircumCircle (a, b, c) then c:=v EndIf EndFor Divide P into P_{E}and P_{D}giving P=P_{E}+{c}+P_{D}TriangulatePseudoPolygon(P_{E}, ac, T) TriangulatePseudoPolygon(P_{D}, cd, T) EndIf If P is not empty then Add triangle with vertices a, b, c into T EndIf EndProc

I will pass the same conclusions about time complexity for edge insertion and for triangulation of
pseudo-polygons as Anglada [8] stated. Let *e* be the number of triangles of CDT cut by
edge *ab* and *m* number of edges. The edges are set in a way, that every edge is
defined by two vertices which are already in CDT. So finding a starting point a is done in constant
time. Number of surrounding triangles can vary between one and maximum, as in the case of the number of edges of convex
polygon with a point in the center. Finding the first triangle then requires *O(m)* time.
Construction of upper and lower pseudo-polygons takes a time complexity of *O(e ^{2})*,
since in each recursive call, total number of points decreases in one unit.
To put all together, procedure

I tested the algorithm on graphs, using different numbers of points and edges. The procedures for insertion of edges and points were measured separately. Table 1 shows some I observed CPU time in relation with number of inserted edges. Usually, greater number of edges means longer processing. But not always. Another important factor is the number of triangles being cut by edges. Tested example d) has less edges than example e) (Table 1), but takes the same time for edge processing. The reason lies in greater number of removed triangles. Removal itself is not that time demanding. More important is the number of points that compose pseudo-polygons, which is proportional to number of removed triangles.

Figure 11 | Edges | Points | Triangles cut | InsertPoint | InsertEdge | Total |
---|---|---|---|---|---|---|

(a) | 368 | 1200 | 988 | 0.01 | 0.01 | 0.02 |

(b) | 4259 | 1856 | 1311 | 0.02 | 0.02 | 0.04 |

(c) | 5453 | 2345 | 1108 | 0.03 | 0.02 | 0.05 |

(d) | 10000 | 10000 | 9462 | 0.08 | 0.05 | 0.13 |

(e) | 13667 | 6132 | 2412 | 0.08 | 0.05 | 0.13 |

(f) | 19044 | 16013 | 9092 | 0.211 | 0.09 | 0.301 |

**Figure 11: **Some tested examples.

**[1]** Aurenhammer, F., Voronoi diagrams-a survey of a fundamental geometric data
structure, ACM Computing Surveys, vol. 23, no. 3, pp. 345-405, 1991.

**[2]** Chew, L. P., Constrained Delaunay triangulations, Proceedings of the 3rd annual
symposium on Computational geometry, ACM Press, pp. 215-222, 1987, Waterloo, Ontario, Canada.

**[3]** Fortune, S., A sweep line algorithm for voronoi diagrams, Algorithmica, vol. 2,
1987, pp. 153-174.

**[4]** Guibas L., Knuth D., Sharir M., Randomized incremental construction of Delaunay and
Voronoi diagrams, Algorithmica, no. 7, pp. 381-413, 1992.

**[5]** Hardwick, J. C., Implementation and evaluation of an efficient parallel Delaunay
triangulation algorithm, Proceedings of the 9th annual ACM symposium on Parallel algorithms and architectures,
ACM Press, pp. 239-248, 1997, Newport, Rhode Island, USA.

**[6]** Ruppert, J., A new and simple algorithm for quality 2-dimensional mesh generation,
Proceedings of the 4th annual ACM-SIAM Symposium on Discrete algorithms, ACM Press, pp. 83-92, 1993, Austin,
Texas, USA.

**[7]** Shewchuk, J. R., Sweep algorithms for constructing higher-dimensional constrained
Delaunay triangulations, Proceedings of the 16th annual symposium on Computational geometry, ACM Press, pp.
350-359, 2000, Clear Water Bay, Kowloon, Hong Kong.

**[8]** Anglada, M. V., An improved incremental algorithm for constructing restricted
Delaunay triangulations, Computers \& Graphics, vol. 21, p. 215-223, 1997.

**[9]** Žalik, B., Kolingerova I., An incremental construction algorithm for Delaunay
triangulation using the nearest-point paradigm, Int. J. Geographical information science, vol. 17, no. 2, pp.
119-138, 2003