xi. Polygon Clipping and Union

Polygon Clipping

Sutherland-Hodgman clipping of any polygon by a convex one, area of intersection of two convex polygons, triangulation of a simple polygon, and 2D affine transforms as 3x3 matrices.

vector<PD> c = clipPolygon(subject, clip);            // clip: convex CCW; result may be empty
double a = convexIntersectionArea(p, q);               // both convex CCW
vector<array<int, 3>> tris = earClipping(poly);        // simple polygon CCW, n - 2 triangles (indices)
Affine t = Affine::rotate(theta) * Affine::translate(dx, dy);   // applied right to left
PD q2 = t.apply(p);   Affine s = Affine::scale(sx, sy), r = t.inverse();

Complexity: clipping $O(nm)$; intersection $O(nm)$; ear clipping $O(n^2)$; affine $O(1)$.

vector<PD> clipPolygon(vector<PD> subject, const vector<P>& clip) {   // repeated single-edge clips
    int m = clip.size();
    for (int i = 0; i < m && !subject.empty(); i++) subject = clipLeft(subject, clip[i], clip[(i + 1) % m]);
    return subject;
}
 
double convexIntersectionArea(const vector<PD>& p, const vector<P>& q) {
    auto c = clipPolygon(p, q);
    if (c.empty()) return 0;
    double a2 = 0;
    for (int i = 0; i < (int)c.size(); i++) a2 += cross(c[i], c[(i + 1) % c.size()]);
    return fabs(a2) / 2;
}

`clipLeft` (one half-plane clip step) does the real work; chaining it over all $m$ edges of a convex clip polygon is Sutherland-Hodgman. Concave clip polygons need Weiler-Atherton or the polygon-union sweep below.

Ear clipping — triangulating a simple polygon

A convex vertex whose triangle holds no other vertex of the polygon is an ear; removing ears one at a time keeps the remainder simple, and every simple polygon with $n\ge 4$ vertices has at least one ear (the two ears theorem).

vector<array<int, 3>> earClipping(const vector<P>& poly) {
    int n = poly.size();
    vector<int> idx(n);
    iota(idx.begin(), idx.end(), 0);
    vector<array<int, 3>> res;
    auto inside = [&](P p, P a, P b, P c) {
        return orient(a, b, p) >= 0 && orient(b, c, p) >= 0 && orient(c, a, p) >= 0;
    };
    while (idx.size() > 3) {
        int k = idx.size();
        bool found = false;
        for (int i = 0; i < k && !found; i++) {
            int pa = idx[(i + k - 1) % k], pb = idx[i], pc = idx[(i + 1) % k];
            if (orient(poly[pa], poly[pb], poly[pc]) <= 0) continue;   // reflex or flat
            bool ear = true;
            for (int j = 0; j < k && ear; j++) {
                int v = idx[j];
                if (v == pa || v == pb || v == pc) continue;
                if (inside(poly[v], poly[pa], poly[pb], poly[pc])) ear = false;
            }
            if (!ear) continue;
            res.push_back({pa, pb, pc});
            idx.erase(idx.begin() + i);
            found = true;
        }
        if (!found) break;   // degenerate input
    }
    if (idx.size() == 3) res.push_back({idx[0], idx[1], idx[2]});
    return res;
}

Complexity: $O(n^2)$ (the ear search is $O(n)$ per removal, over $n-2$ removals); $O(n\log n)$ algorithms exist (sweep-based, or monotone-polygon decomposition) but are rarely needed.

Proof: the two ears theorem and ear clipping's correctness

Two ears theorem. Every simple polygon with $n\ge 4$ vertices has at least two ears (not sharing an edge). Proof sketch by induction: take any convex vertex $v$ (one exists, e.g. the lexicographically smallest point, which is always convex). If triangle $(prev(v),v,next(v))$ contains no other vertex, $v$ is an ear — done. Otherwise, among the vertices strictly inside that triangle, let $w$ be the one farthest from line $prev(v)next(v)$; then the segment from $v$ to $w$ cannot be blocked by any other vertex (by the extremal choice of $w$), so it is a valid diagonal splitting the polygon into two simple sub-polygons, each with fewer vertices; by strong induction each sub-polygon has an ear, and at most one of the two new ears in each piece can be the artificial vertex $v$ or $w$ shared by the split, so at least one genuine ear survives in the original polygon on each side, giving at least two.
Correctness of the algorithm. Removing a true ear produces a strictly smaller simple polygon (the removed triangle lies entirely outside… rather entirely inside the original polygon and does not overlap the rest, since no vertex is inside it and the two new boundary edges replace the two old ones without crossing anything else) — so the invariant “the remaining vertices form a simple polygon” is maintained at every step, and by the two-ears theorem a valid ear always exists to remove until only a triangle remains. Each removal is found in $O(n)$ by scanning for a convex vertex whose triangle contains no other remaining vertex, giving $O(n)$ removals $\times$ $O(n)$ search $=O(n^2)$. $■$

2D affine transforms as matrices

Homogeneous coordinates turn translation into a matrix product, so any chain of rotations, scalings and shifts is one $3\times 3$ matrix, applied right to left.

struct Affine {
    array<array<double, 3>, 3> m{};
    static Affine identity() { Affine a; for (int i = 0; i < 3; i++) a.m[i][i] = 1; return a; }
    static Affine translate(double dx, double dy) { Affine a = identity(); a.m[0][2] = dx; a.m[1][2] = dy; return a; }
    static Affine rotate(double theta) {                 // counter-clockwise around the origin
        Affine a = identity();
        a.m[0][0] = cos(theta); a.m[0][1] = -sin(theta); a.m[1][0] = sin(theta); a.m[1][1] = cos(theta);
        return a;
    }
    static Affine scale(double sx, double sy) { Affine a = identity(); a.m[0][0] = sx; a.m[1][1] = sy; return a; }
    friend Affine operator*(const Affine& a, const Affine& b) {
        Affine c;
        for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) c.m[i][j] += a.m[i][k] * b.m[k][j];
        return c;
    }
    PD apply(PD p) const { return {m[0][0] * p.x + m[0][1] * p.y + m[0][2], m[1][0] * p.x + m[1][1] * p.y + m[1][2]}; }
    Affine inverse() const {                              // of the 2x2 linear part plus the translation
        double det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
        Affine r = identity();
        r.m[0][0] = m[1][1] / det; r.m[0][1] = -m[0][1] / det;
        r.m[1][0] = -m[1][0] / det; r.m[1][1] = m[0][0] / det;
        r.m[0][2] = -(r.m[0][0] * m[0][2] + r.m[0][1] * m[1][2]);
        r.m[1][2] = -(r.m[1][0] * m[0][2] + r.m[1][1] * m[1][2]);
        return r;
    }
};

Proof: homogeneous composition and the inverse

Writing a point as $(x,y,1{)}^T$, a linear map $L$ plus a translation $t$ acts as $p\mapsto Lp+t$, which is exactly $\left(\begin{array}{cc}L& t\\ 0& 1\end{array}\right)\left(\begin{array}{c}p\\ 1\end{array}\right)$ — the given $3\times 3$ form, with the bottom row $(0,0,1)$ implicit (never written to, so it starts and stays that way from identity). Composing two such maps, $p\mapsto L_2(L_{1}p+t_1)+t_2=L_2{L}_{1}p+(L_2{t}_1+t_2)$, matches exactly the ordinary $3\times 3$ matrix product of the two homogeneous matrices — which is why operator* is plain matrix multiplication and transforms compose right-to-left (the rightmost matrix is applied to the point first). The inverse of $\left(\begin{array}{cc}L& t\\ 0& 1\end{array}\right)$ is $\left(\begin{array}{cc}{L}^{-1}& -L^{-1}t\\ 0& 1\end{array}\right)$ (block inverse: apply $L^{-1}$ then undo the translation $t$ in the original, un-rotated frame), matching the code’s $2\times 2$ cofactor inverse for L^{-1} and the two translation lines for $-L^{-1}t$. $■$

Tested: earClipping’s triangles sum to the polygon’s shoelace area and tile it without gaps or overlaps (checked via total area and via no two triangles’ interiors intersecting) on random simple polygons; Affine composition and inverse checked against direct coordinate transformation and round-trip (t.inverse().apply(t.apply(p)) == p) on random transforms; clipPolygon/convexIntersectionArea tested in Half-Plane Intersection’s test suite (same underlying clipLeft).

Polygon Union

Exact-up-to-precision area covered by at least one of several simple polygons.

vector<vector<PD>> polys;                  // each polygon CCW (CW = hole, see below)
double a = polygonUnionArea(polys);

Complexity: $O(N^2\log N)$, $N=$ total number of vertices.

Idea (Green’s theorem): the signed area of a simple polygon is $\frac{1}{2}{\sum }_{\text{edges}}\mathrm{cross}(A,B)$. For a union of possibly-overlapping polygons, apply the same sum but weight each edge only by the fraction of it that is not covered by any other polygon — overlapping regions get credited to exactly one edge, not double-counted.

double polygonUnionArea(const vector<vector<PD>>& poly) {
    auto ratio = [](PD a, PD b) { return fabs(b.x) > 1e-12 ? a.x / b.x : a.y / b.y; };
    double res = 0;
    for (int i = 0; i < (int)poly.size(); i++) {
        for (int v = 0; v < (int)poly[i].size(); v++) {
            PD A = poly[i][v], B = poly[i][(v + 1) % poly[i].size()];
            vector<pair<double, int>> segs = {{0, 0}, {1, 0}};
            for (int j = 0; j < (int)poly.size(); j++) {
                if (i == j) continue;
                for (int u = 0; u < (int)poly[j].size(); u++) {
                    PD C = poly[j][u], D = poly[j][(u + 1) % poly[j].size()];
                    int sc = sgn(cross(A, B, C)), sd = sgn(cross(A, B, D));
                    if (sc != sd) {
                        double sa = cross(C, D, A), sb = cross(C, D, B);
                        if (min(sc, sd) < 0) segs.emplace_back(sa / (sa - sb), sgn(sc - sd));
                    } else if (!sc && !sd && j < i && sgn(dot(B - A, D - C)) > 0) {
                        segs.emplace_back(ratio(C - A, B - A), 1);
                        segs.emplace_back(ratio(D - A, B - A), -1);
                    }
                }
            }
            sort(segs.begin(), segs.end());
            for (auto& s : segs) s.first = min(max(s.first, 0.0), 1.0);
            double sum = 0;
            int cnt = segs[0].second;
            for (int k = 1; k < (int)segs.size(); k++) {
                if (!cnt) sum += segs[k].first - segs[k - 1].first;
                cnt += segs[k].second;
            }
            res += cross(A, B) * sum;
        }
    }
    return res / 2;
}

Clockwise polygons subtract (holes): a CW polygon’s edges contribute negative cross(A,B), so its area is removed from the total — exactly the usual signed-area convention extended to a multi-polygon union with holes.

Proof: polygonUnionArea

Fix one edge $A\to B$ of polygon $i$. Parametrize it by $t\in [0,1]$, $P(t)=A+t(B-A)$. The claim is that the edge should contribute to the total (signed) area sum only over the sub-interval of $t$ where $P(t)$ is not “covered” in the specific sense needed for the sum to equal the true union area — precisely, where the ray from $P(t)$ in the edge’s outward normal direction exits the union through this edge rather than through some other polygon’s boundary lying “in front of” it from the perspective of Green’s theorem’s line integral. For each other polygon $j$, each of its edges $C\to D$ that crosses the line through $A,B$ (strictly, $sc\ne sd$) marks the $t$-value where the sweep along $AB$ transitions into or out of being blocked by polygon $j$‘s interior at that crossing — computed by intersecting line $CD$ with line $AB$ (sa/(sa-sb), the same two-line intersection formula as Lines and Segments), and only recorded (min(sc,sd)<0) when the crossing happens on the side of line $CD$ that matters for blocking (one of the two crossing orientations marks entering coverage, the other leaving — sgn(sc-sd) records which). Collinear overlapping edges (!sc && !sd) are handled separately: when edge $CD$ runs along the same line as $AB$ in the same direction (dot(B-A,D-C)>0), they cover each other’s extent directly (as a position-ratio interval via ratio), and the j<i guard ensures that among several mutually-collinear same-direction duplicate edges across different polygons, only the first polygon’s copy is treated as the “real” one — avoiding both double-subtracting and double-crediting when polygons share a boundary exactly.
After collecting and sorting all these breakpoints along $[0,1]$, cnt tracks the net “coverage depth” entering/leaving as $t$ sweeps — a standard coordinate-compressed interval-coverage sweep (the same pattern as an interval-union length computation) — and sum accumulates exactly the length of the sub-interval where cnt==0, i.e. where this edge is not shadowed by any other polygon. Multiplying by $\mathrm{cross}(A,B)$ and summing over all edges of all polygons gives, by the additivity of Green’s theorem’s line integral over the edges actually forming the union’s true boundary, exactly twice the union’s signed area (holes subtract automatically since a CW polygon’s edges contribute with the opposite sign). $■$

Tested against Monte Carlo area estimation (random sample points in a bounding box, fraction classified inside any polygon via `inPolygon`, scaled by the box area) on random sets of 2-4 overlapping convex and non-convex polygons, including a polygon fully inside another and a CW polygon used as a hole.