Segment Intersection and the Sweep Line
Decide exactly whether two segments meet, in every degenerate case, then find whether any of n segments intersect in O(n log n) with the Shamos-Hoey sweep.
Do any of these segments cross?
A polygon is simple if its edges meet only at shared corners. A road network imported from two map layers should have a junction wherever two roads cross. A circuit board fails its design check if two traces on the same layer touch. A path planner must know whether a straight move crosses a wall. Each of these is the same question about a set of line segments: does any pair intersect, and if so, which?
For two segments the answer comes from four orientation tests, and the only difficulty is the degenerate cases: segments that touch at an endpoint, lie on the same line, or overlap. Get one of those wrong and a polygon with a vertex resting on another edge passes as simple.
For n segments, testing every pair costs n(n-1)/2 tests: about 50 million for 10,000 segments. The sweep line idea, due to Shamos and Hoey, answers the yes-or-no question in O(n log n). It is also the first appearance in this course of a pattern that comes back in later lessons: keep a sorted structure of the objects that cross an imaginary moving line, and only ever compare neighbors.
This lesson uses the integer Point and exact orientation from the previous lesson, with coordinates bounded by |x|, |y| <= 10^9.
Two segments cross when each straddles the other
Take segments s from s.a to s.b and t from t.a to t.b. In the general case, where no three of the four endpoints are collinear, they cross exactly when:
s.aands.blie on opposite sides of the line throught, andt.aandt.blie on opposite sides of the line throughs.
Each condition takes two orientation tests: opposite sides means the two orientations have opposite signs, so their product is negative.
A worked example. Let s run from (0, 0) to (4, 4) and t from (0, 4) to (4, 0), an X. Walking along t, the direction is (4, -4). For s.a, the vector from t.a is (0, -4), and the cross product is 4·(-4) - (-4)·0 = -16: right. For s.b, the vector is (4, 0) and the cross product is 4·0 - (-4)·4 = 16: left. Opposite signs, so s straddles the line through t. The same computation the other way round shows that t straddles s, and the segments cross.
Both conditions are needed. Let s run from (0, 0) to (4, 0) and t from (5, -1) to (5, 1). The endpoints of t are on opposite sides of the line through s (one below, one above), but both endpoints of s are to the left of the vertical line through t. The lines cross at (5, 0), which is beyond the end of s, so the segments do not.
Every way two segments can meet
When an orientation is zero, an endpoint lies on the other segment's line, and the straddle rule says nothing. Each zero needs a second question: is that endpoint actually within the other segment? Because the point is already known to be on the line, checking that it lies in the segment's bounding box is enough.
If all four orientations are zero, the segments lie on one line, and the question becomes one-dimensional: do the two intervals along that line overlap? Lexicographic order (by x, then by y) is a valid order along any line, including a vertical one, so the program sorts each segment's endpoints that way and compares the later start with the earlier end.
| Situation | Orientations | Result |
|---|---|---|
| X: proper crossing | d1·d2 < 0 and d3·d4 < 0 | crossing |
| T, or a shared endpoint | one is 0, that endpoint in the other's box | touching |
| Same line, overlapping | all 0, intervals share more than a point | overlapping |
| Same line, end to end | all 0, intervals share one point | touching |
| Same line, apart | all 0, intervals disjoint | disjoint |
| Anything else | disjoint |
The program classifies one segment against a set that covers every row of the table. Both segments must have distinct endpoints; a segment of length zero has no line through it.
#include <algorithm>
#include <print>
#include <utility>
struct Point { long long x, y; };
struct Segment { Point a, b; };
long long cross(Point a, Point b, Point c) { // twice the signed area of abc, exact
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) { // -1 right turn, 0 collinear, +1 left
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }
// For p collinear with a and b: is p within the segment's bounding box?
bool inBox(Point a, Point b, Point p) {
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);
}
enum class Contact { Disjoint, Crossing, Touching, Overlapping };
Contact classify(Segment s, Segment t) {
int d1 = orientation(t.a, t.b, s.a), d2 = orientation(t.a, t.b, s.b);
int d3 = orientation(s.a, s.b, t.a), d4 = orientation(s.a, s.b, t.b);
if (d1 * d2 < 0 && d3 * d4 < 0) return Contact::Crossing; // the general case
if (d1 == 0 && d2 == 0) { // all four points on one line: compare along it
if (lexLess(s.b, s.a)) std::swap(s.a, s.b);
if (lexLess(t.b, t.a)) std::swap(t.a, t.b);
Point lo = lexLess(s.a, t.a) ? t.a : s.a; // the later start
Point hi = lexLess(s.b, t.b) ? s.b : t.b; // the earlier end
if (lexLess(hi, lo)) return Contact::Disjoint;
return lexLess(lo, hi) ? Contact::Overlapping : Contact::Touching;
}
// One endpoint on the other segment, without the lines being the same.
if ((d1 == 0 && inBox(t.a, t.b, s.a)) || (d2 == 0 && inBox(t.a, t.b, s.b)) ||
(d3 == 0 && inBox(s.a, s.b, t.a)) || (d4 == 0 && inBox(s.a, s.b, t.b)))
return Contact::Touching;
return Contact::Disjoint;
}
const char* name(Contact c) {
switch (c) {
case Contact::Disjoint: return "disjoint";
case Contact::Crossing: return "crossing";
case Contact::Touching: return "touching";
case Contact::Overlapping: return "overlapping";
}
return "?";
}
int main() {
const std::pair<const char*, Segment> others[] = {
{"X shape ", {{2, -2}, {2, 2}}},
{"T shape ", {{2, 0}, {2, 3}}},
{"shared endpoint ", {{4, 0}, {6, 2}}},
{"near miss ", {{5, -1}, {5, 1}}},
{"same line, apart", {{5, 0}, {7, 0}}},
{"end to end ", {{4, 0}, {7, 0}}},
{"overlap ", {{3, 0}, {7, 0}}},
{"parallel ", {{0, 1}, {4, 1}}},
};
const Segment base{{0, 0}, {4, 0}};
for (auto [label, t] : others)
std::println("{} ({},{})-({},{}): {}", label, t.a.x, t.a.y, t.b.x, t.b.y,
name(classify(base, t)));
}Press Run to compile and run this example.
X shape (2,-2)-(2,2): crossing
T shape (2,0)-(2,3): touching
shared endpoint (4,0)-(6,2): touching
near miss (5,-1)-(5,1): disjoint
same line, apart (5,0)-(7,0): disjoint
end to end (4,0)-(7,0): touching
overlap (3,0)-(7,0): overlapping
parallel (0,1)-(4,1): disjointThe near miss is the second worked example. The two collinear cases at the bottom are why a separate one-dimensional test exists: (5,0)-(7,0) and (4,0)-(7,0) give the same four zero orientations against the base segment, and only the interval comparison tells "apart" from "end to end".
Testing the test on every small case
A function with this many branches needs a test that does not share its logic. A good one comes from a different formulation of the same question. Write the points of the two segments as s.a + λ·r and t.a + μ·q, where r and q are the segments' direction vectors and λ and μ run from 0 to 1. Setting them equal and solving with Cramer's rule gives λ = lamNum / den and μ = muNum / den, where den is the cross product of the two directions. The segments share exactly one point when both parameters are in [0, 1]. If den is zero the directions are parallel, and the collinear case becomes an interval test on positions along r. All of this stays in exact integer arithmetic: compare lamNum with 0 and den instead of dividing.
On a 5 by 5 grid there are 300 segments with distinct endpoints, and 90,000 ordered pairs of them. That is small enough to try every one, and it contains every degenerate configuration many times over.
#include <algorithm>
#include <print>
#include <utility>
#include <vector>
struct Point { long long x, y; };
struct Segment { Point a, b; };
long long cross(Point a, Point b, Point c) { // twice the signed area of abc, exact
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) { // -1 right turn, 0 collinear, +1 left
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }
// For p collinear with a and b: is p within the segment's bounding box?
bool inBox(Point a, Point b, Point p) {
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);
}
enum class Contact { Disjoint, Crossing, Touching, Overlapping };
Contact classify(Segment s, Segment t) {
int d1 = orientation(t.a, t.b, s.a), d2 = orientation(t.a, t.b, s.b);
int d3 = orientation(s.a, s.b, t.a), d4 = orientation(s.a, s.b, t.b);
if (d1 * d2 < 0 && d3 * d4 < 0) return Contact::Crossing; // the general case
if (d1 == 0 && d2 == 0) { // all four points on one line: compare along it
if (lexLess(s.b, s.a)) std::swap(s.a, s.b);
if (lexLess(t.b, t.a)) std::swap(t.a, t.b);
Point lo = lexLess(s.a, t.a) ? t.a : s.a; // the later start
Point hi = lexLess(s.b, t.b) ? s.b : t.b; // the earlier end
if (lexLess(hi, lo)) return Contact::Disjoint;
return lexLess(lo, hi) ? Contact::Overlapping : Contact::Touching;
}
// One endpoint on the other segment, without the lines being the same.
if ((d1 == 0 && inBox(t.a, t.b, s.a)) || (d2 == 0 && inBox(t.a, t.b, s.b)) ||
(d3 == 0 && inBox(s.a, s.b, t.a)) || (d4 == 0 && inBox(s.a, s.b, t.b)))
return Contact::Touching;
return Contact::Disjoint;
}
// An independent definition: solve s.a + lam*r = g.a + mu*q with Cramer's rule,
// keeping lam = lamNum/den and mu = muNum/den as exact fractions.
Contact parametric(Segment s, Segment g) {
Point r{s.b.x - s.a.x, s.b.y - s.a.y}, q{g.b.x - g.a.x, g.b.y - g.a.y};
Point w{g.a.x - s.a.x, g.a.y - s.a.y};
long long den = r.x * q.y - r.y * q.x;
long long lamNum = w.x * q.y - w.y * q.x, muNum = w.x * r.y - w.y * r.x;
if (den != 0) {
if (den < 0) { den = -den; lamNum = -lamNum; muNum = -muNum; }
if (lamNum < 0 || lamNum > den || muNum < 0 || muNum > den)
return Contact::Disjoint;
bool inside = 0 < lamNum && lamNum < den && 0 < muNum && muNum < den;
return inside ? Contact::Crossing : Contact::Touching;
}
if (muNum != 0) return Contact::Disjoint; // parallel, on different lines
// Same line: positions of g's endpoints along r, in units of 1/(r.r).
long long rr = r.x * r.x + r.y * r.y;
long long p0 = w.x * r.x + w.y * r.y;
long long p1 = (g.b.x - s.a.x) * r.x + (g.b.y - s.a.y) * r.y;
long long lo = std::max(std::min(p0, p1), 0LL), hi = std::min(std::max(p0, p1), rr);
if (lo > hi) return Contact::Disjoint;
return lo == hi ? Contact::Touching : Contact::Overlapping;
}
int main() {
std::vector<Segment> all; // every segment with distinct endpoints on a 5 by 5 grid
for (int a = 0; a < 25; ++a)
for (int b = a + 1; b < 25; ++b)
all.push_back({{a % 5, a / 5}, {b % 5, b / 5}});
int count[4] = {}, mismatches = 0;
for (Segment s : all)
for (Segment t : all) {
Contact c = classify(s, t);
++count[static_cast<int>(c)];
mismatches += c != parametric(s, t);
}
std::println("segments: {}, ordered pairs: {}", all.size(), all.size() * all.size());
std::println("disjoint {}, crossing {}, touching {}, overlapping {}",
count[0], count[1], count[2], count[3]);
std::println("disagreements with the parametric definition: {}", mismatches);
}Press Run to compile and run this example.
segments: 300, ordered pairs: 90000
disjoint 54672, crossing 14796, touching 19368, overlapping 1164
disagreements with the parametric definition: 0More than a fifth of the pairs are "touching", which random points with large coordinates would almost never produce. A small grid is the cheapest way to reach the cases that matter. Try it yourself, exercise 1, shows how quickly this test catches a missing branch.
Where they cross is a construction, not a predicate
Deciding that two segments cross is a predicate, and it is exact. Computing where they cross is a construction (lesson 2's term for an operation that produces new coordinates), and with integer input its result is rational. Solving the parametric equations gives the point s.a + r·lamNum/den. Over the common denominator den, its coordinates are (s.a.x·den + r.x·lamNum) / den and (s.a.y·den + r.y·lamNum) / den.
How many bits does that take? With coordinates up to 10^9, den and lamNum are each at most 8·10^18 in absolute value and fit in long long. The numerators reach about 2.4·10^28, roughly 94 bits, which does not fit in 64 bits but does fit in __int128. So one honest representation is a fraction with __int128 numerators. The other common choice is double, which is fine for drawing or reporting but must not be fed back into predicates, as the third example shows.
#include <optional>
#include <print>
#include <string>
#include <utility>
struct Point { long long x, y; };
struct Segment { Point a, b; };
using i128 = __int128; // GCC/Clang extension, not standard C++
std::string str(i128 v) { // std::format has no __int128 support
if (v < 0) return "-" + str(-v);
std::string s = v >= 10 ? str(v / 10) : "";
return s + char('0' + int(v % 10));
}
i128 gcd(i128 a, i128 b) { // std::gcd may reject __int128 in strict ISO mode
if (a < 0) a = -a;
if (b < 0) b = -b;
while (b != 0) a = std::exchange(b, a % b);
return a;
}
// The intersection point of the lines through s and t, as exact fractions
// x = xNum / den, y = yNum / den. Requires |coordinates| <= 10^9.
struct RationalPoint { i128 xNum, yNum, den; };
std::optional<RationalPoint> lineIntersection(Segment s, Segment t) {
long long rx = s.b.x - s.a.x, ry = s.b.y - s.a.y;
long long qx = t.b.x - t.a.x, qy = t.b.y - t.a.y;
long long den = rx * qy - ry * qx; // <= 8e18: fits in long long
if (den == 0) return std::nullopt; // parallel or the same line
long long lamNum = (t.a.x - s.a.x) * qy - (t.a.y - s.a.y) * qx;
// s.a + r * lamNum / den, over a common denominator: needs about 94 bits
i128 xNum = i128(s.a.x) * den + i128(rx) * lamNum;
i128 yNum = i128(s.a.y) * den + i128(ry) * lamNum;
if (den < 0) { xNum = -xNum; yNum = -yNum; den = -den; }
i128 g = gcd(gcd(xNum, yNum), den);
return RationalPoint{xNum / g, yNum / g, den / g};
}
int main() {
const Segment pairs[][2] = {
{{{0, 0}, {4, 4}}, {{0, 4}, {4, 0}}},
{{{0, 0}, {7, 3}}, {{0, 2}, {5, 0}}},
{{{0, 0}, {999'999'999, 1}}, {{0, 1}, {1'000'000'000, -999'999'998}}},
};
for (auto& [s, t] : pairs) {
auto p = lineIntersection(s, t);
if (!p) { std::println("parallel"); continue; }
double x = double(p->xNum) / double(p->den), y = double(p->yNum) / double(p->den);
std::println("exact: ({}/{}, {}/{})", str(p->xNum), str(p->den),
str(p->yNum), str(p->den));
std::println(" as doubles: ({:.17g}, {:.17g})", x, y);
// Is the point left of the vertical line x = 1? Ask both versions.
std::println(" x < 1? exact: {}, double: {}", p->xNum < p->den, x < 1);
}
}Press Run to compile and run this example.
exact: (2/1, 2/1)
as doubles: (2, 2)
x < 1? exact: false, double: false
exact: (70/29, 30/29)
as doubles: (2.4137931034482758, 1.0344827586206897)
x < 1? exact: false, double: false
exact: (999999999000000000/999999999000000001, 1000000000/999999999000000001)
as doubles: (1, 1.000000001e-09)
x < 1? exact: true, double: falseThe two lines in the last example cross at an x just below 1, at 1 - 1/999999999000000001. The nearest double is exactly 1. Ask "is the crossing point to the left of the vertical line x = 1?" of the exact fraction and the answer is yes; ask it of the rounded point and the answer is no. This is the predicate-versus-construction distinction from lesson 2 again: decide with the original input and exact predicates, and use constructed coordinates only for output. When an algorithm must reason about constructed points, as Bentley-Ottmann does below, it has to compare the exact fractions, which needs the extra bits.
std::format has no support for __int128, so the program prints it with a small recursive helper. It also uses its own gcd: std::gcd requires a type that the standard library recognizes as an integer, and in strict ISO mode libstdc++ does not count __int128 as one.
Sweeping a line across the plane
Imagine a vertical line moving from left to right across the segments. At any moment some segments cross it, and those segments have an order along the line, from bottom to top. Call the set of crossing segments, kept in that order, the status. The status changes only at events: when the line reaches a left endpoint, a segment joins; when it reaches a right endpoint, a segment leaves.
The key observation: if two segments intersect, then just before the line reaches their leftmost intersection point, they are adjacent in the status. Nothing can be between them, because anything between them would have to pass through the gap, which shrinks to the intersection point itself. So they became neighbors at some earlier event: either one of them was inserted next to the other, or a segment between them left. That gives the algorithm:
- Sort the 2n endpoints.
- At a left endpoint, insert the segment into the status and test it against its new neighbors above and below.
- At a right endpoint, test the segment's two neighbors against each other (they are about to become adjacent), then remove the segment.
- Stop at the first intersecting pair.
Each event costs O(log n) for a balanced search tree plus at most two segment tests, so the whole algorithm runs in O(n log n) time in the worst case. It is not guaranteed to report the leftmost intersection; it reports an intersection no later than the sweep reaches the leftmost one, and stops. If the answer is no, every event is processed, which is the worst case.
The algorithm is from Shamos and Hoey's 1976 paper. It decides whether any intersection exists; reporting all of them is a different problem and comes later in this lesson.
Shamos-Hoey in C++: the comparator is the hard part
The status is a std::set<int> of segment indices, and all the geometry lives in its comparator. The textbook comparator says "compare the y-coordinates of the two segments at the sweep line's current x". Implemented literally, that computes a division, depends on a global sweep position, and changes meaning as the sweep moves, which a std::set comparator is not allowed to do. The implementation below avoids all three problems.
A static comparator. As long as no two active segments intersect, their vertical order never changes while both are in the status. So the order can be decided once, at the moment the later of the two was inserted: is the later segment's first endpoint above or below the earlier segment's line? That is one orientation test, exact, and independent of where the sweep is now. If the endpoint lies exactly on the earlier segment, the two intersect and the algorithm will report it; the comparator still has to return a consistent answer, so it uses the direction of the later segment's other endpoint, and for two collinear overlapping segments it falls back to the indices.
Lexicographic events. Events are sorted by x, then by y. This handles vertical segments without special code: a vertical segment's lower endpoint is its first event and its upper endpoint its last, as if the sweep line were tilted by an infinitesimally small angle, so that among points with the same x the lower one is reached first.
Shared endpoints first. If two segments share an endpoint, they intersect, and the program reports it before sweeping by looking for equal adjacent points in the sorted event list. That removes the most delicate tie from the sweep: an insertion and a removal at the same point. (Try it yourself, exercise 2, shows what happens without it.)
Erase by iterator. The program stores the iterator returned by insert for each segment and erases with it. Erasing by value would make the set search with the comparator again.
The test generator uses tiny coordinate ranges and short segments, a third of them vertical, so that touching, collinear and shared-endpoint cases are common, and compares the sweep with brute force over every pair.
#include <algorithm>
#include <optional>
#include <print>
#include <set>
#include <utility>
#include <vector>
struct Point { long long x, y; };
struct Segment { Point a, b; };
long long cross(Point a, Point b, Point c) { // twice the signed area of abc, exact
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) { // -1 right turn, 0 collinear, +1 left
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }
bool inBox(Point a, Point b, Point p) {
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 intersects(Segment s, Segment t) { // closed segments: touching counts
int d1 = orientation(t.a, t.b, s.a), d2 = orientation(t.a, t.b, s.b);
int d3 = orientation(s.a, s.b, t.a), d4 = orientation(s.a, s.b, t.b);
if (d1 * d2 < 0 && d3 * d4 < 0) return true;
return (d1 == 0 && inBox(t.a, t.b, s.a)) || (d2 == 0 && inBox(t.a, t.b, s.b)) ||
(d3 == 0 && inBox(s.a, s.b, t.a)) || (d4 == 0 && inBox(s.a, s.b, t.b));
}
using Pair = std::optional<std::pair<int, int>>; // indices of an intersecting pair
// Shamos-Hoey: does any pair intersect? Segments must have distinct endpoints.
Pair anyIntersection(std::vector<Segment> s) {
const int n = static_cast<int>(s.size());
for (Segment& g : s)
if (lexLess(g.b, g.a)) std::swap(g.a, g.b); // a is the first end swept
struct Event { Point p; int seg; bool isEnd; };
std::vector<Event> events;
for (int i = 0; i < n; ++i) {
events.push_back({s[i].a, i, false});
events.push_back({s[i].b, i, true});
}
std::ranges::sort(events, lexLess, &Event::p);
for (size_t k = 0; k + 1 < events.size(); ++k) // two segments share an endpoint
if (events[k].p == events[k + 1].p)
return std::pair{events[k].seg, events[k + 1].seg};
// Is segment i below segment j? Decided where the later of the two starts.
auto below = [&](int i, int j) {
if (i == j) return false;
bool iFirst = lexLess(s[i].a, s[j].a);
const Segment& u = iFirst ? s[i] : s[j]; // u started first
const Segment& v = iFirst ? s[j] : s[i];
int side = orientation(u.a, u.b, v.a); // +1: v starts above u
if (side == 0) side = orientation(u.a, u.b, v.b); // v starts on u: use its end
if (side == 0) return i < j; // collinear overlap: any consistent order
return iFirst ? side > 0 : side < 0;
};
std::set<int, decltype(below)> status(below);
std::vector<decltype(status)::iterator> where(n);
for (const Event& e : events) {
if (!e.isEnd) {
auto it = status.insert(e.seg).first;
where[e.seg] = it;
if (it != status.begin() && intersects(s[*std::prev(it)], s[e.seg]))
return std::pair{*std::prev(it), e.seg};
if (std::next(it) != status.end() && intersects(s[*std::next(it)], s[e.seg]))
return std::pair{*std::next(it), e.seg};
} else {
auto it = where[e.seg]; // erase by iterator: no comparisons needed
if (it != status.begin() && std::next(it) != status.end() &&
intersects(s[*std::prev(it)], s[*std::next(it)]))
return std::pair{*std::prev(it), *std::next(it)};
status.erase(it);
}
}
return std::nullopt;
}
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
};
Pair bruteForce(const std::vector<Segment>& s) { // every pair: O(n^2)
for (size_t i = 0; i < s.size(); ++i)
for (size_t j = i + 1; j < s.size(); ++j)
if (intersects(s[i], s[j])) return std::pair{int(i), int(j)};
return std::nullopt;
}
int main() {
int cases = 0, withIntersection = 0, failures = 0;
for (int seed = 0; seed < 20000; ++seed) {
Lcg rng{static_cast<unsigned long long>(seed)};
const int range = rng.range(3, 30), len = rng.range(1, 4), n = rng.range(2, 25);
std::vector<Segment> segs;
while (static_cast<int>(segs.size()) < n) { // short segments on a small grid
Point a{rng.range(0, range), rng.range(0, range)};
Point b{a.x + rng.range(-len, len), a.y + rng.range(-len, len)};
if (rng.range(0, 3) == 0) b.x = a.x; // plenty of vertical segments
if (!(a == b)) segs.push_back({a, b});
}
Pair slow = bruteForce(segs), fast = anyIntersection(segs);
++cases;
withIntersection += slow.has_value();
bool ok = slow.has_value() == fast.has_value() &&
(!fast || intersects(segs[fast->first], segs[fast->second]));
failures += !ok;
}
std::println("random segment sets: {}", cases);
std::println("with an intersection: {}", withIntersection);
std::println("sweep disagrees: {}", failures);
}Press Run to compile and run this example.
random segment sets: 20000
with an intersection: 13355
sweep disagrees: 0The sweep agrees with brute force on all 20,000 sets, of which about a third have no intersection at all. Those are the important ones, because the sweep has to process every event to answer no. Passing tests are not a proof, but a sweep with a broken tie rule fails this test immediately, as exercise 2 shows.
Counting the work
To see where the savings come from, count how many pairs of segments each algorithm tests with intersects. The program builds sets of n segments that do not intersect (each segment inside its own cell of a grid), so neither algorithm can stop early.
The count leaves out the comparisons that std::set makes internally while it keeps the status sorted. There are O(log n) of them per insertion, but their exact number depends on how the standard library implements its tree: the same program makes different numbers of comparator calls with libstdc++ (GCC's library) and libc++ (Clang's on macOS). A count that changes with the library cannot go in an output block, so the program counts only what our own code decides to do.
#include <algorithm>
#include <cmath>
#include <optional>
#include <print>
#include <set>
#include <utility>
#include <vector>
struct Point { long long x, y; };
struct Segment { Point a, b; };
long long cross(Point a, Point b, Point c) { // twice the signed area of abc, exact
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) { // -1 right turn, 0 collinear, +1 left
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }
bool inBox(Point a, Point b, Point p) {
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);
}
long pairTests = 0; // instrumentation: how many pairs of segments were tested?
bool intersects(Segment s, Segment t) { // closed segments: touching counts
++pairTests;
int d1 = orientation(t.a, t.b, s.a), d2 = orientation(t.a, t.b, s.b);
int d3 = orientation(s.a, s.b, t.a), d4 = orientation(s.a, s.b, t.b);
if (d1 * d2 < 0 && d3 * d4 < 0) return true;
return (d1 == 0 && inBox(t.a, t.b, s.a)) || (d2 == 0 && inBox(t.a, t.b, s.b)) ||
(d3 == 0 && inBox(s.a, s.b, t.a)) || (d4 == 0 && inBox(s.a, s.b, t.b));
}
using Pair = std::optional<std::pair<int, int>>; // indices of an intersecting pair
// Shamos-Hoey: does any pair intersect? Segments must have distinct endpoints.
Pair anyIntersection(std::vector<Segment> s) {
const int n = static_cast<int>(s.size());
for (Segment& g : s)
if (lexLess(g.b, g.a)) std::swap(g.a, g.b); // a is the first end swept
struct Event { Point p; int seg; bool isEnd; };
std::vector<Event> events;
for (int i = 0; i < n; ++i) {
events.push_back({s[i].a, i, false});
events.push_back({s[i].b, i, true});
}
std::ranges::sort(events, lexLess, &Event::p);
for (size_t k = 0; k + 1 < events.size(); ++k) // two segments share an endpoint
if (events[k].p == events[k + 1].p)
return std::pair{events[k].seg, events[k + 1].seg};
// Is segment i below segment j? Decided where the later of the two starts.
auto below = [&](int i, int j) {
if (i == j) return false;
bool iFirst = lexLess(s[i].a, s[j].a);
const Segment& u = iFirst ? s[i] : s[j]; // u started first
const Segment& v = iFirst ? s[j] : s[i];
int side = orientation(u.a, u.b, v.a); // +1: v starts above u
if (side == 0) side = orientation(u.a, u.b, v.b); // v starts on u: use its end
if (side == 0) return i < j; // collinear overlap: any consistent order
return iFirst ? side > 0 : side < 0;
};
std::set<int, decltype(below)> status(below);
std::vector<decltype(status)::iterator> where(n);
for (const Event& e : events) {
if (!e.isEnd) {
auto it = status.insert(e.seg).first;
where[e.seg] = it;
if (it != status.begin() && intersects(s[*std::prev(it)], s[e.seg]))
return std::pair{*std::prev(it), e.seg};
if (std::next(it) != status.end() && intersects(s[*std::next(it)], s[e.seg]))
return std::pair{*std::next(it), e.seg};
} else {
auto it = where[e.seg]; // erase by iterator: no comparisons needed
if (it != status.begin() && std::next(it) != status.end() &&
intersects(s[*std::prev(it)], s[*std::next(it)]))
return std::pair{*std::prev(it), *std::next(it)};
status.erase(it);
}
}
return std::nullopt;
}
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
};
Pair bruteForce(const std::vector<Segment>& s) { // every pair: O(n^2)
for (size_t i = 0; i < s.size(); ++i)
for (size_t j = i + 1; j < s.size(); ++j)
if (intersects(s[i], s[j])) return std::pair{int(i), int(j)};
return std::nullopt;
}
// n segments, each inside its own 10 by 10 cell, so no two intersect. "No" is the
// worst case for both algorithms: neither can stop early.
std::vector<Segment> disjointSegments(int n, Lcg& rng) {
const int cols = static_cast<int>(std::sqrt(n)) + 1;
std::vector<Segment> segs;
for (int i = 0; i < n; ++i) {
long long x0 = 10LL * (i % cols), y0 = 10LL * (i / cols);
Point a{x0 + rng.range(0, 8), y0 + rng.range(0, 8)}, b = a;
while (b == a) b = {x0 + rng.range(0, 8), y0 + rng.range(0, 8)};
segs.push_back({a, b});
}
return segs;
}
int main() {
Lcg rng{7};
std::println("{:>5} {:>12} {:>8} {:>8}", "n", "brute force", "sweep", "sweep/n");
for (int n : {125, 250, 500, 1000, 2000}) {
auto segs = disjointSegments(n, rng);
pairTests = 0;
bool slow = bruteForce(segs).has_value();
long slowTests = pairTests;
pairTests = 0;
bool fast = anyIntersection(segs).has_value();
long fastTests = pairTests;
if (slow || fast) std::println("unexpected intersection");
std::println("{:>5} {:>12} {:>8} {:>8.2f}", n, slowTests, fastTests,
double(fastTests) / n);
}
}Press Run to compile and run this example.
n brute force sweep sweep/n
125 7750 214 1.71
250 31125 494 1.98
500 124750 1063 2.13
1000 499500 2225 2.23
2000 1999000 4700 2.35Brute force tests every pair, n(n-1)/2 of them, and quadruples when n doubles. The sweep tests a segment only against its neighbors: at most two tests when it enters and one when it leaves, so never more than 3n. The sweep/n column creeps up towards that limit as larger sets give more segments two neighbors. At n = 2000 the sweep tests 4,700 pairs instead of almost 2 million. The rest of its cost is sorting the 2n events and 2n operations on the std::set, each O(log n), which is where the O(n log n) total comes from.