Polygon Clipping: Sutherland-Hodgman and Boolean Operations
Clip segments and polygons against windows with Liang-Barsky and Sutherland-Hodgman, intersect half-planes, and see why rounded intersection points and touching edges make general Boolean operations hard.
Cutting shapes with other shapes
A map renderer draws only the part of each road and lake that falls inside the screen. A CNC toolpath generator intersects the part's outline with the stock. A GIS answers "how much of this parcel lies in the flood zone?" A game computes what a light or a camera can see. Each of these takes one polygon, the subject, and keeps the part of it inside another, the clip polygon or window. That is clipping, and its general form, union, intersection and difference of arbitrary polygons, is called Boolean operations.
Clipping brings something new to this course: the output contains points that are not in the input. Where a subject edge crosses a window edge, the result has a vertex with rational, not integer, coordinates. Up to now we only asked questions about the input and answered them exactly. Now we have to construct new geometry, decide how to represent it, and live with what that choice does to the next step.
This lesson builds the classic algorithms in order of difficulty: a segment against a rectangle (Liang-Barsky), a polygon against a convex window (Sutherland-Hodgman), and the intersection of many half-planes. Then it looks at what goes wrong in general Boolean operations and how the main libraries deal with it. It uses the segment test from the segment intersection lesson and the shoelace area from the polygons lesson. Polygons are std::vector<Point>, counter-clockwise, without a repeated closing vertex; input coordinates are integers with |x|, |y| <= 10^9.
A segment against a rectangle
Start with one segment and an axis-aligned rectangle. Write the segment as P(t) = a + t * (b - a) for t from 0 to 1. Each side of the rectangle is one inequality on t. For the left side x >= xmin, it reads a.x + t * dx >= xmin, which is -dx * t <= a.x - xmin. All four have the form p * t <= q.
- If
p < 0, the inequality is a lower bound ont: the segment is entering the half-plane att = q / p. - If
p > 0, it is an upper bound: the segment is leaving att = q / p. - If
p = 0, the segment is parallel to that side; it is entirely outside ifq < 0, and the side imposes nothing otherwise.
Keep the largest entering value t0 (start at 0) and the smallest leaving value t1 (start at 1). The visible part is t in [t0, t1], or nothing if t0 > t1. That is the Liang-Barsky algorithm (Liang and Barsky, 1984). The older Cohen-Sutherland algorithm gives each endpoint a 4-bit code saying which sides it is beyond and cuts the segment one side at a time; Liang-Barsky does the same job without ever shortening the segment step by step: it compares four parameter values and computes the endpoints once.
The decision is a comparison of fractions, and fractions with integer numerators and denominators compare exactly: a/b < c/d is a*d < c*b when b and d are positive. So the program decides visibility exactly and computes the endpoints only for printing. The last part checks every segment between points of a 9 by 9 grid against an independent definition: the segment meets the closed rectangle if an endpoint is inside or it intersects one of the four sides.
#include <algorithm>
#include <numeric>
#include <print>
#include <string>
struct Point { long long x, y; };
struct Rect { long long xmin, ymin, xmax, ymax; };
// A parameter t = num / den along the segment, with den > 0.
struct Param { long long num, den; };
bool before(Param a, Param b) { return a.num * b.den < b.num * a.den; } // exact
struct Clipped { bool visible; Param t0, t1; };
// Liang-Barsky: the segment is P(t) = a + t * (b - a) for t in [0, 1]. Each side
// of the rectangle is one inequality p * t <= q. Entering sides raise t0,
// leaving sides lower t1, and the segment is visible if t0 <= t1 at the end.
Clipped liangBarsky(Point a, Point b, Rect r) {
long long dx = b.x - a.x, dy = b.y - a.y;
const long long p[4] = {-dx, dx, -dy, dy};
const long long q[4] = {a.x - r.xmin, r.xmax - a.x, a.y - r.ymin, r.ymax - a.y};
Param t0{0, 1}, t1{1, 1};
for (int k = 0; k < 4; ++k) {
if (p[k] == 0) { // parallel to this side
if (q[k] < 0) return {false, t0, t1}; // and outside it
continue;
}
// t = q / p, written with a positive denominator.
Param t = p[k] > 0 ? Param{q[k], p[k]} : Param{-q[k], -p[k]};
if (p[k] < 0) { if (before(t0, t)) t0 = t; } // entering: keep the latest
else { if (before(t, t1)) t1 = t; } // leaving: keep the earliest
}
return {!before(t1, t0), t0, t1};
}
std::string show(Param t) {
long long g = std::gcd(t.num, t.den);
return t.den == g ? std::to_string(t.num / g)
: std::to_string(t.num / g) + "/" + std::to_string(t.den / g);
}
// The construction: a point on the segment, rounded only when printed.
void printAt(Point a, Point b, Param t) {
double x = double(a.x) + double(b.x - a.x) * double(t.num) / double(t.den);
double y = double(a.y) + double(b.y - a.y) * double(t.num) / double(t.den);
std::print("({:g}, {:g})", x, y);
}
// Independent check: a segment meets a closed rectangle if an endpoint is
// inside it or the segment intersects one of its four sides (lesson 3's test).
int orientation(Point a, Point b, Point c) {
long long v = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
return (v > 0) - (v < 0);
}
bool onSegment(Point a, Point b, Point p) { // p is collinear with a and b
return std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x) &&
std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}
bool segmentsMeet(Point a, Point b, Point c, Point d) {
int o1 = orientation(a, b, c), o2 = orientation(a, b, d);
int o3 = orientation(c, d, a), o4 = orientation(c, d, b);
if (o1 * o2 < 0 && o3 * o4 < 0) return true; // a proper crossing
return (o1 == 0 && onSegment(a, b, c)) || (o2 == 0 && onSegment(a, b, d)) ||
(o3 == 0 && onSegment(c, d, a)) || (o4 == 0 && onSegment(c, d, b));
}
bool meetsRect(Point a, Point b, Rect r) {
auto in = [&](Point p) {
return r.xmin <= p.x && p.x <= r.xmax && r.ymin <= p.y && p.y <= r.ymax;
};
Point c[4] = {{r.xmin, r.ymin}, {r.xmax, r.ymin}, {r.xmax, r.ymax}, {r.xmin, r.ymax}};
if (in(a) || in(b)) return true;
for (int i = 0; i < 4; ++i)
if (segmentsMeet(a, b, c[i], c[(i + 1) % 4])) return true;
return false;
}
int main() {
const Rect window{0, 0, 10, 6};
struct Case { const char* name; Point a, b; };
const Case cases[] = {
{"inside", {2, 2}, {8, 4}},
{"crosses two sides", {-4, 1}, {14, 5}},
{"misses the window", {-2, 4}, {4, 12}},
{"touches a corner", {8, 8}, {12, 4}},
{"runs along an edge", {-3, 6}, {4, 6}},
{"a single point", {5, 5}, {5, 5}},
};
for (const Case& c : cases) {
Clipped r = liangBarsky(c.a, c.b, window);
std::print("{:<19} ", c.name);
if (!r.visible) { std::println("rejected"); continue; }
std::print("t in [{}, {}] ", show(r.t0), show(r.t1));
printAt(c.a, c.b, r.t0);
std::print(" to ");
printAt(c.a, c.b, r.t1);
std::println("");
}
// Every segment between points of a 9 by 9 grid around a 4 by 2 window.
long tested = 0, disagree = 0;
const Rect small{2, 3, 6, 5};
for (int i = 0; i < 81 * 81; ++i) {
Point a{i % 9, i / 9 % 9}, b{i / 81 % 9, i / 729};
++tested;
disagree += liangBarsky(a, b, small).visible != meetsRect(a, b, small);
}
std::println("{} segments checked, {} disagreements", tested, disagree);
}Press Run to compile and run this example.
inside t in [0, 1] (2, 2) to (8, 4)
crosses two sides t in [2/9, 7/9] (0, 1.88889) to (10, 4.11111)
misses the window rejected
touches a corner t in [1/2, 1/2] (10, 6) to (10, 6)
runs along an edge t in [3/7, 1] (0, 6) to (4, 6)
a single point t in [0, 1] (5, 5) to (5, 5)
6561 segments checked, 0 disagreementsLook at the boundary cases. A segment through the corner (10, 6) has t0 = t1 = 1/2: it touches the window in a single point, and the comparison t0 <= t1 keeps it. A segment lying on the top edge is clipped to the part that overlaps the edge. A segment of length zero is visible because its one point is inside. Whether touching counts as visible is a policy; this code says yes, and the 6,561-segment check agrees with the same policy written another way. Changing !before(t1, t0) to before(t0, t1), a plausible "strictly visible" edit, makes 812 of them disagree, so the test also pins down the policy.
Sutherland-Hodgman: one edge at a time
Clipping a polygon is not the same as clipping each of its edges: the result must be a closed polygon, so where the subject leaves the window and comes back, the result needs new edges along the window's boundary. Sutherland and Hodgman (1974) handle that with a simple reduction. A convex window is the intersection of the half-planes to the left of its edges. Clip the whole subject against the first half-plane, clip the result against the second, and so on.
Clipping a polygon against one half-plane is a walk around its vertices. For each edge from prev to cur:
- Both inside: output
cur. - Leaving (prev inside, cur outside): output the crossing point.
- Entering (prev outside, cur inside): output the crossing point, then
cur. - Both outside: output nothing.
The walk connects the crossing points in order, and that connection is the new edge along the window's boundary. The inside test is the orientation predicate, exact as always. The crossing point is a construction, and the code below rounds it to the nearest integer point. Why that is a reasonable choice, and what it costs, gets its own section below. Points exactly on the window's edge count as inside, and the cases are written so that such a point is output once and never produces a second, identical crossing point.
The drawing clips a concave star against a convex hexagon. As always in this course, the helper flips the y-axis because SVG's y grows downward.
#include <format>
#include <print>
#include <string>
#include <vector>
struct Point { long long x, y; };
class Svg {
public:
Svg(double width, double height) : w_(width), h_(height) {}
void line(Point a, Point b, const char* color = "#1f2937", double width = 1.5) {
add("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"{}\" stroke-width=\"{}\"/>",
a.x, fy(a.y), b.x, fy(b.y), color, width);
}
void dot(Point p, const char* color = "#1f2937", double r = 3) {
add("<circle cx=\"{}\" cy=\"{}\" r=\"{}\" fill=\"{}\"/>", p.x, fy(p.y), r, color);
}
void polygon(const std::vector<Point>& pts, const char* stroke = "#1f2937",
const char* fill = "none") {
std::string list;
for (Point p : pts) list += std::format("{},{} ", p.x, fy(p.y));
add("<polygon points=\"{}\" stroke=\"{}\" fill=\"{}\" stroke-width=\"1.5\"/>", list, stroke, fill);
}
void text(Point p, const std::string& s, double size = 8) {
add("<text x=\"{}\" y=\"{}\" font-size=\"{}\">{}</text>", p.x, fy(p.y), size, s);
}
void print() const {
std::println("<svg viewBox=\"0 0 {} {}\">", w_, h_);
std::print("{}", body_);
std::println("</svg>");
}
private:
double fy(double y) const { return h_ - y; } // SVG's y axis points down
template <class... Args>
void add(std::format_string<Args...> fmt, Args&&... args) {
body_ += std::format(fmt, std::forward<Args>(args)...);
body_ += '\n';
}
double w_, h_;
std::string body_;
};
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
// n / d rounded to the nearest integer (halves away from zero).
long long roundDiv(__int128 n, __int128 d) {
if (d < 0) n = -n, d = -d;
__int128 q = n >= 0 ? (2 * n + d) / (2 * d) : -((-2 * n + d) / (2 * d));
return static_cast<long long>(q);
}
// Where segment ab crosses the line through c and d (a and b on opposite
// sides), rounded to the integer grid. The exact point is rational.
Point crossing(Point a, Point b, Point c, Point d) {
__int128 sa = cross(c, d, a), sb = cross(c, d, b);
__int128 den = sa - sb; // non-zero: the signs differ
return {roundDiv(a.x * den + (b.x - a.x) * sa, den),
roundDiv(a.y * den + (b.y - a.y) * sa, den)};
}
// Sutherland-Hodgman: clip a polygon against each edge of a convex,
// counter-clockwise window in turn. Points on an edge count as inside.
std::vector<Point> clip(const std::vector<Point>& subject,
const std::vector<Point>& window) {
std::vector<Point> out = subject;
for (size_t i = 0; i < window.size() && !out.empty(); ++i) {
Point c = window[i], d = window[(i + 1) % window.size()];
std::vector<Point> in = std::move(out);
out.clear();
for (size_t j = 0; j < in.size(); ++j) {
Point prev = in[(j + in.size() - 1) % in.size()], cur = in[j];
int sp = orientation(c, d, prev), sc = orientation(c, d, cur);
if (sc >= 0) {
if (sp < 0 && sc > 0) out.push_back(crossing(prev, cur, c, d));
out.push_back(cur);
} else if (sp > 0) {
out.push_back(crossing(prev, cur, c, d));
}
}
}
return out;
}
int main() {
// Subject: a concave, counter-clockwise star. Window: a convex hexagon.
const std::vector<Point> star{{150, 190}, {125, 130}, {40, 140}, {105, 95},
{60, 15}, {150, 70}, {240, 15}, {195, 95},
{260, 140}, {175, 130}};
const std::vector<Point> window{{90, 40}, {210, 40}, {250, 100}, {210, 160},
{90, 160}, {50, 100}};
std::vector<Point> result = clip(star, window);
Svg svg(300, 200);
svg.polygon(window, "#2563eb", "rgba(37,99,235,0.15)");
svg.polygon(star, "#1f2937");
svg.polygon(result, "#dc2626", "rgba(220,38,38,0.25)");
for (Point p : result) svg.dot(p, "#dc2626", 2.5);
svg.text({8, 188}, std::format("result: {} vertices", result.size()), 10);
svg.print();
}Press Run to compile and run this example.
The result, in red, is the star with its five tips cut off by the hexagon. Each tip contributes two crossing points, joined by a piece of the hexagon's boundary; at the two lower tips that piece includes a corner of the hexagon. Ten crossing points, five inner vertices of the star and two hexagon corners make the 17 vertices. Each half-plane step is linear in the size of its input, so a subject with n vertices and a window with m edges costs O(m * k), where k is the largest intermediate polygon. For a convex subject k <= n + m.
The artifact of a concave subject
Sutherland-Hodgman always outputs one polygon. That is correct when the subject is convex, because the intersection of two convex polygons is a single convex polygon. When the subject is concave, the true intersection can fall apart into several pieces, and the algorithm joins them with edges that run along the window's boundary there and back: zero-width bridges.
The program clips a U shape with a window that keeps only the tops of its two arms. The correct answer is two separate rectangles of area 150 each.
#include <algorithm>
#include <print>
#include <vector>
struct Point { long long x, y; };
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
// n / d rounded to the nearest integer (halves away from zero).
long long roundDiv(__int128 n, __int128 d) {
if (d < 0) n = -n, d = -d;
__int128 q = n >= 0 ? (2 * n + d) / (2 * d) : -((-2 * n + d) / (2 * d));
return static_cast<long long>(q);
}
// Where segment ab crosses the line through c and d (a and b on opposite
// sides), rounded to the integer grid. The exact point is rational.
Point crossing(Point a, Point b, Point c, Point d) {
__int128 sa = cross(c, d, a), sb = cross(c, d, b);
__int128 den = sa - sb; // non-zero: the signs differ
return {roundDiv(a.x * den + (b.x - a.x) * sa, den),
roundDiv(a.y * den + (b.y - a.y) * sa, den)};
}
// Sutherland-Hodgman: clip a polygon against each edge of a convex,
// counter-clockwise window in turn. Points on an edge count as inside.
std::vector<Point> clip(const std::vector<Point>& subject,
const std::vector<Point>& window) {
std::vector<Point> out = subject;
for (size_t i = 0; i < window.size() && !out.empty(); ++i) {
Point c = window[i], d = window[(i + 1) % window.size()];
std::vector<Point> in = std::move(out);
out.clear();
for (size_t j = 0; j < in.size(); ++j) {
Point prev = in[(j + in.size() - 1) % in.size()], cur = in[j];
int sp = orientation(c, d, prev), sc = orientation(c, d, cur);
if (sc >= 0) {
if (sp < 0 && sc > 0) out.push_back(crossing(prev, cur, c, d));
out.push_back(cur);
} else if (sp > 0) {
out.push_back(crossing(prev, cur, c, d));
}
}
}
return out;
}
long long twiceArea(const std::vector<Point>& poly) { // shoelace, exact
long long sum = 0;
for (size_t i = 0; i < poly.size(); ++i) {
Point a = poly[i], b = poly[(i + 1) % poly.size()];
sum += a.x * b.y - a.y * b.x;
}
return sum;
}
// Pairs of edges that lie on one line, point in opposite directions and
// overlap in more than a point: a zero-width "bridge" in the output.
int bridges(const std::vector<Point>& poly) {
int found = 0;
size_t n = poly.size();
for (size_t i = 0; i < n; ++i)
for (size_t j = i + 1; j < n; ++j) {
Point a = poly[i], b = poly[(i + 1) % n], c = poly[j], d = poly[(j + 1) % n];
if (orientation(a, b, c) != 0 || orientation(a, b, d) != 0) continue;
long long dx = b.x - a.x, dy = b.y - a.y;
auto at = [&](Point p) { return (p.x - a.x) * dx + (p.y - a.y) * dy; };
long long lo = std::max(0LL, std::min(at(c), at(d)));
long long hi = std::min(dx * dx + dy * dy, std::max(at(c), at(d)));
bool opposite = (d.x - c.x) * dx + (d.y - c.y) * dy < 0;
if (opposite && lo < hi) {
std::println(" bridge: ({}, {})->({}, {}) and ({}, {})->({}, {})",
a.x, a.y, b.x, b.y, c.x, c.y, d.x, d.y);
++found;
}
}
return found;
}
void report(const char* name, const std::vector<Point>& poly) {
std::print("{}: {} vertices:", name, poly.size());
for (Point p : poly) std::print(" ({}, {})", p.x, p.y);
std::println("\n twice the area: {}", twiceArea(poly));
}
int main() {
// A U shape, counter-clockwise: two arms joined by a base along y = 0..10.
const std::vector<Point> u{{0, 0}, {30, 0}, {30, 30}, {20, 30},
{20, 10}, {10, 10}, {10, 30}, {0, 30}};
// A window that keeps only the top of the arms.
const std::vector<Point> window{{-5, 15}, {35, 15}, {35, 40}, {-5, 40}};
std::vector<Point> result = clip(u, window);
report("result", result);
std::println("{} bridge(s) found", bridges(result));
}Press Run to compile and run this example.
result: 8 vertices: (0, 15) (30, 15) (30, 30) (20, 30) (20, 15) (10, 15) (10, 30) (0, 30)
twice the area: 600
bridge: (0, 15)->(30, 15) and (20, 15)->(10, 15)
1 bridge(s) foundThe area is right: 600 is twice 300. But the output is one polygon in which the edge from (0,15) to (30,15) lies on top of the edge from (20,15) to (10,15), running the opposite way. Filled, it looks right; outlined, it shows a line between the arms; passed to code that assumes a simple polygon, such as a triangulator, it can misbehave, because the polygon touches itself along a segment.
The bridge detector above is one way to notice the problem: look for pairs of collinear edges that run in opposite directions and overlap in more than a point. To remove it, split the polygon at the overlapping parts into separate rings. In practice, if the subject can be concave and you need separate pieces, use an algorithm that produces them, which is the subject of the last part of this lesson.
Testing it against a different algorithm
To check Sutherland-Hodgman for convex polygons, compute the same intersection in an independent way. The intersection of two convex polygons is the convex hull of three sets of points: the vertices of each polygon that lie inside or on the other, and the points where their edges cross. The reference below computes that hull with Andrew's monotone chain and takes its area.
Because the crossing points are constructions, the two areas cannot be compared exactly. They can be compared with a tolerance derived from the representation: every rounded vertex is within about 0.71 units of the true point, so the area can change by at most a small multiple of the perimeter. The test allows twice the perimeter. Shapes are drawn from a small grid and then scaled by 10^6, which keeps the coincidences of the small grid (shared edges, vertices on edges, identical polygons) while making the rounding tiny in comparison.
#include <algorithm>
#include <cmath>
#include <print>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
// n / d rounded to the nearest integer (halves away from zero).
long long roundDiv(__int128 n, __int128 d) {
if (d < 0) n = -n, d = -d;
__int128 q = n >= 0 ? (2 * n + d) / (2 * d) : -((-2 * n + d) / (2 * d));
return static_cast<long long>(q);
}
// Where segment ab crosses the line through c and d (a and b on opposite
// sides), rounded to the integer grid. The exact point is rational.
Point crossing(Point a, Point b, Point c, Point d) {
__int128 sa = cross(c, d, a), sb = cross(c, d, b);
__int128 den = sa - sb; // non-zero: the signs differ
return {roundDiv(a.x * den + (b.x - a.x) * sa, den),
roundDiv(a.y * den + (b.y - a.y) * sa, den)};
}
// Sutherland-Hodgman: clip a polygon against each edge of a convex,
// counter-clockwise window in turn. Points on an edge count as inside.
std::vector<Point> clip(const std::vector<Point>& subject,
const std::vector<Point>& window) {
std::vector<Point> out = subject;
for (size_t i = 0; i < window.size() && !out.empty(); ++i) {
Point c = window[i], d = window[(i + 1) % window.size()];
std::vector<Point> in = std::move(out);
out.clear();
for (size_t j = 0; j < in.size(); ++j) {
Point prev = in[(j + in.size() - 1) % in.size()], cur = in[j];
int sp = orientation(c, d, prev), sc = orientation(c, d, cur);
if (sc >= 0) {
if (sp < 0 && sc > 0) out.push_back(crossing(prev, cur, c, d));
out.push_back(cur);
} else if (sp > 0) {
out.push_back(crossing(prev, cur, c, d));
}
}
}
return out;
}
// Convex hull (Andrew's monotone chain), counter-clockwise, no collinear points.
std::vector<Point> hull(std::vector<Point> p) {
std::ranges::sort(p, [](Point a, Point b) {
return a.x != b.x ? a.x < b.x : a.y < b.y;
});
std::vector<Point> h(2 * p.size());
size_t k = 0;
for (size_t i = 0; i < p.size(); ++i) {
while (k >= 2 && orientation(h[k - 2], h[k - 1], p[i]) <= 0) --k;
h[k++] = p[i];
}
for (size_t i = p.size() - 1, lower = k + 1; i-- > 0;) {
while (k >= lower && orientation(h[k - 2], h[k - 1], p[i]) <= 0) --k;
h[k++] = p[i];
}
h.resize(k > 1 ? k - 1 : k);
return h;
}
bool insideOrOn(Point p, const std::vector<Point>& convex) {
for (size_t i = 0; i < convex.size(); ++i)
if (orientation(convex[i], convex[(i + 1) % convex.size()], p) < 0) return false;
return true;
}
struct Pd { double x, y; };
double area(const std::vector<Pd>& p) {
double s = 0;
for (size_t i = 0; i < p.size(); ++i)
s += p[i].x * p[(i + 1) % p.size()].y - p[i].y * p[(i + 1) % p.size()].x;
return s / 2;
}
// Independent reference: the intersection of two convex polygons is the convex
// hull of the vertices of each inside the other and the edge crossing points.
double referenceArea(const std::vector<Point>& a, const std::vector<Point>& b) {
std::vector<Pd> pts;
for (Point p : a) if (insideOrOn(p, b)) pts.push_back({double(p.x), double(p.y)});
for (Point p : b) if (insideOrOn(p, a)) pts.push_back({double(p.x), double(p.y)});
for (size_t i = 0; i < a.size(); ++i)
for (size_t j = 0; j < b.size(); ++j) {
Point p = a[i], q = a[(i + 1) % a.size()];
Point r = b[j], s = b[(j + 1) % b.size()];
long long d1 = cross(r, s, p), d2 = cross(r, s, q);
if ((d1 < 0 && d2 > 0) || (d1 > 0 && d2 < 0))
if (orientation(p, q, r) * orientation(p, q, s) < 0) {
double t = double(d1) / double(d1 - d2);
pts.push_back({p.x + t * (q.x - p.x), p.y + t * (q.y - p.y)});
}
}
if (pts.size() < 3) return 0;
// Their convex hull, by the monotone chain again, here in doubles.
std::ranges::sort(pts, [](Pd u, Pd v) { return u.x != v.x ? u.x < v.x : u.y < v.y; });
auto turn = [](Pd o, Pd u, Pd v) {
return (u.x - o.x) * (v.y - o.y) - (u.y - o.y) * (v.x - o.x);
};
std::vector<Pd> h(2 * pts.size() + 1);
size_t k = 0;
for (size_t i = 0; i < pts.size(); ++i) {
while (k >= 2 && turn(h[k - 2], h[k - 1], pts[i]) <= 0) --k;
h[k++] = pts[i];
}
for (size_t i = pts.size() - 1, lower = k + 1; i-- > 0;) {
while (k >= lower && turn(h[k - 2], h[k - 1], pts[i]) <= 0) --k;
h[k++] = pts[i];
}
h.resize(k > 1 ? k - 1 : k);
pts = h;
return area(pts);
}
int main() {
Lcg rng{10};
const long long scale = 1'000'000; // shapes from a small grid, then scaled up
long cases = 0, failed = 0, touching = 0;
for (int t = 0; t < 3'000; ++t) {
int side = t % 2 ? 4 : 30; // a 4 by 4 grid makes shared edges common
auto randomConvex = [&] {
std::vector<Point> p;
for (int i = rng.range(3, 8); i > 0; --i)
p.push_back({rng.range(0, side) * scale, rng.range(0, side) * scale});
return hull(p);
};
std::vector<Point> a = randomConvex(), b = randomConvex();
if (a.size() < 3 || b.size() < 3) continue;
++cases;
std::vector<Point> r = clip(a, b);
std::vector<Pd> rd;
for (Point p : r) rd.push_back({double(p.x), double(p.y)});
double got = area(rd), want = referenceArea(a, b);
double perimeter = 0;
for (size_t i = 0; i < a.size(); ++i)
perimeter += std::hypot(double(a[(i + 1) % a.size()].x - a[i].x),
double(a[(i + 1) % a.size()].y - a[i].y));
// Each rounded vertex moves at most 0.71 units: allow 2 units per unit length.
failed += std::abs(got - want) > 2 * perimeter;
touching += !r.empty() && want == 0;
}
std::println("{} pairs of convex polygons, {} failed", cases, failed);
std::println("{} pairs touch without overlapping (zero-area result)", touching);
}Press Run to compile and run this example.
2896 pairs of convex polygons, 0 failed
64 pairs touch without overlapping (zero-area result)Every pair agrees, including 64 pairs that only touch along an edge or at a corner, where the true intersection has no area. The test catches real mistakes. Deleting the line that outputs the crossing point when the subject enters the half-plane makes 2,491 pairs fail. Treating a point on the window's edge as outside when it is the current vertex but as inside when it is the previous one, an inconsistency that is easy to write, loses vertices and makes 808 fail.
Half-planes and the kernel of a polygon
The heart of Sutherland-Hodgman, "clip against the half-plane to the left of a directed line", is useful on its own. Many problems are an intersection of half-planes: the feasible region of a two-variable linear program, the set of points from which a shape is fully visible, the region a robot can reach under several linear constraints. Start with a large square and clip it against each half-plane in turn. That is O(n^2) for n half-planes, because the intermediate polygon can have up to n edges. An O(n log n) algorithm exists (de Berg et al., chapter 4), and two convex polygons can be intersected in O(n + m) time (O'Rourke, Chien, Olson and Naddor, 1982).
There is a pleasant consequence. Our clip function never checks that the window is convex; it just clips against the half-plane of each window edge. Given a concave window, it no longer computes the intersection with the window. It computes the intersection of all the half-planes to the left of the window's edges, which is the window's kernel: the set of points that can see the whole polygon. A polygon with a non-empty kernel is called star-shaped, and a single guard anywhere in the kernel watches the whole room. The program computes the kernel of a concave room by clipping a large square against it.
#include <format>
#include <print>
#include <string>
#include <vector>
struct Point { long long x, y; };
class Svg {
public:
Svg(double width, double height) : w_(width), h_(height) {}
void line(Point a, Point b, const char* color = "#1f2937", double width = 1.5) {
add("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"{}\" stroke-width=\"{}\"/>",
a.x, fy(a.y), b.x, fy(b.y), color, width);
}
void dot(Point p, const char* color = "#1f2937", double r = 3) {
add("<circle cx=\"{}\" cy=\"{}\" r=\"{}\" fill=\"{}\"/>", p.x, fy(p.y), r, color);
}
void polygon(const std::vector<Point>& pts, const char* stroke = "#1f2937",
const char* fill = "none") {
std::string list;
for (Point p : pts) list += std::format("{},{} ", p.x, fy(p.y));
add("<polygon points=\"{}\" stroke=\"{}\" fill=\"{}\" stroke-width=\"1.5\"/>", list, stroke, fill);
}
void text(Point p, const std::string& s, double size = 8) {
add("<text x=\"{}\" y=\"{}\" font-size=\"{}\">{}</text>", p.x, fy(p.y), size, s);
}
void print() const {
std::println("<svg viewBox=\"0 0 {} {}\">", w_, h_);
std::print("{}", body_);
std::println("</svg>");
}
private:
double fy(double y) const { return h_ - y; } // SVG's y axis points down
template <class... Args>
void add(std::format_string<Args...> fmt, Args&&... args) {
body_ += std::format(fmt, std::forward<Args>(args)...);
body_ += '\n';
}
double w_, h_;
std::string body_;
};
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
// n / d rounded to the nearest integer (halves away from zero).
long long roundDiv(__int128 n, __int128 d) {
if (d < 0) n = -n, d = -d;
__int128 q = n >= 0 ? (2 * n + d) / (2 * d) : -((-2 * n + d) / (2 * d));
return static_cast<long long>(q);
}
// Where segment ab crosses the line through c and d (a and b on opposite
// sides), rounded to the integer grid. The exact point is rational.
Point crossing(Point a, Point b, Point c, Point d) {
__int128 sa = cross(c, d, a), sb = cross(c, d, b);
__int128 den = sa - sb; // non-zero: the signs differ
return {roundDiv(a.x * den + (b.x - a.x) * sa, den),
roundDiv(a.y * den + (b.y - a.y) * sa, den)};
}
// Sutherland-Hodgman: clip a polygon against each edge of a convex,
// counter-clockwise window in turn. Points on an edge count as inside.
std::vector<Point> clip(const std::vector<Point>& subject,
const std::vector<Point>& window) {
std::vector<Point> out = subject;
for (size_t i = 0; i < window.size() && !out.empty(); ++i) {
Point c = window[i], d = window[(i + 1) % window.size()];
std::vector<Point> in = std::move(out);
out.clear();
for (size_t j = 0; j < in.size(); ++j) {
Point prev = in[(j + in.size() - 1) % in.size()], cur = in[j];
int sp = orientation(c, d, prev), sc = orientation(c, d, cur);
if (sc >= 0) {
if (sp < 0 && sc > 0) out.push_back(crossing(prev, cur, c, d));
out.push_back(cur);
} else if (sp > 0) {
out.push_back(crossing(prev, cur, c, d));
}
}
}
return out;
}
long long twiceArea(const std::vector<Point>& poly) {
long long sum = 0;
for (size_t i = 0; i < poly.size(); ++i) {
Point a = poly[i], b = poly[(i + 1) % poly.size()];
sum += a.x * b.y - a.y * b.x;
}
return sum;
}
int main() {
// A room shaped like a concave, counter-clockwise polygon.
const std::vector<Point> room{{40, 20}, {150, 55}, {260, 20}, {215, 100},
{280, 180}, {160, 140}, {60, 185}, {95, 100}};
// Clipping a large square by every edge of the room intersects the half-planes
// to the left of all edges: the room's kernel.
const std::vector<Point> big{{0, 0}, {300, 0}, {300, 200}, {0, 200}};
std::vector<Point> kernel = clip(big, room);
Svg svg(300, 200);
svg.polygon(room, "#1f2937", "none");
if (!kernel.empty()) svg.polygon(kernel, "#2563eb", "rgba(37,99,235,0.15)");
svg.text({8, 8}, std::format("kernel: {} vertices, area {}", kernel.size(),
double(twiceArea(kernel)) / 2), 10);
svg.print();
}Press Run to compile and run this example.
The kernel is the shaded octagon in the middle, bounded by the extensions of the room's edges near the four reflex corners. Move one of those corners inward and the kernel shrinks; push it far enough and the kernel disappears, and no single guard can see everything. Computing the kernel this way takes O(n^2); Lee and Preparata (1979) showed it can be done in O(n).
The same observation is a warning. If you pass a concave window to clip by mistake, nothing fails: you get a smaller, convex result that is wrong. A clip function that accepts only convex windows should check that its window is convex, with the orientation test at each vertex, before it runs.
Representing constructed points
Every crossing point is rational. With integer input bounded by 10^9, the crossing of segment ab with the line cd is a + (b - a) * sa / (sa - sb), where sa and sb are the cross products from the orientation test, each up to about 8 * 10^18. The numerator of each coordinate needs about 95 bits and the denominator about 64, which fits in __int128 (a GCC and Clang extension, not standard C++, as discussed in the robustness lesson). So one clip can be exact. The next clip has to compute cross products of these fractions, and clearing their denominators multiplies the number of bits by two or three, more than any built-in integer holds. Repeated exact clipping needs arbitrary-precision rationals, which is what CGAL's exact-construction kernels provide, at a large cost in speed and memory.
The other choice, used throughout this lesson, is to round each crossing point to the nearest integer point. Every predicate stays exact, because it works on integer points again, and the error of each construction is bounded: at most half a unit in each coordinate. Scale the input so that one unit is negligible (a micrometer instead of a meter) and the error is too. The price is that the result is no longer exactly the intersection, and the program measures what that means.
#include <algorithm>
#include <cmath>
#include <print>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
// n / d rounded to the nearest integer (halves away from zero).
long long roundDiv(__int128 n, __int128 d) {
if (d < 0) n = -n, d = -d;
__int128 q = n >= 0 ? (2 * n + d) / (2 * d) : -((-2 * n + d) / (2 * d));
return static_cast<long long>(q);
}
// Where segment ab crosses the line through c and d (a and b on opposite
// sides), rounded to the integer grid. The exact point is rational.
Point crossing(Point a, Point b, Point c, Point d) {
__int128 sa = cross(c, d, a), sb = cross(c, d, b);
__int128 den = sa - sb; // non-zero: the signs differ
return {roundDiv(a.x * den + (b.x - a.x) * sa, den),
roundDiv(a.y * den + (b.y - a.y) * sa, den)};
}
// Sutherland-Hodgman: clip a polygon against each edge of a convex,
// counter-clockwise window in turn. Points on an edge count as inside.
std::vector<Point> clip(const std::vector<Point>& subject,
const std::vector<Point>& window) {
std::vector<Point> out = subject;
for (size_t i = 0; i < window.size() && !out.empty(); ++i) {
Point c = window[i], d = window[(i + 1) % window.size()];
std::vector<Point> in = std::move(out);
out.clear();
for (size_t j = 0; j < in.size(); ++j) {
Point prev = in[(j + in.size() - 1) % in.size()], cur = in[j];
int sp = orientation(c, d, prev), sc = orientation(c, d, cur);
if (sc >= 0) {
if (sp < 0 && sc > 0) out.push_back(crossing(prev, cur, c, d));
out.push_back(cur);
} else if (sp > 0) {
out.push_back(crossing(prev, cur, c, d));
}
}
}
return out;
}
// Convex hull (Andrew's monotone chain), counter-clockwise, no collinear points.
std::vector<Point> hull(std::vector<Point> p) {
std::ranges::sort(p, [](Point a, Point b) {
return a.x != b.x ? a.x < b.x : a.y < b.y;
});
std::vector<Point> h(2 * p.size());
size_t k = 0;
for (size_t i = 0; i < p.size(); ++i) {
while (k >= 2 && orientation(h[k - 2], h[k - 1], p[i]) <= 0) --k;
h[k++] = p[i];
}
for (size_t i = p.size() - 1, lower = k + 1; i-- > 0;) {
while (k >= lower && orientation(h[k - 2], h[k - 1], p[i]) <= 0) --k;
h[k++] = p[i];
}
h.resize(k > 1 ? k - 1 : k);
return h;
}
bool insideOrOn(Point p, const std::vector<Point>& convex) {
for (size_t i = 0; i < convex.size(); ++i)
if (orientation(convex[i], convex[(i + 1) % convex.size()], p) < 0) return false;
return true;
}
int bits(__int128 v) { // bits needed for |v|
int b = 0;
for (v = v < 0 ? -v : v; v > 0; v >>= 1) ++b;
return b;
}
int main() {
// 1. One construction, exactly: segment ab crosses the line cd at a rational point.
Point a{0, 0}, b{999'999, 333'334}, c{700'001, -5}, d{299'999, 1'000'000};
__int128 sa = cross(c, d, a), sb = cross(c, d, b), den = sa - sb;
__int128 nx = a.x * den + (b.x - a.x) * sa;
Point r = crossing(a, b, c, d);
std::println("exact x = n / d with n of {} bits and d of {} bits; x is about {:.4f}",
bits(nx), bits(den), double(nx) / double(den));
std::println("rounded to the grid: ({}, {})", r.x, r.y);
// 2. Rounded constructions, repeated: clip random convex polygons twice.
Lcg rng{3};
long results = 0, escaped = 0, changed = 0;
double worst = 0; // largest distance of a result vertex outside the window
for (int t = 0; t < 2'000; ++t) {
std::vector<Point> subject, window;
for (int i = 0; i < 8; ++i) {
subject.push_back({rng.range(0, 1'000'000), rng.range(0, 1'000'000)});
window.push_back({rng.range(0, 1'000'000), rng.range(0, 1'000'000)});
}
subject = hull(subject);
window = hull(window);
std::vector<Point> once = clip(subject, window);
if (once.empty()) continue;
++results;
bool out = false;
for (Point p : once)
for (size_t i = 0; i < window.size(); ++i) {
Point e0 = window[i], e1 = window[(i + 1) % window.size()];
if (orientation(e0, e1, p) >= 0) continue;
out = true;
double len = std::hypot(double(e1.x - e0.x), double(e1.y - e0.y));
worst = std::max(worst, -double(cross(e0, e1, p)) / len);
}
escaped += out;
std::vector<Point> twice = clip(once, window);
changed += !std::ranges::equal(once, twice, [](Point p, Point q) {
return p.x == q.x && p.y == q.y;
});
}
std::println("{} results: {} have a vertex outside the window, by at most {:.2f}",
results, escaped, worst);
std::println("and {} change when clipped by the same window again", changed);
}Press Run to compile and run this example.
exact x = n / d with n of 60 bits and d of 41 bits; x is about 617645.9585
rounded to the grid: (617646, 205883)
1998 results: 1817 have a vertex outside the window, by at most 0.99
and 1370 change when clipped by the same window againThe first two lines show one crossing exactly and rounded: the exact x-coordinate is a fraction with a 60-bit numerator, and rounding moves the point by less than one unit. The rest of the program clips random convex polygons with random convex windows. In 1,817 of 1,998 results, at least one vertex lies outside the window, by less than one unit. And because those vertices are outside, clipping the result again by the same window changes it in 1,370 cases. A mathematician's clip is idempotent; a rounded one is not.
That is acceptable when you know it and plan for it: test containment with a tolerance of one unit, never clip the same output twice, and do not expect two polygons that should share an edge to share it exactly after rounding. Clipper2 takes exactly this approach. Its documentation says it works on integer coordinates internally, shows intersection points moved by rounding, removes resulting near-zero-area triangles, and warns that tiny self-intersections can remain. A more careful method, snap rounding, moves every segment that passes through a rounded point's grid cell as well, so that rounding cannot change the topology (Hobby, 1999).