cp-algorithms 垂直分解（Vertical Decomposition）算法详解：O(n² log n) 三角形并集面积与凸多边形求交

导读

垂直分解（Vertical Decomposition）是计算几何中的一类经典技术：将平面用若干垂直直线切分成纵向条带（vertical stripe），使每个条带内的问题结构"变好"，再对条带逐一独立求解。本文以本仓库 src/geometry/vertical_decomposition.md 为骨架，完整讲解如何用垂直分解计算 $n$ 个三角形的并集面积（从 $O(n^3\log n)$ 优化到 $O(n^2\log n)$ 时间、$O(n)$ 内存），并介绍其在两个凸多边形线性求交中的应用；同时结合本仓库的完整 C++ 实现与 test/test_vertical_decomposition.cpp 测试用例，说明代码如何落地运行。读完本文，你将掌握条带分解、括号平衡扫描线、事件驱动差分这三种核心思想及其可复用的实现模板。

---

一、垂直分解：核心思想概述

垂直分解（Vertical decomposition）被广泛应用于各类几何问题。其一般思路是：

1. 切分平面：用若干条竖直直线（即垂直于 $Ox$ 轴的线）将平面切成若干"纵向条带"；
2. 条带具备"好性质"：在每个条带内部，几何对象之间的相对关系变得简单、规整；
3. 独立求解：在每条条带上分别解决问题，再汇总结果。

"好性质"因问题而异，但通常包括：条带内没有线段相交、线段之间的上下顺序固定、任意一条横穿条带的曲线恰好把条带内的区域切成若干梯形等。下面以仓库文档中的两个经典问题为例展开。

---

二、应用一：三角形并集面积（Union Area of Triangles）

2.1 问题描述与朴素思路

假设平面上给定 $n$ 个三角形，要求计算它们并集的面积。若三角形互不相交，问题就退化为简单多边形面积求和（可参考 area-of-simple-polygon.md 中的鞋带公式方法）；但一般情形下三角形会彼此交叠，直接求和会重复计算重叠区域。

垂直分解给出的方案是：

1. 过所有三角形的顶点以及任意两条不同三角形边之间的交点，作竖直分割线。每条分割线都是 $x = x_0$ 形式的垂直线。
2. 顶点与交点合计可能有 $O(n^2)$ 个不同的 $x$ 坐标，因此平面被切成 $O(n^2)$ 个纵向条带。
3. 观察单个条带可以发现两个关键性质：
   • 任意一条非竖直线段要么从左到右横穿该条带，要么完全不穿过它；
   • 条带内部不存在两条线段严格相交（因为所有交点都已经被切到了条带边界上）。

这两点意味着：该条带内三角形并集的部分，恰好由若干底边落在条带两侧边界上的互不相交梯形组成。计算这些梯形面积之和，即为整个并集面积。

2.2 条带内扫描线：括号匹配视角

每个横穿条带的线段，要么是"上边界"，要么是"下边界"，取决于对应三角形内部位于线段的哪一侧：

• 若三角形内部在线段上方，则此线段为该三角形的上边界；
• 若三角形内部在线段下方，则此线段为该三角形的下边界。

一个优美的观察是：把每条上边界想象成左括号 (，把每条下边界想象成右括号 )，那么一条条带内的线段在竖直方向上的叠放关系恰好构成一段括号序列。条带内梯形的一个竖直切面被哪些线段包围、如何围成梯形，对应着把这段括号序列拆分成若干正确的（合法）括号子序列的过程。于是：

• 一个梯形的上下边界分别是某对匹配的"上括号/下括号"线段；
• 当竖直扫描线扫过条带时，括号序列的**平衡值（balance）**为 0 的区间，就是封闭梯形的区间。

对每条条带独立做这样的扫描线，总复杂度为 $O(n^3\log n)$（排序 $O(n^2)$ 个条带内的事件，每层再乘上 $O(n)$ 条线段），内存为 $O(n^2)$（需要显式保存所有条带）。这就是文档给出的第一版算法。

朴素版成本分析：条带数量 $O(n^2)$，每个条带内有 $O(n)$ 条线段需要排序，故总时间 $O(n^3\log n)$；因为要逐条带生成并保存梯形，内存 $O(n^2)$。

---

三、优化一：把竖直扫描线变成水平扫描线（O(n² log n)）

朴素算法为每个条带都跑一次竖直扫描线，成本过高。文档给出的第一个优化是不再按条带生成梯形，而是按线段枚举"它何时成为某个梯形的边"。

3.1 关键观察

固定某条三角形边（线段 $s = (s_0, s_1)$），我们关心的是：哪些条带中，$s$ 是某个梯形的边？

从括号平衡的视角看，$s$ 成为梯形边的充要条件是：在它下方（对下边界则是上方）的括号平衡值为 0。也就是说，我们不需要对每个条带分别处理，而只需找出所有会改变"以 $s$ 为参照的括号平衡值"的其他线段片段，然后对它们跑一次水平扫描线。

为表述简便，以下只讨论上边界的情形；下边界的处理完全对称。

3.2 对另一条线段 t 的三分类讨论

考虑任意一条非竖直线段 $t = (t_0, t_1)$，设 $s$ 与 $t$ 在 $Ox$ 轴上投影的交集为 $[x_1, x_2]$：

• 若该交集为空或退化为一个点，则 $s$ 与 $t$ 不可能穿过同一条条带的内部，$t$ 可以直接丢弃。

否则，考察 $s$ 与 $t$ 的几何关系 $I = s \cap t$，分三种情况：

情况 1：$I = \varnothing$（两线段不相交）

此时在 $[x_1, x_2]$ 上，$t$ 要么完全在 $s$ 上方，要么完全在 $s$ 下方：

• 若 $t$ 在 $s$ 上方：它不会影响 $s$ 是否构成某梯形边，可忽略；
• 若 $t$ 在 $s$ 下方：它会对 $[x_1, x_2]$ 内所有条带的括号平衡值产生贡献。贡献为 $+1$ 还是 $-1$，取决于 $t$ 本身是上边界还是下边界。

情况 2：$I$ 恰为单个交点 $p$

把 $[x_1, x_2]$ 拆成 $[x_1, p_x]$ 与 $[p_x, x_2]$ 两个区间，即可归约到情况 1 分别处理。由于 $p$ 本身是 $s$ 与 $t$ 的交点，必被竖直分割线切开，两个子区间内部 $s$ 与 $t$ 不再相交。

情况 3：$I$ 为一段线段 $l$（共线重叠）

这说明 $s$ 与 $t$ 在 $[x_1, x_2]$ 上完全重合：

• 若 $t$ 是下边界，则 $s$ 显然不可能成为梯形边（同一位置既有上边界又有下边界，无法围出梯形）；
• 若 $t$ 也是上边界，则 $s$、$t$ 都有可能被当成梯形边，存在二义性。解决办法是约定只有编号更小的线段才能作为梯形边（这里假定三角形边已被某种方式编号）。于是：
  • 若 $index(s) < index(t)$，忽略该情况（$t$ 不作为边，$s$ 不受影响）；
  • 否则，标记 $s$ 在 $[x_1, x_2]$ 上永远不可能成为边——实现上可以添加一个平衡值为 $-2$ 的事件，使该区间内的平衡值不可能回到 0，从而自然屏蔽 $s$。

下图直观展示了三种情况的几何形态：

3.3 区间加值的事件化与扫描线处理

对 $[x_1, x_2]$ 上"统一加 $w$"（$w$ 为 $+1$、$-1$ 或 $-2$ 等）的操作，不必逐条带更新，而是用差分事件表示：

• 创建事件 $(x_1, w)$（进入区间时加 $w$）；
• 创建事件 $(x_2, -w)$（离开区间时减 $w$）。

然后对所有事件按 $x$ 坐标排序，用一条扫描线顺序处理：维护当前平衡值，当平衡值为 0 时，说明 $s$ 在该区间内构成某个梯形的边，此时用 $s$ 的直线方程累加相应梯形面积即可。因为所有事件只需排序一次，总复杂度降为 $O(n^2\log n)$。

需要说明：事件排序是瓶颈，$O(n)$ 条线段每一条会产生 $O(n)$ 个事件，共 $O(n^2)$ 个事件，排序开销 $O(n^2\log n)$，这与朴素版的 $O(n^3\log n)$ 相比是本质性改进。

---

四、优化二：隐式条带，内存降到 O(n)

第二个优化十分直接：采用上一节的"按线段 + 事件扫描线"策略后，不再需要显式生成并保存所有条带。条带只存在于推理中，代码里只保存线段、类型与事件队列等 $O(n)$ 规模的数据结构。

• 时间：$O(n^2\log n)$；
• 内存：$O(n)$。

这使该算法可以在 $n$ 达到数千量级时依然可行（前提是事件常数较小、排序可控）。

---

五、应用二：两个凸多边形的线性时间求交

垂直分解的另一个经典应用是在 $O(n)$ 时间内求两个凸多边形的交集：

1. 过两个多边形的每个顶点作竖直分割线，把平面切成若干纵向条带；
2. 考虑其中一个凸多边形与任意一个条带的交集——由于凸多边形在每个竖直切片上的表现是"一段连续竖直区间"，这个交集只可能是梯形、三角形或一个点三种形状之一；
3. 在每个条带内，把两个多边形对应的这些简单形状两两求交，再把各条带的交集按顺序合并，即可拼出完整的交集多边形。

整个过程没有排序、没有求交事件，每个顶点只被处理常数次，因此是严格的线性时间。它体现了垂直分解"条带化后局部结构简单"这一核心收益。

---

六、完整实现：union_area 源码逐段解读

以下是仓库 src/geometry/vertical_decomposition.md 中给出的完整代码，实现了 $O(n^2\log n)$ 时间、$O(n)$ 内存的三角形并集面积计算。我们分段解读其关键结构。

6.1 浮点比较与点、线基础类型

typedef double dbl;

const dbl eps = 1e-9;

inline bool eq(dbl x, dbl y){
    return fabs(x - y) < eps;
}

inline bool lt(dbl x, dbl y){
    return x < y - eps;
}

inline bool gt(dbl x, dbl y){
    return x > y + eps;
}

inline bool le(dbl x, dbl y){
    return x < y + eps;
}

inline bool ge(dbl x, dbl y){
    return x > y - eps;
}

使用 eps = 1e-9 的容差比较（eq/lt/gt/le/ge），避免浮点误差导致条带划分与括号判断出错。这是所有几何算法的基础设施。

struct pt{
    dbl x, y;
    inline pt operator - (const pt & p)const{
        return pt{x - p.x, y - p.y};
    }
    inline pt operator + (const pt & p)const{
        return pt{x + p.x, y + p.y};
    }
    inline pt operator * (dbl a)const{
        return pt{x * a, y * a};
    }
    inline dbl cross(const pt & p)const{
        return x * p.y - y * p.x;
    }
    inline dbl dot(const pt & p)const{
        return x * p.x + y * p.y;
    }
    inline bool operator == (const pt & p)const{
        return eq(x, p.x) && eq(y, p.y);
    }
};

struct Line{
    pt p[2];
    Line(){}
    Line(pt a, pt b):p{a, b}{}
    pt vec()const{
        return p[1] - p[0];
    }
    pt& operator [](size_t i){
        return p[i];
    }
};

• pt 提供向量运算：减法、加法、数乘、叉积 cross、点积 dot 以及基于 eq 的相等比较；
• Line 用两个端点 p[0]、p[1] 表示线段，vec() 返回方向向量；
• 这两个类型与仓库 segments-intersection.md 中线段求交所用到的点/线基础结构同源。

6.2 字典序比较与线段求交

inline bool lexComp(const pt & l, const pt & r){
	if(fabs(l.x - r.x) > eps){
		return l.x < r.x;
	}
	else return l.y < r.y;
}

vector<pt> interSegSeg(Line l1, Line l2){
    if(eq(l1.vec().cross(l2.vec()), 0)){
        if(!eq(l1.vec().cross(l2[0] - l1[0]), 0))
            return {};
        if(!lexComp(l1[0], l1[1]))
            swap(l1[0], l1[1]);
        if(!lexComp(l2[0], l2[1]))
            swap(l2[0], l2[1]);
        pt l = lexComp(l1[0], l2[0]) ? l2[0] : l1[0];
        pt r = lexComp(l1[1], l2[1]) ? l1[1] : l2[1];
        if(l == r)
            return {l};
        else return lexComp(l, r) ? vector<pt>{l, r} : vector<pt>();
    }
    else{
        dbl s = (l2[0] - l1[0]).cross(l2.vec()) / l1.vec().cross(l2.vec());
        pt inter = l1[0] + l1.vec() * s;
        if(ge(s, 0) && le(s, 1) && le((l2[0] - inter).dot(l2[1] - inter), 0))
            return {inter};
        else
            return {};
    }
}

• lexComp 按 $(x, y)$ 做带容差的字典序比较，用于统一线段方向（保证线段"从左到右"）以及后续梯形边编号判定；
• interSegSeg 返回两线段的交集：
  • 方向向量叉积为 0 时两线段平行/共线：先检查是否共线（否则返回空），再把两线段各按字典序排序，取"左端点的较大者"与"右端点的较小者"得到重叠部分——这与 segments-intersection.md 中平行情形的处理思路一致；
  • 非平行时用参数 $s$ 求交点，并检查 $s\in[0,1]$ 且交点确实落在线段 2 上（利用点积非正判断）才返回该交点。

6.3 线段类型：上边界 / 下边界 / 竖直

inline char get_segtype(Line segment, pt other_point){
    if(eq(segment[0].x, segment[1].x))
        return 0;
    if(!lexComp(segment[0], segment[1]))
        swap(segment[0], segment[1]);
    return (segment[1] - segment[0]).cross(other_point - segment[0]) > 0 ? 1 : -1;
}

get_segtype 决定一条边的"角色"，返回值三种：

• 0：竖直边（两端点 $x$ 相等）——竖直边横穿条带不构成上/下边界，直接排除；
• 1：上边界（三角形内部位于边的上方，即第三个顶点 other_point 在边的左侧/上方，叉积为正）；
• -1：下边界（三角形内部位于边的下方）。

这个分类正是第二节"括号匹配"中"左括号/右括号"的来源。

6.4 主函数 union_area

dbl union_area(vector<tuple<pt, pt, pt> > triangles){
    vector<Line> segments(3 * triangles.size());
    vector<char> segtype(segments.size());
    for(size_t i = 0; i < triangles.size(); i++){
        pt a, b, c;
        tie(a, b, c) = triangles[i];
        segments[3 * i] = lexComp(a, b) ? Line(a, b) : Line(b, a);
        segtype[3 * i] = get_segtype(segments[3 * i], c);
        segments[3 * i + 1] = lexComp(b, c) ? Line(b, c) : Line(c, b);
        segtype[3 * i + 1] = get_segtype(segments[3 * i + 1], a);
        segments[3 * i + 2] = lexComp(c, a) ? Line(c, a) : Line(a, c);
        segtype[3 * i + 2] = get_segtype(segments[3 * i + 2], b);
    }
    vector<dbl> k(segments.size()), b(segments.size());
    for(size_t i = 0; i < segments.size(); i++){
        if(segtype[i]){
            k[i] = (segments[i][1].y - segments[i][0].y) / (segments[i][1].x - segments[i][0].x);
            b[i] = segments[i][0].y - k[i] * segments[i][0].x;
        }
    }
    dbl ans = 0;
    for(size_t i = 0; i < segments.size(); i++){
        if(!segtype[i])
            continue;
        dbl l = segments[i][0].x, r = segments[i][1].x;
        vector<pair<dbl, int> > evts;
        for(size_t j = 0; j < segments.size(); j++){
            if(!segtype[j] || i == j)
                continue;
            dbl l1 = segments[j][0].x, r1 = segments[j][1].x;
            if(ge(l1, r) || ge(l, r1))
                continue;
            dbl common_l = max(l, l1), common_r = min(r, r1);
            auto pts = interSegSeg(segments[i], segments[j]);
            if(pts.empty()){
                dbl yl1 = k[j] * common_l + b[j];
                dbl yl = k[i] * common_l + b[i];
                if(lt(yl1, yl) == (segtype[i] == 1)){
                    int evt_type = -segtype[i] * segtype[j];
                    evts.emplace_back(common_l, evt_type);
                    evts.emplace_back(common_r, -evt_type);
                }
            }
            else if(pts.size() == 1u){
                dbl yl = k[i] * common_l + b[i], yl1 = k[j] * common_l + b[j];
                int evt_type = -segtype[i] * segtype[j];
                if(lt(yl1, yl) == (segtype[i] == 1)){
                    evts.emplace_back(common_l, evt_type);
                    evts.emplace_back(pts[0].x, -evt_type);
                }
                yl = k[i] * common_r + b[i], yl1 = k[j] * common_r + b[j];
                if(lt(yl1, yl) == (segtype[i] == 1)){
                    evts.emplace_back(pts[0].x, evt_type);
                    evts.emplace_back(common_r, -evt_type);
                }
            }
            else{
                if(segtype[j] != segtype[i] || j > i){
                    evts.emplace_back(common_l, -2);
                    evts.emplace_back(common_r, 2);
                }
            }
        }
        evts.emplace_back(l, 0);
        sort(evts.begin(), evts.end());
        size_t j = 0;
        int balance = 0;
        while(j < evts.size()){
            size_t ptr = j;
            while(ptr < evts.size() && eq(evts[j].first, evts[ptr].first)){
                balance += evts[ptr].second;
                ++ptr;
            }
            if(!balance && !eq(evts[j].first, r)){
                dbl next_x = ptr == evts.size() ? r : evts[ptr].first;
                ans -= segtype[i] * (k[i] * (next_x + evts[j].first) + 2 * b[i]) * (next_x - evts[j].first);
            }
            j = ptr;
        }
    }
    return ans/2;
}

逐段对应文档理论：

1. 线段生成：每个三角形 3 条边，共 $3n$ 条线段；每条边都通过 lexComp 统一为"左端点在前"，并用第三个顶点确定类型（上/下边界或竖直）。
2. 直线方程：对非竖直线段预计算斜率 $k$ 与截距 $b$（$y = kx + b$），供后续快速比较高低与累加梯形面积。
3. 外层循环（固定 $s$）：对每条非竖直线段 $s = segments[i]$，在投影区间 $[l, r]$ 上收集影响其括号平衡的事件：
   • 投影无交集（ge(l1, r) || ge(l, r1)）直接跳过；
   • 情况 1（pts.empty()）：在公共区间端点比较 $t$ 与 $s$ 的高低；只有当 $t$ 位于 $s$ 的"相关一侧"（对 $s$ 为上边界时 $t$ 在其下方，即 lt(yl1, yl) == (segtype[i] == 1)）才产生 $\pm1$ 事件；
   • 情况 2（pts.size() == 1u）：以交点 $x$ 坐标为界拆成两个区间，各自按情况 1 处理（对应文档"把 $[x_1,x_2]$ 拆成 $[x_1,p_x]$ 与 $[p_x,x_2]$"）；
   • 情况 3（共线重叠）：只有不同类或编号较大者才添加 -2/+2 事件，把 $s$ 在该区间"永远不可能成为边"这一事实编码进平衡值（对应文档"添加平衡值 $-2$ 的事件"）。
4. 事件处理：区间加值被差分化为 (x_1, w) 与 (x_2, -w) 两个事件，全部排序后用一条水平扫描线顺序扫描。相同 $x$ 坐标的事件合并处理，避免同一位置重复结算。
5. 面积累加：当平衡值 balance == 0 且未到右端点 $r$ 时，说明 $s$ 在此区间内是某梯形边界，把区间 $[evts[j].first, next_x]$ 上的梯形面积（利用直线方程 $y = kx + b$ 的定积分形式 $\int (kx+b)\,dx = \frac{k}{2}(x_2^2-x_1^2)+b(x_2-x_1)$，代码中按 $(next\_x + evts[j].first)$ 与 $(next\_x - evts[j].first)$ 的乘积形式写出）累加到 ans；上/下边界分别取正负号，因此最终 ans 是"带符号面积的两倍"，除以 2 即得并集面积。

整体来看：外层 $O(n)$ 条线段 × 内层 $O(n)$ 条候选线段 + 每层 $O(n)$ 规模事件的排序，总复杂度 $O(n^2\log n)$，辅助存储只有事件队列与常数级数组，内存 $O(n)$。

---

七、如何运行与验证：仓库测试链路

本仓库为这篇文章的实现提供了完整的自动化测试闭环，可以按以下链路验证正确性：

1. 测试用例：test/test_vertical_decomposition.cpp 通过 #include "triangle_union.h" 引用从文档抽取的代码片段。它把 4 个矩形拆成 8 个三角形（每个矩形按对角线切成两个三角形），再调用 union_area 求并集面积并与解析值比对：

   • test1：4 个边长为 3 的正方形彼此交叠成十字形，期望并集面积为 $1.5$，容差 $10^{-6}$；
   • test2：两个菱形相交，期望并集面积为 $1.1428571429$（即 $\frac{8}{7}$），容差 $10^{-6}$。

   测试先把矩形转为三角形、把"矩形并集面积"映射为"三角形并集面积"（solve 中 sum/union_area(tr) 的比值恒为 1），从而用矩形这种易手工验证的形状间接检验三角形并集算法。

2. 代码抽取：test/extract_snippets.py 扫描 src/ 下所有 Markdown，用正则 ^\s*```\{.cpp\s+file=(\S+)\}$ 匹配形如 {.cpp file=triangle_union} 的代码块，把其中的源码抽取为 triangle_union.h 头文件——这正是文档代码片段能直接参与编译测试的机制。

3. 编译与运行：test/test.sh 依次执行 python extract_snippets.py 抽取代码，再用 g++ -std=c++17 -fsanitize=undefined -fno-sanitize-recover 编译每个 *.cpp 测试并运行；通过 assert 断言，若两个测试全部通过即验证算法实现正确。

本地复现步骤：在仓库 test/ 目录下运行 bash test.sh（或手动执行 python extract_snippets.py && g++ -std=c++17 test_vertical_decomposition.cpp -o test && ./test）。依赖环境为 Python 3 与支持 C++17 的 GCC/Clang。

---

八、复杂度总结与适用边界

阶段	时间	内存	说明
朴素垂直分解	$O(n^3\log n)$	$O(n^2)$	逐条带生成梯形并扫描
优化一（水平扫描线）	$O(n^2\log n)$	$O(n^2)$	事件化区间加值
优化二（隐式条带）	$O(n^2\log n)$	$O(n)$	不再显式保存条带
凸多边形求交	$O(n)$	$O(n)$	逐条带合并简单形状

需要说明的适用前提：

• 算法基于实数坐标与 eps = 1e-9 容差，输入含大量近似共线或退化点（如三点近共线、边近竖直）时需谨慎，竖直边（segtype == 0）会被直接排除；
• 复杂度分析假定三角形边数 $3n$，若输入为任意多边形需先三角剖分；
• 浮点累加在大 $n$ 下存在误差累积，竞赛环境下建议使用 long double 或额外归一化（如需更高精度）。

---

九、相关练习与延伸阅读

垂直分解与三角形并集是经典计算几何题目，常见练习包括：

• Codeforces 62C Inquisition：经典三角形并集面积题，正是本文 union_area 的直接应用；
• Codeforces 107E Darts：掷飞镖与面积统计，可转化为并集面积计算。

仓库内可继续延伸阅读的关联资料：

• segments-intersection.md：本文 interSegSeg 求交思路的完整版（含平行、共线、单点退化处理）；
• area-of-simple-polygon.md：多边形面积计算，与梯形面积累加同源；
• convex-hull.md：凸包构造，可作为"凸多边形求交"的前置知识；
• test/test_vertical_decomposition.cpp 与 test/extract_snippets.py：本文实现的可运行测试与代码抽取机制。

掌握垂直分解 = 掌握"用竖直条带简化局部结构 + 用括号平衡判定梯形 + 用差分事件驱动扫描线"这套组合拳，它同样适用于区间覆盖面积、曲线并集、多面体投影等更广泛的问题。