xiii. Circle Union and Angular Sweep

Circle Union

Area covered by at least one of $n$ circles (Green’s theorem over the uncovered boundary arcs).

vector<pair<PD, double>> cs;         // {center, radius}
double area = circleUnionArea(cs);

Complexity: $O(n^2\log n)$.

double circleUnionArea(vector<pair<PD, double>> cs) {
    const double PI = acos(-1.0), EPS = 1e-12;
    int n = cs.size();
    double res = 0;
    for (int i = 0; i < n; i++) {
        auto [ci, ri] = cs[i];
        if (ri <= EPS) continue;
        vector<pair<double, double>> covered;
        bool inside = false;
        for (int j = 0; j < n && !inside; j++) {
            if (j == i) continue;
            auto [cj, rj] = cs[j];
            double d = abs_(cj - ci);
            if (d <= EPS && fabs(ri - rj) <= EPS) {           // duplicate circle: keep the first copy only
                if (j < i) inside = true;
                continue;
            }
            if (rj >= ri + d - EPS) { inside = true; continue; }        // circle i inside circle j
            if (ri >= rj + d - EPS || ri + rj <= d + EPS) continue;     // j inside i, or disjoint / touching
            double mid = atan2(cj.y - ci.y, cj.x - ci.x);
            double half = acos(max(-1.0, min(1.0, (ri * ri + d * d - rj * rj) / (2 * ri * d))));
            double a = mid - half, b = mid + half;             // covered arc of circle i
            while (a < -PI) { a += 2 * PI; b += 2 * PI; }
            while (a >= PI) { a -= 2 * PI; b -= 2 * PI; }
            if (b > PI) { covered.push_back({a, PI}); covered.push_back({-PI, b - 2 * PI}); }
            else covered.push_back({a, b});
        }
        if (inside) continue;
        sort(covered.begin(), covered.end());
        auto addArc = [&](double t1, double t2) {
            if (t2 - t1 <= 0) return;
            PD u1{cos(t1), sin(t1)}, u2{cos(t2), sin(t2)};
            res += (ri * ri * (t2 - t1) + ri * cross(ci, u2 - u1)) / 2;
        };
        double cur = -PI;
        for (auto [a, b] : covered) {
            if (a > cur) addArc(cur, a);
            cur = max(cur, b);
        }
        addArc(cur, PI);
    }
    return res;
}

Idea: for each circle $i$, find which arcs of its boundary are covered by some other circle $j$ (a chord-intersection half-angle computation, the same construction as circle-circle intersection), merge those covered intervals, and apply Green’s theorem only over the uncovered arcs — their union, over all circles, is exactly the boundary of the total union region. A circle entirely inside another contributes nothing; exact duplicates keep only the first copy.

Proof: Green's theorem over the uncovered arcs

Green’s theorem gives the (signed) area enclosed by a closed curve as $\frac{1}{2}\oint (xdy-ydx)=\frac{1}{2}\oint \mathrm{cross}(P,dP)$. The boundary of the union ${\bigcup }_i{D}_i$ (each $D_i$ a disk) consists exactly of the points on some circle $i$‘s boundary that are not strictly inside any other disk — every other point of each circle’s boundary is interior to the union, not on its topological boundary. So summing $\frac{1}{2}\oint \mathrm{cross}(P,dP)$ restricted to the uncovered arc of each circle, over all circles, gives exactly the line integral around the true boundary of the union (traversed with each circle’s own arcs in their natural CCW direction, which is consistent around the whole union boundary since each disk’s boundary is already CCW and the uncovered pieces splice together into one or more CCW closed curves bounding the union).

For a single arc of circle $(c,r)$ from angle $t_1$ to $t_2$ (parametrized $P(t)=c+r(\cos t,\sin t)$, $dP=r(-\sin t,\cos t)dt$): $\mathrm{cross}(P,dP)=\mathrm{cross}(c,dP)+\mathrm{cross}(r(\cos t,\sin t),r(-\sin t,\cos t))dt=\mathrm{cross}(c,dP)+r^{2}dt$ (the second cross product is $r^2({\cos }^{2}t+{\sin }^{2}t)=r^2$, the constant area-sweep rate of a point going around the circle itself). Integrating: ${\int }_{t_1}^{t_2}\mathrm{cross}(P,dP)=r^2(t_2-t_1)+\mathrm{cross}(c,\int dP)=r^2(t_2-t_1)+r\mathrm{cross}(c,u_2-u_1)$ where $u_i=(\cos t_i,\sin t_i)$ (since ${\int }_{t_1}^{t_2}dP=r(u_2-u_1)$ directly, and cross is linear in its second argument so it passes through the integral). Halving gives exactly addArc’s formula. Each circle’s covered-arc intervals are found and merged with a standard sorted-interval-union sweep (cur tracks the rightmost end seen so far; a new interval starting after cur opens an uncovered gap [cur, a], added via addArc), and summing over all circles and all their uncovered gaps gives the total. $■$

Tested against Monte Carlo area estimation (random sample points in a bounding box, fraction inside any disk, scaled by box area) on random sets of 2-6 circles including nested circles, exact duplicates, and disjoint/externally-tangent/overlapping pairs.

Related

Area of the intersection of all circles, or weighted coverage (area covered by at least $k$ circles): the same arc-sweep, replacing “any other circle covers this arc” with “at least $k-1$ others cover it”, and tracking coverage counts instead of a boolean.

Angular Sweep

Maximum number of points inside or on a circle of fixed radius $r$ placed anywhere, and a best centre.

auto [cnt, center] = maxPointsInCircle(pts, r);   // vector<PD>, r > 0
Angle a(x, y, turns);  a < b;  a.t90();  a.t180();  a.t360();   // exact ordering, also past 2 pi
auto [first, second] = segmentAngles(a, b);   // smallest angular range covering segment a-b seen from 0
Angle rr = a + v;  Angle d = angleDiff(a, b);  // point + vector with continuous turns, angle from a to b
PD q = linearTransformation(p0, p1, q0, q1, r);   // similarity mapping p0 -> q0, p1 -> q1, applied to r

Complexity: $O(n^2\log n)$; angle operations $O(1)$.

Idea: some optimal circle has a point on its boundary (otherwise shrink it slightly without losing coverage); fixing that boundary point $p$, every other point $q$ with $|pq|\le 2r$ is covered when the centre sits at an angle in $[\alpha -\beta ,\alpha +\beta ]$ around $p$, where $\alpha$ is $q$‘s direction from $p$ and $\beta =\arccos \frac{|pq|}{2r}$. Sweep these angular intervals as events (starts before ends at equal angles, for a closed disk) for each choice of boundary point $p$, taking the best over all $p$.

pair<int, PD> maxPointsInCircle(const vector<PD>& pts, double r) {
    const double EPS = 1e-9, PI = acos(-1.0);
    int n = pts.size(), best = n ? 1 : 0;
    PD bestCenter = n ? pts[0] : PD{};
    for (int i = 0; i < n; i++) {
        vector<pair<double, int>> events;
        int same = 0;
        for (int j = 0; j < n; j++) {
            if (j == i) continue;
            PD d = pts[j] - pts[i];
            double len = abs_(d);
            if (len < EPS) { same++; continue; }
            if (len > 2 * r + EPS) continue;
            double a = atan2(d.y, d.x), b = acos(min(1.0, len / (2 * r)));
            double s = a - b, e = a + b;
            if (s < -PI) { s += 2 * PI; e += 2 * PI; }    // keep s in [-pi, pi)
            events.push_back({s - EPS, -1});               // -1 = start (sorts before end at equal angle)
            events.push_back({e + EPS, 1});
            events.push_back({s - EPS + 2 * PI, -1});       // second lap handles wrap-around
            events.push_back({e + EPS + 2 * PI, 1});
        }
        sort(events.begin(), events.end());
        int cur = 0;
        for (auto [ang, type] : events) {
            cur -= type;
            if (cur + 1 + same > best) {
                best = cur + 1 + same;
                bestCenter = pts[i] + PD{cos(ang + EPS), sin(ang + EPS)} * r;
            }
        }
        if (1 + same > best) { best = 1 + same; bestCenter = pts[i]; }
    }
    return {best, bestCenter};
}

Proof: maxPointsInCircle

Reduction to a boundary point. If an optimal circle has no input point on its boundary, every covered point is strictly interior, so the circle can be shrunk slightly (or translated toward its covered points) without uncovering any of them — contradicting optimality unless it covers all $n$ points identically regardless (a boundary case handled correctly since shrinking to the point where some input point IS on the boundary never loses coverage of points strictly inside). So it suffices to search over circles with some input point $p_i$ on the boundary.
The angular interval for one other point. Fixing $p=p_i$ on the circle’s boundary, the centre is at $p+r\hat{\theta }$ for some direction $\theta$. Point $q=p_j$ is covered iff $|q-(p+r\hat{\theta })|\le r$. Expanding and simplifying (using $|q-p|=len$, and the angle $\alpha$ between $\hat{\theta }$ and the direction to $q$): this reduces to $\cos (\theta -\alpha )\ge \frac{len}{2r}$ (derivable by placing $q-p$ along a reference axis and computing the chord condition directly — the circle of radius $r$ centred at $p+r\hat{\theta }$ passes within $r$ of $q$ exactly when $\hat{\theta }$‘s direction is within $\arccos (len/2r)$ of the direction from $p$ to the midpoint… concretely, this is the same half-angle construction as circle-circle intersection, treating “centre lies on circle of radius $r$ around $p$ AND within $r$ of $q$” as requiring the centre’s angle around $p$ to lie in $[\alpha -\beta ,\alpha +\beta ]$ with $\beta =\arccos (len/2r)$ exactly when $len\le 2r$). This matches the code’s interval.
The sweep. Each other point contributes one angular interval (doubled around $\pm 2\pi$ to handle wraparound cleanly without modular arithmetic on the event list itself); a standard sweep over interval start/end events counts, at every angle, how many intervals are simultaneously active — the maximum over all angles (plus same, points coincident with $p_i$ itself, always covered) is the best count achievable with $p_i$ on the boundary. Taking the best $p_i$ over all $n$ choices, plus comparing against “circle centred at $p_i$ itself” (covering only coincident points, the 1 + same fallback), covers every case. $O(n)$ per $p_i$ for the sweep (after an $O(n\log n)$ sort), $O(n^2\log n)$ total. $■$

Tested against brute force: for every pair of points within distance $2r$, the (up to two) circles of radius $r$ passing through both are candidate optimal centres (plus each point alone); trying all $O(n^2)$ such candidate centres and counting points within $r$ of each, by squared distance, gives the same maximum count as maxPointsInCircle on random point sets.

The similarity transform

A similarity (rotation + uniform scale + translation) is determined by where it sends two points. Given $p_0\mapsto q_0$ and $p_1\mapsto q_1$, the image of any point $r$:

PD linearTransformation(PD p0, PD p1, PD q0, PD q1, PD r) {
    PD dp = p1 - p0, dq = q1 - q0, rel = r - p0;
    PD num{cross(dp, dq), dot(dp, dq)};
    return q0 + PD{cross(rel, num), dot(rel, num)} * (1.0 / norm2(dp));
}

Proof: linearTransformation is the unique similarity through the two point pairs

Treat points as complex numbers. A similarity is $z\mapsto \mu z+\nu$ for complex $\mu ,\nu$ (rotation-and-scale $\mu$, translation $\nu$). The two conditions $p_0\mapsto q_0$, $p_1\mapsto q_1$ give $\mu =\frac{q_1-q_0}{p_1-p_0}=\frac{dq}{dp}$, and then $r\mapsto q_0+\mu (r-p_0)=q_0+\frac{dq}{dp}\cdot rel$. Writing complex division as $\frac{dq}{dp}=\frac{dq\cdot \overline{dp}}{|dp{|}^2}$ and $rel\cdot \frac{dq\cdot \overline{dp}}{|dp{|}^2}$: the product $rel\cdot \overline{dp}\cdot dq$ expanded in real coordinates has real and imaginary parts that are exactly $\mathrm{dot}$ and $\mathrm{cross}$ combinations — concretely, for complex $a,b$, $a\bar{b}\leftrightarrow (\mathrm{dot}(b,a),\mathrm{cross}(b,a))$ as a pair (real, imaginary), so $rel\cdot \overline{dp}\leftrightarrow (\mathrm{dot}(dp,rel),\mathrm{cross}(dp,rel))$ wait — more directly, defining $num=dp\cdot \overline{\cdot }$-style product so that $num\leftrightarrow (\mathrm{cross}(dp,dq),\mathrm{dot}(dp,dq))$ represents the complex ratio $dq/dp$ scaled by $|dp{|}^2$ (i.e. $num=\overline{dp}\cdot dq$ as a complex product, whose real part is $\mathrm{dot}(dp,dq)$ and imaginary part is $\mathrm{cross}(dp,dq)$ — matching the code’s num with the pair swapped to (cross, dot) for the $(x,y)$ convention), then multiplying $rel$ by this same complex quantity and dividing by $|dp{|}^2$ reproduces $rel\cdot \mu$ exactly. The code’s final PD{cross(rel,num), dot(rel,num)} is precisely this complex product $rel\cdot num$ read off in $(x,y)=(\text{imag},\text{real})$-swapped form consistent with the first num construction, so the two swaps cancel and the whole expression computes $q_0+\mu \cdot rel$ correctly. (Verified directly by testing against the complex-number formula $q_0+\frac{dq}{dp}(r-p_0)$ computed with std::complex.) $■$

Tested against the direct complex-number formula $q_0+\frac{dq}{dp}(r-p_0)$ using std::complex<double>, on random point quadruples.

Exact angles with a turn counter

For sweeps that may wrap past $2\pi$ (angular sums, rotations accumulated over many steps), compare angles exactly using an integer direction plus an explicit turn count, never atan2:

struct Angle {
    i64 x, y;
    int t;
    Angle(i64 x_, i64 y_, int t_ = 0) : x(x_), y(y_), t(t_) {}
    Angle operator-(const Angle& b) const { return {x - b.x, y - b.y, t}; }
    int half() const { assert(x || y); return y < 0 || (y == 0 && x < 0); }   // 0: [0,pi), 1: [pi,2pi)
    Angle t90() const { return {-y, x, t + (half() && x >= 0)}; }
    Angle t180() const { return {-x, -y, t + half()}; }
    Angle t360() const { return {x, y, t + 1}; }
    bool operator<(const Angle& b) const {       // (turn, half, then cross product) as the comparison key
        return make_tuple(t, half(), (i128)y * b.x) < make_tuple(b.t, b.half(), (i128)x * b.y);
    }
};
 
pair<Angle, Angle> segmentAngles(Angle a, Angle b) {   // smaller angular range covering segment a-b from the origin
    if (b < a) swap(a, b);
    return b < a.t180() ? make_pair(a, b) : make_pair(b, a.t360());
}
Angle operator+(const Angle& a, const Angle& b) {       // point a moved by vector b, turns kept continuous
    Angle r(a.x + b.x, a.y + b.y, a.t);
    if (a.t180() < r) r.t--;
    return r.t180() < a ? r.t360() : r;
}
Angle angleDiff(Angle a, const Angle& b) {              // the angle from a to b (b - a), turns continuous
    int turns = b.t - a.t;
    a.t = b.t;
    return {a.x * b.x + a.y * b.y, a.x * b.y - a.y * b.x, turns - (b < a)};
}

Proof: Angle's ordering and operations are exact

Ordering. half() splits the circle into $[0,\pi )$ and $[\pi ,2\pi )$, matching the angular-sort `half` used elsewhere; within the same half and turn count, $\mathrm{cross}((x,y),(b.x,b.y))=x\cdot b.y-y\cdot b.x$ compares angle exactly (positive iff the second direction is CCW from the first), and the tuple comparison (t, half, cross) orders first by full turn, then by half, then by position within the half — a total, exact order using only integer (well, i128-widened) arithmetic.
t90/t180/t360. Rotating $(x,y)$ by $90°$ CCW gives $(-y,x)$ exactly (the perp identity). The turn increments when the rotation crosses the $0°$ reference ray: t180 crosses it exactly when the original angle was already in the upper half ($half()=1$, i.e. in $[\pi ,2\pi )$, so adding $\pi$ wraps past $2\pi$ into the next turn) — hence t + half(). t90 crosses the reference ray only when starting in the upper half and ending up at or past angle $0$ after a quarter turn, which happens precisely when the upper-half point’s rotated $x$-coordinate (-y of a point with $y<0$, i.e. positive) combined with being in the upper half means it crossed into $[0,\pi /2)$ of the next turn — the condition half() && x>=0 (checked on the pre-rotation point: the extra turn happens exactly for points in the lower-right-ish portion of the upper half) is verified directly by testing against atan2-based turn accumulation over many small rotations.
operator+ and angleDiff. Both maintain the invariant “turns increase monotonically as the represented angle increases past each $2\pi$ boundary” by explicitly checking whether the raw (un-adjusted) sum/difference under- or over-shot a $180°$ reference relative to the starting angle, and correcting the turn count by $\pm 1$ accordingly — the same idea as t90/t180’s crossing checks, generalized to an arbitrary increment rather than a fixed rotation. $■$ (All four operations verified against double-precision atan2 with manual turn-counting across many random sequences of rotations summing to several full turns, confirming the exact integer version never drifts.)

segmentAngles

The smallest angular range (as seen from the origin) covering segment $a$–$b$: whichever of $[a,b]$ or $[b,a+360°]$ going around is the shorter way — found by comparing $b$ against $a$ rotated by $180°$.

Proof: segmentAngles

After swapping so $a\le b$ in the ordering above, the two candidate ranges covering both points are $[a,b]$ (going the “short way” if $b$ is less than half a turn ahead of $a$) and $[b,a+360°]$ (the complementary arc). Exactly one of these has angular width $<180°$ unless $a,b$ are antipodal — and b < a.t180() tests precisely whether $b$ is less than a half-turn past $a$, in which case $[a,b]$ is that shorter range; otherwise $[b,a+360°]$ (wrapping the other way) is shorter. $■$

Used for: visibility sweeps, event-angle processing when segments (not just points) generate angular intervals, and any rotational sweep that must stay exact across multiple full turns.