Segment Trees for Geometry, in depth: coordinate compression, cover counts, rectangle union area, perimeter and deepest overlap

The problems this solves

The geometric problems this family solves share a shape: a set of axis-aligned intervals or rectangles, and an aggregate over the region they cover.

All of them start the same way. Real coordinates are arbitrary numbers, but only the distinct coordinates that appear as edges matter, so the first step is always to compress them.

Coordinate compression and elementary intervals

Collect every distinct y-coordinate from rectangle bottoms and tops and sort them: ys = [y0, y1, ..., y(m-1)]. Between consecutive values lie m-1 elementary intervals, [ys[i], ys[i+1]). No rectangle edge falls strictly inside an elementary interval, so each one is either fully covered by a rectangle or not at all. Those m-1 intervals, not the m points, are the leaves of the tree. A rectangle spanning [y1, y2) covers leaves index[y1] through index[y2] - 1 inclusive.

Two classic bugs come from this step. Building the tree over m points instead of m-1 intervals double counts or drops a boundary. Treating intervals as closed on both ends makes two rectangles that merely touch appear to overlap and adds zero-width slivers. Use half-open intervals throughout and the arithmetic stays consistent: a leaf's length is ys[i+1] - ys[i], and a node spanning leaves lo..hi has length ys[hi+1] - ys[lo].

The cover-count tree

For union problems the tree must answer one question at any moment: how much of the y-axis is covered by at least one active rectangle? The obvious approach, lazily adding +1 and -1 to ranges and counting nonzero leaves, does not work, because a count-of-nonzero is not composable under range add. The standard trick exploits a property of the sweep: every removal exactly matches an earlier insertion of the same range. So a range is always removed at the same canonical nodes where it was added, and nodes never need to push anything down.

Each node keeps cnt, the number of active ranges that cover this node's whole span via canonical decomposition, and length, the covered length within its span. The invariant is: if cnt is positive the node is fully covered; otherwise its covered length is the sum of its children's, or zero for a leaf.

class CoverTree:
    """Covered length over elementary intervals [ys[i], ys[i+1])."""
    def __init__(self, ys):
        self.ys = ys
        self.n = len(ys) - 1                 # leaves = elementary intervals
        self.cnt = [0] * (4 * self.n)
        self.length = [0] * (4 * self.n)

    def _pull(self, node, lo, hi):
        if self.cnt[node] > 0:
            self.length[node] = self.ys[hi + 1] - self.ys[lo]
        elif lo == hi:
            self.length[node] = 0
        else:
            self.length[node] = self.length[2 * node] + self.length[2 * node + 1]

    def update(self, l, r, delta, node=1, lo=0, hi=None):
        """Add delta to the cover count of leaves l..r inclusive."""
        if hi is None:
            hi = self.n - 1
        if r < lo or hi < l:
            return
        if l <= lo and hi <= r:
            self.cnt[node] += delta
        else:
            mid = (lo + hi) // 2
            self.update(l, r, delta, 2 * node, lo, mid)
            self.update(l, r, delta, 2 * node + 1, mid + 1, hi)
        self._pull(node, lo, hi)

    def covered(self):
        return self.length[1]

Note the subtle point in _pull: when a node's count drops back to zero, its length is recomputed from children whose own counts may still be positive from other rectangles. That is why the children's values must always be correct, and why it is enough to fix up nodes on the update path.

Union area with a sweep line

Each rectangle (x1, y1, x2, y2) produces two events: at x1 add +1 over its y-range, at x2 add -1. Sort events by x. Before applying an event, the covered length has been constant since the previous event, so add covered * (x - prev_x) to the area.

def union_area(rects):
    ys = sorted({y for _, y1, _, y2 in rects for y in (y1, y2)})
    index = {y: i for i, y in enumerate(ys)}
    events = []
    for x1, y1, x2, y2 in rects:
        if x1 < x2 and y1 < y2:              # skip degenerate boxes
            events.append((x1, +1, y1, y2))
            events.append((x2, -1, y1, y2))
    events.sort()
    tree = CoverTree(ys)
    area, prev_x = 0, None
    for x, delta, y1, y2 in events:
        if prev_x is not None:
            area += tree.covered() * (x - prev_x)
        tree.update(index[y1], index[y2] - 1, delta)
        prev_x = x
    return area

The order of events that share an x does not matter for area, because the slab between them has zero width. It does matter for perimeter and overlap, as the next sections show. With n rectangles there are 2n events, each costing O(log n), plus the sort: O(n log n) time and O(n) memory.

Worked example

Take A = (0, 0, 4, 3), B = (2, 1, 6, 5) and C = (5, 0, 8, 2). The compressed y-values are 0, 1, 2, 3, 5, giving four leaves: [0,1), [1,2), [2,3) and [3,5). The events in x order and the slab each one closes:

The total is 29. Check it by inclusion and exclusion: the areas are 12, 16 and 6, summing to 34; A and B overlap in a 2 by 2 square and B and C in a 1 by 1 square, so 34 - 4 - 1 = 29. At x = 5, shown in the diagram, the event for C adds +1 to leaves 0 and 1, which the tree covers with one canonical node, [0,2). B already sits on node [2,5) and on leaf [1,2). The root has count zero, so its length is 2 + 3 = 5.

Deepest overlap: range add and max

To find the point covered by the most rectangles, the per-node aggregate changes from covered length to maximum count, and range add composes with max: max(a + d, b + d) = max(a, b) + d. Each node stores its pending add and the maximum within its span including that add, so no push-down is needed either.

class MaxAddTree:
    def __init__(self, n):
        self.n, self.mx, self.add = n, [0] * (4 * n), [0] * (4 * n)

    def update(self, l, r, delta, node=1, lo=0, hi=None):
        if hi is None:
            hi = self.n - 1
        if r < lo or hi < l:
            return
        if l <= lo and hi <= r:
            self.mx[node] += delta
            self.add[node] += delta
            return
        mid = (lo + hi) // 2
        self.update(l, r, delta, 2 * node, lo, mid)
        self.update(l, r, delta, 2 * node + 1, mid + 1, hi)
        self.mx[node] = self.add[node] + max(self.mx[2 * node], self.mx[2 * node + 1])

def max_overlap(rects):
    ys = sorted({y for _, y1, _, y2 in rects for y in (y1, y2)})
    index = {y: i for i, y in enumerate(ys)}
    # kind 0 = leave, 1 = enter: at equal x, leave first, so touching boxes do not overlap
    events = sorted([(x1, 1, y1, y2) for x1, y1, x2, y2 in rects] +
                    [(x2, 0, y1, y2) for x1, y1, x2, y2 in rects])
    tree, best = MaxAddTree(len(ys) - 1), 0
    for x, kind, y1, y2 in events:
        tree.update(index[y1], index[y2] - 1, 1 if kind else -1)
        best = max(best, tree.mx[1])
    return best

For the example the answer is 2. Tie-breaking encodes the semantics: processing leaves before entries treats rectangles as half-open, which is what you want for scheduling (a meeting ending at 10:00 does not clash with one starting at 10:00). Reverse the order if touching counts as overlap.

Union perimeter

The perimeter of the union has two parts. Vertical edges appear at events: each event changes the covered length, and the absolute change is the vertical boundary added there. Horizontal edges run along the slabs: if the cross-section consists of k disjoint covered runs, the slab contributes 2k times its width. The tree therefore also tracks the number of runs per node and whether each end of the span is covered, merging two children by adding their runs and subtracting one when the left child's right end and the right child's left end are both covered. Process entries before exits at equal x so that rectangles sharing an edge merge instead of producing a phantom boundary. For the example the perimeter is 28.

Choosing a structure

The segment tree is not always the right structure. For stabbing queries over a dynamic set of intervals (which intervals contain this point?), an interval tree reports the k answers in O(log n + k) and handles insertions without knowing coordinates in advance. When you only need counts at points and every update is a range add, a Fenwick tree over difference arrays is smaller and faster. For static counting of points in a query rectangle, a merge sort tree (a segment tree over x whose nodes store sorted y-lists) answers in O(log^2 n), and fractional cascading reduces that to O(log n) at the cost of complexity. If you need to query the cross-section as it was at an earlier sweep position, a persistent segment tree keeps every version.

Failure modes

• Floating-point coordinates. Compression by exact equality breaks when 0.1 + 0.2 appears as a coordinate. Snap inputs to an integer grid (for example micrometres or cents) before compressing.
• Integer overflow. Area is a product of two coordinate ranges; with 32-bit languages and 1e9 coordinates, accumulate in 64-bit or larger.
• Recursion depth. Python recursion is fine for depth around 20, but for very large inputs an iterative bottom-up tree is faster by a constant factor of several.
• Zero-width rectangles. They produce events that change nothing for area but can confuse perimeter logic; filter them up front.
• Array size. Allocate 4 times the number of leaves for the recursive layout; 2 times is enough only for the iterative layout.

What to do next

1. Implement CoverTree and union_area and test them against a brute-force grid count on random small inputs.
2. Reproduce the worked example by hand, including which canonical nodes each event touches.
3. Add max_overlap and decide explicitly whether touching boundaries count as overlap.
4. Extend the tree with run counts and end flags and verify the perimeter of the example is 28.
5. Snap any real-valued input to an integer grid before compression.
6. For a production workload, compare against an interval tree or Fenwick tree before committing.