"Inside" is a remarkably subtle concept for a triangle. It becomes even more subtle for more complex polygons.
Introduction
Geometric algorithms are important tools in fields such as graphics, image processing, and pattern recognition, to name a few. This two-part series presents a set of classes that implement some standard and efficient techniques in the field of Computational Geometry. You can use or extend these classes to aid in development of your own geometric-intensive code. In Part 1 of this series, I introduced the mathematical basis of these tools. I also presented class definitions and implementations to handle some of the simpler entities in geometric operations: points, line segments, and triangles. I showed the definition of a Polygon class, but lacked the space to present the implementation. In this second part I present the Polygon implementation and describe some geometric operations related to polygons.
But first, a little review. As I explained last month, most of the algorithms are based on computation of the area of a triangle given the x-y coordinates of its vertices. This computation is efficiently performed using a formula that involves only three real-number multiplications and five additions. (See the section, "The Mathematics of 2-D Operations" in Part 1 of this series.)
This simple operation provides several pieces of information, including the orientation of a triplet of points. That is, for a triplet (p1, p2, p3), this operation tells whether the point p3 is at the left or at the right of the oriented line segment extending from p1 to p2. Finding these orientations of point triplets makes it easy to determine whether two segments intersect, as well as whether a point is inside a triangle. These three fundamental operations (determining orientation, segment intersection, and point inclusion in a triangle) are the basis (almost exclusively) for the algorithms that I present in this article. They form the basis of many other complex geometric algorithms as well.
Using the classes presented in Part 1 (for Point, Segment, Triangle, and Polygon) enables you to write simple expressions such as the following:
Point p, origin; Segment line; Triangle t; Polygon P; if (p.is_inside (t)) // ... if (t.contains (p)) // equivalent to // the previous expression // ... Segment seg (origin, p); if (line.intersects (seg)) // ... Polygon::const_iterator i = P.begin(); while (t.contains (*i) && tti != P.begin()); if (i == P.begin()) // all the vertices are inside // t...From the code above one of the payoffs of these classes ought to be obvious: the ability to do geometric operations using simple, intuitive notation.
Geometric Algorithms for Polygons
Now that I've given a quick review of Part 1, I am prepared to describe the geometric algorithms related to polygons, and how I implemented them in class Polygon (see Listing 1, polygon.cpp). Some of the geometric operations related to polygons are implemented as member functions, since they correspond to operations inherent to a polygon. Some operations are implemented as global functions that use the polygons and their related operations to perform the processing required.
Unless otherwise specified, this article assumes polygons are simple (i.e., no edges intersect, and no two vertices are the same point.). Also, when it is necessary to assume an orientation, this article will assume that the vertices are oriented in counter-clockwise sense. When dealing with polygons, the term orientation refers to the apparent sense of rotation when traversing a polygon's vertices. In particular, if the orientation is counter-clockwise, then the interior of the polygon will always "touch" the edges on their left. This further assumes that we are dealing with edges as oriented segments from one vertex to the next vertex — i.e., from Vn to Vn+1. (Note that in part 1 of this series, when I discussed the properties of the vector product, I used the term orientation to refer to a different concept. When talking about the vector product, orientation refers to where the result vector is pointing. This is not the same kind of orientation as the kind associated with polygons and point triplets.)
Testing for Point Inclusion in a Polygon
To test if a point p is inside a polygon P, we need only find a point that is outside the polygon, and count how many times the line L from p to the point outside intersects the polygon. To find a point outside the polygon, simply find the vertex with highest x-coordinate and add one to that x-coordinate.
Each time the line L intersects the polygon, it passes from inside to outside or vice versa. Thus, if the number of intersections is even, the point in question is outside the polygon. If it is odd, the point is inside.
We must be careful only when the line L passes exactly over a vertex of the polygon. In that case, we must check if L intersects the polygon (Figure 1a) or merely touches it (Figure 1b). We can determine this by checking if the polygon segments connected to the intersecting vertex are on different sides of L. If they are on different sides of L, it intersects the polygon; if they're on the same side, L merely touches the polygon.
Furthermore, if several consecutive vertices fall exactly on the line L, they must be discarded before checking if L intersects or touches the polygon (see Listing 1).
Testing for Polygon Orientation
The orientation of a polygon is often known in advance, since it may be guaranteed by the mechanism used to obtain the polygon. In some situations, though, we may have no more at our disposal than a sequence of points, with the knowledge that it represents a simple polygon. Determining the orientation of such a polygon consists of three steps:
1) Find a triplet of consecutive vertices such that the triangle formed by those three vertices does not contain any other vertex of the polygon. If no other vertex is inside, then the region inside the triangle is either entirely inside the polygon, or entirely outside the polygon. Part 1 of this series presented the techniques required for this purpose.
2) Determine whether the interior of the triangle found in 1) lies inside or outside of the polygon. To do this, take any point inside the triangle and test if it is inside the polygon. If a point inside the triangle is inside the polygon, the interior of the triangle lies inside the polygon as well. If a point inside the triangle is outside of the polygon, the interior of the triangle is outside of the polygon as well. To find a point guaranteed to be inside the triangle, use the formula
pi = (p1 + p2 + p3)/3where p1, p2, and p3 are the vertices of the triangle.
3) Determine the orientation of the triplet of points making up the triangle. (Part 1 of this series showed how to determine the orientation of a triplet.) If the the triplet is oriented counter-clockwise, and the triangle is inside the polygon, then the polygon is counter-clockwise. Also if the triplet is oriented clockwise, and the triangle is outside the polygon, then the polygon is counter-clockwise. The other two combinations correspond to a polygon with a clockwise orientation (See Listing 1).
Finding the Convex Hull of a Simple Polygon
The convex hull of a set R is defined as the smallest convex region (in area) that contains the entire set R, where R can be a set of points in the plane, in a region, in a polygon, in the interior of a polygon, etc. If the set R is a finite set of points or a polygon, it can be easily proven that the convex hull is a convex polygon. Figure 2 shows a simple polygon outlined with solid lines, and the convex hull of that polygon, in dashed lines.
For an arbitrary set of N points, it has been proven that no algorithm for finding the convex hull can have an order of complexity of less than O(N logN). If the points are known to form a simple polygon, the scenario is different. For a simple polygon, several algorithms have been proposed to find the convex hull in linear time (O(N) operations). Most of those algorithms have been proved to be incorrect, including the first O(N) algorithm proposed, known as the Sklansky scan (also known as the three coins algorithm). This algorithm was accepted and used for six years, until a counterexample was published, proving that the algorithm would work incorrectly for most arbitrary polygons. A detailed discussion about the various convex hull algorithms would be beyond the scope of this article. I want only to make it clear that the problem is far more complicated than it looks, despite the simplicity of the algorithm that I present in this article (the Melkman algorithm).
Melkman proposed an extremely elegant and straightforward algorithm to compute the convex hull of a simple polygonal chain in linear time. Clearly, such algorithm can be used for a simple polygon, if we remove the edge that joins the last vertex to the first vertex, we obtain a simple polygonal chain, and the convex hull of that polygonal chain is the same as that of the original polygon. Informally, a simple polygonal chain is a sequence of connected segments that form part of a simple polygon. An example appears in Figure 3. Unlike a simple polygon, a simple polygonal chain is not closed. The key idea in Melkman's algorithm is to maintain an updated version of the convex hull, and for each new point encountered in the input polygonal, to test the current convex hull to be updated.
The output at any step of the algorithm is a convex polygon with counter-clockwise orientation, whose first and last elements are the same point, and which correspond to the last point appended to the output. Thus, the algorithm starts by sending to the output the first three points of the input polygon (with the appropriate orientation).
At the step n of the algorithm, we test if the new point (point number n+1 of the input polygon) is outside the convex hull. If it is, then the convex hull must be updated to account for the new point. If not, the new point is discarded and the output remains the same. Figure 4 shows an example of the output convex hull, showing the first two and the last two points, as well as several possibilities of location for the new point (Pn+1).
Since the algorithm should work in linear time, the test to determine if the new point is inside the polygon must be done in constant time. Thus, we cannot use the function previously described to test for point inclusion in a polygon, since it does not operate in constant time. In this case, given the way we construct the current convex hull, Melkman showed a way to determine in constant time if the point is inside the current convex hull.
If the new point is at the right of the oriented segment P2-P1 and at the left of the oriented segment Pn-1-Pn, then it is certain that the point is inside the convex hull, since the line going from Pn to Pn+1 (which is one edge of the input polygon) would have to intersect the convex hull for the point Pn+1 to be outside. If the line intersects the convex hull, then it intersects the input polygon itself. But since the input polygon is simple, we know there are no edge intersections. Thus, by testing the orientation for these two triplets of points ((P2, P1, Pn+1) and (Pn-1, Pn, Pn+1)), we determine if the new point should be discarded or not. If the new point is inside the currently defined convex hull, we discard it. If the point is outside, then we append it at both ends of the output, and now we must readjust the output polygon by removing any vertices that became concave after the insertion of this new point.
All these updates involve only tests of triplet orientations, since we know that the output polygon is in a counter-clockwise orientation, meaning that a left turn represents a convex vertex, and a right turn represents a concave vertex.
Figure 5 shows an example of this updating process. Notice that P1 and Pn are actually the same physical point. I show them as two separate points in the figure only for illustration purposes.
Despite the fact that one step (after O(n) insertions) might involve O(n) removals, the algorithm does work in linear time (worst-case complexity O(N)). Each point from the input polygon can be inserted at most once into the output, and can be removed at most once. Therefore, we have at most 4*N operations. (Insertions are done at both ends of the output, and removals may be done from both ends.) Since the test to decide if the point should be discarded is done in constant time, the algorithm works in linear time.
Listing 1 shows the implementation details of this algorithm. It is implemented as a member function that returns a polygon without modifying the object that calls the member function.
Triangulating a Polygon
Triangulation is an important and very useful operation in many geometry applications. It provides a representation of a complex figure in simpler pieces — triangles.
The result of a triangulation process is a set of triangles for which each and every vertex must be one of the vertices of the polygon being triangulated. Also, the triangles' internal regions must add up the entire interior of the polygon. (This is not a formal and rigorous definition of a triangulation, but you get the idea.) Figure 6 shows an example of a polygon and a triangulation of it. Notice that for a given polygon, its triangulation is not necessarily unique.
Some applications require the output of a triangulation to be expressed as a graph in which adjacent triangles (triangles that share one edge) are adjacent in the graph. In this article, I restrict the scope of the triangulation algorithm to provide a non-ordered set of triangles.
The algorithm uses a divide-and-conquer recursive approach, which means that its expected complexity is O(N logN). Its worst case complexity, however, is O(N2).
The algorithm requires knowledge of the input polygon's orientation, so it just assumes that the input polygon has a counter-clockwise orientation. This assumption enables the algorithm to identify a convex vertex, that is, when the triplet (previous point, current point, next point) makes a left turn.
The algorithm consists of the following steps:
1) Find a convex vertex in the input polygon.
2) Consider the triangle Vprev-Vcurr-Vnext containing this convex vertex, such that Vprev is the point previous to the convex vertex in the polygon; Vcurr is the convex vertex; and Vnext is the point succeeding the convex vertex in the polygon.
3) Determine if this triangle contains any other vertices of the input polgon.
a) If this triangle contains no other vertices of the input polygon, then the segment Vprev-Vnext does not intersect any edge of the polgon. The triangle is sent to the output and the vertex Vcurr is removed from the polygon.
b) If the triangle does contain other vertices of the input polygon, find one of these vertices such that it can be joined to the vertex Vcurr without intersecting any edges of the polygon. Call this vertex Vs. The segment
Vcurr-Vs divides the input polgon into two simple polygons. Recursively call the triangulation function to triangulate each half of the input polygon (see Figure 7). This recursive process ends when the input polygon has only three vertices — that is, when it is already a triangle.Finding a vertex that can be joined to the vertex Vcurr, as required in step 3b) above, is not as simple as it appears. Many wrong approaches have been used in the past, including the most intuitive, which is finding the nearest vertex to Vcurr. Figure 8a illustrates a situation in which this approach fails, since it produces an output with triangles that intersect the input polygon. The appropriate way to find such vertex is to find the point with highest perpendicular distance to the straight line containing the segment Vprev-Vnext, as illustrated in Figure 8b.
Naturally, to obtain this point the algorithm does not directly compute a single distance (it would be too hard, since it is perpendicular distance). Instead, the algorithm finds the point P for which the area of the triangle Vprev-Vnext-P is maximized. A triangle's area is equal to the base of the triangle times the height divided by two. The base of the triangle is the length of segment Vprev-Vnext; the height is the perpendicular distance to the base segment. As I mentioned in part 1, the member function area actually returns twice the area. This is no problem, since the algorithm just needs to compare areas to find the point that produces the largest one — the actual value of the area is of no interest here.
The prototype for the triangulation function is:
void triangulate (const Polygon &, list<Triangle> &);Listing 2, triangle.cpp, shows the implementation of this algorithm.
Finding the Convex Hull of an Arbitrary Set of Points
If you want to find the convex hull of a set of points, and the points are not guaranteed to be ordered in such a way that they form a simple polygonal chain, the problem is different from that presented in the section "Finding the Convex Hull of a Simple Polygon." Obviously, an approach to do this in O(N logN) would be to sort the points by their x-coordinate, which would allow the building of a simple polygonal chain. Since the points were sorted by their x-coordinates, there could be no edge intersections, and we could then use a linear time algorithm to obtain the convex hull. A much more straightforward approach, known as the Jarvi's March, is the following: first find a point that is on the convex hull — for example, the point with lowest y-coordinate. This can be done in linear time. After that, find the right-most point with respect to that first point. The right-most point with respect to point P is the point X for which no other point of the set is at the right of the oriented segment from P to X. The right-most point is also on the convex hull, so append it to the output. Continue to find the right-most point with respect to the last point added. Keep going until the right-most point found is the first point.
Notice that finding the right-most point requires no angle calculations — it can be done by computing the orientations of point triplets. Listing 3 shows the code for this algorithm.
Conclusion
Geometric algorithms and techniques are important and convenient tools in many applications, including Image Processing, Graphics, and Pattern Recognition.
Efficiency is a very important (maybe the most important) issue, since applications that involve this type of processing in general act on huge amounts of information. Polygons consisting of several hundreds or thousands of vertices are common. Furthermore, this type of application often requires real-time processing, with very demanding timing restrictions.
The approach that I presented in this article is based on very efficient mathematical tools. Most geometric operations require no computation of trigonometric functions, square roots, etc. Instead, the operations require a small number of multiplications and additions per item processed. In some cases, I sacrificed a certain amount of efficiency in favor of "straightforwardness" of the code. I did this only in situations where I considered that the impact on the overall speed of execution would not be considerable.
The object-oriented approach turns out to be quite appropriate for this application. Not only are geometric entities intuitively represented as software objects, the OO approach leads to expressions that are intuitive and readable. Also, the OO approach does not imply inefficient implementation. The C++ compilers I have used provide very efficient machine code for this application. As I mentioned before, efficiency is one of the most important considerations in this type of application.
Bibliography
Bjarne Stroustrup. The C++ Programming Language, 3rd Edition (Addison-Wesley, 1997).
Franco Preparata and Michael I. Shamos. Computational Geometry, An Introduction (Springer-Verlag, 1985).
Avraham A. Melkman. "On-line Construction of the Convex Hull of a Simple Polyline," Information Processing Letters 25, 1987.
Carlos Moreno has a Bachelor's Degree in Electronics Engineering and a Master Engineering diploma in Computer and Electrical Engineering. He has ten years' experience in the development of low-level applications, and currently works as an instructor/trainer in C, C++, Object-Oriented Programming, and Visual Basic, and as an independent consultant and software developer. His main interests are digital signal processing, audio and image processing applications, communications, data encryption and compression, and (of course) computer games development. He can be reached at moreno@mochima.com or moreno@cyberglobe.net, or, on the web at http://www.mochima.com.