An almost linear-time algorithm for trapezoidation of GIS polygons

Abstract

The decomposition of planar polygons into triangles is a well studied area of computer graphics with particular relevance to geographic information systems (GIS). Trapezoidation is often performed as a first step to triangulation. There is a complex linear-time algorithm [Discrete Comput. Geom. 6 (1991) 485] for the decomposition of a simple polygon into triangles. However, it is extremely complicated and in practice O(n log n) algorithms are used. Our motivation in trapezoidation of large GIS polygons is the fast display of such polygons. It is much faster to display simple shapes like triangles or trapezoids on raster graphics devices, compared to complex polygons. Hence, quite often complex polygons are decomposed into triangles or trapezoids for displaying. Since triangulation is usually more difficult compared to trapezoidation, we are interested in trapezoidation of GIS polygons for faster display. We present a very simple algorithm for the trapezoidation of simple polygons without holes. Our algorithm runs in O(n) time in practice with a very small hidden constant. We have extensively tested our algorithm for polygons in a GIS database. Our algorithm is easy to implement compared to existing algorithms and runs extremely fast even for polygons with thousands of vertices.

Introduction

The decomposition of a simple polygon into triangles is a well studied field of computer graphics [1] and computational geometry. In computer graphics and GIS applications, it is comparatively easier to render a simple shape like a triangle than an arbitrary polygon. Similarly, in computational geometry, a large number of data structures related to polygons are based on triangulation of those polygons. For many triangulation algorithms, a trapezoidal decomposition of the polygon is computed as a first step to triangulation. A trapezoid is a four-sided polygon with at least two parallel sides. The key difference between a trapezoidal decomposition and a triangulation is that the vertices defining a particular trapezoid can be additional points on the boundary of the polygon, whereas, in a triangulation it is not allowed to introduce new vertices.

One of the first triangulation algorithms was designed by Garey et al. [5]. Their algorithm is based on the method of dividing the input polygon into monotone pieces and it runs in O(n log n) time. The best known triangulation algorithm is by Chazelle [2] and it runs in O(n) time, where n is the number of vertices of the polygon. However, Chazelle’s algorithm is extremely complicated and is not practical for implementation. Even now, the algorithm by Garey et al. [5] is extensively used for its simplicity and programming ease. Chazelle and Incerpi [3] improved upon this algorithm by introducing the notion of sinuosity which is the number of spiraling and anti-spiraling polygonal chains on the polygon boundary. Their algorithm runs in O(n log s) time, where s is the sinuosity of the polygon. The sinuosity of a polygon is usually a small constant in most practical cases and the algorithm [3] is an almost linear-time algorithm. However, the algorithm by Chazelle and Incerpi [3] is based on a divide-and-conquer strategy and still quite complex to implement.

In this paper, we are interested in the trapezoidal decomposition of a simple polygon. In computational geometry, a triangulation is performed for building efficient data structures based on the decomposition of a polygon into regular pieces. Such data structures are used for many applications including processing the shortest path queries and ray shooting. On the other hand, in almost all applications in computer graphics and geographic information systems (GIS), the main requirement is the fast rendering of a polygon. A trapezoid is as simple a shape as a triangle and it is easy to render like a triangle in raster graphic devices. The fastest known algorithm for trapezoidation of general GIS data sets is by Lorenzetto et al. [6]. Their algorithm improves upon a previously published algorithm by Žalik and Clapworthy [8] and can perform a trapezoidal decomposition of general GIS polygons with holes in O(n log n) time, where n is the total number of vertices defining the polygon and all holes inside it. However, many GIS polygons do not contain holes inside them and it is interesting to investigate whether it is possible to design a simpler and faster algorithm for decomposing such polygons.

In this paper, we present a simple algorithm for computing the trapezoidal decomposition of a simple polygon without holes. In many GIS applications, it is often necessary to render many GIS polygons with a large number of vertices simultaneously. Common graphics APIs are usually slow in rendering such large polygons. Our algorithm can be used for fast rendering of large GIS data sets which quite often consist of polygons with thousands or hundreds of thousands of vertices. Our algorithm is simple and terminates in O(n) time for almost all polygons encountered in GIS applications, where n is the number of vertices of a given polygon. We have implemented our algorithm in C++ and tested it on many GIS data sets which include polygons with extremely large number of vertices. Our experimental results show that the running time is O(n) for almost all polygons encountered in GIS data sets. The hidden constant in big-Oh is very small for all the data sets. Our algorithm always produces a correct trapezoidation and this is ensured through an in-built correctness test in the algorithm. See Fig. 1 for an example polygon divided into trapezoids.

The method of rendering polygons in raster graphics devices is through scan conversion. The main idea behind this method is the setting of pixel colors inside a polygon by checking whether the pixel is inside or outside the polygon. This check is much easier if the polygon is convex. For convex polygons, a sweepline can be moved from the least y-coordinate vertex upwards and the span of the sweepline can be modified according to its intersection with the new polygon edges. It is sufficient to color the pixels on the current sweepline and as the sweepline progresses upwards, all the pixels inside the polygon are colored progressively. Since a convex polygon has only two intersections with a sweepline, the modification of the sweepline is simple for convex polygons. However, an arbitrary non-convex polygon may have many intersections with a sweepline and it is difficult to identify the parts of the sweepline inside the polygons. These are the parts of the sweepline that should be colored for coloring the interior of the polygon. Checking which parts of the sweepline are inside the polygon is a computationally intensive process and hence rendering an arbitrary polygon through scan conversion is infeasible for large polygons. Most graphic APIs adopt the alternative approach of converting an arbitrary polygon first into a collection of convex polygons. Since the simplest convex polygons are triangles and trapezoids, usually an arbitrary polygon is either triangulated or decomposed into trapezoids. In this paper, we adopt the approach of trapezoidal decomposition as it is usually easier to compute compared to triangulation. The trapezoidal decomposition has the added advantage that a single trapezoidal decomposition can be stored in a cache and rendered repeatedly if needed.

The rest of the paper is organized as follows. We discuss some preliminaries and an algorithm for trapezoidal decomposition of monotone polygons in Section 2. We present some background for the trapezoidation algorithm for general polygons in Section 3. Our algorithm for trapezoidation of general polygons is discussed in Section 4. Finally, we discuss the complexity of the algorithm, our experimental results and some concluding remarks in Section 5.