扫描线算法：从线段交到矩形面积并

计算几何的核心困难在于”二维”——点、线、面在平面上纠缠，暴力枚举的复杂度往往是 O(n^2) 甚至更高。扫描线（Sweep Line）是一种将二维问题降维为一维动态问题的通用范式：想象一条竖直的线从左向右匀速扫过整个平面，在扫描过程中维护一个一维数据结构来回答当前切面上的问题。每当扫描线遇到关键的”事件点”（端点、边界、交点），就更新这个数据结构。

这个思想看似简单，却能解决大量看起来毫不相关的问题——线段交点枚举、矩形面积并、矩形周长并、最近点对、甚至 Voronoi 图的构造。本文将从最经典的 Bentley-Ottmann 线段交点算法出发，逐步展开到矩形面积并、周长并、最近点对，最后讨论扫描线在工业界（特别是 EDA）中的实际应用与工程陷阱。

一、扫描线范式：降维的艺术

1.1 核心思想

扫描线的本质可以用一句话概括：

用一条虚拟的线（通常是竖直线 x = t）从左到右扫过平面，将二维静态问题转化为一维动态问题。

具体来说，任何可以用扫描线解决的问题都遵循以下模板：

预处理：将所有几何对象的”关键位置”提取为事件（Event），按扫描方向排序后放入事件队列（Event Queue）。
维护状态：用一个一维数据结构（状态结构，Status Structure）维护扫描线当前位置处与几何对象的交互信息。
处理事件：每当扫描线到达一个事件点，更新状态结构，并可能产生新的事件（如交点事件）插入事件队列。
输出结果：在处理事件的过程中逐步输出答案。

1.2 抽象模板

用伪代码描述这个范式：

SWEEP-LINE(objects):
    Q ← 初始化事件队列（从 objects 提取所有端点/边界）
    T ← 初始化空的状态结构
    while Q 非空:
        event ← Q.extractMin()
        处理 event，更新 T
        可能向 Q 插入新事件
    return 收集到的结果

这里的关键设计决策有三个：

事件的定义：什么构成一个事件？对于不同问题，事件的含义不同。
状态结构的选择：BST、线段树、还是别的？
事件之间的交互：处理一个事件时是否会产生新事件？

1.3 为什么是”从左到右”

扫描方向的选择是任意的——从左到右、从下到上、甚至旋转任意角度都可以。但从左到右是最常见的约定，因为 x 轴通常是”第一坐标”。在某些问题中（如矩形面积并），从下到上更自然。关键是保持一致。

1.4 扫描线与分治的关系

扫描线和分治经常可以解决同一类问题（如最近点对），但思路截然不同：

分治：空间上切分，递归求解子问题，合并时处理跨界情况。
扫描线：时间上推进，增量地维护全局信息。

扫描线通常更容易实现增量式输出，也更容易与在线查询结合。分治在理论上有时能给出更好的复杂度界（如用随机化），但实现起来合并步骤往往更复杂。

二、Bentley-Ottmann 算法：线段交点枚举

2.1 问题定义

给定平面上 n 条线段，找出所有交点。

暴力做法是 O(n^2)——枚举所有线段对，检查是否相交。但如果交点数 k 远小于 n^2（这在大多数实际场景中成立），我们希望复杂度与 k 相关。Bentley 和 Ottmann 在 1979 年提出的算法达到了 O((n+k) log n)。

2.2 核心观察

Bentley-Ottmann 算法基于一个关键观察：

两条线段如果相交，那么在相交之前的某个时刻，它们在扫描线的状态结构中必然是相邻的。

这意味着我们不需要检查所有线段对，只需要在状态结构中检查相邻的线段对是否相交。

2.3 三种事件

算法处理三种事件：

左端点事件（Left Endpoint）：线段开始。将线段插入状态结构 T，检查它与上下邻居是否相交。
右端点事件（Right Endpoint）：线段结束。从 T 中删除线段，检查它的原上下邻居（现在变成相邻的）是否相交。
交点事件（Intersection）：两条线段相交。报告交点，在 T 中交换两条线段的顺序，然后检查各自的新邻居是否相交。

2.4 数据结构

事件队列 Q：优先队列（通常用平衡 BST 实现），按 x 坐标排序。需要支持插入和删除（因为可能插入重复的交点事件后需要去重）。
状态结构 T：平衡 BST，按线段在当前扫描线位置的 y 坐标排序。注意：比较函数是动态的——同一对线段在不同 x 位置的顺序可能不同。

2.5 算法流程

BENTLEY-OTTMANN(segments):
    Q ← 所有线段端点，按 x 坐标排序
    T ← 空 BST
    intersections ← []

    while Q 非空:
        p ← Q.extractMin()
        sweep_x ← p.x    // 更新扫描线位置

        if p 是左端点:
            s ← p 对应的线段
            T.insert(s)
            above ← T.successor(s)
            below ← T.predecessor(s)
            if above 存在 且 s 与 above 相交于扫描线右侧:
                Q.insert(交点事件)
            if below 存在 且 s 与 below 相交于扫描线右侧:
                Q.insert(交点事件)

        else if p 是右端点:
            s ← p 对应的线段
            above ← T.successor(s)
            below ← T.predecessor(s)
            T.delete(s)
            if above 存在 且 below 存在 且 above 与 below 相交于扫描线右侧:
                Q.insert(交点事件)

        else:  // p 是交点事件
            s1, s2 ← p 对应的两条线段（s1 在上，s2 在下）
            intersections.append(p)
            T.swap(s1, s2)    // 交换顺序
            // s1 现在在下面，检查 s1 的新下邻居
            new_below ← T.predecessor(s1)
            if new_below 存在 且 s1 与 new_below 相交于右侧:
                Q.insert(交点事件)
            // s2 现在在上面，检查 s2 的新上邻居
            new_above ← T.successor(s2)
            if new_above 存在 且 s2 与 new_above 相交于右侧:
                Q.insert(交点事件)

    return intersections

2.6 复杂度分析

初始事件数：2n（每条线段两个端点）。
交点事件数：k（最多 O(n^2)，但通常远小于此）。
每个事件的处理：O(log n)（BST 操作）。
总复杂度：O((n+k) log n)。
空间复杂度：O(n+k)。

一个重要的事实：事件队列中同时存在的事件数不超过 O(n)——这是因为在任何时刻，每对相邻线段最多产生一个待处理的交点事件。

2.7 正确性论证

算法的正确性依赖于以下不变量：

在扫描线到达任何交点之前，该交点对应的两条线段在 T 中已经成为邻居，因此该交点已经被检测到并加入了 Q。

证明思路：考虑交点 p = s_i 与 s_j 的交点。在扫描线从左向右推进的过程中，s_i 和 s_j 的相对顺序最终会在 p 处翻转。在翻转之前，它们之间可能有其他线段，但这些”阻隔”线段的右端点会在 p 之前被处理（从 T 中删除），或者它们与 s_i / s_j 的交点会在 p 之前被处理（导致顺序交换）。无论哪种情况，s_i 和 s_j 最终会在 p 之前成为邻居。

三、Bentley-Ottmann 的 C 实现

下面给出一个完整但精简的 C 实现。为了聚焦核心逻辑，使用了简单的数组排序代替平衡 BST，并对事件队列使用了最小堆。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

#define MAX_SEGMENTS 1024
#define MAX_EVENTS   8192
#define EPS          1e-9

/* -------- 基本几何类型 -------- */

typedef struct {
    double x, y;
} Point;

typedef struct {
    Point left, right;  /* left.x <= right.x */
    int id;
} Segment;

typedef enum {
    EVT_LEFT,       /* 左端点 */
    EVT_RIGHT,      /* 右端点 */
    EVT_INTERSECT   /* 交点 */
} EventType;

typedef struct {
    Point    pt;
    EventType type;
    int      seg1;   /* 相关线段索引 */
    int      seg2;   /* 仅交点事件使用 */
} Event;

/* -------- 全局状态 -------- */

static Segment segments[MAX_SEGMENTS];
static int     n_segments;

/* 事件最小堆 */
static Event event_heap[MAX_EVENTS];
static int   heap_size;

/* 活跃线段表（简化：有序数组） */
static int   active[MAX_SEGMENTS];
static int   n_active;

/* 结果 */
static Point results[MAX_EVENTS];
static int   n_results;

/* 当前扫描线 x 坐标 */
static double sweep_x;

/* -------- 几何工具函数 -------- */

static double cross2d(double ax, double ay, double bx, double by)
{
    return ax * by - ay * bx;
}

static int segments_intersect(const Segment *a, const Segment *b, Point *out)
{
    double dx1 = a->right.x - a->left.x;
    double dy1 = a->right.y - a->left.y;
    double dx2 = b->right.x - b->left.x;
    double dy2 = b->right.y - b->left.y;
    double denom = cross2d(dx1, dy1, dx2, dy2);

    if (fabs(denom) < EPS) return 0;  /* 平行或共线 */

    double t = cross2d(b->left.x - a->left.x, b->left.y - a->left.y,
                        dx2, dy2) / denom;
    double u = cross2d(b->left.x - a->left.x, b->left.y - a->left.y,
                        dx1, dy1) / denom;

    if (t < -EPS || t > 1.0 + EPS || u < -EPS || u > 1.0 + EPS)
        return 0;

    out->x = a->left.x + t * dx1;
    out->y = a->left.y + t * dy1;
    return 1;
}

/* 求线段在 x = sweep_x 处的 y 值 */
static double y_at_sweep(int seg_idx)
{
    const Segment *s = &segments[seg_idx];
    if (fabs(s->right.x - s->left.x) < EPS)
        return (s->left.y + s->right.y) / 2.0;
    double t = (sweep_x - s->left.x) / (s->right.x - s->left.x);
    return s->left.y + t * (s->right.y - s->left.y);
}

/* -------- 事件堆操作 -------- */

static int event_less(const Event *a, const Event *b)
{
    if (fabs(a->pt.x - b->pt.x) > EPS)
        return a->pt.x < b->pt.x;
    return a->pt.y < b->pt.y;
}

static void heap_swap(int i, int j)
{
    Event tmp = event_heap[i];
    event_heap[i] = event_heap[j];
    event_heap[j] = tmp;
}

static void heap_push(Event e)
{
    int i = heap_size++;
    event_heap[i] = e;
    while (i > 0) {
        int parent = (i - 1) / 2;
        if (event_less(&event_heap[i], &event_heap[parent])) {
            heap_swap(i, parent);
            i = parent;
        } else break;
    }
}

static Event heap_pop(void)
{
    Event top = event_heap[0];
    event_heap[0] = event_heap[--heap_size];
    int i = 0;
    for (;;) {
        int smallest = i;
        int left = 2 * i + 1, right = 2 * i + 2;
        if (left < heap_size && event_less(&event_heap[left], &event_heap[smallest]))
            smallest = left;
        if (right < heap_size && event_less(&event_heap[right], &event_heap[smallest]))
            smallest = right;
        if (smallest == i) break;
        heap_swap(i, smallest);
        i = smallest;
    }
    return top;
}

/* -------- 活跃表操作 -------- */

static int active_find(int seg_idx)
{
    for (int i = 0; i < n_active; i++)
        if (active[i] == seg_idx) return i;
    return -1;
}

static void active_insert(int seg_idx)
{
    double y = y_at_sweep(seg_idx);
    int pos = n_active;
    for (int i = 0; i < n_active; i++) {
        if (y < y_at_sweep(active[i]) - EPS) {
            pos = i;
            break;
        }
    }
    for (int i = n_active; i > pos; i--)
        active[i] = active[i - 1];
    active[pos] = seg_idx;
    n_active++;
}

static void active_remove(int seg_idx)
{
    int pos = active_find(seg_idx);
    if (pos < 0) return;
    for (int i = pos; i < n_active - 1; i++)
        active[i] = active[i + 1];
    n_active--;
}

static void active_swap_adj(int seg1, int seg2)
{
    int p1 = active_find(seg1);
    int p2 = active_find(seg2);
    if (p1 >= 0 && p2 >= 0) {
        active[p1] = seg2;
        active[p2] = seg1;
    }
}

/* 检查两条活跃线段是否相交于扫描线右侧，若是则加入事件堆 */
static void check_intersection(int seg1, int seg2)
{
    if (seg1 < 0 || seg2 < 0) return;
    Point ip;
    if (segments_intersect(&segments[seg1], &segments[seg2], &ip)) {
        if (ip.x > sweep_x + EPS) {
            Event e;
            e.pt = ip;
            e.type = EVT_INTERSECT;
            e.seg1 = seg1;
            e.seg2 = seg2;
            heap_push(e);
        }
    }
}

/* -------- 主算法 -------- */

void bentley_ottmann(void)
{
    heap_size = 0;
    n_active = 0;
    n_results = 0;

    /* 初始化事件：所有线段的左右端点 */
    for (int i = 0; i < n_segments; i++) {
        /* 确保 left.x <= right.x */
        if (segments[i].left.x > segments[i].right.x) {
            Point tmp = segments[i].left;
            segments[i].left = segments[i].right;
            segments[i].right = tmp;
        }
        segments[i].id = i;

        Event el = { .pt = segments[i].left,  .type = EVT_LEFT,  .seg1 = i, .seg2 = -1 };
        Event er = { .pt = segments[i].right, .type = EVT_RIGHT, .seg1 = i, .seg2 = -1 };
        heap_push(el);
        heap_push(er);
    }

    while (heap_size > 0) {
        Event ev = heap_pop();
        sweep_x = ev.pt.x;

        switch (ev.type) {
        case EVT_LEFT: {
            int s = ev.seg1;
            active_insert(s);
            int pos = active_find(s);
            if (pos > 0)
                check_intersection(active[pos - 1], s);
            if (pos < n_active - 1)
                check_intersection(s, active[pos + 1]);
            break;
        }
        case EVT_RIGHT: {
            int s = ev.seg1;
            int pos = active_find(s);
            int below = (pos > 0) ? active[pos - 1] : -1;
            int above = (pos < n_active - 1) ? active[pos + 1] : -1;
            active_remove(s);
            if (below >= 0 && above >= 0)
                check_intersection(below, above);
            break;
        }
        case EVT_INTERSECT: {
            int s1 = ev.seg1, s2 = ev.seg2;
            results[n_results++] = ev.pt;
            active_swap_adj(s1, s2);
            int p1 = active_find(s1);
            int p2 = active_find(s2);
            if (p1 > 0)
                check_intersection(active[p1 - 1], s1);
            if (p2 < n_active - 1)
                check_intersection(s2, active[p2 + 1]);
            break;
        }
        }
    }
}

/* -------- 测试 -------- */

int main(void)
{
    n_segments = 3;
    segments[0] = (Segment){ {1, 1}, {5, 5}, 0 };
    segments[1] = (Segment){ {1, 5}, {5, 1}, 0 };
    segments[2] = (Segment){ {0, 3}, {6, 3}, 0 };

    bentley_ottmann();

    printf("Found %d intersection(s):\n", n_results);
    for (int i = 0; i < n_results; i++)
        printf("  (%.4f, %.4f)\n", results[i].x, results[i].y);

    return 0;
}

上面的实现为了可读性使用了数组而非平衡 BST。在生产环境中，活跃表应当用红黑树或 AVL 树实现，以保证 O(log n) 的插入、删除和邻居查询。

四、矩形面积并

4.1 问题定义

给定 n 个轴对齐矩形，求它们的面积并（即所有矩形覆盖区域的总面积）。

这是扫描线最经典的应用之一，也是竞赛中的常见题目。

4.2 思路

将每个矩形拆成两条竖直边：左边（x = x1，开始事件）和右边（x = x2，结束事件）。
将所有事件按 x 坐标排序。
用一个数据结构维护当前扫描线上被覆盖的 y 区间总长度。
每处理一个事件，面积增量 = 前后两个事件的 x 间距 * 当前被覆盖的 y 总长度。

4.3 关键：线段树维护覆盖

“当前被覆盖的 y 区间总长度”可以用一棵线段树高效维护。这棵线段树有以下特点：

叶节点对应离散化后的 y 坐标区间。
每个节点维护两个值：cnt（被完整覆盖的次数）和 len（被覆盖的总长度）。
区间修改：对于一条”开始”边，将对应 y 区间的 cnt 加一；对于”结束”边，减一。
查询：根节点的 len 就是当前总覆盖长度。

关键细节：这棵线段树不需要懒标记下推。因为 cnt 的语义是”整个区间被完整覆盖的次数”，当 cnt > 0 时 len 就是区间全长；当 cnt == 0 时 len 是左右子节点 len 之和。这使得实现比一般的线段树更简单。

4.4 离散化

y 坐标可能是实数，需要先离散化。收集所有矩形的 y1 和 y2，排序去重后得到 m 个不同的 y 值，形成 m-1 个基本区间。线段树建立在这些基本区间上。

4.5 完整 C 实现

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAXN 100005

/* -------- 事件与离散化 -------- */

typedef struct {
    double x, y1, y2;
    int type;  /* +1 = 左边（进入），-1 = 右边（离开） */
} SweepEvent;

static SweepEvent events[MAXN * 2];
static int n_events;

static double ys[MAXN * 2];  /* 离散化用 */
static int    n_ys;

static int event_cmp(const void *a, const void *b)
{
    const SweepEvent *ea = (const SweepEvent *)a;
    const SweepEvent *eb = (const SweepEvent *)b;
    if (ea->x < eb->x) return -1;
    if (ea->x > eb->x) return  1;
    return ea->type - eb->type;  /* 先处理进入 */
}

static int double_cmp(const void *a, const void *b)
{
    double da = *(const double *)a;
    double db = *(const double *)b;
    if (da < db) return -1;
    if (da > db) return  1;
    return 0;
}

static int y_index(double val)
{
    int lo = 0, hi = n_ys - 1;
    while (lo <= hi) {
        int mid = (lo + hi) / 2;
        if (ys[mid] < val - 1e-9)      lo = mid + 1;
        else if (ys[mid] > val + 1e-9) hi = mid - 1;
        else return mid;
    }
    return -1;
}

/* -------- 线段树（区间覆盖） -------- */

typedef struct {
    int    cnt;     /* 被完整覆盖的次数 */
    double len;     /* 被覆盖的总长度 */
} Node;

static Node tree[MAXN * 8];

static void build(int nd, int l, int r)
{
    tree[nd].cnt = 0;
    tree[nd].len = 0.0;
    if (l + 1 >= r) return;
    int mid = (l + r) / 2;
    build(nd * 2, l, mid);
    build(nd * 2 + 1, mid, r);
}

static void push_up(int nd, int l, int r)
{
    if (tree[nd].cnt > 0) {
        tree[nd].len = ys[r] - ys[l];
    } else if (l + 1 >= r) {
        tree[nd].len = 0.0;
    } else {
        tree[nd].len = tree[nd * 2].len + tree[nd * 2 + 1].len;
    }
}

static void update(int nd, int l, int r, int ql, int qr, int val)
{
    if (ql >= r || qr <= l) return;
    if (ql <= l && r <= qr) {
        tree[nd].cnt += val;
        push_up(nd, l, r);
        return;
    }
    int mid = (l + r) / 2;
    update(nd * 2, l, mid, ql, qr, val);
    update(nd * 2 + 1, mid, r, ql, qr, val);
    push_up(nd, l, r);
}

/* -------- 主算法 -------- */

double rectangle_area_union(int rects[][4], int n)
{
    n_events = 0;
    n_ys = 0;

    for (int i = 0; i < n; i++) {
        double x1 = rects[i][0], y1 = rects[i][1];
        double x2 = rects[i][2], y2 = rects[i][3];

        events[n_events++] = (SweepEvent){ x1, y1, y2, +1 };
        events[n_events++] = (SweepEvent){ x2, y1, y2, -1 };

        ys[n_ys++] = y1;
        ys[n_ys++] = y2;
    }

    /* 排序事件 */
    qsort(events, n_events, sizeof(SweepEvent), event_cmp);

    /* 离散化 y 坐标 */
    qsort(ys, n_ys, sizeof(double), double_cmp);
    int unique = 1;
    for (int i = 1; i < n_ys; i++)
        if (ys[i] - ys[i - 1] > 1e-9)
            ys[unique++] = ys[i];
    n_ys = unique;

    /* 建线段树 */
    build(1, 0, n_ys - 1);

    double area = 0.0;
    for (int i = 0; i < n_events; i++) {
        if (i > 0) {
            area += tree[1].len * (events[i].x - events[i - 1].x);
        }
        int lo = y_index(events[i].y1);
        int hi = y_index(events[i].y2);
        update(1, 0, n_ys - 1, lo, hi, events[i].type);
    }

    return area;
}

int main(void)
{
    /* 测试：两个重叠矩形 */
    int rects[][4] = {
        {1, 1, 4, 4},
        {2, 2, 5, 5}
    };
    int n = 2;
    double area = rectangle_area_union(rects, n);
    printf("Area union = %.2f\n", area);  /* 期望 15.00 */
    return 0;
}

4.6 复杂度分析

事件数：2n。
离散化后 y 坐标数：O(n)。
每个事件处理：线段树更新 O(log n)。
总复杂度：O(n log n)。
空间：O(n)。

4.7 线段树不需要 lazy 的原因

普通线段树做区间修改需要懒标记（lazy propagation），但矩形面积并的线段树不需要，原因在于：

cnt 字段表示”该区间被完整覆盖的次数”，而非”该区间内所有点的覆盖次数”。
修改操作总是成对出现：每个 +1 都有对应的 -1，且修改的区间完全相同。
查询只发生在根节点——我们只需要 tree[1].len。

这使得 push_up 可以正确地根据 cnt 和子节点信息计算 len，无需下推。

五、矩形周长并

5.1 问题转化

矩形周长并比面积并复杂，但仍然可以用扫描线解决。关键观察：周长并的水平边贡献和竖直边贡献可以分开计算。

对于水平边的贡献，扫描线从下到上扫描：

每遇到一个事件（矩形的下边或上边），水平方向的周长增量 = 当前覆盖的 y 区间被覆盖段数的变化量的绝对值 * 相关的 x 区间长度。

更具体地说：

水平周长贡献 = |当前覆盖的段数变化| × 2
竖直周长贡献 = 覆盖的总长度变化的绝对值 × 2

5.2 线段树需要额外维护的信息

与面积并相比，周长并的线段树需要额外维护：

num_seg：当前被覆盖区间的连续段数。
left_covered：区间最左端是否被覆盖。
right_covered：区间最右端是否被覆盖。

合并时：

parent.num_seg = left_child.num_seg + right_child.num_seg
    - (left_child.right_covered && right_child.left_covered ? 1 : 0)

5.3 算法流程

RECTANGLE-PERIMETER-UNION(rectangles):
    // 计算竖直边的贡献
    events ← 所有矩形的下边（+1）和上边（-1），按 y 排序
    prev_len ← 0
    perimeter ← 0

    for each event:
        更新线段树
        curr_len ← tree[1].len
        perimeter += |curr_len - prev_len|  // 竖直边贡献
        prev_len = curr_len

        // 水平边贡献 = num_seg * 2 * (下一个事件的 y - 当前 y)
        perimeter += tree[1].num_seg * 2 * (next_y - curr_y)

    // 类似地处理水平方向（或直接在同一次扫描中处理）
    return perimeter

5.4 复杂度

同样是 O(n log n)，但常数因子比面积并大一些，因为线段树节点需要维护更多信息。

六、最近点对的扫描线解法

6.1 经典分治 vs 扫描线

最近点对问题最著名的解法是分治算法，复杂度 O(n log n)。扫描线给出了另一种同样优美的 O(n log n) 解法。

6.2 算法思路

将所有点按 x 坐标排序。
维护一个集合 S（用平衡 BST 按 y 坐标排序），表示”活跃点”。
维护当前已知的最近距离 d（初始为正无穷）。
从左到右扫描每个点 p：
- 从 S 中删除所有 x 坐标与 p.x 相差超过 d 的点（它们不可能再贡献更近的距离）。
- 在 S 中查询 y 坐标在 [p.y - d, p.y + d] 范围内的所有点，计算与 p 的距离，更新 d。
- 将 p 插入 S。

6.3 关键性质

看起来第 4 步中”查询范围内的所有点”可能很多，但实际上：

在 d * 2d 的矩形区域内，最多只有 O(1) 个点（具体来说不超过 6 个）。

这是因为这些点两两距离至少为 d（否则 d 应该更小），而在 d * 2d 的矩形中最多放 6 个两两距离至少为 d 的点。

因此每次查询处理的点数是 O(1)，总复杂度 O(n log n)。

6.4 实现要点

/* 扫描线求最近点对——核心逻辑 */

typedef struct {
    double x, y;
} Point;

/* 假设 points 已按 x 排序 */
/* 使用一个按 y 排序的集合（此处用排序数组模拟） */

double closest_pair_sweep(Point *pts, int n)
{
    /* 按 x 排序 */
    /* ... (省略排序代码) ... */

    double d = 1e18;
    int left = 0;  /* 活跃窗口左边界 */

    /* active_set: 按 y 排序的活跃点集合 */
    /* 此处用简化实现 */

    for (int i = 0; i < n; i++) {
        /* 移除不再活跃的点 */
        while (pts[i].x - pts[left].x > d)
            left++;

        /* 在活跃集合中查找 y 范围内的点 */
        for (int j = left; j < i; j++) {
            double dy = pts[i].y - pts[j].y;
            if (dy > d) continue;
            if (dy < -d) continue;
            double dx = pts[i].x - pts[j].x;
            double dist = sqrt(dx * dx + dy * dy);
            if (dist < d) d = dist;
        }
    }
    return d;
}

注意：上面的简化实现在最坏情况下仍然是 O(n^2)。要达到严格的 O(n log n)，需要将活跃集合用平衡 BST 实现，支持 O(log n) 的范围查询。

七、数值稳定性与退化情况

扫描线算法在理论上非常优美，但在实现中充满了陷阱。数值问题是最大的敌人。

7.1 浮点精度问题

Bentley-Ottmann 算法中，线段在状态结构中的比较依赖于它们在当前扫描线位置的 y 值。这意味着：

比较函数使用浮点运算，可能因舍入误差导致不一致的比较结果。
BST 的不变量可能被破坏（如 a < b 且 b < c 但 a > c）。
交点计算的误差可能导致事件顺序错误。

/* 危险：直接比较浮点数 */
if (y_at_sweep(s1) < y_at_sweep(s2))  /* 可能因舍入而不稳定 */

/* 稍好：使用 epsilon */
if (y_at_sweep(s1) < y_at_sweep(s2) - EPS)  /* 但 EPS 取多少？ */

7.2 退化情况

以下退化情况需要特别处理：

多条线段共享端点：多个事件具有相同的 x 坐标。必须定义处理顺序（通常先处理左端点，再处理交点，最后处理右端点）。
竖直线段：x 坐标范围为零。需要特殊处理，通常将竖直线段微旋转一个极小角度。
多条线段交于一点：一个交点涉及多于两条线段。Bentley-Ottmann 原始算法假设交点只涉及两条线段；处理多重交点需要修改算法。
重叠线段：两条线段部分或完全重叠。需要作为特殊情况处理。
端点恰好在另一条线段上：T 型交叉。

7.3 解决方案

问题	解决方案	代价
浮点精度	使用精确算术（如 GMP/MPFR）	性能下降 10-100 倍
浮点精度	使用自适应精度（如 Shewchuk predicates）	实现复杂，但性能接近原生浮点
浮点精度	使用整数坐标 + 精确叉积	限制输入范围
共享端点	定义严格的事件处理优先级	代码复杂度增加
竖直线段	符号扰动（Symbolic perturbation）	理论上优雅，实现有技巧
多重交点	将事件扩展为处理一组线段	算法复杂度增加
重叠线段	预处理合并重叠段	需要额外的预处理步骤

7.4 符号扰动

符号扰动（Symbolic Perturbation，简称 SoS）是处理退化情况的一种通用技术：

假设每个坐标都加上了一个无穷小的扰动 epsilon_i，使得所有点都处于”一般位置”（general position）——没有三点共线、没有两点共 x 坐标等等。

在实现中，不需要真的加上扰动，只需要在比较函数中定义当原始值相等时的打破平局规则。例如，当两个点的 x 坐标相等时，比较它们的 y 坐标；当 y 坐标也相等时，比较它们的索引。

/* 符号扰动的比较函数 */
int compare_events(const Event *a, const Event *b)
{
    if (fabs(a->pt.x - b->pt.x) > EPS)
        return (a->pt.x < b->pt.x) ? -1 : 1;
    if (fabs(a->pt.y - b->pt.y) > EPS)
        return (a->pt.y < b->pt.y) ? -1 : 1;
    /* 打破平局：左端点 < 交点 < 右端点 */
    if (a->type != b->type)
        return a->type - b->type;
    /* 最终：按线段索引 */
    return a->seg1 - b->seg1;
}

八、扫描线在 EDA 中的应用

8.1 EDA 与几何运算

电子设计自动化（Electronic Design Automation，EDA）是扫描线算法最重要的工业应用领域。现代芯片设计中，一个芯片上可能有数十亿个几何图形（矩形、多边形），需要进行各种几何运算：

设计规则检查（DRC）：检查相邻导线间距是否满足工艺要求。
布尔运算：计算两层金属的交集、并集、差集。
面积计算：统计各层的覆盖面积。
连通性提取：确定哪些几何图形电气相连。

这些运算的输入规模通常在百万到十亿级别，暴力算法完全不可行。

8.2 EDA 中的扫描线特点

EDA 中的扫描线与学术版本有以下不同：

坐标是整数：芯片设计使用离散网格（通常以纳米为单位），所有坐标都是整数。这消除了浮点精度问题，是一个巨大的优势。
图形以矩形为主：虽然也有斜线（45 度或任意角度），但大部分图形是轴对齐矩形。这使得专用的矩形扫描线可以大幅优化。
规模极大：需要处理数十亿个图形。内存和缓存效率至关重要。
分区并行：将芯片分成若干区域，每个区域独立进行扫描线处理，最后合并边界结果。

8.3 DRC 中的间距检查

设计规则检查中最核心的操作是间距检查：

给定同一层上的所有矩形，检查是否存在两个不同矩形之间的最短距离小于规定的最小间距 d。

这可以用扫描线高效解决：

将每个矩形沿各方向膨胀 d/2，得到”膨胀矩形”。
对膨胀矩形做面积并（或交点检测）。
如果两个膨胀矩形相交但原矩形不相交，则存在间距违规。

更高效的做法是直接在扫描过程中维护邻近关系，避免显式膨胀。

8.4 布尔运算

两层几何图形的布尔运算（AND、OR、XOR、NOT）可以用统一的扫描线框架实现：

BOOLEAN-OP(polygons_A, polygons_B, op):
    events ← 合并 A 和 B 的所有边
    按 x 排序

    for each event:
        更新两棵状态结构（分别维护 A 和 B 的覆盖情况）
        根据 op 计算结果区域的边界
        输出结果多边形的边

在工业实现中，这通常称为”polygon clipping”或”polygon overlay”，是 EDA 工具链的核心算法之一。

8.5 性能优化技巧

工业级 EDA 扫描线的性能优化包括：

事件压缩：将同一 x 坐标的所有事件合并处理，减少状态切换次数。
分层处理：先粗粒度扫描确定热区，再细粒度扫描处理热区。
内存池：预分配内存，避免频繁的 malloc/free。
SIMD 加速：对大量矩形的坐标比较使用向量指令。
增量更新：设计修改后只重新扫描受影响的局部区域。

九、工程陷阱速查表

在实际实现扫描线算法时，以下是常见的错误和注意事项：

陷阱	症状	解决方案
浮点比较不一致	BST 损坏，段错误或无限循环	使用 epsilon 比较，或改用整数坐标
事件处理顺序错误	遗漏交点或重复报告	严格定义同 x 坐标事件的处理优先级
忘记处理竖直线段	除以零	特判或符号扰动
线段树区间端点理解错误	面积偏大或偏小	画图确认：节点 [l,r] 表示 ys[l] 到 ys[r] 的区间
离散化后忘记用原始坐标计算长度	面积完全错误	线段树的 `len` 使用原始浮点坐标而非索引
活跃表在交点处未正确交换	遗漏后续交点	仔细检查 swap 逻辑，确保新邻居对被检查
同一交点被多次插入事件队列	重复报告	在事件队列中去重（用集合而非堆）
内存分配不足	数组越界	交点数最多 O(n^2)，预留足够空间
未处理端点共线	产生虚假交点	检查交点是否在两条线段的公共端点上
数值误差累积	后期事件处理严重错误	使用自适应精度或定期重新计算

十、扫描线的变体与扩展

10.1 旋转扫描线

标准扫描线使用竖直线从左到右扫描。旋转扫描线（Rotational Sweep）则使用一条从某个定点出发的射线，绕定点旋转。

典型应用：

可见性问题：从一个点出发，哪些线段可见？
角度排序问题：按极角排序处理几何对象。

10.2 径向扫描线

径向扫描线使用一个不断扩大的圆，从某个中心点向外扩展。这在最近邻查询和 Voronoi 图构造中有应用。

10.3 平面扫描（Plane Sweep）的推广

扫描线可以推广到更高维度：

三维：用一个平面从左到右扫过三维空间，将三维问题降维为二维动态问题。但二维动态问题本身可能仍然很困难。
时间维度：将时间作为扫描维度，处理动态几何问题。

10.4 离线与在线

标准扫描线是离线算法——需要预先知道所有几何对象。在线版本需要支持动态插入/删除几何对象，通常使用动态线段树或持久化数据结构实现。

十一、完整应用示例：矩形面积并与周长并的综合实现

下面给出一个综合示例，同时计算矩形面积并和周长并。这展示了如何在一次扫描中提取多种信息。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXN 50005

/* -------- 数据结构 -------- */

typedef struct {
    double x;
    double y1, y2;
    int type;  /* +1 进入, -1 离开 */
} Event;

static Event events[MAXN * 2];
static int   n_events;
static double ys[MAXN * 2];
static int    n_ys;

/* 线段树节点 */
typedef struct {
    int    cnt;            /* 完整覆盖次数 */
    double len;            /* 被覆盖的总长度 */
    int    num_seg;        /* 不连续的被覆盖段数 */
    int    left_covered;   /* 左端点是否被覆盖 */
    int    right_covered;  /* 右端点是否被覆盖 */
} STNode;

static STNode tree[MAXN * 8];

/* -------- 排序与离散化 -------- */

static int event_cmp(const void *a, const void *b)
{
    const Event *ea = (const Event *)a;
    const Event *eb = (const Event *)b;
    if (ea->x < eb->x) return -1;
    if (ea->x > eb->x) return  1;
    return ea->type - eb->type;
}

static int dcmp(const void *a, const void *b)
{
    double d = *(const double *)a - *(const double *)b;
    return (d > 1e-9) - (d < -1e-9);
}

static int find_y(double v)
{
    int lo = 0, hi = n_ys - 1;
    while (lo <= hi) {
        int mid = (lo + hi) / 2;
        if (ys[mid] < v - 1e-9)      lo = mid + 1;
        else if (ys[mid] > v + 1e-9) hi = mid - 1;
        else return mid;
    }
    return -1;
}

/* -------- 线段树 -------- */

static void st_build(int nd, int l, int r)
{
    tree[nd] = (STNode){0, 0.0, 0, 0, 0};
    if (l + 1 >= r) return;
    int mid = (l + r) / 2;
    st_build(nd * 2, l, mid);
    st_build(nd * 2 + 1, mid, r);
}

static void st_push_up(int nd, int l, int r)
{
    if (tree[nd].cnt > 0) {
        tree[nd].len = ys[r] - ys[l];
        tree[nd].num_seg = 1;
        tree[nd].left_covered = 1;
        tree[nd].right_covered = 1;
    } else if (l + 1 >= r) {
        tree[nd].len = 0.0;
        tree[nd].num_seg = 0;
        tree[nd].left_covered = 0;
        tree[nd].right_covered = 0;
    } else {
        tree[nd].len = tree[nd * 2].len + tree[nd * 2 + 1].len;
        tree[nd].num_seg = tree[nd * 2].num_seg + tree[nd * 2 + 1].num_seg
            - (tree[nd * 2].right_covered && tree[nd * 2 + 1].left_covered ? 1 : 0);
        tree[nd].left_covered = tree[nd * 2].left_covered;
        tree[nd].right_covered = tree[nd * 2 + 1].right_covered;
    }
}

static void st_update(int nd, int l, int r, int ql, int qr, int val)
{
    if (ql >= r || qr <= l) return;
    if (ql <= l && r <= qr) {
        tree[nd].cnt += val;
        st_push_up(nd, l, r);
        return;
    }
    int mid = (l + r) / 2;
    st_update(nd * 2, l, mid, ql, qr, val);
    st_update(nd * 2 + 1, mid, r, ql, qr, val);
    st_push_up(nd, l, r);
}

/* -------- 主函数 -------- */

int main(void)
{
    int n;
    printf("输入矩形数量: ");
    scanf("%d", &n);

    n_events = 0;
    n_ys = 0;

    for (int i = 0; i < n; i++) {
        double x1, y1, x2, y2;
        scanf("%lf %lf %lf %lf", &x1, &y1, &x2, &y2);
        events[n_events++] = (Event){x1, y1, y2, +1};
        events[n_events++] = (Event){x2, y1, y2, -1};
        ys[n_ys++] = y1;
        ys[n_ys++] = y2;
    }

    qsort(events, n_events, sizeof(Event), event_cmp);
    qsort(ys, n_ys, sizeof(double), dcmp);
    int unique = 1;
    for (int i = 1; i < n_ys; i++)
        if (fabs(ys[i] - ys[i - 1]) > 1e-9)
            ys[unique++] = ys[i];
    n_ys = unique;

    st_build(1, 0, n_ys - 1);

    double area = 0.0;
    double perimeter = 0.0;
    double prev_len = 0.0;
    int    prev_seg = 0;

    for (int i = 0; i < n_events; i++) {
        if (i > 0) {
            double dx = events[i].x - events[i - 1].x;
            area += tree[1].len * dx;
            perimeter += prev_seg * 2.0 * dx;  /* 水平边贡献 */
        }
        int lo = find_y(events[i].y1);
        int hi = find_y(events[i].y2);
        st_update(1, 0, n_ys - 1, lo, hi, events[i].type);
        perimeter += fabs(tree[1].len - prev_len);  /* 竖直边贡献 */
        prev_len = tree[1].len;
        prev_seg = tree[1].num_seg;
    }

    printf("面积并 = %.2f\n", area);
    printf("周长并 = %.2f\n", perimeter);

    return 0;
}

这段代码在一次从左到右的扫描中同时计算了面积并和周长并。周长的计算分为两部分： - 竖直边贡献：每次事件更新后，覆盖长度的变化量就是新增的竖直边。 - 水平边贡献：两个相邻事件之间，每个被覆盖的连续段产生两条水平边。

十二、个人思考与总结

12.1 扫描线的思维模型

我个人认为，扫描线是计算几何中最值得深入理解的算法范式——不是因为它最难，而是因为它揭示了一种普遍的问题解决思路：当面对高维问题时，引入一个有序的扫描过程，将问题降维。 这个思路不仅适用于几何，也适用于数据库查询优化（sort-merge join）、事件驱动模拟、甚至编译器的寄存器分配（线性扫描分配器）。

12.2 理论与实践的鸿沟

Bentley-Ottmann 算法在理论上非常优美，O((n+k) log n) 的复杂度是最优的（在比较模型下）。但在实践中，它的实现难度和数值问题使得很多工程师宁愿使用更简单但更鲁棒的方法——比如基于网格的空间散列加暴力检查，或者使用 CGAL 这样经过工业验证的几何库。

我的建议是：

理解算法思想：即使不自己实现 Bentley-Ottmann，理解它的事件驱动和邻居检查机制对思考其他问题很有帮助。
使用成熟的库：在生产代码中，优先使用 CGAL、Boost.Geometry、Clipper 等经过充分测试的库。
注意数值问题：如果必须自己实现，整数坐标 + 精确谓词是最务实的选择。

12.3 竞赛 vs 工业

在算法竞赛中，矩形面积并和周长并是扫描线的经典题型，通常坐标范围有限（整数），不需要担心太多数值问题。重点是快速正确地实现线段树。

在工业应用（特别是 EDA）中，扫描线面对的是完全不同的挑战：超大规模数据、严格的正确性要求、复杂的退化情况。这要求算法实现不仅正确，还要高效、鲁棒、可维护。

12.4 扫描线的局限

扫描线并不万能。它的主要局限包括：

不适合在线问题：标准扫描线需要离线处理所有数据。
不适合曲线：对于非线性几何对象（圆弧、样条曲线），交点计算和排序变得复杂得多。
常数因子较大：由于需要维护复杂的数据结构，扫描线的实际运行时间可能不如简单的分治或随机化方法。
并行化困难：扫描线本质上是串行的——每一步依赖前一步的状态。虽然可以通过空间分区实现一定程度的并行，但不如分治算法那样天然适合并行。

12.5 推荐资源

de Berg 等《Computational Geometry: Algorithms and Applications》：第二章详细讲解扫描线和 Bentley-Ottmann，配图精美，证明严谨。
Preparata & Shamos《Computational Geometry: An Introduction》：经典教材，对扫描线范式有深入讨论。
CGAL 文档：工业级实现的参考。了解 CGAL 如何处理退化情况和数值问题非常有教育意义。
Shewchuk 的精确谓词论文：“Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates”——如果你需要自己处理数值问题，这是必读文献。

参考文献

Bentley, J.L., Ottmann, T.A. “Algorithms for Reporting and Counting Geometric Intersections.” IEEE Transactions on Computers, C-28(9), 1979.
Shamos, M.I., Hoey, D. “Geometric Intersection Problems.” 17th Annual Symposium on Foundations of Computer Science, 1976.
de Berg, M., Cheong, O., van Kreveld, M., Overmars, M. “Computational Geometry: Algorithms and Applications.” 3rd Edition, Springer, 2008.
Preparata, F.P., Shamos, M.I. “Computational Geometry: An Introduction.” Springer, 1985.
Shewchuk, J.R. “Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates.” Discrete & Computational Geometry, 18(3):305-363, 1997.
Nievergelt, J., Preparata, F.P. “Plane-Sweep Algorithms for Intersecting Geometric Figures.” Communications of the ACM, 25(10):739-747, 1982.
Greiner, G., Hormann, K. “Efficient Clipping of Arbitrary Polygons.” ACM Transactions on Graphics, 17(2):71-83, 1998.
CGAL: The Computational Geometry Algorithms Library. https://www.cgal.org/

算法系列导航：上一篇：凸包算法 | 下一篇：Voronoi 图与 Delaunay 三角剖分

相关阅读：线段树与树状数组 | R-tree 与空间索引

目录