The Problem
OpenGL and other low-level rendering APIs are limited to rendering convex polygons. This causes a complications to arise when trying to render concave primitives. In my case, I wanted to be able to render arbitrary 2D "ground," and not succumb to manually breaking apart the polygon into OpenGL-compatible primitives. Thus the search for a triangulation algorithm arose. All I could find were algorithm names without accompanying implementations or explanations, and long mathematical abstracts on computational geometry. Any implementations I did find seemed to be extremely over-engineered and over-complicated for easy comprehension and adaptation. I finally settled on writing my own implementation of the ear-clipping algorithm to facilitate triangulation.What We Will Accomplish
 |
 |
| 139-vertex polygon | Triangulated |
What is "Ear Clipping"?
Ear clipping is a triangulation algorithm most easily understood recursively. An "ear" is described as a vertex that forms a triangle with its 2 adjacent vertices that contains no other vertices of the polygon within it. These ears are found and removed from a polygon recursively. I found a high-level implementation of this algorithm by David Eberly in his paper on ear clipping here. There are only a few steps:1. Iterate through vertices
a. If the vertex is concave, skip
b. If the triangle formed by the vertex has other vertices in it, skip
b. Otherwise, ear found
2. "Cut off" the ear from the polygon and store it
3. Go to 1
If you'd like to see a functioning version of what we're going to be doing, and you have Python/Pygame installed, you can run this script to see the algorithm in action. Click the left mouse button to create vertices, right click to finalize a shape, and hit "T" to triangulate the polygon.
Algorithmic Complexity
As described by Eberly, a straightforward implementation of the algorithm is theoretically O(n3), but in practice the average case is Θ(n2). It is possible to optimize the algorithm to be O(n2) in the worst case, but is somewhat more complicated and less intuitive than this tutorial is intended to be. We will be implementing the O(n3) implementation for clarity's sake. We will include several optimizations to reduce iterations. Since we stop looking as soon as we find an ear, it's hardly ever the case where this takes more than a few iterations.In testing, a shape with 10 vertices took 93 iterations to triangulate, 30 took 770, and an incredibly complex polygon with 80 vertices took ~6100 iterations. Thus we see that although the recursive algorithm we employ is O(n3) in theory, the average case in practice is Θ(n2). In addition, purely convex shapes find an ear on the first iteration every time, which guarantees O(n2) complexity.
Finding Ears
Ears have several limitations that allow for us to process and find them faster.- Ears cannot be reflex vertices (θ > 180°).
- The triangle formed by the ear cannot contain any other vertices in the polygon.
- The triangle formed by the ear must have the same orientation as the polygon.
There are multiple ways to find ears, some of which are more computationally expensive than others.
The first way is to determine the angle inside of the polygon that the vertex forms with its two adjacent vertices. This is not easily done, because simply finding the angle between two vectors will always return an angle <= 180° due to the way arc-tangent works.
Let's call the vector formed by the vertex n and n+1 A, and the vector formed by the vertex n and n-1 B. If you compare the angle of the normal vector of A to that of B, they will be <= 90° if the angle is convex, and the angle is concave (and thus reflex) if the angle is > 90°. The diagram below outlines why this process works.
Pseudo-code for the reflexive test could look something like this:
// Assume we are at the n-th index.
POINT curr = points[n];
POINT next = points[n+1];
POINT prev = points[n-1];
VECTOR One(prev.x - curr.x, prev.y - curr.y);
VECTOR Two(next.x - curr.x, next.y - curr.y);
VECTOR Nml(One.y, -One.x);
float magnitude = Nml.Magnitude() * Two.Magnitude();
float angl = acos(Nml.Dot(Two) / magnitude);
angl *= 180.0 / 3.14159265359; // Convert to degrees.
// If angl <= 90.0, the n-th vertex is not a reflex vertex.
A limitation of this method is that it only works as-is in one orientation; it assumes the polygon is formed in a clockwise fashion. A much simpler and faster way of determining whether an ear vertex is reflex is by checking the orientation of its respective triangle, and comparing it to the orientation of the polygon.
To get an accurate polygon orientation, we find the topleft-most vertex in the polygon, and calculate the orientation of its respective triangle. This requires a single pass through the whole polygon prior to clipping. The algorithm is fairly straightforward, calculating the cross-product of two sides of the triangle. Pseudo-code looks something like this:
// Input
VECTOR a, b, c; // the parts of our triangle.
// if (b.x - a.x) * (c.y - b.y) - (c.x - b.x) * (b.y - a.y) > 0
// we are counter-clockwise
Before clipping, we run through the polygon like so:
// Input
POLYGON p;
int vcount;
// We assume len(shape) > 0
int index = 0;
VERTEX leftmost = p[0];
for i = 0 to vcount - 1
if p[i].x < leftmost.x or (p[i].x == leftmost.x and p[i].y < leftmost.y):
leftmost = p[i];
index = i
// Make a triangle from the vertex.
POLYGON tri = {
p[index - 1 if index - 1 > 0 else vcount - 1],
p[index],
p[index + 1 if index + 1 < vcount else 0]
};
bool orient = is_ccw(tri[0], tri[1], tri[2]);
Finally, we need a function to test if a given vertex lies within a triangle. This is done by calculating barycentric coordinates. A sample algorithm was given on StackOverflow, for which a pseudo-code adaptation is given below.
// Input
POINT p;
TRIANGLE t; // A set of three POINTs.
float denom = (t[1].y - t[2].y) * (t[0].x - t[2].x) +
(t[2].x - t[1].x) * (t[0].y - t[2].y);
// Avoid division by zero.
if denom == 0: return true;
denom = 1.0 / denom;
float alpha = denom * ((tri[1].y - tri[2].y) * (p.x - tri[2].x) +
(tri[2].x - tri[1].x) * (p.y - tri[2].y));
float beta = denom * ((tri[2].y - tri[0].y) * (p.x - tri[2].x) +
(tri[0].x - tri[2].x) * (p.y - tri[2].y));
float gamma = 1.0 - alpha - beta;
// if alpha < 0, the point is NOT in the triangle.
// if beta < 0, the point is NOT in the triangle.
// if gamma < 0, the point is NOT in the triangle.
// else, the point IS in the triangle (go figure).
Armed with our slew of helper functions, we can begin the actual ear clipping algorithm. A pseudo-code implementation is provided below, and here is a link to a full C++ implementation. You can also view a Python implementation in the sample script I provided above.
And there you have it. An algorithm for creating completely triangulated, OpenGL-compatible polygons that you can then use to draw arbitrary landscapes, create collision maps, or anything else you may need concave polygons for.
A comparison of Ear Clipping and a New Triangulation Algorithm - Ran LiuSources / Further Reading
Triangulation By Ear Clipping - David Eberly
Triangle orientation - GameDev.net
A True O(n2) Triangulation Algorithm in C++ - urraka
Great post!
ReplyDeleteBTW, the link to the Python/Pygame script is dead; is there an alternative link to it?
Hi brouwer, I've updated the links! Thanks for noticing; I haven't done much posting in a while. :)
Delete