"Inside" is an easy predicate for people to determine, but rather harder for computers.
Introduction
Geometric algorithms and techniques are important tools in graphics, image processing, and computer vision applications. They also offer convenient and efficient solutions to pattern recognition problems such as the Nearest Neighbor rules, clustering, and image or polygon preprocessing.
In this two-part series, I present a C++ implementation of some standard and very efficient techniques in the field of Computational Geometry. This implementation includes fundamental operations such as orientation, test for point inclusion in a triangle, and segment intersection. These fundamental operations provide a basis for more complex operations, such as polygon orientation, point inclusion in a polygon, triangulation, convex hull computation for a set of points and for a simple polygon, and other standard algorithms. In the first part of this series, I present the mathematical basis and the fundamental operations, as well as the class definitions to implement these operations. The second part will describe some geometric algorithms and operations related to polygons.
The tools are based on the implementation of four classes, to represent points, segments, triangles, and polygons, all of them in the plane (two-dimensional).
The Mathematics of 2-D Operations
The basic operation at the heart of the algorithms presented is computation of a triangle's area given its three vertices (i.e., given the x-y coordinates of those vertices). This computation is efficiently performed based on the properties of the vector product.
Given two vectors v1, v2, their vector product v1 X v2 is perpendicular to both vectors; its magnitude is equal to the product of the magnitudes of v1 and v2 times the sine of the angle between them; and its orientation can be obtained using the Right Hand Rule.
If the two vectors are in the x-y plane (i.e., their z-coordinates are zero), their vector product will be parallel to the z-axis, and its z-coordinate will be positive if vector v2 is "at the left" of vector v1, zero if v1 and v2 are parallel, and negative if v2 is "at the right" of v1. (The term orientation refers to whether this z-coordinate is positive or negative.)
Given three points p1, p2, p3 in the x-y plane, the z-coordinate of the vector product between p2 - p1 (p2 minus p1) and p3 - p2 is given by the following formula:
z = x1(y2 - y3) + x2(y3 - y1) + x3(y1 - y2)Given the properties just discussed for the vector product, the above formula yields the following information (see Figure 1):
- The magnitude of z is twice the area of the triangle p1, p2, p3.
- The sign of z tells whether the triplet p1, p2, p3 represents a right turn or left turn (that is, if the point p3 is at the right or at the left of the oriented segment from p1 to p2).
Notice that calculating this value requires only three real-number multiplications and five additions. No trigonometric function calculations are involved.
Surprisingly enough, this simple operation forms the basis of most of the algorithms I present in this series, as well as many other extremely complex algorithms. Even more surprising is the fact that what is (almost exclusively) used in the algorithms is the sign of z (i.e., the orientation of triplets of points), rather than the magnitude or the value of z itself.
In addition to this fundamental operation and the basic arithmetic (including the scalar product), two other operations complete the basis for all the algorithms to be presented: a test for segment intersection, and a test for point inclusion in a triangle. These two operations are also based exclusively on the orientation of triplets of points.
Testing for Segment Intersection
Given two segments (p1, p2) and (q1, q2), they intersect if and only if the orientation of the triplet (p1, p2, q1) is different from the orientation of the triplet (p1, p2, q2) and the orientation of the triplet (q1, q2, p1) is different from the orientation of the triplet (q1, q2, p2). The first condition means that q1 is on one side of the segment (p1, p2), and q2 is on the other side. The second condition means that p1 is on one side of the segment (q1, q2), and p2 is on the other side. Clearly, the segments intersect if and only if both conditions are met, as illustrated in Figures 2a and 2b.
Testing for Point Inclusion in a Triangle
Given a triangle (v1, v2, v3) and a point p, the test for p's inclusion in the triangle is explained as follows. If we traverse the points v1, v2, v3 and the point is inside, we will always see the point on the same side of the segment we are visiting. If v1, v2, v3 are arranged in a counter-clockwise sense, the points inside it are always at the left of the segments. If the point is outside, for at least one of the segments the point will be at the right (see Figure 3). If the vertices are arranged in a clockwise sense, the reasoning is identical, except that a point that is inside the triangle will always be at the right of the segment we are visiting.
Thus, to determine if the point p is inside the triangle (v1, v2, v3), we need only obtain the directions of rotation along the triplets (v1, v2, p), (v2, v3, p), and (v3, v1, p). The point is inside if and only if the three directions are equal.
A Class to Represent Points
The representation for a point includes a pair of real numbers to store the x and y coordinates of the point (I use double), as well as member functions to read each coordinate (get_x and get_y). I provide the basic arithmetic operations in the form of overloaded operators. The arithmetic operations include addition, subtraction, and multiplication (scalar product). Also, I provide multiplication and division by real numbers in the form of overloaded operators.
The basic arithmetic operations allow you, for example, to obtain a point that is inside a triangle given its vertices:
p = (v1 + v2 + v3)/3;Also, an overloaded member function is provided to test for inclusion in a triangle and in a polygon. These functions allow for clear, concise expressions, as shown below:
Point p; Triangle t; Polygon P; if (p.is_inside (t)) // ... if (p.is_inside (P)) // ...Similarly, I provide a function to test if a point is exactly on a segment. This is determined by checking if the segment endpoints (p1, p2) and the point p are collinear, and if the segments p-p1 and p-p2 have opposite directions (which means that p is between p1 and p2).
Figure 4 shows the code for the class definition and implementation of class Point.
A Class to Represent Segments
The class Segment contains, as data members, two Point objects that represent a non-oriented segment. (However, most algorithms deal with a Segment as an oriented segment from p1 to p2.) I provide no operators for this class, since there are no natural arithmetic operations between segments.
Class Segment includes a member function to test for segment intersection. This member function allows expressions such as the following:
Segment s, line; if (s.intersects (line)) // ...Figure 5 shows the code for the class interface and implementation of class Segment.
Both the Point and Segment classes provide draw member functions. These member functions are platform dependent. I provide a simplified implementation for a Win32 platform (without scaling or other considerations). If you want to use these tools in an application that requires graphical display, you should provide the implementation for the draw member functions for Point and Segment.
A Class to Represent Triangles
The class Triangle includes three Point objects representing the vertices of the triangle. It provides member functions to compute the area and the orientation. (Class Triangle also provides a private member function to compute the area with sign, which enables other member functions to obtain the area and the orientation of the triangle.)
Since we are never interested in the exact value of the triangle's area (at least not in these types of applications), the member function actually computes twice the area of the triangle, for efficiency reasons. This may seem a little strange, but the reasons for computing twice the area should become clear in the second part of this series.
Class Triangle also provides a member function to test if a point is inside the triangle. This member function allows expressions such as the following:
Triangle t; Point p; if (t.contains (p)) // ...Note that the Point class member function to test for inclusion is equivalent to this Triangle class operation. The Point's member function is obviously implemented in terms of this Triangle member function (see Figure 4).
Class Triangle also includes a member function draw. In this case (as with the Polygon class, yet to be presented), the draw function is platform independent, since it is implemented in terms of the Segment class draw member function.
Figure 6 shows the code for the class definition and implementation of class Triangle.
None of the above-mentioned classes contains any pointers to dynamically allocated data or any allocated resources. Therefore, I provide no destructors, copy constructors, or overloaded assignment operators for either of these classes.
A Class to Represent Polygons
Representing a polygon requires a circular, doubly-linked list of Points. The list must be doubly-linked because typical operations on polygons require insertions and removals of points in the sequence. Some of these operations must traverse the sequence in reverse, visiting more than one point in succession. The linked list must be circular because polygons are closed sequences of vertices, and operations on pairs or triplets of consecutive vertices must not be limited by a vertex being the last or the first in the chain.
The implementation of class Polygon uses the STL (Standard Template Library) list container. For convenience of use, I provide an iterator for this class, as the standard containers do. Thus, polygons can be manipulated in the same manner as most of the standard containers, as shown below:
Polygon P; Polygon::iterator current = P.begin(); do { // ... } while (++current != P.begin());(notice that there is no end member function for the Polygon container, since it is a circular list.)
Also, the typical Polygon operations are similar to their equivalents for the standard containers (e.g., push_back, push_front, insert, remove), which should make polygon manipulation intuitive to C++ programmers.
No destructor, constructors, or assignment operator are required for this class, since the container included as a data member encapsulates those operations.
Figure 7 shows the definition of class Polygon and its iterator. Notice that the iterator implementation is based on the list iterator. Its operations are almost a direct map to the operations of the list iterator, except for the ++ and -- operators, since the Polygon is a circular list of points.
In the second part of this series, I will describe several geometric algorithms and polygon operations, most of them implemented as member functions of the class Polygon. I will also present the class implementation of Polygon in the second part of this series.
Summary
In Part 1 of this series, I have introduced the mathematical basis and some fundamental tools used in two-dimensional geometric algorithms. I have also presented part of the implementation of these tools in C++.
These fundamental operations are based on a simple and efficient computation of the orientations of point triplets, based on the properties of the vector product.
As Part 2 of this series will show, many other complex geometric algorithms are based almost exclusively on this operation.
Bibliography
Bjarne Stroustrup. The C++ Programming Language, Third 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 at http://www.mochima.com.