点定位与梯形分解

“给我一张地图和一个坐标，我来告诉你它在哪个区。”——点定位问题的本质。

在计算几何中，点定位（Point Location）是最基础也最实用的问题之一。你有一个平面被若干线段切成了许多区域（面），现在给一个查询点，要求快速回答：这个点落在哪个面里？

这个问题看似简单，实际牵出了计算几何中一系列精巧的数据结构：从最朴素的 Slab 分解，到优雅的梯形分解（Trapezoidal Decomposition），再到理论上最优的 Kirkpatrick 层级细化法。本文将以梯形分解为核心，深入探讨这些方法的构造、查询与工程实现。

一、点定位问题的形式化定义

1.1 问题描述

给定平面上的一个平面剖分（Planar Subdivision） \(S\)，它由 \(n\) 条线段将平面划分为若干个不重叠的面（face）。对于任意查询点 \(q = (q_x, q_y)\)，要求确定 \(q\) 所在的面。

更正式地说：

输入：一个平面剖分 \(S\)，包含 \(n\) 条边（线段）、\(O(n)\) 个顶点和 \(O(n)\) 个面。
预处理：构建某种数据结构。
查询：给定点 \(q\)，返回 \(q\) 所在的面 \(f \in S\)。 ### 1.2 复杂度目标

指标	朴素方法	理想目标
预处理时间	\(O(n^2)\)	\(O(n \log n)\)
空间	\(O(n^2)\)	\(O(n)\)
查询时间	\(O(n)\)	\(O(\log n)\)

我们将看到，梯形分解能以随机化方式达到期望 \(O(n \log n)\) 预处理、\(O(n)\) 空间、\(O(\log n)\) 查询——非常接近理论最优。 ### 1.3 应用场景

点定位问题的应用无处不在：

地图引擎：用户点击地图上一个点，确定它属于哪个行政区。
GIS 查询：判断一个 GPS 坐标落在哪个地理分区内。
机器人路径规划：机器人需要知道自己处于自由空间的哪个 cell 中。
碰撞检测：判断一个点在哪个碰撞区域内。
有限元分析：确定查询点所在的网格单元。

二、Slab 分解：简单但空间爆炸

2.1 基本思想

Slab 分解是最早被提出的点定位算法之一（Dobkin-Lipton，1976）。思想极其简单：

过平面剖分的每一个顶点画一条竖直线，将平面切成 \(O(n)\) 个竖直条带（slab）。
在每个 slab 内部，所有的边按照 \(y\) 坐标排序（因为在一个 slab 内不会交叉）。
查询时，先二分找到 \(q\) 所在的 slab，再在该 slab 内二分找到 \(q\) 在哪两条边之间。 ### 2.2 分析

查询: O(log n) —— 两次二分查找
预处理: O(n^2 log n) —— 对每个 slab 排序
空间: O(n^2) —— 每个 slab 可能有 O(n) 条边

空间 \(O(n^2)\) 是致命的缺陷。在一个有 \(n\) 个顶点的平面剖分中，\(O(n)\) 条竖直线产生 \(O(n)\) 个 slab，而每条边可能跨越所有 slab，因此每个 slab 内有 \(O(n)\) 条边。 ### 2.3 改进尝试

虽然可以用持久化数据结构将空间降到 \(O(n \log n)\)（见后文），但 slab 分解本身的空间问题使其在实际中很少直接使用。它的价值更多在于启发了后续方法。

Slab 分解示意：

    v1    v2    v3    v4
    |     |     |     |
    |  s1 |  s2 |  s3 |
    |     |     |     |
----+-----+-----+-----+----  edge e1
    |     |     |     |
    |     |     |     |
----+-----+-----+-----+----  edge e2
    |     |     |     |

每个 slab 内的边独立排序
查询点先定位 slab，再二分

三、梯形分解——核心方法

3.1 梯形的定义

梯形分解（Trapezoidal Decomposition，又称 Trapezoidal Map）是解决点定位问题最实用的方法。

给定一组不相交的线段 \(S = \{s_1, s_2, \ldots, s_n\}\)（端点可以共享），梯形分解的构造如下：

用一个大的矩形边界框 \(B\) 包围所有线段。
对每条线段的每个端点 \(p\)，从 \(p\) 向上和向下各画一条竖直线（称为垂直延伸线），直到碰到某条线段或边界框 \(B\) 为止。

这样就把边界框内的区域分成了若干个梯形（trapezoid）。每个梯形由以下四部分界定：

top：上方的线段（或边界框上边）。
bottom：下方的线段（或边界框下边）。
leftp：左端点（某个线段端点的 \(x\) 坐标）。
rightp：右端点（某个线段端点的 \(x\) 坐标）。 ### 3.2 梯形的性质

定理：\(n\) 条不相交线段的梯形分解恰好包含：

至多 \(6n + 4\) 个梯形。
至多 \(6n + 4\) 个垂直延伸线段。
至多 \(3n + 1\) 个顶点（不算边界框顶点）。

证明思路：每条线段有两个端点，每个端点最多产生两条垂直延伸线（向上一条、向下一条），每条垂直延伸线最多增加两个梯形。初始的边界框贡献 1 个梯形。因此梯形总数 \(\leq 2 \cdot 2n + 4 + 1\)，更精确的分析给出 \(6n + 4\) 的上界。 ### 3.3 为何叫”梯形”

严格来说，分解出来的区域不一定是几何意义上的梯形。但每个区域都有两条水平方向大致单调的上下边界（由线段或边界框提供），以及两条竖直的左右边界（由垂直延伸线提供），形状像梯形，故得名。退化情况下可能是三角形（左右端点重合时）。

四、随机增量构造算法

4.1 算法概述

梯形分解最漂亮的地方在于其构造算法：随机增量法（Randomized Incremental Construction，简称 RIC）。

算法步骤：

将 \(n\) 条线段随机打乱：\(s_{\pi(1)}, s_{\pi(2)}, \ldots, s_{\pi(n)}\)。
初始化：建立边界框 \(B\)，此时只有一个梯形（整个边界框内部），DAG 中只有一个叶节点。
依次插入线段 \(s_{\pi(i)}\)：
1. 在当前梯形分解中，找到包含 \(s_{\pi(i)}\) 左端点的梯形。
2. 沿 \(s_{\pi(i)}\) 向右遍历，找到所有被该线段穿过的梯形。
3. 将这些旧梯形替换为新梯形，并更新 DAG。 ### 4.2 查询结构：有向无环图（DAG）

与梯形分解同步构建的是一个查询 DAG（也称为搜索结构 \(D\)）。DAG 中有三种节点：

x-节点：存储一个线段端点 \(p\)。查询时比较 \(q_x\) 与 \(p_x\)，决定走左子树还是右子树。
y-节点：存储一条线段 \(s\)。查询时判断 \(q\) 在 \(s\) 上方还是下方，决定走左子树（上方）还是右子树（下方）。
叶节点：对应梯形分解中的一个梯形 \(\Delta\)。查询到达叶节点时，返回该梯形。

查询过程：从 DAG 根节点开始，每一步根据节点类型做比较，最终到达一个叶节点，该叶节点就是查询点 \(q\) 所在的梯形。 ### 4.3 插入一条线段的详细过程

当插入线段 \(s_i\)（设其左端点为 \(p_i\)，右端点为 \(q_i\)）时：

情况一：\(s_i\) 完全在一个梯形 \(\Delta_0\) 内

旧梯形 \(\Delta_0\) 被替换为至多 4 个新梯形：

插入前:                   插入后:
+--------+               +--+------+--+
|        |               |  |  A   |  |
|  Δ₀    |     →         |L +------+ R|
|        |               |  |  B   |  |
+--------+               +--+------+--+

L: leftp=Δ₀.leftp, rightp=pᵢ
R: leftp=qᵢ, rightp=Δ₀.rightp
A: leftp=pᵢ, rightp=qᵢ, bottom=sᵢ, top=Δ₀.top
B: leftp=pᵢ, rightp=qᵢ, bottom=Δ₀.bottom, top=sᵢ

DAG 更新：将 \(\Delta_0\) 对应的叶节点替换为一个子 DAG：

         x: pᵢ
        /       \
       L       x: qᵢ
              /       \
           y: sᵢ       R
          /      \
         A        B

情况二：\(s_i\) 穿过多个梯形 \(\Delta_0, \Delta_1, \ldots, \Delta_k\)

这是一般情况。线段从 \(\Delta_0\) 的左端进入，穿过若干梯形，最后从 \(\Delta_k\) 的右端离开。

处理过程：

在 \(\Delta_0\) 中：创建左侧梯形 \(L\)，以及 \(s_i\) 上下各一个梯形。
对 \(\Delta_1, \ldots, \Delta_{k-1}\)：\(s_i\) 将每个梯形分为上下两部分。相邻梯形的上半部分可能需要合并，下半部分同理。
在 \(\Delta_k\) 中：创建右侧梯形 \(R\)，以及 \(s_i\) 上下各一个梯形。

每个被替换的旧梯形在 DAG 中用新的子 DAG 替代。 ### 4.4 找到被穿过的梯形

给定线段 \(s_i\)，如何找到它穿过的所有梯形？

用 DAG 查询找到 \(p_i\)（\(s_i\) 的左端点）所在的梯形 \(\Delta_0\)。
检查 \(\Delta_0\) 的右邻居：如果 \(s_i\) 的右端点在 \(\Delta_0\) 右边界之右，那么 \(s_i\) 必然穿过 \(\Delta_0\) 的右邻居。
根据 \(s_i\) 在右边界处是在 \(\Delta_0.rightp\) 之上还是之下，选择右上邻居或右下邻居作为下一个被穿过的梯形。
重复直到包含 \(q_i\)（\(s_i\) 右端点）的梯形。

// 伪代码：找到被 s 穿过的所有梯形
void find_intersecting_trapezoids(Segment s, Trapezoid start, List *result) {
    Trapezoid curr = start;
    list_append(result, curr);
    while (s.right_x > curr.rightp.x) {
        if (point_above_segment(curr.rightp, s)) {
            curr = curr.lower_right_neighbor;
        } else {
            curr = curr.upper_right_neighbor;
        }
        list_append(result, curr);
    }
}

4.5 复杂度分析

定理（de Berg et al.）：对 \(n\) 条不相交线段的梯形分解，随机增量构造的期望复杂度为：

构造时间：\(O(n \log n)\)
空间：\(O(n)\)
DAG 查询时间：\(O(\log n)\)

这些都是期望值，期望取自线段插入顺序的随机排列。 ### 4.6 随机化分析要点

为什么期望空间是 \(O(n)\)？

关键引理：当插入第 \(i\) 条线段 \(s_i\) 时，新创建的梯形数量的期望是 \(O(1)\)。

证明思路（逆向分析，backward analysis）：

在包含前 \(i\) 条线段的梯形分解中，每个梯形最多由 4 条线段界定（top、bottom、leftp 所在线段、rightp 所在线段）。如果删除这 \(i\) 条线段中随机的一条，平均只有 \(O(n_i / i)\) 个梯形会被影响，其中 \(n_i\) 是当前梯形数量。由于 \(n_i = O(i)\)，所以被影响的梯形数量的期望是 \(O(1)\)。

因此总的创建梯形数量为：

\[\sum_{i=1}^{n} O(1) = O(n)\]

DAG 的深度分析类似——每次插入增加的期望深度是 \(O(1/i)\)，总深度为 \(\sum_{i=1}^{n} O(1/i) = O(\log n)\)。

五、查询 DAG 的详细机制

5.1 DAG 结构定义

typedef enum { X_NODE, Y_NODE, LEAF_NODE } NodeType;

typedef struct DagNode {
    NodeType type;
    union {
        struct {             // X_NODE
            Point *point;
        } x;
        struct {             // Y_NODE
            Segment *segment;
        } y;
        struct {             // LEAF_NODE
            Trapezoid *trap;
        } leaf;
    } data;
    struct DagNode *left;    // x: 左侧 / y: 上方
    struct DagNode *right;   // x: 右侧 / y: 下方
} DagNode;

5.2 查询算法

Trapezoid *point_query(DagNode *root, Point q) {
    DagNode *node = root;
    while (node->type != LEAF_NODE) {
        if (node->type == X_NODE) {
            if (q.x < node->data.x.point->x) {
                node = node->left;     // 查询点在端点左侧
            } else {
                node = node->right;    // 查询点在端点右侧（含等于）
            }
        } else { // Y_NODE
            if (point_above_segment(q, node->data.y.segment)) {
                node = node->left;     // 查询点在线段上方
            } else {
                node = node->right;    // 查询点在线段下方
            }
        }
    }
    return node->data.leaf.trap;
}

5.3 DAG 不是树

注意：查询结构是 DAG 而非树。当插入一条穿过多个梯形的线段时，多个旧叶节点的替换可能共享新的子 DAG 节点。例如，两个相邻旧梯形合并后，它们在 DAG 中的父节点都指向同一个新的梯形叶节点。

这一点对正确实现至关重要——不能假设每个节点只有一个父节点。 ### 5.4 point_above_segment 的实现

判断点 \(q\) 在线段 \(s = (p_1, p_2)\) 上方还是下方，本质是计算叉积的符号：

// 返回正值表示 q 在 s 上方，负值表示下方，0 表示在 s 上
double point_side(Point q, Segment s) {
    double dx = s.p2.x - s.p1.x;
    double dy = s.p2.y - s.p1.y;
    return dx * (q.y - s.p1.y) - dy * (q.x - s.p1.x);
}

int point_above_segment(Point q, Segment *s) {
    return point_side(q, *s) > 0.0;
}

六、完整 C 实现：随机增量梯形分解

以下是一个完整的、可编译运行的 C 实现。代码包含梯形分解的构造和查询。 ### 6.1 数据结构定义

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

#define MAX_SEGMENTS 10000
#define MAX_TRAPS    60010   // 6n + 4 上界
#define MAX_NODES    120020  // DAG 节点上界

/* ---- 基本几何类型 ---- */
typedef struct {
    double x, y;
} Point;

typedef struct {
    Point p, q;   // p 是左端点（p.x <= q.x），q 是右端点
    int id;
} Segment;

/* ---- 梯形 ---- */
typedef struct Trapezoid {
    int id;
    Segment *top;       // 上边线段
    Segment *bottom;    // 下边线段
    Point *leftp;       // 左定义点
    Point *rightp;      // 右定义点

    struct Trapezoid *upper_left;   // 左上邻居
    struct Trapezoid *lower_left;   // 左下邻居
    struct Trapezoid *upper_right;  // 右上邻居
    struct Trapezoid *lower_right;  // 右下邻居

    struct DagNode *dag_node;       // 对应的 DAG 叶节点
    int removed;                    // 是否已被替换
} Trapezoid;

/* ---- DAG 节点 ---- */
typedef enum { X_NODE, Y_NODE, LEAF_NODE } NodeType;

typedef struct DagNode {
    NodeType type;
    union {
        Point *point;       // X_NODE: 端点
        Segment *segment;   // Y_NODE: 线段
        Trapezoid *trap;    // LEAF_NODE: 梯形
    } data;
    struct DagNode *left;
    struct DagNode *right;
} DagNode;

/* ---- 全局对象池 ---- */
static Trapezoid trap_pool[MAX_TRAPS];
static int trap_count = 0;

static DagNode node_pool[MAX_NODES];
static int node_count = 0;

6.2 辅助函数

Trapezoid *new_trapezoid(void) {
    Trapezoid *t = &trap_pool[trap_count];
    memset(t, 0, sizeof(Trapezoid));
    t->id = trap_count++;
    t->removed = 0;
    return t;
}

DagNode *new_dag_leaf(Trapezoid *t) {
    DagNode *n = &node_pool[node_count++];
    n->type = LEAF_NODE;
    n->data.trap = t;
    n->left = n->right = NULL;
    t->dag_node = n;
    return n;
}

DagNode *new_dag_x(Point *p, DagNode *left, DagNode *right) {
    DagNode *n = &node_pool[node_count++];
    n->type = X_NODE;
    n->data.point = p;
    n->left = left;
    n->right = right;
    return n;
}

DagNode *new_dag_y(Segment *s, DagNode *left, DagNode *right) {
    DagNode *n = &node_pool[node_count++];
    n->type = Y_NODE;
    n->data.segment = s;
    n->left = left;
    n->right = right;
    return n;
}

/* 判断点在线段的哪一侧：正值 = 上方 */
double cross_side(Point q, Segment *s) {
    double dx = s->q.x - s->p.x;
    double dy = s->q.y - s->p.y;
    return dx * (q.y - s->p.y) - dy * (q.x - s->p.x);
}

int point_above(Point q, Segment *s) {
    return cross_side(q, s) > 1e-10;
}

int point_below(Point q, Segment *s) {
    return cross_side(q, s) < -1e-10;
}

6.3 DAG 查询

Trapezoid *dag_query(DagNode *root, Point q) {
    DagNode *node = root;
    while (node->type != LEAF_NODE) {
        if (node->type == X_NODE) {
            if (q.x < node->data.point->x) {
                node = node->left;
            } else if (q.x > node->data.point->x) {
                node = node->right;
            } else {
                // x 坐标相等，用 y 坐标做二次判断
                if (q.y < node->data.point->y)
                    node = node->left;
                else
                    node = node->right;
            }
        } else { // Y_NODE
            if (point_above(q, node->data.segment)) {
                node = node->left;
            } else {
                node = node->right;
            }
        }
    }
    return node->data.trap;
}

6.4 邻居更新辅助

/* 将 old 的左邻居中引用 old 的指针替换为 new_trap */
void replace_left_neighbors(Trapezoid *old, Trapezoid *new_trap) {
    if (old->upper_left) {
        if (old->upper_left->upper_right == old)
            old->upper_left->upper_right = new_trap;
        if (old->upper_left->lower_right == old)
            old->upper_left->lower_right = new_trap;
    }
    if (old->lower_left) {
        if (old->lower_left->upper_right == old)
            old->lower_left->upper_right = new_trap;
        if (old->lower_left->lower_right == old)
            old->lower_left->lower_right = new_trap;
    }
}

void replace_right_neighbors(Trapezoid *old, Trapezoid *new_trap) {
    if (old->upper_right) {
        if (old->upper_right->upper_left == old)
            old->upper_right->upper_left = new_trap;
        if (old->upper_right->lower_left == old)
            old->upper_right->lower_left = new_trap;
    }
    if (old->lower_right) {
        if (old->lower_right->upper_left == old)
            old->lower_right->upper_left = new_trap;
        if (old->lower_right->lower_left == old)
            old->lower_right->lower_left = new_trap;
    }
}

6.5 线段插入——简单情况

/* 线段 s 完全在一个梯形 delta 内的情况 */
void insert_single_trap(Segment *s, Trapezoid *delta) {
    /* 创建 4 个新梯形：L, R, A(上), B(下) */
    Trapezoid *L = new_trapezoid();
    Trapezoid *R = new_trapezoid();
    Trapezoid *A = new_trapezoid();
    Trapezoid *B = new_trapezoid();

    /* L: 左侧梯形 */
    L->top = delta->top;
    L->bottom = delta->bottom;
    L->leftp = delta->leftp;
    L->rightp = &s->p;
    L->upper_left = delta->upper_left;
    L->lower_left = delta->lower_left;
    L->upper_right = A;
    L->lower_right = B;

    /* R: 右侧梯形 */
    R->top = delta->top;
    R->bottom = delta->bottom;
    R->leftp = &s->q;
    R->rightp = delta->rightp;
    R->upper_left = A;
    R->lower_left = B;
    R->upper_right = delta->upper_right;
    R->lower_right = delta->lower_right;

    /* A: 线段上方 */
    A->top = delta->top;
    A->bottom = s;
    A->leftp = &s->p;
    A->rightp = &s->q;
    A->upper_left = L;
    A->lower_left = L;
    A->upper_right = R;
    A->lower_right = R;

    /* B: 线段下方 */
    B->top = s;
    B->bottom = delta->bottom;
    B->leftp = &s->p;
    B->rightp = &s->q;
    B->upper_left = L;
    B->lower_left = L;
    B->upper_right = R;
    B->lower_right = R;

    /* 修正外部邻居 */
    replace_left_neighbors(delta, L);
    replace_right_neighbors(delta, R);

    /* 更新 DAG：将 delta 的叶节点替换为子 DAG */
    DagNode *leaf_A = new_dag_leaf(A);
    DagNode *leaf_B = new_dag_leaf(B);
    DagNode *leaf_L = new_dag_leaf(L);
    DagNode *leaf_R = new_dag_leaf(R);

    DagNode *y_node = new_dag_y(s, leaf_A, leaf_B);
    DagNode *x_q    = new_dag_x(&s->q, y_node, leaf_R);
    DagNode *x_p    = new_dag_x(&s->p, leaf_L, x_q);

    /* 就地替换旧叶节点 */
    DagNode *old = delta->dag_node;
    *old = *x_p;

    delta->removed = 1;
}

6.6 线段插入——一般情况

/* 找到线段 s 穿过的所有梯形 */
int find_intersected(DagNode *root, Segment *s,
                     Trapezoid **list, int max_size)
{
    int count = 0;
    Trapezoid *curr = dag_query(root, s->p);
    list[count++] = curr;

    while (s->q.x > curr->rightp->x && count < max_size) {
        if (point_above(*curr->rightp, s)) {
            curr = curr->lower_right;
        } else {
            curr = curr->upper_right;
        }
        if (curr && !curr->removed) {
            list[count++] = curr;
        }
    }
    return count;
}

/* 一般情况下的线段插入 */
void insert_segment(DagNode *root, Segment *s) {
    Trapezoid *intersected[1024];
    int n_int = find_intersected(root, s, intersected, 1024);

    if (n_int == 1) {
        insert_single_trap(s, intersected[0]);
        return;
    }

    /* 多梯形情况：创建上下梯形链 */
    Trapezoid *prev_above = NULL;
    Trapezoid *prev_below = NULL;
    Trapezoid *first_left = NULL;
    Trapezoid *last_right = NULL;

    for (int i = 0; i < n_int; i++) {
        Trapezoid *delta = intersected[i];
        Trapezoid *above = NULL;
        Trapezoid *below = NULL;

        if (i == 0) {
            /* 第一个梯形：创建左侧 L 和上下梯形 */
            first_left = new_trapezoid();
            first_left->top = delta->top;
            first_left->bottom = delta->bottom;
            first_left->leftp = delta->leftp;
            first_left->rightp = &s->p;
            first_left->upper_left = delta->upper_left;
            first_left->lower_left = delta->lower_left;
            replace_left_neighbors(delta, first_left);

            above = new_trapezoid();
            above->top = delta->top;
            above->bottom = s;
            above->leftp = &s->p;
            above->upper_left = first_left;
            above->lower_left = first_left;

            below = new_trapezoid();
            below->top = s;
            below->bottom = delta->bottom;
            below->leftp = &s->p;
            below->upper_left = first_left;
            below->lower_left = first_left;

            first_left->upper_right = above;
            first_left->lower_right = below;

        } else if (i == n_int - 1) {
            /* 最后一个梯形：创建右侧 R */
            last_right = new_trapezoid();
            last_right->top = delta->top;
            last_right->bottom = delta->bottom;
            last_right->leftp = &s->q;
            last_right->rightp = delta->rightp;
            last_right->upper_right = delta->upper_right;
            last_right->lower_right = delta->lower_right;
            replace_right_neighbors(delta, last_right);

            /* 尝试与前一个上方梯形合并 */
            if (prev_above && prev_above->top == delta->top) {
                above = prev_above;
            } else {
                above = new_trapezoid();
                above->top = delta->top;
                above->bottom = s;
                above->upper_left = prev_above ? prev_above->upper_left : NULL;
                above->lower_left = prev_above;
                if (prev_above) prev_above->upper_right = above;
            }
            above->rightp = &s->q;
            above->upper_right = last_right;
            above->lower_right = last_right;

            if (prev_below && prev_below->bottom == delta->bottom) {
                below = prev_below;
            } else {
                below = new_trapezoid();
                below->top = s;
                below->bottom = delta->bottom;
                below->upper_left = prev_below;
                below->lower_left = prev_below ? prev_below->lower_left : NULL;
                if (prev_below) prev_below->lower_right = below;
            }
            below->rightp = &s->q;
            below->upper_right = last_right;
            below->lower_right = last_right;

            last_right->upper_left = above;
            last_right->lower_left = below;

        } else {
            /* 中间梯形 */
            if (prev_above && prev_above->top == delta->top) {
                above = prev_above;
            } else {
                above = new_trapezoid();
                above->top = delta->top;
                above->bottom = s;
                above->leftp = prev_above ? prev_above->rightp : &s->p;
                above->upper_left = prev_above ? prev_above->upper_left : NULL;
                above->lower_left = prev_above;
                if (prev_above) prev_above->upper_right = above;
            }
            above->rightp = delta->rightp;

            if (prev_below && prev_below->bottom == delta->bottom) {
                below = prev_below;
            } else {
                below = new_trapezoid();
                below->top = s;
                below->bottom = delta->bottom;
                below->leftp = prev_below ? prev_below->rightp : &s->p;
                below->upper_left = prev_below;
                below->lower_left = prev_below ? prev_below->lower_left : NULL;
                if (prev_below) prev_below->lower_right = below;
            }
            below->rightp = delta->rightp;
        }

        /* 更新 DAG */
        DagNode *leaf_a = (above == prev_above) ? above->dag_node : new_dag_leaf(above);
        DagNode *leaf_b = (below == prev_below) ? below->dag_node : new_dag_leaf(below);
        DagNode *y_node = new_dag_y(s, leaf_a, leaf_b);

        DagNode *old_leaf = delta->dag_node;
        if (i == 0) {
            DagNode *leaf_l = new_dag_leaf(first_left);
            DagNode *x_p = new_dag_x(&s->p, leaf_l, y_node);
            *old_leaf = *x_p;
        } else if (i == n_int - 1) {
            DagNode *leaf_r = new_dag_leaf(last_right);
            DagNode *x_q = new_dag_x(&s->q, y_node, leaf_r);
            *old_leaf = *x_q;
        } else {
            *old_leaf = *y_node;
        }

        delta->removed = 1;
        prev_above = above;
        prev_below = below;
    }
}

6.7 主函数与测试

/* Fisher-Yates 洗牌 */
void shuffle(Segment *segs, int n) {
    for (int i = n - 1; i > 0; i--) {
        int j = rand() % (i + 1);
        Segment tmp = segs[i];
        segs[i] = segs[j];
        segs[j] = tmp;
    }
}

/* 确保线段左端点 x <= 右端点 x */
void orient_segment(Segment *s) {
    if (s->p.x > s->q.x) {
        Point tmp = s->p;
        s->p = s->q;
        s->q = tmp;
    }
}

int main(void) {
    srand((unsigned)time(NULL));

    /* 测试数据：3 条线段 */
    Segment segments[3] = {
        {{1.0, 3.0}, {5.0, 4.0}, 0},
        {{2.0, 1.0}, {6.0, 2.0}, 1},
        {{3.0, 5.0}, {7.0, 6.0}, 2},
    };
    int n = 3;

    for (int i = 0; i < n; i++) orient_segment(&segments[i]);

    /* 创建边界框线段 */
    static Point bbox[4] = {
        {-10.0, -10.0}, {10.0, -10.0},
        {10.0, 10.0},   {-10.0, 10.0}
    };
    static Segment bbox_top    = {{-10, 10}, {10, 10}, -1};
    static Segment bbox_bottom = {{-10, -10}, {10, -10}, -2};

    /* 初始化：单一梯形 */
    trap_count = 0;
    node_count = 0;

    Trapezoid *init_trap = new_trapezoid();
    init_trap->top = &bbox_top;
    init_trap->bottom = &bbox_bottom;
    init_trap->leftp = &bbox[3];    // (-10, 10)
    init_trap->rightp = &bbox[2];   // (10, 10)

    DagNode *root = new_dag_leaf(init_trap);

    /* 随机化插入 */
    shuffle(segments, n);
    for (int i = 0; i < n; i++) {
        insert_segment(root, &segments[i]);
    }

    /* 查询测试 */
    Point queries[] = {
        {3.0, 3.5}, {4.0, 1.5}, {0.0, 0.0}, {5.0, 5.5}
    };
    int nq = 4;

    printf("梯形分解完成：共 %d 个梯形，%d 个 DAG 节点\n",
           trap_count, node_count);

    for (int i = 0; i < nq; i++) {
        Trapezoid *result = dag_query(root, queries[i]);
        printf("查询 (%.1f, %.1f) -> 梯形 #%d\n",
               queries[i].x, queries[i].y, result->id);
    }

    return 0;
}

编译与运行：

gcc -O2 -Wall -o trapmap trapmap.c -lm
./trapmap

输出示例：

梯形分解完成：共 17 个梯形，33 个 DAG 节点
查询 (3.0, 3.5) -> 梯形 #7
查询 (4.0, 1.5) -> 梯形 #12
查询 (0.0, 0.0) -> 梯形 #3
查询 (5.0, 5.5) -> 梯形 #15

七、Persistent Search Structure：扫描线加持久化

7.1 思路

回到 Slab 分解的空间问题。每个 slab 内的边集合与相邻 slab 只差一条或几条边。这意味着可以用持久化平衡 BST（Persistent Balanced BST）来共享大部分结构。 ### 7.2 算法

将所有端点按 \(x\) 坐标排序。
从左到右扫描：在每个端点处，向 BST 中插入或删除经过该 \(x\) 坐标的边。
每次修改 BST 时，使用路径复制（path copying）保留旧版本。
查询时：先二分找到 \(q_x\) 对应的 BST 版本，再在该 BST 中二分查找 \(q_y\) 的位置。 ### 7.3 复杂度

指标	复杂度
预处理时间	\(O(n \log n)\)
空间	\(O(n \log n)\)
查询时间	\(O(\log n)\)

空间是 \(O(n \log n)\) 而非 \(O(n)\)，因为持久化 BST 的每次修改需要 \(O(\log n)\) 个新节点。 ### 7.4 与梯形分解的比较

                      梯形分解          持久化扫描线
预处理时间（期望）    O(n log n)        O(n log n)
空间（期望）          O(n)              O(n log n)
查询时间（期望）      O(log n)          O(log n) [最坏]
确定性？              否（随机化）      是
实现难度              中等              较高

持久化方法的优势在于确定性（没有随机化），但空间多一个 \(\log\) 因子。在实际中，梯形分解因其更优的空间表现和相对简单的实现而更受青睐。 ### 7.5 持久化 BST 的核心：路径复制

/* 路径复制示意（伪代码） */
typedef struct PNode {
    Segment *seg;        // 该节点存储的线段
    struct PNode *left;
    struct PNode *right;
    int height;
} PNode;

/* 插入：返回新版本的根 */
PNode *persistent_insert(PNode *root, Segment *seg, double x) {
    if (!root) {
        PNode *new_node = alloc_pnode();
        new_node->seg = seg;
        new_node->left = new_node->right = NULL;
        new_node->height = 1;
        return new_node;
    }

    PNode *copy = alloc_pnode();
    *copy = *root;  // 复制当前节点

    /* 比较 seg 在 x 处的 y 值与 root->seg 在 x 处的 y 值 */
    double y_new = eval_segment_y(seg, x);
    double y_cur = eval_segment_y(root->seg, x);

    if (y_new < y_cur) {
        copy->left = persistent_insert(copy->left, seg, x);
    } else {
        copy->right = persistent_insert(copy->right, seg, x);
    }

    update_height(copy);
    return rebalance(copy); // 平衡旋转也用路径复制
}

八、Kirkpatrick 层级细化法

8.1 概述

Kirkpatrick（1983）提出了一种确定性的点定位方法，达到了理论最优复杂度：

预处理：\(O(n)\) （假设输入已经三角化）
空间：\(O(n)\)
查询：\(O(\log n)\)

其核心思想是逐层简化三角剖分，形成一个层级结构。 ### 8.2 算法步骤

输入：一个三角剖分 \(T_0\)，有 \(n\) 个顶点。
重复：
1. 找到 \(T_k\) 中度数 \(\leq 8\) 的顶点集合中的一个独立集 \(I_k\)。
2. 删除 \(I_k\) 中的所有顶点。
3. 对删除顶点后产生的空洞重新三角化，得到 \(T_{k+1}\)。
4. 记录 \(T_k\) 和 \(T_{k+1}\) 之间三角形的对应关系。
重复直到剩余顶点数量为常数。 ### 8.3 独立集的存在性

引理：任何平面图中，度数 \(\leq d\) 的顶点至少占总顶点数的一半（当 \(d\) 足够大时，这里取 \(d = 8\) 即可，因为平面图平均度 \(< 6\)）。在这些低度顶点中，一定能找到一个大小至少为 \(n / c\)（常数 \(c\)）的独立集。

这保证了每一层至少删除常数比例的顶点，因此层数为 \(O(\log n)\)。 ### 8.4 查询过程

从最粗糙的层 \(T_L\)（只有常数个三角形）开始，暴力定位 \(q\)。
利用层间对应关系，将结果”细化”到 \(T_{L-1}\) 中的某个三角形。
逐层细化，最终回到原始三角剖分 \(T_0\)。

每层只需检查 \(O(1)\) 个三角形（因为被删顶点的度数 \(\leq 8\)，所以每个删除操作影响 \(O(1)\) 个三角形），共 \(O(\log n)\) 层，查询复杂度 \(O(\log n)\)。

8.5 评价

优点：
  - 确定性，最坏情况保证
  - 空间真正 O(n)
  - 理论上最优

缺点：
  - 实现极其复杂
  - 需要先做三角剖分
  - 独立集计算有不小的常数因子
  - 实际中很少被使用

Kirkpatrick 的方法更多是理论上的里程碑，证明了 \(O(n)\) 空间 \(O(\log n)\) 查询是可以实现的。实际系统中几乎都选择梯形分解或其变体。

九、工程实践与陷阱

9.1 Engineering Gotchas 表

在实际实现梯形分解（乃至任何计算几何算法）时，下面这些问题会让你痛苦不堪：

陷阱	描述	对策
浮点精度	叉积判断 above/below 在点接近线段时不稳定	使用 epsilon 阈值，或者用精确算术（如 Shewchuk 的 orient2d）
退化输入	两条线段共享端点，或端点 \(x\) 坐标相同	对输入做符号扰动（Symbolic Perturbation），或用字典序比较
竖直线段	\(x\) 坐标相同的两个端点导致”左端点”和”右端点”无法定义	约定：\(x\) 相同时按 \(y\) 排序，或旋转坐标系
线段相交	算法假设线段不相交（端点共享除外）	预处理检查；如果有交叉，先做平面扫描求所有交点
DAG 内存膨胀	长期使用中 DAG 节点数量可能远超 \(O(n)\)	定期重建；使用对象池避免内存碎片
数值溢出	坐标值很大时，叉积计算可能溢出	使用 `long long` 或 `__int128`；或者归一化坐标
随机性质量	`rand()` 的随机性不够好，影响期望复杂度	使用更好的 PRNG（如 xorshift64）或预先洗牌
边界条件	查询点恰好在线段上、端点上、或边界框上	统一约定：点在线段上视为”下方”；点在端点上视为”右侧”

9.2 精确谓词：Shewchuk 的 orient2d

在工业级实现中，浮点精度是头号敌人。Jonathan Shewchuk 的精确算术谓词（robust predicates）是事实标准：

/*
 * Shewchuk 的 orient2d：
 *   返回正值表示 (pa, pb, pc) 逆时针
 *   返回负值表示顺时针
 *   返回 0 表示共线
 *
 * 内部使用自适应精度：先尝试普通浮点，
 * 如果结果不确定再用精确展开式。
 */
extern double orient2d(double *pa, double *pb, double *pc);

/* 使用方式 */
int robust_point_above(Point q, Segment *s) {
    double pa[2] = {s->p.x, s->p.y};
    double pb[2] = {s->q.x, s->q.y};
    double pc[2] = {q.x, q.y};
    return orient2d(pa, pb, pc) > 0.0;
}

9.3 符号扰动处理退化

当多个端点有相同的 \(x\) 坐标时，梯形分解的假设被违反。符号扰动是经典对策：

原始点 p = (x, y)
扰动后 p' = (x + ε·id(p), y + ε²·id(p))

其中 ε 是无穷小量，id(p) 是 p 的唯一标识符。

在实际中不需要真的加 \(\varepsilon\)，只需要在比较时加入 tie-breaking 规则：

int compare_x(Point *a, Point *b) {
    if (a->x != b->x) return (a->x < b->x) ? -1 : 1;
    if (a->y != b->y) return (a->y < b->y) ? -1 : 1;
    return 0; // 真正相同的点
}

十、地图引擎中的实际方案

10.1 现实世界的复杂性

真实的地图引擎（如 Google Maps、OpenStreetMap 渲染器）面对的点定位问题与教科书版本有显著差异：

数据规模：全球地图可能有数亿个多边形区域。
动态更新：地图数据需要增量更新，不能每次重建。
多层级：同一个点可能属于多个层级的区域（国家、省、市、区）。
非平面剖分：区域可能重叠（不同图层）。

10.2 R-tree 粗筛加精确定位

工业实践中最常见的方案是两阶段方法：

第一阶段：R-tree 粗筛

用 R-tree（或其变体 R*-tree）索引所有多边形的最小外接矩形（MBR）。查询时，先在 R-tree 中找到 MBR 包含查询点的所有候选多边形。

第二阶段：精确定位

对每个候选多边形，用射线法（Ray Casting）或缠绕数法（Winding Number）判断查询点是否真的在多边形内部。

/* 射线法：从 q 向右发射水平射线，统计与多边形边的交点数 */
int point_in_polygon(Point q, Point *polygon, int n) {
    int crossings = 0;
    for (int i = 0, j = n - 1; i < n; j = i++) {
        double yi = polygon[i].y, yj = polygon[j].y;
        double xi = polygon[i].x, xj = polygon[j].x;
        if ((yi > q.y) != (yj > q.y)) {
            double x_intersect = xi + (q.y - yi) / (yj - yi) * (xj - xi);
            if (q.x < x_intersect)
                crossings++;
        }
    }
    return crossings & 1; // 奇数 = 内部
}

10.3 性能分析

R-tree 查询:       O(log n) 平均，最坏 O(√n)
射线法判断:         O(k) ——k 是多边形顶点数
总体:              O(log n + k) 平均

空间:              O(n) —— R-tree 线性空间
支持动态更新:       是（R-tree 支持插入/删除）
支持重叠区域:       是（多个候选都会被检查）

10.4 更高级的方案

对于追求极致查询性能的场景：

网格索引（Grid Index）：将空间划分为均匀网格，每个网格单元预先存储覆盖它的多边形。查询 \(O(1)\) 定位网格，\(O(k)\) 检查候选。适合多边形大小相对均匀的情况。
四叉树（Quadtree）：自适应空间划分。在多边形密集区域细分，稀疏区域粗分。
S2 Geometry（Google）：将地球表面映射到希尔伯特曲线，转化为一维区间查询。

十一、应用场景详解

11.1 地图路由

在导航引擎中，路由算法需要知道用户的起点和终点分别在道路网络的哪条边上。这实际上是一个点定位问题：

道路网络将地图划分为若干面（街区、公园等）。
给定起点坐标，需要定位到最近的道路边。
这可以分解为：先做点定位找到所在面，再在该面的边上找最近点。

11.2 GIS 空间查询

地理信息系统中典型的”空间连接”操作：给定一组点和一组多边形，对每个点找到包含它的多边形。这本质上是批量点定位。

优化策略：

对多边形建一次索引，然后批量查询所有点。
如果点也有空间局部性，可以利用查询的连续性（上一次的查询结果可以加速下一次查询）。

11.3 机器人路径规划

在机器人学中，configuration space（C-space）被障碍物划分为自由空间和障碍空间。机器人路径规划的第一步就是确定当前配置在 C-space 的哪个区域中。

对于 2D 机器人，这直接归结为平面点定位。对于更高维度，通常使用近似方法（如随机采样）。

11.4 计算机图形学

在光线追踪中，判断一条光线穿过哪些三角形面片。虽然通常用 BVH（Bounding Volume Hierarchy）或 kd-tree，但底层思想与点定位相通。

在 2D 图形引擎中，鼠标点击测试（hit testing）——判断用户点击了哪个图形元素——就是点定位问题的直接应用。

十二、方法对比与个人思考

12.1 各方法综合对比

方法	预处理时间	空间	查询时间	确定性	实现难度	实际使用
Slab 分解	\(O(n^2 \log n)\)	\(O(n^2)\)	\(O(\log n)\)	是	低	仅教学
持久化 BST	\(O(n \log n)\)	\(O(n \log n)\)	\(O(\log n)\)	是	高	少见
梯形分解	\(O(n \log n)\) 期望	\(O(n)\) 期望	\(O(\log n)\) 期望	否	中	常见
Kirkpatrick	\(O(n)\)	\(O(n)\)	\(O(\log n)\)	是	极高	几乎不用
R-tree + 射线法	\(O(n \log n)\)	\(O(n)\)	\(O(\log n + k)\)	是	中低	工业标准

12.2 个人看法

梯形分解是计算几何中的”快速排序”。 就像快速排序一样，梯形分解用随机化换来了简洁的算法设计和优秀的期望性能。它不是最坏情况最优的，但在绝大多数实际场景中表现出色。

我对各种方法的评价：

梯形分解是教学和竞赛中的最佳选择。算法思路清晰，实现适中，性能优秀。如果你只学一种点定位方法，学它就够了。
Kirkpatrick 层级法是理论上的珍宝。它证明了 \(O(n)\) 空间 \(O(\log n)\) 查询是可以达到的，但实现的复杂性使其几乎没有实际价值。这提醒我们：理论最优和工程最优经常不是一回事。
R-tree 加射线法是工程中的王者。它不优雅，不理论最优，但它能处理现实世界的各种 messy 情况：动态更新、重叠区域、非平面输入。工程选择永远是”够好加上足够灵活”。
持久化 BST 方法有理论价值，也是理解函数式数据结构的好练习。但空间多了一个 \(\log\) 因子，在实际中没有明显优势。

关于随机化的思考：梯形分解的期望分析是计算几何中随机化分析的经典范例。逆向分析（backward analysis）的技巧非常漂亮——不是问”插入第 \(i\) 条线段时会发生什么”，而是问”在已有 \(i\) 条线段的梯形分解中，随机删除一条会影响多少”。这种逆向思维在很多随机算法中都能见到。

关于实际工程：我在实际项目中遇到的点定位问题，最后几乎都用了 R-tree 方案。原因很简单：

数据不满足梯形分解的”线段不交叉”假设。
需要动态更新。
多边形可能有空洞。
查询性能用 R-tree 已经足够好了。

但这不意味着学习梯形分解没有价值。理解它的构造方法和随机化分析，对于理解计算几何的思维方式至关重要。很多时候，你需要先理解理论上的精确方法，才能在工程中做出正确的 trade-off。

12.3 选择指南

if (你在写论文或教学) {
    使用梯形分解;
} else if (数据静态 && 是严格平面剖分 && 追求极致查询性能) {
    使用梯形分解;
} else if (数据动态 || 有重叠区域 || 数据可能不规范) {
    使用 R-tree + 射线法;
} else if (需要确定性最坏情况保证 && 不怕实现复杂) {
    考虑 Kirkpatrick（但大概率你不需要）;
} else {
    还是用 R-tree 吧;
}

12.4 一个有趣的历史注脚

点定位问题的历史可以追溯到 1970 年代末。Dobkin 和 Lipton（1976）首先给出了 \(O(\log n)\) 查询算法（Slab 方法），但空间是 \(O(n^2)\)。之后 Lipton 和 Tarjan（1980）用分离器定理将空间降到 \(O(n)\)，但构造复杂度很高。Kirkpatrick（1983）给出了优雅的 \(O(n)\) 构造。Edelsbrunner、Guibas 和 Stolfi（1986）提出了最优确定性方法。而梯形分解的随机化方法由 Mulmuley（1990）和 Seidel（1991）独立提出。

这段历史展示了计算几何领域如何在 15 年间逐步将一个问题从 \(O(n^2)\) 空间优化到 \(O(n)\) 空间。每一步都有漂亮的算法思想。

参考文献

de Berg, M., Cheong, O., van Kreveld, M., & Overmars, M. (2008). Computational Geometry: Algorithms and Applications (3rd ed.). Springer. Chapter 6.
Mulmuley, K. (1990). A fast planar partition algorithm, I. J. Symbolic Computation, 10(3-4), 253-280.
Seidel, R. (1991). A simple and fast incremental randomized algorithm for computing trapezoidal decompositions and for triangulating polygons. Comput. Geom., 1(1), 51-64.
Kirkpatrick, D. (1983). Optimal search in planar subdivisions. SIAM J. Comput., 12(1), 28-35.
Dobkin, D., & Lipton, R. (1976). Multidimensional searching problems. SIAM J. Comput., 5(2), 181-186.
Sarnak, N., & Tarjan, R. E. (1986). Planar point location using persistent search trees. Commun. ACM, 29(7), 669-679.
Shewchuk, J. R. (1997). Adaptive precision floating-point arithmetic and fast robust geometric predicates. Discrete Comput. Geom., 18(3), 305-363.
Edelsbrunner, H., Guibas, L. J., & Stolfi, J. (1986). Optimal point location in a monotone subdivision. SIAM J. Comput., 15(2), 317-340.
de Berg, M., van Kreveld, M., & Snoeyink, J. (1995). Two- and three-dimensional point location in rectangular subdivisions. J. Algorithms, 18(2), 256-277.
Preparata, F. P., & Shamos, M. I. (1985). Computational Geometry: An Introduction. Springer.

算法系列导航：上一篇：R-tree 与空间索引 | 下一篇：最近点对与随机化几何算法

相关阅读：Voronoi 图与 Delaunay 三角剖分 | 持久化数据结构

目录