Delaunay Triangulation of Polygons

Motivation

For realtime rendering purposes (and similar), it may be desirable to convert some polygon meshes into a list of triangles. Such conversion is not easy however: many simplistic approaches could prove useless if the polygon to convert is star-shaped, for instance. A nice and very stable solution is Delaunay triangulation. (If you do not know it yet, check out the Wikipedia article; an image says everything.) This text is to relieve you of the eventual fear from mathematical theory that might have come upon you.

Algorithm

• Loop around vertices of the polygon.
• Circumscribe a circle to the current vertex and its two neighbors along the polygon.
• Check that no vertex of the polygon lies inside this circle.
• If not, output the three vertices as a triangle and remove the current vertex from the polygon.
• If yes, just continue.
• If you did not remove anything while making the last complete loop, remove any vertex you please and output it as a triangle along with its two neighbors.
• If the polygon has only three vertices, output it and exit.

The circumscribed circle need not be actually calculated. The function in_circle in code below checks if a point lies within the circle circumscribed to three other points perhaps in the most concise way possible.

Defense

According to the definition of Delaunay triangulation, we should check all triplets of vertices. How come that trying each vertex with its two neighbors is enough? It turns out that any triangulation of a polygon contains such a triangle of three consecutive vertices. You can prove that easily by counting the half-edges of a triangulation. The rest of the proof of correctness comes easily by induction.

The total complexity is obviously N³ in the worst case. For really large polygons, this might be a disadvantage. It can especially hurt when somebody tells you that this problem is solvable in N · log(N). Hopefully the pain subsides when you look at a classical implementation.

Example code

This piece of code should be quite easy to tailor to your specific needs. It uses STL containers only. There are several C++11 constructs but they should be rather easy to get rid of if necessary. The Point struct used here is 2D only; if you want to apply this to 3D polygons, it may be very helpful for you to use Newell's technique of calculating a surface normal.

Use this code at your free will.

#include <vector>
#include <algorithm>
#include <cassert>

struct Point
{
    float x, y;
    bool operator==(const Point & other) const { return x == other.x and y == other.y; }
};

using Polygon = std::vector<Point>;

inline float pow2(float x)
{
    return x * x;
}

float cross_product(const Point & a, const Point & b, const Point & c)
{
    return (b.y - a.y) * (c.x - b.x) - (b.x - a.x) * (c.y - b.y);
}

template<typename F>
inline F determinant(F aa, F ab, F ac, F ba, F bb, F bc, F ca, F cb, F cc)
{
    return aa * bb * cc + ba * cb * ac + ca * ab * bc - ac * bb * ca - bc * cb * aa - cc * ab * ba;
}

bool in_circle(Point a, Point b, Point c, Point p)
{
    a = {a.x - p.x, a.y - p.y};
    b = Point{b.x - p.x, b.y - p.y};
    c = Point{c.x - p.x, c.y - p.y};
    float inside = determinant(
        a.x, a.y, pow2(a.x) + pow2(a.y),
        b.x, b.y, pow2(b.x) + pow2(b.y),
        c.x, c.y, pow2(c.x) + pow2(c.y));
    return inside < 0;
}

inline bool nothing_in_circle(Point a, Point b, Point c, const Polygon & poly)
{
    if (cross_product(a, b, c) < 0) {
        return false;
    }
    return std::none_of(poly.begin(), poly.end(), [a, b, c](const Point & p) {
        if (a == p or b == p or c == p) {
            return false;
        }
        return in_circle(a, b, c, p);
    });
}

Polygon extract_delaunay_triangle(Polygon & poly)
{
    assert (poly.size() >= 3);
    Point prev = *(poly.end() - 2);
    auto here = poly.end() - 1;
    for (auto next=poly.begin(); next != poly.end(); prev = *here, here = next, next++) {
        if (nothing_in_circle(prev, *here, *next, poly)) {
            Point tmp_here = *here;
            Point tmp_next = *next;
            poly.erase(here);
            return Polygon({prev, tmp_here, tmp_next});
        }
    }
    Point tmp_here = *here;
    poly.erase(here);
    return Polygon({prev, tmp_here, poly.front()});
}

std::vector<Polygon> delaunay_triangulation(Polygon input)
{
    std::vector<Polygon> result;
    while (input.size() > 3) {
        result.emplace_back(extract_delaunay_triangle(input));
    }
    return result;
}