- Email: [email protected]
A universal trapezoidation algorithm for planar polygons
Computers & Graphics 23 (1999) 353}363
Technical Section
A universal trapezoidation algorithm for planar polygons Borut Z[ alik *, Gordon J. Clapworthy University of Maribor, Faculty of Electrical Engineering and Computer Science, Smetanova 17, 2000 Maribor, Slovenia De Montfort University, Department of Computer and Information Sciences, Hammerwood Gate, Kerts Hill, Milton Keynes, MK7 6HP, UK
Abstract A new algorithm for decomposing non-monotone planar polygons, which may contain holes, into trapezoids is described. The holes may be nested and may have common edges. In the "rst part of the paper, the main idea is explained for non-monotone polygons without holes, and the algorithm is then extended to polygons containing holes; the holes can also be decomposed into trapezoids, if desired. Finally, it is shown that the algorithm performs trapezoidation of general polygons in O(n log n) time, where n is the common number of polygon vertices. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Polygon decomposition; Trapezoidation; Computational geometry
1. Introduction Techniques for decomposing complex geometric shapes into simpler components are of great importance in computational geometry and its applications. Depending upon the particular problem, there are di!erent ideas about how this should be performed and the choice of element into which the shape should be decomposed. For example, planar polygons can be decomposed into triangles, trapezoids, or even star-shaped polygons [1]. While triangulation is, without doubt, the most popular and the most investigated technique, splitting a complex polygon into trapezoids is an alternative, and a generalised method for performing trapezoidation on an arbitrary polygon with any number of holes, nested to any level, is considered in the paper. The motivation for the current work arose from the display of land use within a given geographical area. The available software could colour triangles or trapezoids to display di!erent uses but was more e$cient with trap-
* Corresponding author. Tel.: #368-62-220-7471; fax: #386-62-211-178. E-mail address: [email protected] (B. Z[ alik)
ezoids because of the reduced number of primitives involved. The individual areas to be displayed could be nested, and their boundaries could coincide with each other or the boundary of the overall area. However, there is no reason why the application should not be used more broadly, particularly to display planar variation of physical quantities. It can also be applied to "nd, e$ciently, the intersection and/or union of polygons. Trapezoidation was "rst considered by Chazelle and Incerpi [2] and Fournier and Montuno [3]. A randomised algorithm for trapezoidation was later suggested by Seidel [4] as the key for triangulation and this was recently used by Narkhede and Manocha [5] for the same purpose. However, none of these authors considered general planar polygons containing holes; indeed, to the authors' knowledge, no previously published method for trapezoidation has successfully addressed polygons with holes, nor the trapezoidation within holes. A trapezoid is a four-sided polygon in which two of the edges are parallel. Note that a triangle can be considered as a degenerate trapezoid in which one of the parallel edges has zero length; hence, some triangles might result from applying any trapezoidation algorithm. A horizontal trapezoidation of an n-sided polygon involves drawing horizontal lines (we shall call them scan lines) through every vertex of the polygon (see Fig. 1). The
0097-8493/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 9 3 ( 9 9 ) 0 0 0 4 4 - 8
354
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
is swept vertically down the plane and an initial data structure is updated at discrete `eventsa during this process. In this algorithm, an event is the intersection of the sweep line with one or more polygon vertices. By maintaining the vertices in a dynamic ordered list, only the leading element has to be checked for intersection at any time. The precise action taken when a new vertex is intersected depends upon the vertex type } these types are identi"ed in Section 2.1. The "rst step in the algorithm considers all intersections of the current scan line with the polygon vertices and edges, sorting them into x-increasing order. Those which are relevant to trapezoidation between the current scan line and the previous scan line are retained, the others are discarded. The second step identi"es the points on the current and previous scan lines which are to be joined to create the trapezoids. Fig. 1. Horizontal trapezoidation of a monotone polygon.
whole operation can be completed in O(n log n) time if the polygon is monotone [4] } in the context of horizontal trapezoidation, a monotone polygon is one in which the boundary consists of two vertically monotone chains. In practice, polygons are frequently non-monotone and may enclose a number of holes. This paper describes a universal algorithm to perform the trapezoidation of a non-monotone planar polygon containing holes. In the "rst part of the paper, the main idea of the algorithm is explained for non-monotone polygons without holes. In the following section, the algorithm is extended to polygons containing holes; the holes can themselves be trapezoided by the algorithm, if required. A straightforward application of the algorithm results in more trapezoids than strictly necessary, and Section 4 indicates how this number can be reduced. Section 5 provides suggestions for an e$cient implementation, Section 6 contains an estimation of the time complexity of the algorithm, and Section 7 describes tests made with both arti"cially constructed polygons and real data from the geographical application.
2. Polygons without holes The algorithm uses the well-known greedy approach [6,7]. In the "rst step, the algorithm prepares the input data such that in the second step the problem can be solved in one pass. Adopting vocabulary from topology [8], the sequence of vertices forming the border of a polygon is named a loop, and a sequence of vertices forming the border of an inner closed polygonal region is called a ring. Initially, the polygon vertices are sorted by y-coordinate. A `plane sweepa is employed, which is a technique used in many geometrical applications } a horizontal line
2.1. Vertex classixcation The solution presented is based on a simple analysis of each polygon vertex, a brief description of which is given by O'Rourke [9], without implementation details. However, following Seidel [4], O'Rourke explicitly ignores the possibility that a scan line may intersect more than one vertex; if this should occur, one of the vertices is notionally moved to an additional scan line. For clarity, our algorithm includes the possibility of multiple vertices on a scan line, so the range of possible vertex types is larger. The following types of relationship can exist between a polygon vertex p and its neighbouring vertices p J and p : P E Intersection: when one of neighbouring vertices has a greater and the other a lower y-coordinate than p (Fig. 2a) p is handled as if it were an intersection of the scan line with a polygon edge, so it is categorised as INT. E Local minimum: when both neighbouring vertices have greater y-coordinates than p (see Fig. 2b), p is categorised as MIN. E Local maximum: when both neighbouring vertices have lower y-coordinates than p (see Fig. 2c), p is categorised as MAX. E Horizontal line seument: when one of the neighbours has the same y-coordinate as p within tolerance e, then p is an end point of a horizontal line segment; by considering the position of the second neighbour, vertex p is given the code E 䡩 HMAX if the second neighbour has a lower y-coordinate (Fig. 2d) or E 䡩 HMIN if the second neighbour has a higher y-coordinate (Fig. 2e). E If both neighbours have the same y coordinate as p within tolerance e (Fig. 2f ), then p will not form part of any horizontal trapezoid generated later and it can be erased from the list of polygon points at this stage.
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
355
Fig. 2. The types of local behaviour of vertex p.
Fig. 3. A non-monotone polygon.
2.2. Initialisation
E if p is of type MAX, it is inserted twice into G array OldPoints
Consider the polygon in Fig. 3 in which the vertices are denoted by p and indicated by circular markers. The G polygon is non-monotone and contains vertices belonging to all of the classi"cations illustrated in Fig. 2, apart from the uninteresting case (f ) in which the vertex would be trivially discarded at the outset. The vertex coordinates are stored sequentially in an array named InputPoints. The vertices of the polygon are sorted with respect to their y-coordinates and stored in an array called SortedVerticesY. It is possible, of course, to have more than one vertex with the same y coordinate. The vertices lying on the "rst scan line (line s in Fig. 3) are then considered; clearly, they cannot be of types INT, MIN or HMIN. The vertices are sorted into x-increasing order and stored in an array named OldPoints using the following procedure:
E if p is of type HMAX, it is inserted once into G OldPoints,
(P1)
E vertices within a horizontal line segment are ignored. For the polygon in Fig. 3, OldPoints now contains p p p p . 2.3. Processing a scan line In the main part of the algorithm, the scan lines s (k'0) are considered sequentially, in descending orI der of y-coordinate. If there are n polygon vertices, there will be, at most, n scan lines and, hence, n!1 steps in this part of the algorithm.
356
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
In Figs. 3, 4 and 6, intersections between scan lines and polygon edges, qH, are identi"ed by square markers; here, G the subscript denotes the ith scan line and the superscript indicates the jth intersection on the line. To create an ordered list of intersections of the current scan line with the polygon, the following actions are needed at s : I E determine the number and the types of polygon vertices lying on s , I E calculate all intersections between s and polygon I edges } these points are of type INT; for s in Fig. 3, we obtain two intersection points q and q, E insert the vertices and the intersections into an array named Intersections in order of increasing x-coordinate; this array will be used twice, E create an array CurrentPoints from Intersections using the procedure: E MIN: add point twice to CurrentPoints, E INT, HMIN: add point once to Current Points,
(P2)
E MAX, HMAX: do not add point. In Fig. 3, s passes through only one vertex } p which is of type MIN. After sorting, the array CurrentPoints contains the points q p p q. Because of the way in which the algorithm is constructed, the arrays OldPoints and CurrentPoints will now contain equal (even) numbers of points, and trapezoids can be generated by joining corresponding pairs of points from the two arrays. Between s and s , two trapezoids are created, determined by the points: p qp p and p p qp . It can be seen that the "rst trapezoid has degenerated into a triangle. 2.4. Updating for the next scan line Before moving to the next scan line, we update OldPoints from the points stored in Intersections according to their types, as follows: E MAX: add point twice to OldPoints, E INT, HMAX: add point once to OldPoints,
(P3)
E MIN, HMIN: do not add point. For s , Intersections contained q p q, which have types INT, MIN, INT, respectively, so OldPoints becomes q q. Scan lines continue to be processed in this way until the polygon is exhausted. The basic algorithm described above, involving processes familiar from previous decomposition algorithms and from scan-line "ll algorithms, is generalised in Section 3 to cope with complex geometries in which the
polygon contains several rings, including the possibility of multiply nested rings and partially coincident rings.
3. Polygons with holes This section considers the generalisation of the trapezoidation algorithm for all planar polygons which do not self-intersect. For simplicity, we deal initially with a situation in which the loop and all rings are totally disjoint and all areas within the loop are to be trapezoided, including the interiors of the rings. This is later generalised to allow parts of the rings to coincide with each other and with the loop, and to allow di!erent areas to be omitted from the trapezoidation. 3.1. Simple nested rings The initialisation takes place as previously, with the array SortedVerticesY being generated and the array OldPoints being initialised as before (because the point with the highest y-coordinate of the polygon cannot belong to any of the rings). The changes to the algorithm are introduced in the subsequent stages, described previously in Sections 2.3 and 2.4; these will be illustrated using the polygon in Fig. 4. It can be seen that the polygon of Fig. 4 contains another polygon, de"ning the ring R . The algorithm requires that two new arrays, Ring-j}OldPoints and Ringj}CurrentPoints, are introduced for each ring, R . H Prior to reaching scan line s in Fig. 4, the algorithm behaves as described previously, and the array OldPoints contains q q. At the next scan line, s , the "rst intersection with a ring is encountered. Array CurrentPoints is generated using (P2) from Section 2.3, and it contains, as previously, q q, so the trapezoid q q q q is again produced. Note that all ring vertices on s are of types MAX or HMAX and are thus ignored within procedure (P2). The "rst change comes when updating before processing the next scan line, s . Intersections with the loop are again considered for the updating of OldPoints using (P3) from Section 2.4, but intersections of the scan line with a ring R are now considered for two arrays: H E the array Ring-j}OldPoints belonging to R and H E the array OldPoints for the loop if R is immediately H within the loop, or the array Ring-k}OldPoints if R is H nested immediately within another ring, R . I Similarly, at subsequent scan lines, in addition to the array CurrentPoints, for each ring R an array RingH j}CurrentPoints is created using (P2), so that ring intersections with the current scan line also occur in two arrays.
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
357
Fig. 4. A non-monotone polygon containing a hole.
Whether the current ring is immediately within the loop or another ring can be determined without additional calculations simply by observing the neighbouring points of intersection with the scan line. This is discussed further in Section 5.2. In the example, in moving from s to s , the arrays are as follows: OldPoints: q p p p p p p q, Ring-1}OldPoints: p p p p , CurrentPoints: q q q q p p q q, Ring-1}CurrentPoints: q p p q. Trapezoids are constructed by using old and new intersections as described in Section 2.3, producing six trapezoids: q q q p , p q q p , p p p p , p q q q, p q p p , p p q p , of which the "rst four belong to the loop and the last two belong to ring R . The third and the sixth trapezoids have degenerated into triangles. Note that the trapezoid p q q p , which lies outside the polygon, has been correctly omitted. The process continues in this way until all scan lines have been processed. If it is decided that the interior of any ring, R , is not to be trapezoided, the arrays RingH j}OldPoints and Ring-j}CurrentPoints can simply be discarded. 3.2. Touching and coincident rings When the boundaries of the holes are touching or coincident, the algorithm described in the previous sec-
Fig. 5. Partially coincident rings.
tion has to be adjusted to ensure that inconsistencies are avoided. It should be noted that just because parts of two rings are coincident, they do not necessarily share polygon vertices, as illustrated in Fig. 5, where p is a vertex of ring R , but not of R , for which it will be an intersection point with the scan line. Likewise, a vertex which lies on both rings may have di!erent types on the two rings, for example, in Fig. 5, p is of type INT in ring R but is of type MIN in ring R (see Fig. 2). Thus, to avoid confusion, a vertex which is common to two (or more) rings should be numbered di!erently in each ring, that is, it is regarded as a set of coincident vertices, rather than as a multiple vertex. The algorithm proceeds as previously, and a notional trapezoid is generated `betweena two coincident rings. At
358
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
Fig. 6. A ring partially coincident with the loop.
the "nal stage, these null trapezoids are identi"ed by checking the pairs of points on the previous and current scan lines } if the points on both lines are coincident, the trapezoid is discarded. To illustrate this, we amend the polygon illustrated in Fig. 4 so that the ring R becomes partially coincident with the loop, as shown in Fig. 6. In adherence to the convention indicated above, there is no di!erence in the numbering of the vertices and intersection points in the two cases. As previously, at scan line s the algorithm produces the six trapezoids: q q q p , p q q p , p p p p , p q q q, p q p p , p p q p . However, we note that p "p and q"q, so the second trapezoid is null and can be discarded. A simple test of the x-coordinates of the trapezoid vertices on the previous and the current scan lines is su$cient to identify such cases. Particular care has to be taken in this situation when trying to identify the parent of a ring. Numerical inaccuracies may mean that the intersection of a nested ring with a scan line is calculated to be exterior to its coincident parent. To avoid this, once the intersections with each scan line are sorted in x-increasing order, they are scanned to "nd inconsistencies, looking particularly at coincident points. The inconsistencies are removed by suitable interchanges of the relevant points. If a ring is intersected by the "rst scan line because it coincides with the loop there (explicitly excluded from consideration in Section 3.1), the initialisation step of Section 2.2 must include the creation of an additional array Ring-j}OldPoints for each intersected ring, using procedure (P1), as before.
The minor changes described above do not change the basic structure of the algorithm. However, their inclusion considerably broadens its range of applicability, and removes the necessity for manual pre-processing of the data to avoid the possibility of rings touching each other or the loop. Thus, provided that: E vertices belonging to rings are numbered separately on each ring and E intersections with scan lines are ordered correctly, the algorithm can deal with any polygon containing any number of holes to any depth of nesting, where it is assumed that polygon edges can touch or overlap, but not intersect.
4. Reducing the number of trapezoids In general, the algorithm presented generates a greater number of trapezoids than necessary. This is because, in common with all other scan-line algorithms, it works on a line-by-line basis and does not look ahead to identify simpli"cations. Applications which subsequently use the trapezoidation created (for example, to "nd the intersection of two planar polygons) will work more e$ciently if the number of trapezoids is minimised. This can be achieved by the use of a simple post-processing algorithm. The costs of this post-processing are considerably outweighed by the economies which result from the simpli"ed data structures and processing methods associated with scan-line approach. The scale of the reduction tends to increase with the complexity of the polygon; a 70%
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
reduction in the number of trapezoids is not uncommon for complex polygons. The manner in which trapezoids are created automatically generates a sequence which is sorted "rstly by y-coordinate and secondly by x-coordinate, thus making the identi"cation of the trapezoids which are suitable for amalgamation much more straightforward. Trapezoids at successive scan lines (such as ¹ and H ¹ in Fig. 7) can be joined only if all of the following I conditions are true: E E E E E
¹ and ¹ both belong to the same ring or to the loop, H I the x-coordinates of points p and p are equal,
I the x-coordinates of points p and p are equal, H I vertices p , p , p are collinear, H H I vertices p , p , p are collinear. H H I The tests above are written in order of their ease of application; they can be applied rapidly and often reduce the trapezoid count substantially. It should be noted that amalgamation takes place only between two trapezoids at any time, although the combined trapezoid may, itself, be further combined with another trapezoid at the following scan line. Fig. 8 illustrates the simpli"cations which are produced in the polygon of Fig. 4.
359
5. Hints for implementation 5.1. Which polygon edges intersect a scan line? The critical part of the algorithm described in Section 3 is the determination of which polygon edges intersect the current scan line, and this can be performed economically by using appropriate data structures } see Narkhede and Manocha [5]. In our implementation, associated with each vertex in the array SortedVerticesY are two pointers, i and i , to the neighbouring vertices and a #ag Used which is initially FALSE and then set to TRUE when the sweep reaches the scan line on which the vertex sits. At each scan line, we wish to identify the intersected polygon edges. After the relevant vertex #ags have been set, the array SortedVerticesY is searched, starting from the "rst element and "nishing just before the "rst vertex which resides on the current scan line is encountered. For a vertex p within this set, if the #ag Used of a neighbouring vertex p is FALSE, the edge joining p to p intersects G G the current scan line. If the Used #ags of a vertex and both of its neighbours are TRUE, then the vertex can be removed from the search. If the neighbours are initially ordered so that i points to the neighbour with the larger y coordinate, only the vertex and the i -neighbour have to be tested. 5.2. Quick inclusion test
Fig. 7. A con"guration where trapezoids can be joined.
Given an intersection p of the current scan line, where p belongs to ring R, we are seeking to identify the polygon which immediately encloses R. This can be the loop or another ring, and its identity can be determined by the use of a few simple tests (Fig. 9). The "rst and last points on the scan line must clearly belong to the loop, so p cannot be at these positions. Also, if p is the second or
Fig. 8. Simpli"cation produced by post-processing.
360
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
Fig. 9. What encloses point p?
penultimate point on the scan line, there is no need to test as R must be immediately within the loop. For p in other positions, we move along the scan line from p towards the nearer end of the line (in terms of the number of points), searching for the "rst point q on the scan line which does not belong to ring R. If q belongs to the loop, then R is immediately within the loop; if q belongs to ring Q, we count the number of points of Q which lie on the scan line in the direction of movement: E if this number is odd, then p lies immediately inside Q and we can exit, E if the number is even, then searching continues until a ring having an odd number of intersection points is identi"ed (in which case R is within that ring) or until a loop point is encountered (in which case R is within the loop). It should be noted that, because we do not allow intersecting holes, each ring in the hierarchy has a single parent (though, of course, it may have several children). Thus, once the parent of a ring is established at a scan line, that information can be utilised at all subsequent scan lines, thereby reducing the need to search along these lines. To support this, introduce an integer array RingParent in which the elements relate to the rings. The elements of the array are initialised to !1 and are subsequently set to the ring number of the parent, when identi"ed (the loop is regarded as ring 0). Thus, for p belonging to ring R, test the sign of the corresponding element in RingParent: E if the element is negative, search the scan line as above, E otherwise, P is immediately inside the ring/loop indicated by the array element.
6. Evaluation of the time complexity of the algorithm We now present an analysis of the time complexity of the algorithm. We consider the worst case throughout; an analysis of an average case is di$cult because, while computation time obviously depends upon the physical form of the polygon, it can also be highly dependent upon the orientation of the polygon with respect to the sweep direction. For example, consider the simple E-shape in Fig. 10. In (a), there are 6 scan lines and 5 trapezoids, but the rotation to case (b) increases these numbers to 12 scan lines and 23 trapezoids. A further small rotation to (c) reduces the numbers to 3 scan lines and 4 trapezoids. Clearly, the greater the number of scan lines and trapezoids, the longer will be the computation time involved. The three parts of the algorithm } initialisation, trapezoidation and post-processing } will be treated separately. We assume that there are n vertices in total and that there are m di!erent scan lines passing through them (m)n). 1. Initialisation: as described in Sections 2.1 and 2.2, this can be decomposed into three steps: 1. 䡩 determining the vertex neighbours requires one pass only and is completed in linear time, 1. 䡩 sorting the polygon vertices regarding their y coordinates to create SortedVerticesY; using Quicksort, it can be performed in ¹ (n)"O(n log n), 1. 䡩 initialising the array OldPoints; this depends upon k, the number of polygon vertices which have the maximal y-coordinate } sorting these vertices into x-increasing order and inserting them into the array will take O(k log k). In the worst case (a triangle in
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
361
Fig. 10. Trapezoidation in di!erent orientations.
2.
1.
1.
1.
which all vertices except one lie on the horizontal upper side), there are n!1 such vertices, but in general, k is small and almost independent of n. The time complexity of initialisation is therefore ¹ (n)"O(n log n). Processinu at the scan lines: the main loop of the algorithm contains m!1 steps, each of which involves the following tasks: 䡩 determining k, the number of polygon vertices lying on the scan line; because the input points are sorted this requires k comparisons and the process is thus linear in k (k(n), 䡩 calculating the intersections between the polygon sides and the scan line; as explained in Section 5.1, this involves searching a subset of SortedVerticesY, testing #ags and calculating the j intersection points; these operations are all linear in time, 䡩 sorting the s"j#k points on the scan line into x-increasing order } these will probably be partly sorted already, though scan-line coherence is less useful here than in, for example, a polygon-"ll algorithm; at worst, a Quicksort would require ¹
(s)"O(s log s), s(n.
1. 䡩 if there are rings, identifying the parent ring } for each ring this involves a single search of half of the points on a single scan line (Section 5.2); for n large this is negligible, 1. 䡩 creating trapezoids between the current and previous scan lines } if there are no rings, the number of trapezoids is, at most, s!1; if rings are present, the ring intersections are included in arrays twice, so in the worst case the number of trapezoids is 2s!1; the time for the process will be linear in s (s(n). Thus, the most critical part is the x-sorting of intersections, and taking into account all scan lines, ¹ "(m!1)¹ (s), m, s(n. Given that it is easy
to construct pathological cases in which m&n and s&n, this gives ¹ (n)"O(n log n), though, in
general, it is unlikely that the time dependency will reach such an order. 3. Post-processinu: as explained in Section 4, the trapezoids are presented to the post-processing stage already sorted, "rstly with respect to y-coordinates, and secondly with respect to x-coordinates. If there are p trapezoids above the ith scan line, and q trapezoids below it, we test each set against the other to see whether trapezoids can be joined. At each test, one of the trapezoids is removed from consideration if they cannot be joined; if joining is possible, both are removed from consideration. Thus, the number of tests at each scan line is, at most, p#q. The tests take place on m!2 scan lines. As p#q(2n and m(n, ¹ (n)"O(n). Thus, consid ering the algorithm as a whole, (n)#¹ (n)#¹ (n)"O(n log n)
#O(nlog n)#O(n)"O(n log n). However, in everyday cases, it is likely that the algorithm will perform rather better than this expression suggests.
¹(n)"¹
7. Tests A number of tests were run on a Pentium 133 MHz, using Visual C##5.0, to determine the performance of the algorithm, both with arti"cially constructed data and data from the geographical application. Tables 1 and 2 show data for two polygons without holes; the "rst is almost monotone and the second has many concavities. Times are presented for a single run, although because of the short duration, the times were actually measured over 100 runs. The times for the trapezoidation phase are, by far, the most signi"cant, and it can be seen that, as predicted in Section 6, they do not increase as rapidly as a behaviour O(n log n) ; in fact, the increase is only marginally greater than O(n log n).
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
362
Table 1 Simple polygon with almost vertically monotone chain of vertices Number of vertices Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation CPU time for initialisation (ms) CPU time for trapezoidation (ms) CPU time for optimisation (ms) Total CPU time (ms)
50 49 4 45
100 92 3 89
150 124 5 119
200 186 8 178
0.2 3.3 0.7 4.2
0.6 9.1 0.6 10.3
1.0 15.4 0.7 17.1
1.3 28.4 1.0 30.7
Table 2 Polygon with many concave parts Number of vertices Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation CPU time for initialisation (ms) CPU time for trapezoidation (ms) CPU time for optimisation (ms) Total CPU time (ms)
50 147 43 104
100 297 92 205
150 417 130 287
200 581 193 388
0.3 5.3 1.0 6.6
0.4 12.8 2.1 15.3
0.7 20.2 3.0 23.9
1.0 36.8 4.7 42.5
Table 3 Polygon with many concave parts and holes Number of holes Number of vertices Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation Total CPU time (ms)
Table 4 Examples from geographical data Total number of vertices Total number of loops Number of loops with at least one ring Loops with k vertices, k(5 Loops with k vertices, 5(k(10 Loops with k vertices, 10(k(20 Loops with k vertices, 20(k(100 Loops with k vertices, 100(k(200 Loops with k vertices, 200(k Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation Time needed for trapezoidation (s)
37,310 2556 57 534 671 915 413 18 5 42,130 16,487 25,643
11,876 1086 32 159 529 282 98 18 0 12,410 4586 7824
3.78
1.03
Tables 3 and 4 relate to geographical data, the "rst is urban, the second semi-urban. The data consists of many loops, which have varying numbers of vertices; a number of the loops enclose rings. The method is seen to deal with large numbers of complex polygons e!ectively and economically. In both cases, the optimisation phase reduced the number of trapezoids by approximately 38%. 8. Conclusions
0 100 326 234
1 110 396 286
2 120 462 336
5 150 677 502
10 200 966 704
92
110
126
175
262
16.2
18.6
21.1
32.9
45.9
The greater complexity of the second polygon results in many more trapezoids being created, and thus the postprocessing becomes more signi"cant; the number of trapezoids joined in Table 2 increases almost linearly, as does the associated time, which is again much less than the worst case considered in Section 6. Further tests were applied on a complex polygon with many concavities and a varying number of holes. The loop was the same in each case and contained 100 vertices; each ring introduced contained 10 vertices. The fragmentation of the interior of the loop caused by presence of rings is clearly re#ected in the scale of the optimisation which took place. It is not appropriate to compare times for the di!erent phases of the algorithm as they relate to topologically di!erent cases; however, the overall processing time is included, for reference.
The paper describes an algorithm which decomposes planar polygons into trapezoids. The polygons can be nonmonotone and may contain a number of holes, including nested holes. The holes can be trapezoided, if desired. The algorithm works in three steps: the initial data is prepared, the trapezoids are generated, and the trapezoids produced are optimised. Evaluation of the time complexity shows that the algorithm completes the task in, at worst, O(n log n) time. The classi"cation of the local role of each polygon vertex and the use of the inclusion test proposed in Section 5 ensure that the algorithm is not sensitive to rounding errors, so the implementation is very stable. Seidel [4] shows that, if desired, the triangulation of a polygon can be easily generated from trapezoids in linear time. The algorithm was implemented in C## and was successfully employed in the geographical application alluded to in Section 1. It can be found at: http://www.uni-mb.si/&cgm. References [1] Etzion M, Rappoport A. On compatible start decompositions of simple polygons. IEEE Transactions on Visualization and Computer Graphics 1997;3:87}95.
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363 [2] Chazelle B, Incerpi J. Triangulation and shape complexity. ACM Transactions on Graphics 1984;3:135}52. [3] Fournier A, Montuno DY. Triangulating simple polygons and equivalent problems. ACM Transactions on Graphics 1984;3:153}74. [4] Seidel R. A simple and fast incremental randomized algorithm for computing trapezoidal decompositions and for triangulating polygons. Computational Geometry: Theory and Applications 1991;1:51}64. [5] Narkhede A, Manocha D. Fast polygon triangulation based on Seidel's algorithm. In: Paeth AW,
[6] [7] [8] [9]
363
editor. Graphics gems V. Boston: Academic Press, 1995. p. 394}7. Preparata FP, Shamos MI. Computational Geometry: An Introduction. Berlin: Springer, 1985. Manber U. Introduction to algorithms: a creative approach. Reading, MA: Addison-Wesley, 1989. Mortenson ME. Geometric modeling. New York: Wiley, 1985. O'Rourke J. Computational Geometry in C. Cambridge: Cambridge University Press, 1993.
Technical Section
A universal trapezoidation algorithm for planar polygons Borut Z[ alik *, Gordon J. Clapworthy University of Maribor, Faculty of Electrical Engineering and Computer Science, Smetanova 17, 2000 Maribor, Slovenia De Montfort University, Department of Computer and Information Sciences, Hammerwood Gate, Kerts Hill, Milton Keynes, MK7 6HP, UK
Abstract A new algorithm for decomposing non-monotone planar polygons, which may contain holes, into trapezoids is described. The holes may be nested and may have common edges. In the "rst part of the paper, the main idea is explained for non-monotone polygons without holes, and the algorithm is then extended to polygons containing holes; the holes can also be decomposed into trapezoids, if desired. Finally, it is shown that the algorithm performs trapezoidation of general polygons in O(n log n) time, where n is the common number of polygon vertices. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Polygon decomposition; Trapezoidation; Computational geometry
1. Introduction Techniques for decomposing complex geometric shapes into simpler components are of great importance in computational geometry and its applications. Depending upon the particular problem, there are di!erent ideas about how this should be performed and the choice of element into which the shape should be decomposed. For example, planar polygons can be decomposed into triangles, trapezoids, or even star-shaped polygons [1]. While triangulation is, without doubt, the most popular and the most investigated technique, splitting a complex polygon into trapezoids is an alternative, and a generalised method for performing trapezoidation on an arbitrary polygon with any number of holes, nested to any level, is considered in the paper. The motivation for the current work arose from the display of land use within a given geographical area. The available software could colour triangles or trapezoids to display di!erent uses but was more e$cient with trap-
* Corresponding author. Tel.: #368-62-220-7471; fax: #386-62-211-178. E-mail address: [email protected] (B. Z[ alik)
ezoids because of the reduced number of primitives involved. The individual areas to be displayed could be nested, and their boundaries could coincide with each other or the boundary of the overall area. However, there is no reason why the application should not be used more broadly, particularly to display planar variation of physical quantities. It can also be applied to "nd, e$ciently, the intersection and/or union of polygons. Trapezoidation was "rst considered by Chazelle and Incerpi [2] and Fournier and Montuno [3]. A randomised algorithm for trapezoidation was later suggested by Seidel [4] as the key for triangulation and this was recently used by Narkhede and Manocha [5] for the same purpose. However, none of these authors considered general planar polygons containing holes; indeed, to the authors' knowledge, no previously published method for trapezoidation has successfully addressed polygons with holes, nor the trapezoidation within holes. A trapezoid is a four-sided polygon in which two of the edges are parallel. Note that a triangle can be considered as a degenerate trapezoid in which one of the parallel edges has zero length; hence, some triangles might result from applying any trapezoidation algorithm. A horizontal trapezoidation of an n-sided polygon involves drawing horizontal lines (we shall call them scan lines) through every vertex of the polygon (see Fig. 1). The
0097-8493/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 9 3 ( 9 9 ) 0 0 0 4 4 - 8
354
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
is swept vertically down the plane and an initial data structure is updated at discrete `eventsa during this process. In this algorithm, an event is the intersection of the sweep line with one or more polygon vertices. By maintaining the vertices in a dynamic ordered list, only the leading element has to be checked for intersection at any time. The precise action taken when a new vertex is intersected depends upon the vertex type } these types are identi"ed in Section 2.1. The "rst step in the algorithm considers all intersections of the current scan line with the polygon vertices and edges, sorting them into x-increasing order. Those which are relevant to trapezoidation between the current scan line and the previous scan line are retained, the others are discarded. The second step identi"es the points on the current and previous scan lines which are to be joined to create the trapezoids. Fig. 1. Horizontal trapezoidation of a monotone polygon.
whole operation can be completed in O(n log n) time if the polygon is monotone [4] } in the context of horizontal trapezoidation, a monotone polygon is one in which the boundary consists of two vertically monotone chains. In practice, polygons are frequently non-monotone and may enclose a number of holes. This paper describes a universal algorithm to perform the trapezoidation of a non-monotone planar polygon containing holes. In the "rst part of the paper, the main idea of the algorithm is explained for non-monotone polygons without holes. In the following section, the algorithm is extended to polygons containing holes; the holes can themselves be trapezoided by the algorithm, if required. A straightforward application of the algorithm results in more trapezoids than strictly necessary, and Section 4 indicates how this number can be reduced. Section 5 provides suggestions for an e$cient implementation, Section 6 contains an estimation of the time complexity of the algorithm, and Section 7 describes tests made with both arti"cially constructed polygons and real data from the geographical application.
2. Polygons without holes The algorithm uses the well-known greedy approach [6,7]. In the "rst step, the algorithm prepares the input data such that in the second step the problem can be solved in one pass. Adopting vocabulary from topology [8], the sequence of vertices forming the border of a polygon is named a loop, and a sequence of vertices forming the border of an inner closed polygonal region is called a ring. Initially, the polygon vertices are sorted by y-coordinate. A `plane sweepa is employed, which is a technique used in many geometrical applications } a horizontal line
2.1. Vertex classixcation The solution presented is based on a simple analysis of each polygon vertex, a brief description of which is given by O'Rourke [9], without implementation details. However, following Seidel [4], O'Rourke explicitly ignores the possibility that a scan line may intersect more than one vertex; if this should occur, one of the vertices is notionally moved to an additional scan line. For clarity, our algorithm includes the possibility of multiple vertices on a scan line, so the range of possible vertex types is larger. The following types of relationship can exist between a polygon vertex p and its neighbouring vertices p J and p : P E Intersection: when one of neighbouring vertices has a greater and the other a lower y-coordinate than p (Fig. 2a) p is handled as if it were an intersection of the scan line with a polygon edge, so it is categorised as INT. E Local minimum: when both neighbouring vertices have greater y-coordinates than p (see Fig. 2b), p is categorised as MIN. E Local maximum: when both neighbouring vertices have lower y-coordinates than p (see Fig. 2c), p is categorised as MAX. E Horizontal line seument: when one of the neighbours has the same y-coordinate as p within tolerance e, then p is an end point of a horizontal line segment; by considering the position of the second neighbour, vertex p is given the code E 䡩 HMAX if the second neighbour has a lower y-coordinate (Fig. 2d) or E 䡩 HMIN if the second neighbour has a higher y-coordinate (Fig. 2e). E If both neighbours have the same y coordinate as p within tolerance e (Fig. 2f ), then p will not form part of any horizontal trapezoid generated later and it can be erased from the list of polygon points at this stage.
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
355
Fig. 2. The types of local behaviour of vertex p.
Fig. 3. A non-monotone polygon.
2.2. Initialisation
E if p is of type MAX, it is inserted twice into G array OldPoints
Consider the polygon in Fig. 3 in which the vertices are denoted by p and indicated by circular markers. The G polygon is non-monotone and contains vertices belonging to all of the classi"cations illustrated in Fig. 2, apart from the uninteresting case (f ) in which the vertex would be trivially discarded at the outset. The vertex coordinates are stored sequentially in an array named InputPoints. The vertices of the polygon are sorted with respect to their y-coordinates and stored in an array called SortedVerticesY. It is possible, of course, to have more than one vertex with the same y coordinate. The vertices lying on the "rst scan line (line s in Fig. 3) are then considered; clearly, they cannot be of types INT, MIN or HMIN. The vertices are sorted into x-increasing order and stored in an array named OldPoints using the following procedure:
E if p is of type HMAX, it is inserted once into G OldPoints,
(P1)
E vertices within a horizontal line segment are ignored. For the polygon in Fig. 3, OldPoints now contains p p p p . 2.3. Processing a scan line In the main part of the algorithm, the scan lines s (k'0) are considered sequentially, in descending orI der of y-coordinate. If there are n polygon vertices, there will be, at most, n scan lines and, hence, n!1 steps in this part of the algorithm.
356
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
In Figs. 3, 4 and 6, intersections between scan lines and polygon edges, qH, are identi"ed by square markers; here, G the subscript denotes the ith scan line and the superscript indicates the jth intersection on the line. To create an ordered list of intersections of the current scan line with the polygon, the following actions are needed at s : I E determine the number and the types of polygon vertices lying on s , I E calculate all intersections between s and polygon I edges } these points are of type INT; for s in Fig. 3, we obtain two intersection points q and q, E insert the vertices and the intersections into an array named Intersections in order of increasing x-coordinate; this array will be used twice, E create an array CurrentPoints from Intersections using the procedure: E MIN: add point twice to CurrentPoints, E INT, HMIN: add point once to Current Points,
(P2)
E MAX, HMAX: do not add point. In Fig. 3, s passes through only one vertex } p which is of type MIN. After sorting, the array CurrentPoints contains the points q p p q. Because of the way in which the algorithm is constructed, the arrays OldPoints and CurrentPoints will now contain equal (even) numbers of points, and trapezoids can be generated by joining corresponding pairs of points from the two arrays. Between s and s , two trapezoids are created, determined by the points: p qp p and p p qp . It can be seen that the "rst trapezoid has degenerated into a triangle. 2.4. Updating for the next scan line Before moving to the next scan line, we update OldPoints from the points stored in Intersections according to their types, as follows: E MAX: add point twice to OldPoints, E INT, HMAX: add point once to OldPoints,
(P3)
E MIN, HMIN: do not add point. For s , Intersections contained q p q, which have types INT, MIN, INT, respectively, so OldPoints becomes q q. Scan lines continue to be processed in this way until the polygon is exhausted. The basic algorithm described above, involving processes familiar from previous decomposition algorithms and from scan-line "ll algorithms, is generalised in Section 3 to cope with complex geometries in which the
polygon contains several rings, including the possibility of multiply nested rings and partially coincident rings.
3. Polygons with holes This section considers the generalisation of the trapezoidation algorithm for all planar polygons which do not self-intersect. For simplicity, we deal initially with a situation in which the loop and all rings are totally disjoint and all areas within the loop are to be trapezoided, including the interiors of the rings. This is later generalised to allow parts of the rings to coincide with each other and with the loop, and to allow di!erent areas to be omitted from the trapezoidation. 3.1. Simple nested rings The initialisation takes place as previously, with the array SortedVerticesY being generated and the array OldPoints being initialised as before (because the point with the highest y-coordinate of the polygon cannot belong to any of the rings). The changes to the algorithm are introduced in the subsequent stages, described previously in Sections 2.3 and 2.4; these will be illustrated using the polygon in Fig. 4. It can be seen that the polygon of Fig. 4 contains another polygon, de"ning the ring R . The algorithm requires that two new arrays, Ring-j}OldPoints and Ringj}CurrentPoints, are introduced for each ring, R . H Prior to reaching scan line s in Fig. 4, the algorithm behaves as described previously, and the array OldPoints contains q q. At the next scan line, s , the "rst intersection with a ring is encountered. Array CurrentPoints is generated using (P2) from Section 2.3, and it contains, as previously, q q, so the trapezoid q q q q is again produced. Note that all ring vertices on s are of types MAX or HMAX and are thus ignored within procedure (P2). The "rst change comes when updating before processing the next scan line, s . Intersections with the loop are again considered for the updating of OldPoints using (P3) from Section 2.4, but intersections of the scan line with a ring R are now considered for two arrays: H E the array Ring-j}OldPoints belonging to R and H E the array OldPoints for the loop if R is immediately H within the loop, or the array Ring-k}OldPoints if R is H nested immediately within another ring, R . I Similarly, at subsequent scan lines, in addition to the array CurrentPoints, for each ring R an array RingH j}CurrentPoints is created using (P2), so that ring intersections with the current scan line also occur in two arrays.
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
357
Fig. 4. A non-monotone polygon containing a hole.
Whether the current ring is immediately within the loop or another ring can be determined without additional calculations simply by observing the neighbouring points of intersection with the scan line. This is discussed further in Section 5.2. In the example, in moving from s to s , the arrays are as follows: OldPoints: q p p p p p p q, Ring-1}OldPoints: p p p p , CurrentPoints: q q q q p p q q, Ring-1}CurrentPoints: q p p q. Trapezoids are constructed by using old and new intersections as described in Section 2.3, producing six trapezoids: q q q p , p q q p , p p p p , p q q q, p q p p , p p q p , of which the "rst four belong to the loop and the last two belong to ring R . The third and the sixth trapezoids have degenerated into triangles. Note that the trapezoid p q q p , which lies outside the polygon, has been correctly omitted. The process continues in this way until all scan lines have been processed. If it is decided that the interior of any ring, R , is not to be trapezoided, the arrays RingH j}OldPoints and Ring-j}CurrentPoints can simply be discarded. 3.2. Touching and coincident rings When the boundaries of the holes are touching or coincident, the algorithm described in the previous sec-
Fig. 5. Partially coincident rings.
tion has to be adjusted to ensure that inconsistencies are avoided. It should be noted that just because parts of two rings are coincident, they do not necessarily share polygon vertices, as illustrated in Fig. 5, where p is a vertex of ring R , but not of R , for which it will be an intersection point with the scan line. Likewise, a vertex which lies on both rings may have di!erent types on the two rings, for example, in Fig. 5, p is of type INT in ring R but is of type MIN in ring R (see Fig. 2). Thus, to avoid confusion, a vertex which is common to two (or more) rings should be numbered di!erently in each ring, that is, it is regarded as a set of coincident vertices, rather than as a multiple vertex. The algorithm proceeds as previously, and a notional trapezoid is generated `betweena two coincident rings. At
358
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
Fig. 6. A ring partially coincident with the loop.
the "nal stage, these null trapezoids are identi"ed by checking the pairs of points on the previous and current scan lines } if the points on both lines are coincident, the trapezoid is discarded. To illustrate this, we amend the polygon illustrated in Fig. 4 so that the ring R becomes partially coincident with the loop, as shown in Fig. 6. In adherence to the convention indicated above, there is no di!erence in the numbering of the vertices and intersection points in the two cases. As previously, at scan line s the algorithm produces the six trapezoids: q q q p , p q q p , p p p p , p q q q, p q p p , p p q p . However, we note that p "p and q"q, so the second trapezoid is null and can be discarded. A simple test of the x-coordinates of the trapezoid vertices on the previous and the current scan lines is su$cient to identify such cases. Particular care has to be taken in this situation when trying to identify the parent of a ring. Numerical inaccuracies may mean that the intersection of a nested ring with a scan line is calculated to be exterior to its coincident parent. To avoid this, once the intersections with each scan line are sorted in x-increasing order, they are scanned to "nd inconsistencies, looking particularly at coincident points. The inconsistencies are removed by suitable interchanges of the relevant points. If a ring is intersected by the "rst scan line because it coincides with the loop there (explicitly excluded from consideration in Section 3.1), the initialisation step of Section 2.2 must include the creation of an additional array Ring-j}OldPoints for each intersected ring, using procedure (P1), as before.
The minor changes described above do not change the basic structure of the algorithm. However, their inclusion considerably broadens its range of applicability, and removes the necessity for manual pre-processing of the data to avoid the possibility of rings touching each other or the loop. Thus, provided that: E vertices belonging to rings are numbered separately on each ring and E intersections with scan lines are ordered correctly, the algorithm can deal with any polygon containing any number of holes to any depth of nesting, where it is assumed that polygon edges can touch or overlap, but not intersect.
4. Reducing the number of trapezoids In general, the algorithm presented generates a greater number of trapezoids than necessary. This is because, in common with all other scan-line algorithms, it works on a line-by-line basis and does not look ahead to identify simpli"cations. Applications which subsequently use the trapezoidation created (for example, to "nd the intersection of two planar polygons) will work more e$ciently if the number of trapezoids is minimised. This can be achieved by the use of a simple post-processing algorithm. The costs of this post-processing are considerably outweighed by the economies which result from the simpli"ed data structures and processing methods associated with scan-line approach. The scale of the reduction tends to increase with the complexity of the polygon; a 70%
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
reduction in the number of trapezoids is not uncommon for complex polygons. The manner in which trapezoids are created automatically generates a sequence which is sorted "rstly by y-coordinate and secondly by x-coordinate, thus making the identi"cation of the trapezoids which are suitable for amalgamation much more straightforward. Trapezoids at successive scan lines (such as ¹ and H ¹ in Fig. 7) can be joined only if all of the following I conditions are true: E E E E E
¹ and ¹ both belong to the same ring or to the loop, H I the x-coordinates of points p and p are equal,
I the x-coordinates of points p and p are equal, H I vertices p , p , p are collinear, H H I vertices p , p , p are collinear. H H I The tests above are written in order of their ease of application; they can be applied rapidly and often reduce the trapezoid count substantially. It should be noted that amalgamation takes place only between two trapezoids at any time, although the combined trapezoid may, itself, be further combined with another trapezoid at the following scan line. Fig. 8 illustrates the simpli"cations which are produced in the polygon of Fig. 4.
359
5. Hints for implementation 5.1. Which polygon edges intersect a scan line? The critical part of the algorithm described in Section 3 is the determination of which polygon edges intersect the current scan line, and this can be performed economically by using appropriate data structures } see Narkhede and Manocha [5]. In our implementation, associated with each vertex in the array SortedVerticesY are two pointers, i and i , to the neighbouring vertices and a #ag Used which is initially FALSE and then set to TRUE when the sweep reaches the scan line on which the vertex sits. At each scan line, we wish to identify the intersected polygon edges. After the relevant vertex #ags have been set, the array SortedVerticesY is searched, starting from the "rst element and "nishing just before the "rst vertex which resides on the current scan line is encountered. For a vertex p within this set, if the #ag Used of a neighbouring vertex p is FALSE, the edge joining p to p intersects G G the current scan line. If the Used #ags of a vertex and both of its neighbours are TRUE, then the vertex can be removed from the search. If the neighbours are initially ordered so that i points to the neighbour with the larger y coordinate, only the vertex and the i -neighbour have to be tested. 5.2. Quick inclusion test
Fig. 7. A con"guration where trapezoids can be joined.
Given an intersection p of the current scan line, where p belongs to ring R, we are seeking to identify the polygon which immediately encloses R. This can be the loop or another ring, and its identity can be determined by the use of a few simple tests (Fig. 9). The "rst and last points on the scan line must clearly belong to the loop, so p cannot be at these positions. Also, if p is the second or
Fig. 8. Simpli"cation produced by post-processing.
360
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
Fig. 9. What encloses point p?
penultimate point on the scan line, there is no need to test as R must be immediately within the loop. For p in other positions, we move along the scan line from p towards the nearer end of the line (in terms of the number of points), searching for the "rst point q on the scan line which does not belong to ring R. If q belongs to the loop, then R is immediately within the loop; if q belongs to ring Q, we count the number of points of Q which lie on the scan line in the direction of movement: E if this number is odd, then p lies immediately inside Q and we can exit, E if the number is even, then searching continues until a ring having an odd number of intersection points is identi"ed (in which case R is within that ring) or until a loop point is encountered (in which case R is within the loop). It should be noted that, because we do not allow intersecting holes, each ring in the hierarchy has a single parent (though, of course, it may have several children). Thus, once the parent of a ring is established at a scan line, that information can be utilised at all subsequent scan lines, thereby reducing the need to search along these lines. To support this, introduce an integer array RingParent in which the elements relate to the rings. The elements of the array are initialised to !1 and are subsequently set to the ring number of the parent, when identi"ed (the loop is regarded as ring 0). Thus, for p belonging to ring R, test the sign of the corresponding element in RingParent: E if the element is negative, search the scan line as above, E otherwise, P is immediately inside the ring/loop indicated by the array element.
6. Evaluation of the time complexity of the algorithm We now present an analysis of the time complexity of the algorithm. We consider the worst case throughout; an analysis of an average case is di$cult because, while computation time obviously depends upon the physical form of the polygon, it can also be highly dependent upon the orientation of the polygon with respect to the sweep direction. For example, consider the simple E-shape in Fig. 10. In (a), there are 6 scan lines and 5 trapezoids, but the rotation to case (b) increases these numbers to 12 scan lines and 23 trapezoids. A further small rotation to (c) reduces the numbers to 3 scan lines and 4 trapezoids. Clearly, the greater the number of scan lines and trapezoids, the longer will be the computation time involved. The three parts of the algorithm } initialisation, trapezoidation and post-processing } will be treated separately. We assume that there are n vertices in total and that there are m di!erent scan lines passing through them (m)n). 1. Initialisation: as described in Sections 2.1 and 2.2, this can be decomposed into three steps: 1. 䡩 determining the vertex neighbours requires one pass only and is completed in linear time, 1. 䡩 sorting the polygon vertices regarding their y coordinates to create SortedVerticesY; using Quicksort, it can be performed in ¹ (n)"O(n log n), 1. 䡩 initialising the array OldPoints; this depends upon k, the number of polygon vertices which have the maximal y-coordinate } sorting these vertices into x-increasing order and inserting them into the array will take O(k log k). In the worst case (a triangle in
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
361
Fig. 10. Trapezoidation in di!erent orientations.
2.
1.
1.
1.
which all vertices except one lie on the horizontal upper side), there are n!1 such vertices, but in general, k is small and almost independent of n. The time complexity of initialisation is therefore ¹ (n)"O(n log n). Processinu at the scan lines: the main loop of the algorithm contains m!1 steps, each of which involves the following tasks: 䡩 determining k, the number of polygon vertices lying on the scan line; because the input points are sorted this requires k comparisons and the process is thus linear in k (k(n), 䡩 calculating the intersections between the polygon sides and the scan line; as explained in Section 5.1, this involves searching a subset of SortedVerticesY, testing #ags and calculating the j intersection points; these operations are all linear in time, 䡩 sorting the s"j#k points on the scan line into x-increasing order } these will probably be partly sorted already, though scan-line coherence is less useful here than in, for example, a polygon-"ll algorithm; at worst, a Quicksort would require ¹
(s)"O(s log s), s(n.
1. 䡩 if there are rings, identifying the parent ring } for each ring this involves a single search of half of the points on a single scan line (Section 5.2); for n large this is negligible, 1. 䡩 creating trapezoids between the current and previous scan lines } if there are no rings, the number of trapezoids is, at most, s!1; if rings are present, the ring intersections are included in arrays twice, so in the worst case the number of trapezoids is 2s!1; the time for the process will be linear in s (s(n). Thus, the most critical part is the x-sorting of intersections, and taking into account all scan lines, ¹ "(m!1)¹ (s), m, s(n. Given that it is easy
to construct pathological cases in which m&n and s&n, this gives ¹ (n)"O(n log n), though, in
general, it is unlikely that the time dependency will reach such an order. 3. Post-processinu: as explained in Section 4, the trapezoids are presented to the post-processing stage already sorted, "rstly with respect to y-coordinates, and secondly with respect to x-coordinates. If there are p trapezoids above the ith scan line, and q trapezoids below it, we test each set against the other to see whether trapezoids can be joined. At each test, one of the trapezoids is removed from consideration if they cannot be joined; if joining is possible, both are removed from consideration. Thus, the number of tests at each scan line is, at most, p#q. The tests take place on m!2 scan lines. As p#q(2n and m(n, ¹ (n)"O(n). Thus, consid ering the algorithm as a whole, (n)#¹ (n)#¹ (n)"O(n log n)
#O(nlog n)#O(n)"O(n log n). However, in everyday cases, it is likely that the algorithm will perform rather better than this expression suggests.
¹(n)"¹
7. Tests A number of tests were run on a Pentium 133 MHz, using Visual C##5.0, to determine the performance of the algorithm, both with arti"cially constructed data and data from the geographical application. Tables 1 and 2 show data for two polygons without holes; the "rst is almost monotone and the second has many concavities. Times are presented for a single run, although because of the short duration, the times were actually measured over 100 runs. The times for the trapezoidation phase are, by far, the most signi"cant, and it can be seen that, as predicted in Section 6, they do not increase as rapidly as a behaviour O(n log n) ; in fact, the increase is only marginally greater than O(n log n).
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363
362
Table 1 Simple polygon with almost vertically monotone chain of vertices Number of vertices Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation CPU time for initialisation (ms) CPU time for trapezoidation (ms) CPU time for optimisation (ms) Total CPU time (ms)
50 49 4 45
100 92 3 89
150 124 5 119
200 186 8 178
0.2 3.3 0.7 4.2
0.6 9.1 0.6 10.3
1.0 15.4 0.7 17.1
1.3 28.4 1.0 30.7
Table 2 Polygon with many concave parts Number of vertices Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation CPU time for initialisation (ms) CPU time for trapezoidation (ms) CPU time for optimisation (ms) Total CPU time (ms)
50 147 43 104
100 297 92 205
150 417 130 287
200 581 193 388
0.3 5.3 1.0 6.6
0.4 12.8 2.1 15.3
0.7 20.2 3.0 23.9
1.0 36.8 4.7 42.5
Table 3 Polygon with many concave parts and holes Number of holes Number of vertices Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation Total CPU time (ms)
Table 4 Examples from geographical data Total number of vertices Total number of loops Number of loops with at least one ring Loops with k vertices, k(5 Loops with k vertices, 5(k(10 Loops with k vertices, 10(k(20 Loops with k vertices, 20(k(100 Loops with k vertices, 100(k(200 Loops with k vertices, 200(k Number of trapezoids Number of trapezoids joined Number of trapezoids after optimisation Time needed for trapezoidation (s)
37,310 2556 57 534 671 915 413 18 5 42,130 16,487 25,643
11,876 1086 32 159 529 282 98 18 0 12,410 4586 7824
3.78
1.03
Tables 3 and 4 relate to geographical data, the "rst is urban, the second semi-urban. The data consists of many loops, which have varying numbers of vertices; a number of the loops enclose rings. The method is seen to deal with large numbers of complex polygons e!ectively and economically. In both cases, the optimisation phase reduced the number of trapezoids by approximately 38%. 8. Conclusions
0 100 326 234
1 110 396 286
2 120 462 336
5 150 677 502
10 200 966 704
92
110
126
175
262
16.2
18.6
21.1
32.9
45.9
The greater complexity of the second polygon results in many more trapezoids being created, and thus the postprocessing becomes more signi"cant; the number of trapezoids joined in Table 2 increases almost linearly, as does the associated time, which is again much less than the worst case considered in Section 6. Further tests were applied on a complex polygon with many concavities and a varying number of holes. The loop was the same in each case and contained 100 vertices; each ring introduced contained 10 vertices. The fragmentation of the interior of the loop caused by presence of rings is clearly re#ected in the scale of the optimisation which took place. It is not appropriate to compare times for the di!erent phases of the algorithm as they relate to topologically di!erent cases; however, the overall processing time is included, for reference.
The paper describes an algorithm which decomposes planar polygons into trapezoids. The polygons can be nonmonotone and may contain a number of holes, including nested holes. The holes can be trapezoided, if desired. The algorithm works in three steps: the initial data is prepared, the trapezoids are generated, and the trapezoids produced are optimised. Evaluation of the time complexity shows that the algorithm completes the task in, at worst, O(n log n) time. The classi"cation of the local role of each polygon vertex and the use of the inclusion test proposed in Section 5 ensure that the algorithm is not sensitive to rounding errors, so the implementation is very stable. Seidel [4] shows that, if desired, the triangulation of a polygon can be easily generated from trapezoids in linear time. The algorithm was implemented in C## and was successfully employed in the geographical application alluded to in Section 1. It can be found at: http://www.uni-mb.si/&cgm. References [1] Etzion M, Rappoport A. On compatible start decompositions of simple polygons. IEEE Transactions on Visualization and Computer Graphics 1997;3:87}95.
B. Z[ alik, G.J. Clapworthy / Computers & Graphics 23 (1999) 353}363 [2] Chazelle B, Incerpi J. Triangulation and shape complexity. ACM Transactions on Graphics 1984;3:135}52. [3] Fournier A, Montuno DY. Triangulating simple polygons and equivalent problems. ACM Transactions on Graphics 1984;3:153}74. [4] Seidel R. A simple and fast incremental randomized algorithm for computing trapezoidal decompositions and for triangulating polygons. Computational Geometry: Theory and Applications 1991;1:51}64. [5] Narkhede A, Manocha D. Fast polygon triangulation based on Seidel's algorithm. In: Paeth AW,
[6] [7] [8] [9]
363
editor. Graphics gems V. Boston: Academic Press, 1995. p. 394}7. Preparata FP, Shamos MI. Computational Geometry: An Introduction. Berlin: Springer, 1985. Manber U. Introduction to algorithms: a creative approach. Reading, MA: Addison-Wesley, 1989. Mortenson ME. Geometric modeling. New York: Wiley, 1985. O'Rourke J. Computational Geometry in C. Cambridge: Cambridge University Press, 1993.