Computational Geometry

Algorithms and Applications, Third Edition

Note

This html generated by Sphinx corresponds to file ComputationalGeometryAA3rd.rst, NOT ComputationalGeometryAA3rd.md.
This is because Sphinx didn't support math symbols/equations in Latex very well, so first edit the markdown file, and then export it to corresponding rst file, and use that rst for Sphinx to create html pages.

Notes by Pyrad, 2023-01-03

3 ways to generate corresponding reStructuredText from this markdown text

  1. Use pandoc

pandoc -s ComputationalGeometryAA3rd.md -f markdown -t rst -o ComputationalGeometryAA3rd.rst

  1. In Typora, use Export to function, and choose export as reStructuredText format.

  2. In Obsidian, use Pandoc Plugin: Export as reStructured Text (RST) function after pandoc plugin is installed.

Notes by Pyrad, 2023-01-27

Since now I switch to use MyST-Parser as the enhanced markdown parser, and I add mathjax support for it, so after I enabled the dollar math syntax extension of MyST-Parser, now the html files generated by a markdown file can display Latex symbols/maths correctly.

Thus I don’t need to convert this file (markdown) to .rst file anymore.

Refere to myst_enable_extensions variable set in file conf.py, and the link Dollar delimited math from MyST-Parser.

Notes by Pyrad, 2023-02-03

About the book

Third Edition (March 2008)

Mark de Berg, TU Eindhoven (the Netherlands) Otfried Cheong, KAIST (Korea) Marc van Kreveld, Mark Overmars, Utrecht University (the Netherlands)

Website URL

Published by Springer

Errata Page for 1st Edition

Errata Page for 2nd Edition

Resources

Handbook of Discrete and Computational Geometry, 3rd Edition, Online Page

计算几何有何经典书籍? - 王boy的回答 - 知乎

Words

asymptote /ˈæsɪmˌtoʊt/ n. [数] 渐近线

asymptotically /,æsimp’tɔtikli,-kəli/ adj. 渐近线的

planar /ˈpleɪnər/ adj. 平面的；二维的；平坦的

polygonal /pəˈlɪɡənl/ adj. [数] 多边形的；[数] 多角形的

lexicographic /ˌleksɪkəˈɡræfɪk/ adj. 词典编辑的；字典式的

dent /dent/ v. 使产生凹痕；损害，削弱；n. 凹痕；削减

clutter /ˈklʌtər/ v. 乱堆，塞满；使（脑子里）塞满（乱七八糟的事）；n. 杂乱的东西；杂乱

radiosity /ˌreɪdɪˈəʊsɪtɪ/ n. 光能传递；热辐射

kinematic /ˌkɪnəˈmætɪk/ adj. [力] 运动学上的，[力] 运动学的

vegetation /ˌvedʒəˈteɪʃn/ n. （总称）植物，植被；（植物的）生长；呆板单调的生活

excavation /ˌekskəˈveɪʃ(ə)n/ n. （对古物的）发掘，挖掘；发掘现场； 挖洞，开凿

interpolate /ɪnˈtɜːrpəleɪt/ vt. 篡改；插入新语句；vi. 插入；篡改

interpolating /ɪnˈtɜːrpəleɪtɪŋ/ n. 插值；内插；v. 窜改；加入（额外的事）

thematic /θɪˈmætɪk/ adj. 主题的，主旋律的；题目的；语干的

precipitation /prɪˌsɪpɪˈteɪʃ(ə)n/ n. 降水（如雨，雪，冰雹）；沉淀，淀析；仓促，鲁莽，轻率；坠落

grizzly /ˈɡrɪzli/ adj. 灰色的；n. 灰熊

lemma /ˈlemə/ n. 引理；辅助定理；论点；膜

recap /ˈriːkæp/ v. 扼要重述，摘要说明；翻新胎面；n. 扼要的重述，概述；翻新的轮胎

treat /triːt/ v. 处理，探讨，论述；

deciduous /dɪˈsɪdʒuəs/ adj. 落叶性的，脱落性的；非永久性的

cyclic /ˈsaɪklɪk/ adj. 环的；循环的；周期的

corollary /ˈkɔːrəleri/ n. 推论；必然的结果

paradigm /ˈpærədaɪm/ n. 典范，范例；样板，范式；词形变化表；纵聚合关系语言项

triangulation /traɪˌæŋɡjuˈleɪʃn/ n. [测] 三角测量；三角形划分

overkill /ˈoʊvərkɪl/ n. 过犹不及

dual graph 【数】对偶图

prong /prɔːŋ/ n. 尖齿；方向；（俚语）阴茎；耙子；（多种工具的）尖头；（作战）分支；v. 刺，贯穿

combinatorial /kɒmˌbaɪnəˈtɔːrɪəl/ adj. 组合的

monotone /ˈmɒnətəʊn/ n. 单调；单音调；adj. 单调的

quadrilateral /ˌkwɒdrɪˈlætərəl/ n. 四边形；adj. 四边形的

funnel /ˈfʌn(ə)l/ n. 漏斗；漏斗状物；（船、蒸汽机车上的）烟囱；*v.*通过漏斗，穿过狭窄通道；输送，传送（金钱、货物或信息）；

reflex /ˈriːfleks/ n. （对刺激的）本能反应；反射（作用）；反映物；反映形式；反射光；adj. 本能反应的；（角）大于180度的；（光）被反射的；反折的；反省的

nevertheless /ˌnevəðəˈles/ adv. 然而，不过

colorability /ˌkʌlərəˈbɪlɪti/ 可着色性

quadrilaterals /ˌkwɒdrɪˈlætərəl/ n. 四边形；adj. 四边形的

trapezoid /ˈtræpəzɔɪd/ n. <英>不规则四边形；<美>梯形；（腕部近食指根底处的）小多角骨；adj. 梯形的，不规则四边形的

polytope /’pɔlitəup*/* n. [数] 多胞形，多面体；可剖分空间

tetrahedra /ˌtetrəˈhiːdrə/ n. 四面体（tetrahedron 的变形）

cutlery /ˈkʌtləri/ n. 餐具（刀、叉和匙）；刀具

casting /ˈkɑːstɪŋ/ n. 角色分配，演员挑选；铸件，铸造物

solidify /səˈlɪdɪfaɪ/ v. （使）凝固，变硬；（使）变可靠，（使）变稳固

polyhedral /ˌpɒliˈhiːdrəl; ˌpɒliˈhedrəl/ adj. [数] 多面的；[数] 多面体的

facet /ˈfæsɪt/ n. 部分，方面；（宝石的）琢面，刻面；（构成昆虫或甲壳动物复眼的）小眼面；v. 在……上琢面

cavity /ˈkævəti/ n. 洞，腔；(牙齿的) 龋洞

latex /‘la:tek/ 一种电子排版系统

coplanar /kəʊˈpleɪnə(r)/ adj. [数] 共面的（be coplanar with）

circumflex /ˈsɜːkəmfleks/ n. 音调符号；adj. 弯曲的；有声调符号的；v. 标有抑扬音符；弯曲

inequality /ˌɪnɪˈkwɒləti/ n. （大小、程度、情况等的）不同，不平等；（数学）（两式间的）不等；（数学）不等式

inequation /ˌɪnɪˈkweɪʒən/ n. 不等式

wedge /wedʒ/ n. 楔子，三角木；楔形物，三角形物（尤指食物）；v. （将……）楔入，插入，挤入；

parameterize /pəˈræmɪt(ə)raɪz/ vt. 用参数表示；确定……的参数

coin /kɔɪn/ v.创造新词，首次使用；铸币，造币；

intriguing /ɪnˈtriːɡɪŋ/ adj.非常有趣的，引人入胜的（intrigue现在分词）

mishap /ˈmɪshæp/ n. 灾祸；不幸事故；晦气

formulate /ˈfɔːmjuleɪt/ v. 制定，规划；确切表达，认真阐述；用公式表示

subsume /səbˈsjuːm/ vt. 把……归入；把……包括在内

canonical /kəˈnɒnɪk(ə)l/ adj. 根据教规的，标准的，典范的；（canonical subset 正则子集）

disjoint union 不相交并集

such and such 某某；如此这般的

chronometer /krəˈnɒmɪtə(r)/ n. 精密记时表；高度精确的钟表；航海经线仪

thematic /θɪˈmætɪk/ adj. 主题的，主旋律的；题目的；语干的

Usage

thought experiment

an elastic rubber band 橡皮筋

direct the line through \(p\) and \(q\)

to this end 为了这个目的（formal : as a way of dealing with or doing something）

rule out 排除，除去

windy river 弯曲的河流（不是多风的河流）

coinciding point 共点

in a sense 某种意义上

incident to 由…产生（这里incident是 adj.）

hold for 适用

mass produce 批量生产

give rise to 使发生，引起

necessary and sufficient conditions 充要条件

dot product 点积（注意，不是production）

inequality 不相等

subexponential /sʌbˌekspəˈnenʃl/ adj.（增长）越来越快的；指数的，含有指数的

Define the y-interval of a segment to be its orthogonal projection onto the y-axis.

把一条线段在 y 轴上的正交投影，叫做它的 y-interval

they are far apart in the y-direction

它们在y方向上相距足够远

We denote the event queue by Q

我们把event queue记作\(\mathcal{Q}\)

This horizontal sweeping line is sloping just a tiny bit upward

这条横向的扫描线翘起来一点点

We need an operation that removes the next event that will occur from Q, and returns it so that it can be treated.

需要一个从队列Q里面删除下个event（point）的操作，并且返回它，以便（对它进行）处理。

Therefore we model a gallery as a polygonal region in the plane.

我们把画廊当做一个二维平面上的多边形

different orientations of the object give rise to different molds.

Take the plane spanned by the vectors (we assume both vectors are rooted at the origin)

由两个向量展开的平面，假设这两个向量从原点出发。

is/are dual to 对偶于，对应于，相当于

The problem is dual to the computation of the convex hull of points in the plane.

这个问题相当于在平面上计算凸包的点。

intensive care station 重症监护站

an axis-parallel rectangle 和坐标轴平行的矩形

We want to report all the points inside an axis-parallel query rectangle

Names

  • Graham’s scan

  • output-sensitive algorithm

  • planar graph

  • planar subdivisions

  • Chapter 2

    • output-sensitive

  • Chapter 3

    • Art Gallery Problem (combinatorial geometry)

    • dual graph

    • Monotone polygon

  • Chapter 4

    • Linear optimization (linear programming，就是线性规划)

    • Simplex Algorithm (运筹学中的单纯形算法，in operations research area)

    • Low-dimensional linear programming problems

    • expected running time 期望运行时间（即n个运行时间的数学期望）

    • Linearity of expectation

    • Backward analysis

    • Smallest Enclosing Disc 最小圆覆盖（最小包围圆）

  • Chapter 5

    • Rectangular Range Query （Orthogonal Range Query）矩形区域查询，或正交区域查询

    • Canonical subset 正则子集

    • Fractional cascading 分散层叠

    • General Sets of Points 一般性点集

    • exact match queries 完全匹配查询

    • partial match queries 部分匹配查询

Maths

\(e'\) ：\(e\) prime（或 \(e\) dash）

\(e''\) ：\(e\) double prime（或 \(e\) double dash）

\(\hat{f}\) ：\(f\) hat（或者 \(f\) roof），caret符号更多的用于音调符号（circumflex）

Contents

  • Preface

  • 1 Computational Geometry (Introduction)

  • 2 Line Segment Intersection (Thematic Map Overlay)

  • 3 Polygon Triangulation (Guarding an Art Gallery)

  • 4 Linear Programming (Manufacturing with Molds)

  • 5 Orthogonal Range Searching (Querying a Database)

  • 6 Point Location (Knowing Where You Are)

  • 7 Voronoi Diagrams (The Post Office Problem)

  • 8 Arrangements and Duality (Supersampling in Ray Tracing)

  • 9 Delaunay Triangulations (Height Interpolation)

  • 10 More Geometric Data Structures (Windowing)

  • 11 Convex Hulls (Mixing Things)

  • 12 Binary Space Partitions (The Painter’s Algorithm)

  • 13 Robot Motion Planning (Getting Where You Want to Be)

  • 14 Quadtrees (Non-Uniform Mesh Generation)

  • 15 Visibility Graphs (Finding the Shortest Route)

  • 16 Simplex Range Searching (Windowing Revisited)

  • Bibliography

  • Index

Content in Chinese

图书目录

第1章 计算几何：导言
  1.1 凸包的例子
  1.2 退化及稳健性
  1.3 应用领域
  1.4 注释及评论
  1.5 习题
第2章 线段求交：专题图叠合
  2.1 线段求交
  2.2 双向链接边表
  2.3 计算子区域划分的叠合
  2.4 布尔运算
  2.5 注释及评论
  2.6 习题
第3章 多边形三角剖分：画廊看守
  3.1 覆盖与三角剖分
  3.2 多边形的单调块划分
  3.3 单调多边形的三角剖分
  3.4 注释及评论
  3.5 习题
第4章 线性规划：铸模制造
  4.1 铸造中的几何
  4.2 半平面求交
  4.3 递增式线性规划
  4.4 随机线性规划
  4.5 无界线性规划问题
  *4.6 高维空间中的线性规划
  *4.7 最小包围圆
  4.8 注释及评论
  4.9 习题
第5章 正交区域查找：数据库查询
  5.1 一维区域查找
  5.2 kd树
  5.3 区域树
  5.4 高维区域树
  5.5 一般性点集
  *5.6 分散层叠
  5.7 注释及评论
  5.8 习题
第6章 点定位：找到自己的位置
  6.1 点定位及梯形图
  6.2 随机增量式算法
  6.3 退化情况的处理
  *6.4 尾分析
  6.5 注释及评论
  6.6 习题
第7章 Voronoi图：邮局问题
  7.1 定义及基本性质
  7.2 构造Voronoi图
  7.3 注释及评论
  7.4 习题
第8章 排列与对偶：光线跟踪超采样
  8.1 差异值的计算
  8.2 对偶变换
  8.3 直线的排列
  8.4 层阶与偏差
  8.5 注释及评论
  8.6 习题
第9章 Delaunay三角剖分：高度插值
  9.1 平面点集的三角剖分
  9.2 Delaunay三角剖分
  9.3  构造Delaunay三角剖分
  9.4 分析
  *9.5 随机算法框架
  9.6 注释及评论
  9.7 习题
第10章 更多几何数据结构：截窗
  10.1 区间树
  10.2 优先查找树
  10.3 线段树
  10.4 注释及评论
  10.5 习题
第11章 凸包： 混合物
  11.1 三维凸包的复杂度
  11.2 构造三维凸包
  *11. 3分析
  *11.4 凸包与半空间求交
  *11.5 再论Voronoi图
  11.6 注释及评论
  11.7 习题
第12章 空间二分：画家算法
  12.1 BSP树的定义
  12.2 BSP树及画家算法
  12.3 构造BSP树
  *12.4 三维BSP树的规模
  12.5 注释及评论
  12.6 习题
第13章 机器人运动规划：随意所之
  13.1 工作空间与C空间
  13.2 点机器人
  13.3 Minkowski和
  13.4 平移式运动规划
  *13.5 允许旋转的运动规划
  13.6 注释及评论
  13.7 习题
第14章 四叉树：非均匀网格生成
  14.1 均匀及非均匀网格
  14.2 点集的四叉树
  14.3 从四叉树到网格
  14.4 注释及评论
  14.5 习题
第15章 可见性图：求最短路径
  15.1 点机器人的最短路径
  15.2 构造可见性图
  15.3 平移运动多边形机器人的最短路径
  15.4 注释及评论
  15.5 习题
第16章 单纯形区域查找：再论截窗
  16.1 划分树
  16.2 多层划分树
  16.3 切分树
  16.4 注释及评论
  16.5 习题
参考文献
关键词索引

Preface

序言要点

  • 计算几何兴起与20世纪70年代（1970s），应用于计算机图形（CG）、地理信息系统（GIS）、机器人（robotics）等领域。

  • 本书每章节基本独立，但初学者可以按顺序阅读前10章。

  • 每章节只举例了容易理解和实现的算法（解决方案），并不是所有，而且提供的是高层次的论述，并不深入每个细节。

  • 带星号（*）的章节作为扩展阅读，以及叫做 Notes and Comments的小节，可以通过其了解更多。

  • 不需要应用领域的知识，只需要基本的数据结构和算法知识储备。

  • 有网页可以找到Errata Page以及其他可用资源。

1 Computational Geometry Introduction

校园中寻找最近电话亭（Voronoi diagram，第7章）

避障最短路径（motion planning，第13,15章）

多张地图定位问题（overlay map，第2章）

实际的应用

  • Robotics

  • Computer graphics

  • CAD/CAM

  • Geographic Information System

1.1 An Example: Convex Hulls

好的几何算法问题解决方案，本质上有两方面

  • 理解几何问题的特性

  • 对算法和数据结构的合理使用

本节举例，介绍了两种计算二维平面上凸体的轮廓的算法（二维平面凸体，planar convex hulls）

1.1.1 第一种算法

第一种算法是时间复杂度较高的算法，文中称为 SlowConvexHull 算法。

输入：平面上点的集合 \(\mathcal{P}\) 。

输出：一个点的序列 \(\mathcal{L}\)，表示点集合 \(\mathcal{P}\) 的Convex Hull，点序是顺时针方向。

算法简述：

从集合 \(\mathcal{P}\) 中取任意不同两点 \(p\) 和 \(q\)，组成一有向线段 \(\overrightarrow{pq}\) ，检查集合 \(\mathcal{P}\) 中剩余的任意一点 \(r\)，如果任意一点 \(r\) 都位于有向线段 \(p \rightarrow q\) 的右侧，说明有向线段 \(p \rightarrow q\) 就是最终轮廓上的其中一条线段，将其加入集合 \(\mathcal{E}\) 中。

穷举集合 \(\mathcal{P}\) 中这样两个点 \(p\) 和 \(q\) 的组合，重复上述检查，直至最终遍历完成，得到一个线段集合 \(\mathcal{E}\)。

最后，从集合 \(\mathcal{E}\) 中找出依次连接的线段，并组成一个点列表，按照顺时针方向排序。

算法复杂度：\(O(n^3)\)

对于伪代码中的几个的说明

  • 诸如判断一个点在一条直线（线段）的左边或右边的操作，默认已经有现成的实现可以使用

  • 从集合 \(\mathcal{E}\) 中找出依次连接的有向线段的步骤是，首先从 \(\mathcal{E}\) 中取任一有向线段 \(e_1\)，以其头点（即 \(p \rightarrow q\) 线段的 \(q\) 点）为目标，从集合 \(\mathcal{E}\) 中找出第二条有向线段 \(e_2\)，其尾点（即 \(p \rightarrow q\) 线段的 \(p\) 点）为 \(e_1\) 的头点，然后再以 \(e_2\) 的头点搜索下一条有向线段，直到搜索到的这些线段 \(e_1\), \(e_2\), \(e_3\), …, \(e_N\) 构成一条闭合的折线段。

  • 关于算法复杂度是\(O(n^3)\)。从 \(n\) 个点中取两个点的组合是 \(\frac{n!}{2!(n-2)!}\)，所以是\(O(n^2)\)，对每一条有向线段，查看剩余 \(n-2\) 个点是否在其右侧，这样时间复杂度就达到了\(O(n^3)\)。最后一步依次找出有向线段并按顺序连接，时间复杂度是 \(O(n^2)\) ，所以最终时间复杂度就是\(O(n^3)\)。

关于 degenerate case 或者叫做 degeneracy

在判断一个点 k 是否在有向线段 \(p \rightarrow q\) 右侧时，点 k 是可能落在有向线段 \(p \rightarrow q\) 上的，针对这种退化情况，可以把它也当做是在有向线段右侧的一种（退化）情况。

关于 rounding error 导致的程序健壮性问题（robustness）

在实际情况中，因为使用的是浮点数计算，那么仍然是在判断一个点 \(k\) 是否在有向线段 \(p \rightarrow q\) 右侧时，可能产生微小的误差（rounding errors），导致最终计算出来的convex hull的点集合 \(\mathcal{E}\) 有三种情况：

  1. 要么不是真正意义上的convex hull（但仍然是非常接近实际情况的）

  2. 要么最终的集合 \(\mathcal{E}\) 中的有向线段不是一个闭合的折线段

  3. 要么最终的集合 \(\mathcal{E}\) 中的有向线段除了可以组成一个闭合的折线段外，还有额外剩余的几条有向线段

正是由于这种robust的问题，迫使我们需要寻找一种更为健壮和正确的算法。

1.1.2 第二种算法

第二种算法是时间复杂度比第一种算法低，采用了所谓的 incremental algorithm 的方法，文中名为ConvexHull算法。

这种算法的总体思路是，将所有的点按照 \(x\) 坐标由大到小排序为 \(p_1\), \(p_2\), \(p_3\), …, \(p_N\)，因为前提是凸多边形，所以先按照从左向右的方向，找到这convex hull的上半部分边界 upper hull，即 \(p_1\), \(u_0\), \(u_1\), …, \(p_N\)（其中\(u_0\), \(u_1\), … 都是集合中的点），再找到convex hull的下半部分边界 lower hull，即\(p_1\), \(v_0\), \(v_1\), …, \(p_N\)（其中\(v_0\), \(v_1\), … 也都是集合中的点）。

这个所谓的 incremental algorithm 方法的关键步骤在于，如何在向已有的但不完整的upper/lower hull 添加一个点之后，更新这个不完整的upper/lower hull，使得其向左或向右延伸一段（最终到达最右或最左的点）。

换句话说，假如现已有upper hull的点是 \(p_1\), \(p_2\), …, \(p_{i-1}\), 如何找到下一个点 pi，使得 \(p_1\), \(p_2\), …, \(p_i\) 是最终 upper hull 的一部分。

因为我们约定是按照顺时针方向来标记最终的convex hull的，所以，沿着convex hull的边界行走，一定是**“右转”**的。因此，可以按照此方法来确定如何加入上面提到的\(pi\)点，从而生成一条新的convex hull的一部分。

假设我们现在计算的是 upper hull，那么我们遍历的点一定是按照 \(x\) 坐标有小到大的顺序，那么当加入点 \(p_i\) 时，点 \(p_i\) 的 \(x\) 坐标就是目前的convex hull 点 \(p_1\), \(p_2\), …, \(p_{i-1}\) 里面 \(x\) 坐标最大的。

加入点 \(p_i\) 后，此时点列为 \(p_1\), \(p_2\), …, \(p_{i-1}\), \(p_i\)。此时我们检查最后三个点 \(p_{i-2}\), \(p_{i-1}\), \(p_i\)。

  • 如果这三个点是**“右转”**的，那么新加入的点 \(p_i\)，就是最终 upper convex hull的一部分（但有可能在加入之后的点以后，继续做调整从而删除点 \(p_i\)）。

  • 如果这三个点是**“左转”**的，那么因为目前 \(p_i\) 的 \(x\) 坐标最大，它就一定是在目前遍历过的convex hull上，所以我们就需要从 \(p_{i-1}\) 开始向后检查，每次删除最后3个点的中间的点（即每次的倒数第二个点），做重新调整。

    先删除 \(p_{i-1}\) 这个点，然后检查此时的最后三个点，\(p_{i-3}\), \(p_{i-2}\), \(p_i\)，如果它们组成了**“右转”的折线，那么本次调整到此结束，然后继续加入下一个点 \(p_{i+1}\)；如果它们组成了“左转”的折线，那么就需要再次删除中间点，即\(p_{i-2}\)，然后继续检查时的最后三个点，\(p_{i-4}\), \(p_{i-3}\), \(p_i\)，并重复上述步骤，直到最后三个点组成“右转”**的折线（或者直到剩下最后2个点），本次调整才到此结束，然后继续加入下一个点 \(p_{i+1}\)。

当针对上述两种情况做完调整之后，此时继续加入下一个点 \(p_{i+1}\)，并重复上述步骤，直到加入最右边的点 \(p_N\)，此时就得到了 upper hull。

寻找 lower hull 的incremental的步骤和上述类似。

第二种算法简述

（这个算法实际上是Andrew对Graham’s scan的一种改进算法）

输入：平面上点的集合 \(\mathcal{P}\)。

输出：一个点的序列 \(\mathcal{L}\)，表示点集合 \(\mathcal{P}\) 的Convex Hull，点序是顺时针方向。

算法简述：

  • 将集合 \(P\) 按 x 坐标排序为 \(p_1\), \(p_2\), …, \(p_N\)

  • 把\(p_1\), \(p_2\) 放入序列 \(\mathcal{L}\)，并且 \(p_1\) 是第一个点， \(p_2\) 是第二个点

  • 变量 \(i\)，值从 \(3\) 到 \(N\)，依次遍历加入序列 \(L_1\)，每次加入点 \(p_i\)，检查最后三个点是否组成**“右转”的折线段。如果是，继续遍历下一个值，否则删除当前序列 \(L_1\) 的倒数第二个点，并继续检查最后三个点是否组成“右转”的折线段，以此类推，直到当前序列 \(L_1\) 的最后三个点组成“右转”**的折线段，才继续遍历下一个 \(i\) 值。

  • 当变量 \(i\) 遍历完成时，就得到了convex hull的上半部分 upper hull的点序列是 \(L_1\) 。

  • 把\(p_N\), \(p_{N-1}\) 放入序列 \(L_2\)，并且 \(p_N\) 是第一个点， \(p_{N-1}\) 是第二个点

  • 变量 \(j\)，值从N-2到1，依次加入序列 \(L_2\)，和上面寻找upper hull的办法类似，仍然是确保每加入一点后，调整序列 \(L_2\)，使得其最后三点组成**“右转”**的折线段，然后才继续遍历下一个 \(j\) 值。

  • 把序列 \(L_2\) 的第一个和最后一个点去掉，避免重复点。

  • 把序列 \(L_1\) 和序列 \(L_2\) 合并，即得到最终的点序列 \(\mathcal{L}\)。

时间复杂度：\(O(nlogn)\)

对于该算法的几点说明

  • 在排序时，如果 \(x\) 坐标相同，可以按照 lexicographic 的办法排序，即先按照 \(x\) 坐标排序，如果 \(x\) 坐标相同，就再按照 \(y\) 坐标排序（仅对 \(x\) 坐标相同的点的情况下）。

  • 在上面判断最后三点是否组成**“右转”的折线段时，如果这三点共线，仍然把这种情况归为”左转“**的情况，从而触发删除三点里面中间点的操作处理。

  • 因为使用的是floating point calculation，并且依然存在rounding error，所以最后的点列表，有一定概率并不是实际上真正的convex hull的点列表（比如有三个点靠的很近以至于是一个左转的折线段，但被计算为右转了），但这种结果是可以接受的。

1.1.3 计算convex hull的时间复杂度

Theorem 1.1 The convex hull of a set of n points in the plane can be computed in O(nlogn) time.

关于第二种算法正确性的证明，文中采用了数学归纳法。

以upper hull为例，假如现已有点列 {\(p_1\), \(p_2\), …, \(p_{i-1}\) }，准备加入点 \(p_i\)。根据算法，点列{\(p_1\), \(p_2\), …, \(p_{i-1}\)}中最后三点一定是组成**“右转”**的折线段（即除了这些点，到目前最大的点为止，其他点都在这些点的下方）。我们把此时的upper hull点列叫做 old chain。

在加入点 \(p_i\) 之后，按照字典序（lexicographic），最小的点是 \(p_1\)，最大的点是 \(p_i\)，经过调整，此时的upper hull我们叫做new chain（而且new chain的最后一个点一定是 \(p_i\)）。

可以断言的是old chain一定是在new chain的下方（有可能点pi就是old chain的延伸，但是在算法中，这种共线的情况被当做是左转而被排除掉了）。

按照算法，我们需要证明的是，到目前为止，除了{\(p_1\), \(p_2\), …, \(p_i\)}，所有的点都在new chain的下方。

假如有一个点位于new chain的上方，那么这个点就必须介于 p(i-1) 和 pi之间，因为在加入 \(p_i\) 之前，所有的点都位于old chain的下方。但这又是矛盾的，因为 \(p_{i-1}\) 和 \(p_i\) 之间没有其它点，因为所有点已经是按照字典序排列过了的。

因此归纳出来，到目前为止，除了{\(p_1\), \(p_2\), …, \(p_i\)}，所有的点都在new chain的下方。算法正确性得到证明。

关于时间复杂度的证明。

对于upper hull，按字典序排序，时间复杂度是\(O(nlogn)\)。

for循环是线性的，关键在于其里面用于检查右转折线段和删除中间点的while循环的执行次数。

这个while循环首先可以肯定至少执行一次（检查右转折线段），而额外执行的次数，是为了删除每次得到的序列最后三点的中间点，而因为所有点只会被加入序列一次，所以，每个点最多也只会被删除一次，那么这个for循环里面的while循环执行的上限就是\(O(n)\)。

所以，带有while循环的这个for循环，时间复杂度是\(O(n)\)，而不是\(O(n^2)\)。

因此计算upper hull的时间复杂度就是\(O(nlogn)\)。

对于lower hull也是类似的。所以加起来，整个算法的时间复杂度就是\(O(nlogn)\)。

1.2 Degeneracies and Robustness

提出算法的三个步骤（阶段）

  • 首先，排除次要因素的干扰，因为这些因素是细节问题，不影响算法的整体思路。

  • 其次，再考虑前面可能出现的退化情况（边界条件，特殊和极端情况等问题），调整算法细节以便处理。

  • 最后，实现细节。比如原子操作，如何遍历等等。

比如，在convex hull的算法中，我们可以先假设没有三个共线的点，没有两个点的 x 坐标是相同的。

symbolic perturbation schemes指在设计和实现阶段忽略了special case，但在实际应用过程当中算法仍然正确的方法。

在实现细节的阶段，使用实数（浮点数）计算可能导致假设在某种情况下失效的问题，这是算法健壮性的体现。就像前面第二种算法中提到的，最终的output也许不是真正意义上的真实结果，但也是十分接近真实的结果，在这种情况下，需要预期这种情况可能的后果，并避免有次可能产生的crash问题等等。

使用现有的arithmetic library是其中一种办法，如果不能达到我们所需要的要求，就需要自己实现一些特定情况下的处理。

1.3 Application Domains

这一节主要介绍了Computational Geometry的几种应用领域，已经每个领域要解决的问题。

  • Computer graphics

  • Robotics

  • Geographic information systems

  • CAD/CAM

  • Other applications domains （比如 molecular modeling，pattern recognition等）

1.4 Notes and Comments

本节主要是对本章内容的一些延伸以及参考书籍资料等出处说明，提到了本章算法的来源，其发展的简要历史，以及相似算法的研究情况。

比如，本章所讨论的convex hull问题是Computational Geometry的经典问题，而本章第二种算法，其实是Graham’s scan算法，是Andrew基于最早的Graham提出的算法的改进。

还有其他的一些算法，时间复杂度也是\(O(nlogn)\)。

2 Line Segment Intersection Thematic Map Overlay

引言部分，以旅游为例，讲述了在实际当中，可能需要查看包含不同信息类型的地图，从而找到所需的信息。

在GIS领域中，layer 是指包含某一种信息的地图（map），而需要将多种类型的地图进行交叉引用的合并结果，叫做 overlay。

比如，一个layer（map）只包含城市名的信息，另一个layer（map）只包含河流的信息，还有一个layer（map）只包含了铁路轨道的信息，诸如此类等等。

当查看了城市信息的layer（map）之后，想要得知如何前往，就需要和另一个包含道路信息的layer（map）重叠查看，就是overlay。

GIS中，在overlay上，不同信息有交叉的地方（比如查看河流和道路的重叠情况），有时是一个交叉点，有时是一个交叉的区域。

2.1 Line Segment Intersection

本节要解决的问题是，给定二维平面上一个有 n 个线段的集合，找出所有的交点。

given a set S of n closed segments in the plane, report all intersection points among the segments in S.

其中线段的端点碰到其他的线段，也算作交点。

Brute-forced algorithm的时间复杂度是\(O(n^2)\)，但实际情况，有可能只有很少的一些线段相交，并不必计算每个线段和其他线段的交点。

即，我们希望算法的复杂度依赖的不仅是输入点的个数，而且也是输出的交点的个数，这样的算法叫做output-sensitive algorithm。

可以利用的观察几何结果是：靠的比较近的线段是可能有交点的候选计算对象，而相离较远的线段是不需要计算交点的。

所以思路是，把所有线段向y轴做投影，得到投影线段有重叠的那些线段，就是需要计算交点的候选线段。

为什么没有投影重叠的线段就一定没有交点？这可以通过反证法得出，如果没有投影重叠的线段有交点，那么这个交点的y坐标值一定是介于两个线段的4个端点的y值之间，而这又说明这两条线段是有投影重叠的，因此矛盾，从而的证。

使用到的技术叫做：plane sweep algorithm。

sweep line：一条水平无限长的假想虚线

status：sweep line的“状态”指的是和它当前相交的线段的集合（segments）

event point：sweep line沿着垂直方向从上向下移动，但不是连续移动的，而是离散的，移动到的这些位置的点，叫做event point。这些event point，一部分是每条线段的upper end point（y值较大的点）和lower end point（y值较小的点），另一部分是线段之间的交点。

只有当sweep line移动到这些event point上的时候，算法才做相应的计算或调整，即更新sweep line的status，并测试线段之间有无交点（如有，就计算交点）。

  • 如果event point是一条线段的upper point，那么这条线段就是和sweep line相交，并且应该加入到status里面，同时要计算这条segment和status里面其他segments的交点（后面会提到，只计算当前线段相邻的左右两条segments的交点，而不是计算和status里面所有线段的交点），而且这个交点（如果有）要放入到event point集合的适当位置，以便sweep line依次向下扫描时可以遍历到它。

  • 如果event point是一条线段的lower point，那么这条线段就和sweep line不再相交（即变为相离），就应该从status里面删除。而且这也会导致status里面原先不直接相邻的两条线段，现在变成了直接相邻了，那就要计算这两条相邻线段之间有无交点（如果有，依然要放入event point集合里面去）

  • 如果event point是两条线段的intersection point（这个intersection point是前面计算得到加入进来的），那么在该点之后，相邻的adjacent neighbor就会发生改变，所以就要测试（计算）这两条segments和它们各自左右相邻的segment的交点。

Lemma 2.1 Let si and sj be two non-horizontal segments whose interiors intersect in a single point \(p\), and assume there is no third segment passing through \(p\). Then there is an event point above \(p\) where si and sj become adjacent and are tested for intersection.

因为根据前面遇到的event point是一条线段的upper point时的操作（计算adjacent segment之间的intersection point），这个引理主要想说明，如果两条都不是水平（也不共线）的线段，如果有交点，那么在这个交点的上方，一定有一个event point，在那个event point的时候，这两条线段变成adjacent，并且会被检查（计算）是否有交点。

这里暂时忽略了三种特殊情况：两条线段可能共线（重合），可能有水平的情况，以及有第三天线段穿过交点。

所以，简要叙述，line sweeping algorithm的大体思路如下

Let’s briefly recap the overall approach. We imagine moving a horizontal sweep line ℓ downwards over the plane. The sweep line halts at certain event points; in our case these are the endpoints of the segments, which we know beforehand, and the intersection points, which are computed on the fly.

While the sweep line moves we maintain the ordered sequence of segments intersected by it. When the sweep line halts at an event point the sequence of segments changes and, depending on the type of event point, we have to take several actions to update the status and detect intersections.

假设有一条水平扫描线，从上而下移动，每次移动到一个特殊的点（event point）。这样的event point有两种，一种是每条线段的upper point（end point），另一种是某两条线段的交点（intersection point）。前一种在计算之前就已知，而后一种是在扫描线移动过程中计算得出。

当扫描线移动时，维护一个有序的线段列表，列表中的每个线段是和扫描线相交的。当扫描线移动到下一个event point的时候，更新线段列表使其保持有序，同时根据event point的类型，更新状态（它是和扫描线相交的线段集合，每次操作有可能添加或删除一条线段）并检查某两条线段是否有交点。

sweep line遇到三种不同event point时对应的操作

  • 如果event point是一条线段的upper point（end point），就要检查这个upper point所在的线段，和它左右两个相邻的线段是否有交点，如果有交点，那么这个交点就是一个新的event point。当然，upper point所在的线段要放入status中去。

    因为sweep line上方的event point都是已知的或已经计算过的，所以关注的是sweep line下方的交点。

  • 如果event point是某两条线段的交点（intersection point），那么这两条线段在所维护的有序线段列表（status）里面的位置就要交换，同时因为位置变化，它们各自相邻的线段也发生了变化（但只变化了一个，因为另一个仍然是它们自己中的一个），所以也要检查它们和各自新邻近的线段之间是否有交点，如果有并且是之前没有的event point，那么就有发现了一个或两个新的event point。

  • 如果event point是一条线段的lower point（end point），那么这条线段原先左右两条线段就变成了直接相邻的线段，就要检查（计算）这两条线段是否有交点，同样的，如果有，就是新的event point。当然，这个lower point所在的线段要从status里面移除出去。

算法当中需要的两个数据结构

  • event queue（记作 \(\mathcal{Q}\)）

    需要支持删除一个点（event point）的操作，并返回这个点以便对其处理。

    （如果两个点有相通的y坐标，返回x坐标较小的一个。这个实际上说明，如果一个线段是水平时，当水平的sweep line扫描到这条线段时，upper point是其左边的点，lower point右边的点，即sweep line先遇到的event point是左边的点。）

    需要支持插入一个点（event point）的操作，因为新的intersection point是在sweep line移动过程中计算得出。

    同时，允许两个event point是共点的（coincide，比如两条线段的upper point可能是同一个点），但把它们当做是同一个点，所以需要支持查看一个event point是否在\(\mathcal{Q}\)中已经存在。

    根据上述特点，采用平衡二叉搜索树（Balanced Binary Search Tree，BST），并定义点（event point） \(p\) < \(q\) 的“小于”操作符（<）为

    （1）如果 \(p\) 和 \(q\) 的y坐标相同，那么 \(p\) 的x坐标小于\(q\) 的x坐标

    （2）如果 \(p\) 和 \(q\) 的y坐标不相同，那么 \(p\) 的y坐标小于\(q\) 的y坐标

    需要删除一个点的操作的原因是，sweep line向下移动时，需要event point的顺序，移动到下一个event point上，而这是二叉树删除一个节点并返回的操作（同时二叉搜索树会重新平衡并排序）

    需要插入一个点的操作的原因是，当sweep line移动到不是intersection point的event point的时候，要计算相邻两条线段之间的intersection point，如果有就要插入BST，所以这是BST的插入节点的操作。

  • status（记作 J）

    这个所谓的状态，是指当前和水平的sweep line相交的线段的有序集合。

    对于给定的一条线段，为了计算它和相邻线段的相交情况，它必须是可以动态调整的，即：

    （1）当sweep line遇到一条线段的upper end point的时候，该线段需要放入status，并且需要查看此时它和左右相邻的两条线段的相交情况，如果有交点就需要计算出来，并放入event queue里

    （2）当sweep line遇到一条线段的lower end point的时候，该线段需要从status中移除，同时它原先左右相邻的两条线段现在变为直接相邻，那么也要再次查看并计算这两条线段是否有交点，如果有，同样放入event queue里

    （3）当sweep line遇到的event point是intersection point的时候，那么就需要交换这两条相交的线段在status中的位置，同时在status中，它们各自分别有一条相邻的线段发生了变化，同样需要再查看并计算交点，如果有交点，同样放入event queue里

    同样根据上面的特点，也采用平衡二叉搜索树（Balanced Binary Search Tree，BST），但这里的BST里面，只有叶子节点是存储了线段的信息，而树中间的每个节点（interior nodes），存储的都是其左子树里面最右边（叶）节点的线段信息。

    虽然中间的节点也可以存储线段信息，但为了方便陈述算法，所以中间节点都是用来引导寻找最终叶节点的导引信息（values to guide the search），而不是最终的线段数据信息（data item）。

FindIntersections算法简述

输入：平面上线段的集合 S。

输出：交点的集合 \(\mathcal{L}\)，这些交点都在集合 S 中的某些线段上，同时每个交点还有其对应的线段信息，表示该交点位于哪（几）条线段上。

算法简述：

首先，初始化一个空的event queue，记作 \(\mathcal{Q}\)。然后把集合 S 里线段的end points都插入到 \(\mathcal{Q}\) 中，当一个end point是线段的upper point时，要同时带上其所在的线段的信息（属于那条线段）。

然后，初始化一个空的status 数据结构，记作 J。

之后，依次遍历 \(\mathcal{Q}\)，每次从 \(\mathcal{Q}\) 中返回下一个event point \(p\)（同时 \(p\) 从 \(\mathcal{Q}\) 中被移除），然后根据event point \(p\)，调用对 \(p\) 的处理函数HandleEventPoint(p)（如下）。这个遍历的终止条件是 \(\mathcal{Q}\) 为空。（即这是一个while循环，而在遍历过程中可能有新event point加入 \(\mathcal{Q}\) ）

算法复杂度：

O((n+k)logn)，其中，n是输入线段个数，k是输出个数

或者更具体地，O((n+I)logn)，其中，n是输入线段个数，I是交点个数

HandleEventPoint(p) 步骤简述

  • 输入是点 \(p\)

  • 记 upper end point为 \(p\) 的线段集合为 U(p)，这些线段是和点 \(p\) 对应存储的。如果线段是水平的，它的upper end point是左边的端点。

  • 在status J 中找到所有包含点 \(p\) 的线段，它们都是相邻的，记 L(p) 是lower endpoint为 \(p\) 的线段集合，记 C(p) 是线段中间包含点 \(p\) 的线段集合（即点 \(p\) 是它们之间某两条或几条线段的交点）。

  • 如果 L(p) ∪ U(p) ∪C(p) 至少有一条线段，就说明点 \(p\) 是一个交点

    • 报告这个结构，并同时报告它所在的线段（在L(p)， U(p) 和 C(p) 中）

  • 从status J 中删除L(p) ∪C(p) （即它们的并集）

  • 向status J 中添加U(p) ∪C(p) （即它们的并集），并且插入的这些线段的顺序是，按照它们和sweep line在 \(p\) 稍下方一点位置相交的顺序。如果有线段是水平的，那么它要排在其他线段的最后面。

  • 从上面的两个步骤可以得到，删除了C(p) 又添加了C(p) ，那么C(p) 中的线段在status J 中的顺序逆序了。

  • 如果U(p) ∪C(p) （即它们的并集）为空集

    • 把 status J 中，在 \(p\) 点左右两边的线段记为 sl 和 sr，调用寻找event point的函数FindNewEvent(sl, sr, p)

      如果 sl 或 sr 不存在，就忽略此步骤。

  • 如果U(p) ∪C(p) （即它们的并集）不是空集

    • 把既在 U(p) ∪C(p) 中又在status J 中，最左边的线段记作 \(s_1\)，把在status J 中 \(s_1\) 左边的线段记作 sl，然后调用寻找event point的函数FindNewEvent(sl, s1, p)。

      如果 sl 不存在，就忽略此步骤。

    • 把既在 U(p) ∪C(p) 中又在status J 中，最右边的线段记作 \(s_2\)，把在status J 中 \(s_1\) 右边的线段记作 sr，然后调用寻找event point的函数FindNewEvent(s2, sr, p)

      如果 sr 不存在，就忽略此步骤。

FindNewEvent(sl, sr, p) 步骤简述

  • 如果线段 sl 和 sr 在sweep line的下方相交，或者就在sweep line上相交并且在当前event point \(p\) 的右边，那么这个新的交点就是在 \(\mathcal{Q}\) 中还没出现的新的event point

    • 把这个新的交点加入到 \(\mathcal{Q}\) 中

Lemma 2.2 和Lemma 2.3 分别是这个算法的正确性，以及算法的时间复杂度的证明。

根据这两个引理，得出Theorem 2.4。

Lemma 2.2 Algorithm FINDINTERSECTIONS computes all intersection points and the segments that contain it correctly.

Lemma 2.3 The running time of Algorithm FINDINTERSECTIONS for a set S of n line segments in the plane is O(nlogn+I logn), where I is the number of intersection points of segments in S.

Theorem 2.4 Let S be a set of n line segments in the plane. All intersection points in S, with for each intersection point the segments involved in it, can be reported in O(nlogn+I logn) time and O(n) space, where I is the number of intersection points.

2.2 The Doubly-Connected Edge List

引出了可平面图（planar graph, or planar embedding graph）的概念，引出可平面图的点（vertex）、线（edge）、面（face）。

同时引出了我们需要的应用，即确定哪个面（face）是包含所给定的一个点（given point）的。

引出了数据结构 doubly-connected edge list，即 doubly-connected edge list 包含了一个平面细分（subdivision）上的face，edge，vertex的记录（record），并且除了几何和拓扑信息外，可能还有一些其他额外的信息，这个额外的信息叫做 information attribute （例如，一个face可能代表的是一种植被的覆盖，那么这个植被的种类就可以是这个额外的信息）。

这个 doubly-connected edge list 数据结构上的几何与拓扑信息，需要允许我们支持以下的一些操作

  • 逆时针遍历这些face的edges，同时也能容易地反方向（顺时针）遍历。（这就要求edge直接有指向前一个和后一个的指针）

  • 因为一个edge是两个face的边界，所以edge上需要有两个指针来指向这两个face

  • 为了更方便表示当前描述的edge是哪个face的edge，可以把一条edge拆解为两条 half-edge

    • 这两条half-edge是不同face的，而且每个half-edge都有唯一的指向前一个half-edge和执行后一个half-edge的指针

    • 而这同样意味着，一条half-edge只属于同一个face

    • 对于同一条edge的两条half-edge，我们把它们叫做 twins

    • 我们把half-edge定义为有方向的，沿着half-edge走，face就在它的左边，所以这个方向是逆时针

    • 把half-edge定义为一个向量，origin（起点）是v，终点（destination）是w。所以它的twin half-edge的起点就是w，而终点是v。

    • 根据上面的定义，为了访问face的边界，可以只存储一个指向half-edge的指针，这样就可以沿着逆时针方向遍历这个face的所有half-edge了。

  • 为了在表示洞（hole）时，仍然有沿着half-edge走时，face还在它的左边，就把洞的half-edge的方向定义为顺时针。

    而且，为了表示洞，需要需要有两个指向half-edge的指针，一个逆时针表示包含洞的face的边界，一个顺时针表示洞本身。

  • 还可以存储多个half-edge的指针，而且这些指针沿着这些edge遍历起来的时候，没有重复的edge，这就是isolated island的形式（为了简化期间，书中暂时不作讨论）

总结起来，doubly-connected edge list数据结构有三种记录数据（record）

  • vertex record

    它用来记录每个vertex（记作v）的坐标Coordinate(v)，并且它还有一个指针\(IncidentEdge(v)\)指向一条half-edge，而且这条half-edge的起点就是v

  • face record

    一个face（记作\(f\)）

    • 存储一个指针\(OuterComponent(f)\)，指向的是outer boundary的half-edge。（如果face是unbound，即open edges的话，这个指针就是空？）

    • 还存储一个指针\(InnerComponent(f)\)，指向的是inner boundary的half-edge，这是用来表示洞的

  • half-edge record

    一个half-edge（记作\(e\)）

    • 存储一个指针\(Origin(e)\)指向它的起点（origin）

    • 存储一个指针\(Twin(e)\)指向它的twin half-edge

    • 存储一个指针\(IncidentFace(e)\)，表示它绑定（bound）的face

    • 存储一个指向它前面half-edge的指针\(Prev(e)\)

    • 存储一个指向它后面half-edge的指针\(Next(e)\)

    没有必要存储它的终点（destination），因为可以通过\(Origin(Twin(e))\)得到。

本节还画了vertex，edge，half-edge，face以及上面提到的各种record的示意图，如下。

（图暂时省略，图位于第32页，页码是41）

这里也提到了，有时候有些record在一些应用中不是必须的（比如river和road构成的face，在某些应用中没有太多意义），所以在实现的时候可以适当忽略，以便在算法实现中更方便地调整其他数据。

2.3 Computing the Overlay of Two Subdivisions

简而言之，计算两个subdivision的overlay，就是根据两个subdivision的doubly-connected edge list（记作S1和S2），计算出一个新的doubly-connected edge list表示的subdivision（记作\(O(S_1, S_2)\) ）。

（此处的图为，Figure 2.4，Overlaying two subdivisions）

这个overlay，可以看做是S1的edges被S2的edges所切割，而S1中的大部分edge其实可以在新生成的doubly-connected edge list中来复用，仅那些被S2的edges所真正切割到的S1的edges，才需要在新生成的\(O(S_1, S_2)\) 被更新。

为了计算overlay结果，要把两个doubly-connected edge list（S1和S2），拷贝到一个新的doubly-connected edge list中去。拷贝的结果当然不是一个合法的doubly-connected edge list，因为它不能代表一个平面的细分（subdivision）。overlay算法的任务就是，把这个不合法的doubly-connected edge list，通过计算两个network edges之间的交点，并把两个doubly-connected edge list的部分区域连接起来，从而最终得到一个合法的doubly-connected edge list，即结果\(O(S_1, S_2)\) 。

下面首先讨论的是，最终的overlay结果\(O(S_1, S_2)\) 中的vertex和half-edge records，是如何被计算出来的。（关于新生成的face record，因为比较复杂，稍后再讨论）

计算\(O(S_1, S_2)\) 的办法，利用了前面提到的计算line segments交点的plane sweep algorithm。算法操作的对象是，包含了S1和S2中所有line segment的线段集合（一个新的线段集合拷贝）。

在plane sweep algorithm中，需要两个数据结构，分别是event point的集合Q，以及status structure J。

Q是用来存储event point的（BST实现），而J是用来存储和sweep line相交的那些line segment的集合的（是有序的，在plane上是从左向右依次和sweep line相交的，也是BST实现的）。

除了这个两个数据结构之外，还需要维护一个doubly-connected edge list的数据结构\(D\)，它的初始值是从S1和S2拷贝而来，也就是说它的初始值是包含了S1和S2的所有line segment的集合。而随着sweep line的向下移动，\(D\)会随之而更新，最终变成一个合理的doubly-connected edge list。

如果一个\(D\)中的edge和sweep line相交而要被放入status J中时，我们需要用指针把放入J中的edge和它来自于\(D\)中的哪个half-edge record联系起来，这样当遇到一个intersection point时，我们就能够方便地找到\(D\)中的哪一个half-edge record（或哪一部分）需要被更新和调整。

在sweep line向下扫描的过程中，sweep line上面是已经计算好的最终overlay结果的一部分，是不再变化的。

当遇到一个event point时候的处理：当event point是来自原先同一个subdivision的edges时，那么这个event point是可以被复用的；但如果event point是来自原先两个subdivision的不同edges时，那么我们就需要更新数据结构\(D\)，更新（加入或删除）某些edges，以便把两个subdivision通过新的intersection point而连接起来。

这里通过举例，说明了一个subdivision中的一条edge，是如何和另一个subdivision中的其他几个edge相交，然后做处理的。这个过程比较tedious，但是不难（difficult）

（图为Figure 2.5，图位于第35页，页码是44）

这里主要结合图形，说明了在新生成了两条edge（对应的是两队half-edge pair）之后，如何调整它们以及周围的edge的Next()和Prev()指针。

值得说明的是，这个例子中，一条edge恰好经过的是另一个subdivision的一个vertex，因此，在调整新产生的edge的prev和next的时候，是按照clockwise的转向，找到第一个相邻的edge作为Next()指针所指向的edge，而按照anti-clockwise的转向，找到其第一个相邻的edge作为Prev()指针所指向的edge。这个可以结合图的说明清晰容易地看到。

除了更新生成的新half-edge pair，还要找到\(O(S_1, S_2)\) 中每个face \(f\) 的 \(OuterComponent(f)\) （指向一个表示outer boundary的half-edge）和\(InnerComponent(f)\) （指向一个或几个half-edge的指针，表示一个或多个洞）。还要给每个edge的\(IncidentFace()\)设定合理的指针指向face record。最后，每个\(face\)还要用原先两个subdivision中包含这个\(face\)的face name来给它做label。

如何判断一个half-edges组成的boundary是outer boundary，还是表示hole的inner boundary？

选定leftmost的vertex（in case of ties，choose lowest of leftmost），因为沿着half-edge的走向是clockwise的就是outer boundary，所以计算这个vertex前后两个相邻的（有序的）half-edge的夹角，如果是小于\(90°\)，那么就是outer vertex的half-edge，如果是大于\(90°\)，就是inner boundary的half-edge。这个特性仅适用于leftmost（或lowest of leftmost if ties）的vertex。

（这里的图位于第36页，页码是45）

通过一个图的例子，说明了如何确定一个face \(f\)是由一个或几个cycle组成的。如果是多个cycles组成，一般有几个洞的cycle（half-edges是顺时针的）和一个outer cycle（for outer boundary）组成，而且一个洞要通过对应的数据结构（比如class上的成员变量）连接到另一个洞或outer boundary上，这样才能表明这些cycles组成的是同一个face \(f\)。

Lemma 2.5 Each connected component of the graph \(G\) corresponds exactly to the set of cycles incident to one face.

关于这个lemma的证明，没看懂。

总之，他想说明的是，一个face上的洞，是和这同一个face上的其他洞相连的，或者是和这个face对应的\(OuterComponent(f)\) 相连接，而这些相连接的洞（实际上就是\(InnerComponent(f)\) ？）和\(OuterComponent(f)\) 就组成了这个face \(f\)。

如果构建graph \(G\)？

构建graph \(G\)，实际上是把这些\(InnerComponent(f)\)（即洞）和\(OuterComponent(f)\) 直接合理地用书中所谓的“arc”连接起来。

对于每个表示洞的cycle的leftmost的vertex \(v\)，如果有一条half-edge \(e\)，是这个vertex \(v\) 左边第一个邻近的half-edge，那么就在这两个node直接就用一条arc连接起来。

为了快速（有效）地找到这些node，每个half-edge的record上有指针指向这些node，表示这些node在这个graph \(G\) 的哪个cycle上。

而找到一条vertex左边的、相邻的第一个half-edge，是在plane sweep algorithm中sweep line向下扫描时得出的，而且这个相邻的左边第一个half-edge，是位于另外一个cycle上的。

（这里用来说明的图，位于第37页，页码是46）

最后一件事情是，在overlay结果 \(O(S_1, S_2)\) 中，每个face \(f\) 都要找到它原先分别在 \(S_1\) 和 \(S_2\) 中的label。

假如一个vertex \(v\) 是来自 \(S_1\) 的一条edge \(e_1\) 和 \(S_2\) 的一条edge \(e_2\) 的相交得到新的点，那么可以从edge \(e_1\) 和 \(e_2\) 的 \(IncidentFace(f)\) 得到各自在原先 \(S_1\) 和 \(S_2\) 中的label name。

但如果vertex \(v\) 本身就是来自 \(S_1\) 的一个点（或者\(S_2\) 的一个点），那么我们首先能得知它来自 \(S_1\) 的哪个face（因为能从 \(v\) 对应的half-edge的 \(IncidentFace(f)\) 上得到。其次，就需要找到在 \(S_2\) 上的哪个face包含这个vertex \(v\)。

书中在此处没有展开解释，只说明了仍然使用本章介绍共的plane sweep algorithm就可以找到，而且也不用再次调用这个plane sweep algorithm，而是在原先扫描的过程中，就可以找到。

MapOverlay算法简述

输入：二维平面上的两个平面细分（subdivision）\(S_1\) 和 \(S_2\)，它们都是以doubly-connected edge list表示。

输出：\(S_1\) 和 \(S_2\) 的overlay \(D\)，并且也是以doubly-connected edge list表示。

算法简述：

  1. 新建一个doubly-connected edge list \(D\)，并把两个原始输入\(S_1\) 和 \(S_2\) 拷贝到 \(D\) 中

  2. 通过第2.1节中提到的plane sweep algorithm，计算\(S_1\) 和 \(S_2\) 中每个edge的交点，除了在每个event point时更新 \(J\) （status）和 \(Q\) （event point），还需要处理

    • 更新步骤1中建立的的doubly-connected edge list \(D\)（前面的叙述中有举例如果\(S_1\)的一条edge穿过了\(S_2\)的一个vertex时，如何生成新的half-edge pair，以及复用原先的half-edge pair并调整相应的record指针）

    • 在处理\(D\)中的每个event point之后，记录每个event point左边第一个half-edge的信息

  3. 经过步骤2，\(D\)已经是\(S_1\) 和 \(S_2\) 的overlay结果\(O(S_1, S_2)\)，但是每个face \(f\) 的信息还没有计算出来

  4. 遍历\(D\)，确定\(O(S_1, S_2)\)中的boundary cycles

  5. 构建graph \(G\)。这样的 \(G\) 是一个或多个component组成。每个component由一个或几个表示boundary的cycle(s)组成，如果一个boundary cycle表示的是洞，那么它的leftmost的vertex就要通过一个所谓的”arc”连接到另外一个表示洞的boundary cycle（或者最终连接到一个表示非洞的boundary cycle上）

  6. 对于步骤5中建立的graph \(G\) 的每个component：

    假设 \(C\) 是这个component中唯一的outer boundary cycle，并用 \(f\) 表示由这个cycle所包含的face。

    创建 \(f\) 的face record，设定指针 \(OuterComponent(f)\) 指向 \(C\) 中的某一个half-edge即可；设定指针数组（或列表）\(InnerComponent(f)\) ，它是这个component中每个洞上的某一个half-edge的指针集合；把这个component中每条half-edge所指向face的指针\(IncidentFace(e)\) 设置为指向 \(f\) 的face record。

  7. 结果\(O(S_1, S_2)\)中的每个face，都用\(S_1\) 和 \(S_2\) 中对应的face名字做标记（label）

算法时间复杂度：\(O(nlogn + klogn)\)

Theorem 2.6 Let \(S_1\) be a planar subdivision of complexity \(n_1\), let \(S_2\) be a subdivision of complexity \(n_2\), and let \(n := n_1 +n_2\). The overlay of \(S_1\) and \(S_2\) can be constructed in \(O(nlogn + klogn)\) time, where k is the complexity of the overlay.

算法复杂度的计算

  • 步骤1中，拷贝两个doubly-connected edge list 到 \(D\) 中，算法复杂度是 \(O(n)\)

  • 步骤2中，plane sweep algorithm的时间复杂度是 \(O(nlogn + klogn)\)

  • 步骤4-6中，用来填写face record的时间复杂度是和 \(O(S_1, S_2)\) 线性相关的

  • 步骤7中，把结果中的每个face用\(S_1\) 和 \(S_2\) 中对应的face名字做标记的时间复杂度是 \(O(nlogn + klogn)\)

2.4 Boolean Operations

Map overlay算法最为常见的应用之一，就是polygon的Boolean操作，即 与（\(AND\)， ∩），或（\(OR\)，∪），非（\(NOT\)，\）。

（这里用来说明的图，位于第39页，页码是30）

按照前面所述，把两个polygon看做是两个平面细分（subdivision），记作 \(P_1\) 和 \(P_2\)，那么map overlay的结果\(O(P_1, P_2)\) 是一个新的平面细分，并且也用一个doubly-connected edge list所表示。这里最重要的是，作为结果的平面细分的每个face record \(f\)，都是用原来两个平面细分 \(P_1\) 和 \(P_2\) 共同标识的。

所以，Boolean操作的求解转换为：

  • 如果计算的是\(P_1\) 和 \(P_2\)的交集（ \(P_1 ∩ P_2\)），我们就从overlay结果中找到那些同时带有\(P_1\) 和 \(P_2\) label的face。

  • 如果计算的是\(P_1\) 和 \(P_2\)的并集（ \(P_1 ∪ P_2\)），我们就从overlay结果中找到那些带有\(P_1\) 或 \(P_2\) 或同时带有\(P_1\) 和 \(P_2\)label 的face。

  • 如果计算的是\(P_1\) 和 \(P_2\)的差集（ \(P_1 \ P_2\)），我们就从overlay结果中找到那些只带有\(P_1\) 、不带有 \(P_2\) label的face。

2.5 Notes and Comments

line segment intersection problem是计算几何中最为基础的问题之一。

本章提到的 \(O(nlogn + klogn)\) 时间复杂度的算法是1979年 Bentley 和 Ottmann 给出的。

求得所有线段交点的时间复杂度的下限是 \(\Omega(nlogn + k)\) ，当 \(k\) 值较大时，这样的算法不是最优的。多位研究者研究后，Clarkson&Shor这两人，和Mulmuley分别给出了randomized incremental algorithms，时间复杂度是\(O(nlogn + k)\)，而空间复杂度分别是\(O(n)\)和\(O(k)\)，而且这两种种randomized algorithms也可以用来计算curve。Balaban后来给出了第一种 deterministic algorithm，时间和空间复杂度分别是\(O(nlogn + k)\)和\(O(n)\)。

有一种叫做red-blue line segment intersection problem的问题，是line segment intersection problem的特殊情况。它是指两个line segment的集合（red segments和blue segments），每个segment集合内部两两segment之间没有交点，那么求解这两个集合之间的segment intersection，Mairson 和 Stolfi给出的算法时间和空间复杂度分别是\(O(nlogn + k)\)和\(O(n)\)。

实际上，red-blue line segment intersection problem的问题就是network overlay problem。

line segment intersection counting problem是计算线段交点个数的问题（而不是报告所有交点坐标），所以它的输出就是一个整型数，不依赖于交点个数算法的时间复杂度是\(O(n^\frac{4}{3}log^cn)\)，其中 \(c\) 是某个小值常数。

Plane sweep是设计几何算法中最为重要的范式之一。第3章plane sweep用它来处理polygon triangulation problem，第7章用它来计算Voronoi diagram（维诺图） of a set of points。本章提到的sweep line是一条水平的（虚拟）直线，在某些情况下，sweep line可能是其他的形式，比如第15章提到的可能是rotating line。plane sweep也可以用于更高维度的空间，这时叫做space sweep algorithms。

本章提到的用来存储平面细分（subdivision）的数据结构是the doubly-connected edge list，Muller 和 Preparata描述了这种数据结构。此外，还有Baumgart的the winged edge structure，Guibas 和Stolfi的the quad edge structure等。这些数据结构的差异总体上不大。

2.6 References

4 Linear Programming Manufacturing with Molds

本章引言部分通过塑料或铁器（合金）的铸造，引出了本章需要探讨的话题：对一个给定的铸件（casting），是否存在这样的模具（mold）使得铸件能够从中移出？

为了方便本章的讨论， 我们把要讨论的情况做了以下的几项限定。

  • 我们假定物件（铸件）是多面体（polyhedral）

  • 我们只考虑模具只有一个整体部件组成，不是由多个部件组成。

  • 我们只允许通过一次单独的变换，把铸件从模具中移出（不是旋转出来？）

4.1 The Geometry of Casting

在本章讨论的情况中，铸造的物体有一个水平的顶部平面（top facet），而且它是唯一一面不和模具接触的。

如果一个铸件至少能从模具的一个方向上移除，就把这个铸件称作可铸造的（castable）。

记 \(\mathcal{P}\) 是一个多面体，并假设模具是一个中间有空间可以和 \(\mathcal{P}\) 紧紧贴合的方块。它的顶部平面是和多面体 \(\mathcal{P}\) 共面。我们把多面体 \(\mathcal{P}\) 中不是顶部平面的其他面叫做普通平面（ordinary facet），每一个普通平面 \(f\) 在模具中都有和它对应的一面，记作 \(\hat{f}\) 。

我们需要做的是，决定多面体 \(\mathcal{P}\) 是否可以通过一次转换（translation）就可以从模具中移出。换句话说，我们要找出来是否存在一个向量 \(\overrightarrow{d}\) ，使得多面体 \(\mathcal{P}\) 沿着它可以移动到无限远而不碰到（intersecting）模具。

因为 \(f\) 是铸件的平面，同时它也是和模具平面 \(\hat{f}\) 贴合的平面。如果 \(\overrightarrow{d}\) 和平面 \(f\) 向外的法向量（记作 \(\eta(f)\)）的夹角小于90°，那么在移除是铸件就会被阻挡。所以必要的条件就是 \(\overrightarrow{d}\) 要和铸件 \(\mathcal{P}\) 所有平面的向外的法向量夹角小于90°。下面的定理4.1也说明了这个必要条件也是充分条件（即充要条件）。

定理4.1 如果一个向量 \(\overrightarrow{d}\) 和多面体 \(\mathcal{P}\) 所有面向外的法向量的夹角小于90°，那么这个多面体 \(\mathcal{P}\) 就能沿着方向 \(\overrightarrow{d}\) 从模具中移出。

Lemma 4.1 The polyhedron \(\mathcal{P}\) can be removed from its mold by a translation in direction \(\overrightarrow{d}\) if and only if \(\overrightarrow{d}\) makes an angle of at least 90° with the outward normal of all ordinary facets of \(\mathcal{P}\) .

证明见第72页，页码为64页。

此外，从定理4.1还可以得出一个有意思的结论：如果多面体 \(\mathcal{P}\) 可以通过一系列短的向量移出模具，那么它也能一次通过一个向量移出模具。

为了方便叙述，显然可以把所有讨论的向量的起始点移动到三维空间的原点。因为模具的顶部平面是向上开口的，所以我们要寻找的向量 \(\overrightarrow{d}\) 其实是正 \(z\) 轴的方向（想象一下要把铸件从开口处移除），那么可以取 \(z = 1\) 这个三维空间的平面，这个平面上面的任意一个点 \((x, y, 1)\) 就能表示一个起点是原点的向量。

假设铸件的每个面向外的法向量是 \(\overrightarrow\eta = ( \overrightarrow\eta_x + \overrightarrow\eta_y + \overrightarrow\eta_z )\)，那么一个向量 \(\overrightarrow{d}\) （\(\vec{d} = (x, y, 1)\)）和这个法向量夹角小于90°的充要条件是，它们的点积小于等于0（非正）。因此有下面的不等式（inequality）：

\[ \overrightarrow\eta_xd_x + \overrightarrow\eta_yd_y + \overrightarrow\eta_z \le 0 \]

这个不等式实际上表示了 \(z = 1\) 这个平面上的半个平面，即一条直线的左边平面，或者右边平面。（当然这个叙述不严谨，因为对于铸件面是水平的面时，有 \(\overrightarrow\eta_x = \overrightarrow\eta_y = 0\) 。但在本章讨论的例子中，这个条件不可能满足或者一开始就永远满足，所以可以不做讨论）。

换句话说，铸件的每个面代表了 \(z = 1\) 这个平面上的半个平面（half-plane），而这些平面相交的共同部分的任意一个点，代表了一个向量 \(\overrightarrow{d}\) ，沿着这个向量，铸件就可以从模具中移出。这就是说，我们把这个铸造问题转化为了求解一个平面上多个半平面交集部分的几何问题。如果一个铸件有 \(n\) 个面，那么就需要求解 \(n-1\) 个半平面交集部分（因为顶部平面不需要参与求解）。

Given a set of half-planes, find a point in their common intersection or decide that the common intersection is empty.

如果对于当前的模具，我们转化为的几何求解问题无解，那么可以尝试选择其他面作为顶部平面（top facet），从而再测试铸件（多面体 \(\mathcal{P}\) ）是否可以从模具中移出。

定理4.2 记 \(\mathcal{P}\) 是有 \(n\) 个面的多面体（铸件），以时间复杂度\(O(n^2)\)和空间复杂度\(O(n)\)，可以确定该多面体 \(\mathcal{P}\) 是否可铸造。此外，如果 \(\mathcal{P}\) 是可铸造的，那么可以以相同的时间复杂度和空间复杂度计算出一个模具和一个可以把 \(\mathcal{P}\) 从模具中移出的合法向量（方向）。

Theorem 4.2 Let \(\mathcal{P}\) be a polyhedron with \(n\) facets. In \(O(n^2)\) expected time and using \(O(n)\) storage it can be decided whether \(\mathcal{P}\) is castable. Moreover, if \(\mathcal{P}\) is castable, a mold and a valid direction for removing \(\mathcal{P}\) from it can be computed in the same amount of time.

4.2 Half-Plane Intersection

本章要讨论的问题，可以看做是对（二元一次）不等式方程组的求解

\[\begin{split} \left\{\begin{matrix} a_1x + b_1y \le c_1 \\ a_2x + b_2y \le c_2 \\ \vdots \\ a_nx + b_ny \le c_n \\ \end{matrix}\right. \end{split}\]

方程组中每个不等式就是一个线性约束（linear constraint），记作 \(h_i\)，它实际上代表了平面上的由一条直线分割的半个平面。其中 \((x,y) \in \mathbb{R}^2\)，即\((x,y)\)是二维实数域上的点。

记这样的一组线性约束为 \(H = \{h_1, h_2, \dots, h_n\}\)，它代表了 \(n\) 个半平面的交集。我们要求解的结果，就是属于这个交集的所有的点\((x,y)\)。

因为每个半平面是一个convex（凸包），而convex之间的交集还是convex，所以这些半平面的交集也是平面上的一个convex（凸包）。这个交集区域边界上的点也是某个半平面边界直线上的点，交集区域的每条边也是属于某个半平面边界直线。因为是convex，所以每条半平面边界直线顶多形成交集区域的一条边，所以最后交集的边至多有 \(n\) 条。

下面图中展示了几种交集的情况。

==（举例的图在74页，页码是67）==

首先给出分治算法来计算 \(n\) 个半平面的交集。

算法：\(IntersectHalfPlanes(H)\)

输入：二维平面上 \(n\) 个半平面的一组约束 \(H\) （\(H = \{h_1, h_2, \dots, h_n\}\)，其中 \(h_i = \{(x,y) \in \mathbb{R}^2 | a_ix + b_iy \le c_i \}\)）。

输出：一个凸包多边形区域 \(C:= \textstyle \bigcap_{h \in H} h\)。

算法步骤：

  • 如果 \(H\) 只包含有一个约束条件（即 \(n=1\)），那么 \(C\) 就是这唯一的一个半平面

  • 如果 \(H\) 包含大于一个约束条件（\(n > 1\)），就把 \(H\) 划分为两个子约束集合 \(H_1\) 和 \(H_2\)，每个集合有 \(\lceil n/2 \rceil\) 个约束条件。

  • 调用函数 \(IntersectHalfPlanes(H_1)\)，得到结果 \(C_1\)。

  • 调用函数 \(IntersectHalfPlanes(H_2)\)，得到结果 \(C_2\)。

  • 调用函数 \(IntersectConvexRegions(C_1, C_2)\)，得到结果 \(C\)。

这里函数 \(IntersectConvexRegions\) 实际上就是第二章里面的求解polygon交集的算法，时间复杂度是 \(O(nlogn + klogn)\)，\(n\) 是两个polygon的顶点个数。需要注意的一点是，polygon（区域）可能退化为一条线段或者一个点，或者这个polygon（区域）就是没有边界的。

假设通过递归调用已经计算出来了 \(C_1\) 和 \(C_2\)，因为它们各自都是由 \(n/2 + 1\) 个半平面计算得到，所以它们各自至多有 \(n/2 + 1\) 条edge，根据第二章计算overlay的算法，时间复杂度就是 \(O((n+k)logn)\)，\(k\) 是 \(C_1\) 和 \(C_2\) edge的交点个数。对于这样任意一个交点 \(v\)，它是来自于 \(C_1\) 中的一条edge \(e_1\) 和 \(C_2\) 中的一条edge \(e_2\) 相交得到的交点，因此这个交点也必然属于 \(C_1 \bigcap C_2\)。而 \(C_1 \bigcap C_2\) 是 \(n\) 个半平面相交得到的区域，因此它至多有 \(n\) 条边（也至多有 \(n\) 个顶点），故 \(k \le n\) 。因此时间复杂度就是 \(O(nlogn)\) 。

==（举例的图在75页，页码是68，从上到下第一个图）==

所以算法的复杂度可以用以下的式子表示

\[\begin{split} T(n) = \left\{\begin{matrix} O(1), if \space\space n = 1 \\ O(nlogn) + 2T(n/2), if \space\space n > 1 \end{matrix}\right. \end{split}\]

这个递归表达的最终结果就是 \(T(n) = O(nlog^2n)\)。

因为\(IntersectHalfPlanes(H)\)算法中用来计算两个polygon交集的程序默认两个polygon是任意的（即可以是凸的，也可以是凹的），但其实在这个算法中能保证两个polygon是凸的，这样就可以根据这个特性对算法加速。

为了方便叙述，需要精确描述一个凸的多边形区域 \(C\)（convex polygonal region，可以是无边界的）。把这个凸多边形区域 \(C\) 左侧和右侧的边，分别以排序的半平面列表来分开记录（存储）。在排序的列表中，每个半平面所对应的边界的出现顺序，就是从上到下遍历时依次出现的顺序。把左边的边界列表记作 \(\mathcal{L}_{left}(C)\)，把右边的边界列表记作 \(\mathcal{L}_{right}(C)\)。这里不显式记录顶点，它们可以根据两条连续出现的边的交点而计算得出。

==（举例的图在75页，页码是68，从上到下第二个图）==

下面为了简化叙述，我们假设没有水平的edge。但可以把顶部的水平的edge归到左侧的边界列表中，把底部的水平的edge归到右侧的边界列表中。

经过改进的（两个convex polygon region）算法也是plane sweep algorithm。从上到下移动sweep line，并且维护和sweep line相交的 \(C_1\) 与 \(C_2\) 中的边。因为 \(C_1\) 和 \(C_2\) 都是convex，所以和sweep line相交的边最多只有4条，因此不需要维护一个复杂的数据结构（BST），而是直接使用4个指针指向它们，\(left\_edge\_C1\)， \(right\_edge\_C1\)，\(left\_edge\_C2\) 和 \(right\_edge\_ C2\)。如果和左侧或右侧边界没有相交的edge，指针就指向空NULL。

初始时，记 \(C_1\) 最高的顶点的 \(y\) 坐标值是 \(y_1\)（如果 \(C_1\) 的上方是无边界的，就令 \(y = \infty\)）。类似的 \(C_2\) 最高的顶点的 \(y\) 坐标值是 \(y_2\)。取 \(y_{start} = min(y_1, y_2)\)，所以 \(C_1\) 和 \(C_2\) 所有的顶点都在 \(y_{start}\) 下方，因此sweep line从 \(y_{start}\) 处开始移动，\(left\_edge\_C1\)， \(right\_edge\_C1\)，\(left\_edge\_C2\) 和 \(right\_edge\_ C2\) 的初始值就是sweep line在 \(y_{start}\) 处时和 \(C_1\) 与 \(C_2\) 相交的边。

==（举例的图在76页，页码是69）==

根据plane sweep algorithm，通常需要一个queue来存储event。在本算法中，就是 \(C_1\) 和 \(C_2\) 中的边和sweep line开始（或终止）相交的点，也就是\(left\_edge\_C1\)， \(right\_edge\_C1\)，\(left\_edge\_C2\) 和 \(right\_edge\_ C2\) 的下端的端点中 \(y\) 坐标最高的那个点。（如果有两个 \(y\) 坐标相同的点，就按照从左向右的顺序）

当sweep line到达某个event point的时候，边界上就会出现一条新的边 \(e\)。处理的办法是，首先检查这条新的边 \(e\) 是属于 \(C_1\) 还是 \(C_2\)，再检查它是 \(C_1\)（或 \(C_2\) ）的左侧还是右侧边界，然后调用相应的处理程序。下面以这条边 \(e\) 是 \(C_1\) 左侧边界上的一条边为例，说明处理过程。

记端点 \(p\) 是这条边 \(e\) 的上方的端点，那么当sweep line移动到 \(p\) 时，最终的 \(C\) （\(C \in C_1 \bigcap C_2\)）在 \(p\) 处所在的边可能有三种情况：

  • 这条边 \(e\) 本身，并且 \(p\) 是这条边上方的点

  • \(e \bigcap left\_edge\_C2\) （即\(e\) 和 \(C_2\) 左侧与sweep line相交的边的交集所形成的边），并且 \(p\) 是这条边上方的点。

  • \(e \bigcap right\_edge\_C2\) （即\(e\) 和 \(C_2\) 右侧与sweep line相交的边的交集所形成的边），并且 \(p\) 是这条边上方的点。

判断情况之后，再执行以下处理：

  • 首先检查，如果 \(p\) 位于 \(left\_edge\_C2\) 和 \(right\_edge\_C2\) 之间，那么以 \(p\) 点为upper point的这条edge \(e\) 就是最终 \(C\) 中的一条边，把包含这条 \(e\) 的boundary line所在的半平面（表达式）添加到 \(\mathcal{L}_{left}(C)\) 中。

  • 其次检查，这条边 \(e\) 是否和 \(right\_edge\_C2\) 相交。如果相交，那么交点就是最终 \(C\) 中的顶点。并且有

    • 如果 \(p\) 位于 \(right\_edge\_C2\) 右侧（沿着\(right\_edge\_C2\) 走，\(C_2\) 在其左侧，以此为参考的 \(p\) 位于其 \(right\_edge\_C2\) 的右侧），那么 \(e\) 和 \(right\_edge\_C2\) 从交点开始的各自边的部分，就是 \(C\) 中的两条edge（或edge的一部分）。 这种情况下，就把 \(e\) 所在的半平面（表达式）添加到\(\mathcal{L}_{left}(C)\)中，把\(right\_edge\_C2\) 所在的半平面（表达式）添加到\(\mathcal{L}_{right}(C)\)中。

    • 如果 \(p\) 位于 \(right\_edge\_C2\) 左侧（沿着\(right\_edge\_C2\) 走，\(C_2\) 在其左侧，以此为参考的 \(p\) 位于其 \(right\_edge\_C2\) 的左侧），那么 \(e\) 和 \(right\_edge\_C2\) 各自的边到交点结束的部分，就是 \(C\) 中的两条edge（或edge的一部分）。 这种情况下，不需要把 \(e\) 或 \(right\_edge\_C2\)所在的半平面添加到\(\mathcal{L}_{left}(C)\)或\(\mathcal{L}_{right}(C)\)中，因为它们会在其他时机被添加进去。

  • 最后再检查，这条边 \(e\) 是否和 \(left\_edge\_C2\) 相交。如果相交，那么交点就也是最终 \(C\) 中的顶点。并且有

    • 如果 \(p\) 位于 \(left\_edge\_C2\) 右侧（沿着\(left\_edge\_C2\) 走，\(C_2\) 在其左侧，以此为参考的 \(p\) 位于其 \(left\_edge\_C2\) 的右侧），那么以交点开始的边 \(e\) 的部分就是 \(C\) 中的边，并把对应的half-plane添加到 \(\mathcal{L}_{left}(C)\) 中。

    • 如果 \(p\) 位于 \(left\_edge\_C2\) 左侧（沿着\(left\_edge\_C2\) 走，\(C_2\) 在其左侧，以此为参考的 \(p\) 位于其 \(left\_edge\_C2\) 的左侧），那么以交点开始的边 \(left\_edge\_C2\) 的部分就是 \(C\) 中的边，并把对应的half-plane添加到 \(\mathcal{L}_{left}(C)\) 中。

==（上面三个处理中对应的3个图在77页，页码是70）==

需要注意的是，我们可能向 \(\mathcal{L}_{left}(C)\) 中添加 \(e\) 所在的half-plan，或 \(left\_edge\_C2\) 所在的half-plane。如果以 \(left\_edge\_C2\) 和 \(e\) 的交点开始的 \(left\_edge\_C2\) 的部分edge是最终 \(C\) 中的edge的话， \(left\_edge\_C2\) 才能加到 \(\mathcal{L}_{left}(C)\) 中。而准备添加 \(e\) 时，要确保最终组成 \(C\) 中的edge，要么是 \(e\) 从其upper endpoint开始的，要么是 \(e\) 从它和 \(right\_edge\_C2\) 的交点开始的部分，而且这两种情况下都要先添加 \(e\) 。

从上可知，处理每条edge花费的是常数时间，那么计算两个convex polygon相交的时间复杂度就是 \(O(n)\)。

为了证明算法是正确的，需要证明我们能找到最终 \(C\) 中的所有的边（edge），并且顺序是正确的。 考虑 \(C\) 中的一条edge \(e_C\)，记它的upper endpoint为 \(p\)。

  • 要么 \(p\) 是 \(C_1\) 或 \(C_2\) 中某条edge的upper endpoint。 这种情况下，当sweep line到达 \(p\) 点的时候，我们就能找到 \(C\) 中的这条edge \(e_C\)。

  • 要么它是 \(C_1\) 中的一条edge \(e\) 和 \(C_2\) 中的一条edge \(e'\) 的交点。 这种情况下，当sweep line到达 \(e\) 和 \(e'\) 的两个upper endpoint的较低（lower）的那个时，我们同样也能找到 \(C\) 中的这条edge \(e_C\)。

因此，我们就能找到最终 \(C\) 中的所有的边。而证明其顺序正确性不难（？？）。

定理4.3 计算平面上两个凸多边形区域的交集的时间复杂度是\(O(n)\)。

Theorem 4.3 The intersection of two convex polygonal regions in the plane can be computed in \(O(n)\) time.

根据这个结论，回到之前的算法 \(IntersectHalfPlanes(H)\)，之前我们分析的时间复杂度是

\[\begin{split} T(n) = \left\{\begin{matrix} O(1), if \space\space n = 1 \\ O(nlogn) + 2T(n/2), if \space\space n > 1 \end{matrix}\right. \end{split}\]

其中 \(O(nlogn)\) 是计算两个任意多边形区域交集的时间复杂度。现在我们证明了计算两个凸多边形区域的交集的时间复杂度是\(O(n)\)，那么就有了

\[\begin{split} T(n) = \left\{\begin{matrix} O(1), if \space\space n = 1 \\ O(n) + 2T(n/2), if \space\space n > 1 \end{matrix}\right. \end{split}\]

这个式子最终的结果是 \(T(nlogn)\)。 所以有下面的推论：

推论4.4 计算平面上 \(n\) 个半平面交集的时间复杂度是 \(O(nlogn)\)，并且空间复杂度是线性的。

Corollary 4.4 The common intersection of a set of n half-planes in the plane can be computed in \(O(nlogn)\) time and linear storage.

计算半平面交集的问题和计算凸包（convex hull）的问题紧密相关，并且有一种算法和第一章中的\(ConvexHull\)算法几乎一致。这两种计算问题的关系在第8.2节和第11.4节有更详细的讨论。

4.3 Incremental Linear Programming

上一节讲述的是，计算 \(n\) 个半平面相交的所有解的算法，即计算 \(n\) 个线性约束所有的解。类似排序算法，这种算法的渐近下界是 \(\Omega(nlogn)\)，即最坏的情况下的时间复杂度。而对于我们的铸件移除问题，不需要得出所有的解，只需要其中的一个解即可，这就使得我们可以找到更快的算法。

求解一组线性约束的一个解的问题，和运筹学（Operation Research，OR）领域的一类问题紧密相关，它就是线性规划（或叫线性优化）。英文是linear programming（linear optimization），这里的programming指的是给计算机发送指令（giving instructions to a computer）。

This term was coined before “programming” came to mean “giving instructions to a computer”.

线性规划问题，是求解一组线性约束的一个解，并且这个解使得一个线性目标函数的值最大化。下面是更精确的描述。 \(\mathbb{R}^d\) （\(d\) 维实数空间）上，有一个要使得其最大化的目标函数， Maximize：

\[ c_1x_1 + c_2x_2 + \cdots + c_dx_d \]

Subject to：

\[\begin{split} \left\{\begin{matrix} a_{1,1}x_1 + a_{1,2}x_2 + \cdots + a_{1,d}x_d \le b_1 \\ a_{2,1}x_1 + a_{2,2}x_2 + \cdots + a_{2,d}x_d \le b_2 \\ \vdots \\ a_{n,1}x_1 + a_{n,2}x_2 + \cdots + a_{n,d}x_d \le b_n \end{matrix}\right. \end{split}\]

\(a_{i,j}\)，\(b_i\) 和 \(c_i\) 都是实数，它们是问题的输入参数。 需要最大化的函数叫做目标函数（objective function）。 目标函数和这一组线性约束一起，叫做线性程序（a linear program）。 变量的个数 \(d\)，是这个线性程序的维度（dimension）。 通过上节可知，一组线性约束可看做是\(\mathbb{R}^d\)空间上的 \(n\) 个半平面，而一个线性程序的解，实际上就是这 \(n\) 个半平面的交集区域（区域中的每个点都是线性程序的一个解）。我们把这个交集区域叫作 feasible region，这个区域里面的点（即线性程序的一个解）叫作feasible，而这个区域外面的点，叫做infeasible。 这个交集区域可能是无界限的，也可能是空的。当这个交集区域是空时，我们把这个线性程序也叫做infeasible。 目标函数（objective function）实际上可以看做是\(\mathbb{R}^d\)空间上的一个向量，所以它最大化，可以转换为求一个点 \(p\) = \((x_1, x_2, \dots, x_d)\) 在 \(\mathbb{R}^d\) 空间（\(d\) 维实数空间）上沿方向 \(\vec{c} = (c_1, c_2, \dots, c_d)\) 的极值。 因此，一个线性程序的解，就是求解在feasible region上的一个点，使得其在 \(\vec{c} = (c_1, c_2, \dots, c_d)\) 上取得极值。 我们用符号\(f_{\vec{c}}\) 来表示由方向 \(\vec{c}\) 定义的目标函数（ \(c_1x_1 + c_2x_2 + \cdots + c_dx_d\) ）。

let \(f_{\vec{c}}\) denote the objective function defined by \(\vec{c}\) .

==（上面提到的线性程序的解的示意图在79页，页码是72，位于最上面的第一个图）==

运筹学（Operation Research，OR）中的许多问题可以用线性程序来描述，并由此产生了许多算法，单纯形算法（simplex algorithm）就是实际中最常用最有效的算法之一。

运筹学中解决的一般都是有多个变量的多线性约束问题，即变量个数和约束条件很多。而本章讨论的是两个变量在 \(n\) 个线性约束条件下目标函数的求解，运筹学中传统的线性程序方法在这种低维度的线性程序求解中并不是非常有效率，因此在计算几何中就有其他更有效的算法。

我们以 \(H\) 表示2维平面上一组线性约束， \(\vec{c} = (c_x, c_y)\) 表示目标函数的方向，因此目标函数就是 \(f_{\vec{c}} (p) = c_xp_x + c_yp_y\) 。我们的目的就是，找到 \(\mathbb{R}^2\) 上的一个点（\(p \in \mathbb{R}^2\)），使得 \(p \in \bigcap H\) （所有半平面的交集区域），并且 \(f_{\vec{c}} (p)\) 取得最大值。我们用 \((H, \vec{c})\) 表示这个线性程序，并用 \(C\) 表示feasible region。

从几何上来看，线性程序 \((H, \vec{c})\) 的解有以下四种情况

  • 情况1，无解。即所有半平面的交集是空。这种情况下，线性程序是infeasible（不可行），线性约束就没有解。

  • 情况2，有解。交集区域是无限的。此时，feasible region（可行域）沿着 \(\vec{c}\) 方向是无限的，而且有一条射线 \(\rho\) 是完全被包含在feasible region \(C\) 当中，我们要求的解，就是这样一条射线 \(\rho\)。

  • 情况3，有解。交集区域是有限的。 此时，feasible region（可行域）的一条edge上的向外的法向的点，位于 \(\vec{c}\) 方向上（实际的意思是，这条edge \(e\) 是垂直于 \(\vec{c}\) ）。这时候，线性程序有解，但解不唯一，这条edge \(e\) 上任意一个点就是feasible point，它都能使得目标函数\(f_{\vec{c}} (p)\) 取得最大值。

  • 情况4，有解。交集区域退化为一个点（属于某即个半平面边界直线共同的交点）。如果不属于之前三种情况中的一种，那么就有一个唯一解，这个解就是 feasible region \(C\) 中的点 \(v\)，并且使得目标函数在 \(\vec{c}\) 方向上取得极值。

==（上面提到的线性程序的解的示意图在79页，页码是72，位于中间的四个一排的图）==

我们的2维线性程序的算法是增量式的，它逐一增加约束条件，并且维护每个中间线性程序的解。然而，我们的算法要求每个中间线性程序的解是唯一且明确定义的。换句话说，它要求每个中间线性程序的解，就是像前面提到的情况4一样，feasible region上有一个唯一的optimal vertex。 为了满足这个要求，给线性程序添加两个额外的约束条件，以便保证交集的区域是有限的。比如，当 \(c_x > 0\) 并且 \(c_y > 0\) 时，取一足够大的值 \(M (M \in \mathbb{R})\)，添加两个约束条件 \(p_x \le M\) 和 \(p_y \le M\)。这里 \(M\) 的选取应该足够大，并且使得添加的两个约束条件不影响最终的最优解。

对于这样额外的界定约束条件，在实际的线性程序应用场景中，都有自然的现实限制。在我们讨论的铸造问题当中，机械原理的限制，使得我们不能把多面体沿着几乎水平的方向移出模具（比如，不能以和\(xy\)平面小于\(1^{\circ}\)的方向移出）。这样，\(p_x\)和\(p_y\)的绝对值就显然有上限（即存在一个足够大的值 \(M (M \in \mathbb{R})\)）， 满足\(p_x \le M\) 和 \(p_y \le M\)。

在4.5节会讨论unbounded linear program，以及如何不使用手动添加约束来解决bounded linear program。

为精确期间，把这两个额外的新约束做如下定义。

\[\begin{split} \begin{split} m_1 := \left\{ \begin{matrix} {p_x \le M, \ if\ c_x \gt 0} \\ {-p_x \le M, \ otherwise} \\ \end{matrix} \right.\end{split} \end{split}\]
\[\begin{split} \begin{split} m_2 := \left\{ \begin{matrix} {p_y \le M, \ if\ c_y \gt 0} \\ {-p_y \le M, \ otherwise} \\ \end{matrix} \right.\end{split} \end{split}\]

这里 \(m_1, m_2\) 的选取只和目标函数中的参数 \(\vec{c}\) 有关，和线性约束 \(H\)（即所有的半平面）无关。而 \(C_0 = m_1 \bigcap m_2\) 这个feasible region是一个正交的楔形。 另外如果是前面提到的情况3，有多个不唯一的解，约定俗成地，我们按照字典序，取最小的那个点（解）。概念上，这相当于把 \(\vec{c}\) 轻微旋转了一个小角度，这样它就和任何一个半平面对应的直线就不垂直了。 需要注意的是，哪怕一个有界的线性程序，按照字典序可能也找不到最小的点（解）。我们需要选择约束条件 \(m_1\) 和 \(m_2\) ，使得这样的情况不会发生。 在选择了这样两个额外的约束条件之后，任意有解的线性程序就只有唯一的解，就是feasible region上的一个vertex，我们把这样的vertex叫做 optimal vertex。

记 \((H, \vec{c})\) 是一个线性程序，把 \(n\) 个半平面依次记为 \(h_1, h_2, \dots, h_n\)。 记 \(H_i\) 是前 \(i\) 个线性约束和两个额外约束\(m_1, m_2\)的组合，记 \(C_i\) 是这样一组线性约束的feasible region（可行域）。

\[ H_i := \{m_1, m_2, h_1, h_2, \dots, h_n\} \]
\[ C_i := {m_1 \cap m_2 \cap h_1 \cap h_2 \cap \dots \cap h_n} \]

加上 \(C_0\)，我们的每个feasible region \(C_i\) （\(i \le n\)）都有一个唯一的optimal vertex，记作 \(v_i\)。显然有

\[ C_0 \supseteq C_1 \supseteq C_2 \supseteq \cdots \supseteq C_n = C \]

这也就是说，对于任何一个 \(i\)（\(0 \le i \le n\)），如果 \(C_i = \emptyset\)，那么对于任意一个 \(j\)（\(i \le j \le n\)）就有 \(C_j = \emptyset\)，线性程序就是不可行的（infeasible），我们的算法也就能提前退出。

下面的引理4.5，说明了每当我们增加一个线性约束（半平面）时，optimal vertex如何改变，而这是我们算法的基础。

引理4.5 对\(C_i\) 和 \(v_i\) （\(1 \le i \le n\)），有：

  1. 如果 \(v_{i-1} \in h_i\)，那么 \(v_i = v_{i-1}\)。

  2. 如果 \(v_{i-1} \notin h_i\)，那么要么 \(C_i = \emptyset\)，要么 \(v_i \in \ell_i\)（\(\ell_i\) 是半平面 \(h_i\) 的边界直线）

证明：

令 \(v_{i-1} \in h_i\)。（注意这里的\(h_i\)指的是半平面区域，不是半平面的分界线）。因为 \(C_i = C_{i-1} \cap h_i\)，并且 \(v_{i-1} \in C_{i-1}\)，这就意味着 \(v_{i-1} \in C_i\)（这是显然的）。进一步，因为 \(C_i \subseteq C_{i-1}\)（即 \(C_i\) 属于 \(C_{i-1}\)），所以 \(C_i\) 的optimal vertex不可能比 \(C_{i-1}\) 更好，所以，\(v_{i-1}\) 同样也是 \(C_i\) 的optimal vertex。

令 \(v_{i-1} \notin h_i\)。采用反证法，假设 \(C_i\) 非空，并且 \(v_i\) 不在 \(\ell_i\) 上。考虑线段 \(\overline{v_{i-1}v_i}\) ， \(v_{i-1} \in C_{i-1}\)（因为 \(v_{i-1}\)是\(C_{i-1}\)的optimal vertex），并且因为\(C_i \subset C_{i-1}\)，所以也有 \(v_i \in C_{i-1}\)。又因为 \(C_{i-1}\)是凸的，所以线段 \(\overline{v_{i-1}v_i}\) 完全被包含在 \(C_{i-1}\)中。因为 \(v_{i-1}\)是\(C_{i-1}\)的optimal vertex，而且目标函数 \(f_{\vec{c}}\)是线性的，所以沿着线段 \(\overline{v_{i-1}v_i}\) ，点 \(p\) 从 \(v_i\) 移动到 \(v_{i-1}\)时，\(f_{\vec{c}}(p)\)是单调递增的。因为 \(v_{i-1} \notin h_i\)并且\(C_i \subset C_{i-1}\)，所以线段 \(\overline{v_{i-1}v_i}\) 和 \(\ell_i\)必然有交点，记作 \(q\)。而又因为线段 \(\overline{v_{i-1}v_i}\) 完全被包含在 \(C_{i-1}\)中，所以交点 \(q\) 也必然在 \(C_i\)中（因为\(\ell_i\)是\(C_i\)的边界）。因为 \(q\) 在线段 \(\overline{v_{i-1}v_i}\) 上，所以根据前面的结论（沿着线段 \(\overline{v_{i-1}v_i}\) ，点 \(p\) 从 \(v_i\) 移动到 \(v_{i-1}\)时，\(f_{\vec{c}}(p)\)是单调递增的），就有 \(f_{\vec{c}}(q) \gt f_{\vec{c}}(v_i)\)。这就是说在\(C_i\)中存在另外一个不同于optimal vertex \(v_i\) 的点 \(q\)，它使得\(f_{\vec{c}}\)的取值比 \(v_i\)还要大，但这是和\(v_i\)的定义是矛盾的，因此假设不成立。故的证。

==（用来说明线段 \(\overline{v_{i-1}v_i}\) 的示意图在81页，页码是74，位于上面的一个图）==

书中接着以两个图，说明了添加了一个半平面（\(h_k\)）之后，optimal vertex的变化情况。 ==（两个示意图在81页，页码是74，位于下面的两个图）==

通俗地讲，引理4.5的意思是，如果上一个线性程序 \((H_{i-1}, \vec{c})\) 的可行域里面的optimal vertex位于当前半平面 \(h_i\) 内部，那么当前线性程序 \((H_i, \vec{c})\) 的可行域里面的optimal vertex，就和上一个线性程序可行域里面的optimal vertex相同；如果上一个线性程序 \((H_{i-1}, \vec{c})\) 的可行域里面的optimal vertex不在当前半平面 \(h_i\) 内部，要么当前线性程序不可行（infeasible，即 \(C_i = \emptyset\)），要么当前线性程序 \((H_i, \vec{c})\) 的可行域里面的optimal vertex，就在当前半平面 \(h_i\) 的分界线 \(\ell_i\)上。

虽然引理4.5告诉了我们当添加一个新的半平面到当前的约束组时，optimal vertex的变化情况，但它没有告诉我们如何找到这个optimal vertex。

如果是引理4.5中说明的第一种情况，那么下一个要求解的optimal vertex \(v_i\) 就是当前的optimal vertex \(v_{i-1}\)。

如果是引理4.5中说明的第二种情况，那么问题就转化为：找到\(\ell_i\)（\(\ell_i\) 是半平面 \(h_i\) 的边界直线）上的一个点\(p\)，使得目标函数\(f_{\vec{c}}(p)\)最大化，并且点\(p \in h\)（\(h \in H_{i-1}\)），即点\(p\)是前\(i-1\)个半平面交集（区域）上的点。

为了简化术语，我们假定\(\ell_i\)（\(\ell_i\) 是半平面 \(h_i\) 的边界直线）不是垂直的，这样我们就能仅用\(x\)的坐标来参数化它。定义一个函数 \(\overline{f_{\vec{c}}}:\mathbb{R} \mapsto \mathbb{R}\)，使得对任意一点 \(p \in \ell_i\)，有 \(f_{\vec{c}}(p) = \overline{f_{\vec{c}}(p_x)}\) （还没想明白如何理解？如何找到这样的一个映射函数？好像要解微分方程？）。对一个半平面\(h(h \in H_{i-1})\)，记 \(\sigma(h,\ell_i)\) 为 \(\ell_i\) 和 \(h\) 的边界线的交点的\(x\)坐标。如果没有交点，要么是\(\ell_i\)上的每个点都满足条件约束\(h\)（也即半平面 \(h_i\)的分界线 \(\ell_i\) 完全包含于半平面\(h\)之内，这里\(h \in H_{i-1}\)），要么\(\ell_i\)上没有点满足条件约束\(h\)（也即半平面 \(h_i\)的分界线 \(\ell_i\) 完全位于半平面\(h\)之外，这里\(h \in H_{i-1}\)）。对于前者，我们忽略这个条件约束，对于后者，我们就报告这个线性程序不可行（infeasible）。

根据 \(\ell_i \cap h\) 的交集区域是被限定在左侧还是右侧，可以得到对解的\(x\)坐标的约束形式，即\(x \geqslant \sigma(h, \ell_i)\)或者\(x \leqslant \sigma(h, \ell_i)\)。因此需要求解的问题可以重新表述如下， 为最大化的目标函数为

\[ \overline{f_{\vec{c}}(x)} \]

它服从于

\[ x \geqslant \sigma(h, \ell_i), \space h \in H_{i-1} \space and \space \ell_i \cap h \space is \space bounded \space to \space the \space left \]

以及

\[ x \leqslant \sigma(h, \ell_i), \space h \in H_{i-1} \space and \space \ell_i \cap h \space is \space bounded \space to \space the \space right \]

上面两个约束条件的意思是，在\(H_{i-1}\)（前\(i-1\)个半平面）中，（1）有一部分半平面和第\(i\)个半平面的分界线相交之后，\(\ell_i\)与这些半平面相交的部分位于这些半平面的右侧，所以它是有左侧下限（bounded to the left）的；（2）而剩下的另一部分半平面和第\(i\)个半平面的分界线相交之后，\(\ell_i\)与这些半平面相交的部分位于这些半平面的左侧，所以它是有右侧上限（bounded to the right）的。那么，这个1维线性程序的可行域就是这两者的交集。

这个一维线性程序的解是容易的，令

\[ x_{left} = \max_{h \in H_{i-1} } \space \{ \sigma(h,\ell_i) : \ell_i \cap h \space is \space bounded \space to \space \space the \space left \} \]
\[ x_{right} = \min_{h \in H_{i-1} } \space \{ \sigma(h,\ell_i) : \ell_i \cap h \space is \space bounded \space to \space \space the \space left \} \]

这里 \(x_{left}\)的意思是，\(H_{i-1}\)中的有些半平面和 \(\ell_i\) 相交之后的部分位于这些半平面的右侧，那么它们就在 \(\ell_i\) 上有左侧下限的\(x\)坐标，取这些\(x\)坐标值中的最大值，就是 \(x_{left}\)。而 \(x_{right}\)的意思是，\(H_{i-1}\)中的有些半平面和 \(\ell_i\) 相交之后的部分位于这些半平面的左侧，那么它们就在 \(\ell_i\) 上有右侧上限的\(x\)坐标，取这些\(x\)坐标值中的最小值，就是 \(x_{right}\)。

这样，区间\([x_{left}, x_{right}]\)就是这个一维线性程序的解（feasible region）。因此如果 \(x_{left} > x_{right}\)，那就说明这个线性程序无解（不可行，infeasible），否则根据目标函数，当\(\ell_i\)上的点\(p\) 在\(x_{left}\)或\(x_{right}\)处取得极值。

需要注意的是，因为我们添加了额外的约束\(m_1\)和\(m_2\)，所以一维线性程序（的解）不会是无界的。

引理4.6 可以在线性时间里解出一维线性程序。因此，如果是引理4.5中的第二种情况时，我们能以\(O(i)\)时间复杂度计算出来新的optimal vertex \(v_i\)，或者判定这个线性程序无解（不可行，infeasible）。

下面是更详细的线性规划算法描述。

算法：\(2DBoundedLP(H, \vec{c}, m_1, m_2)\)

输入：线性程序 \((H \cup {m_1, m_2}, \vec{c})\) ，\(H\)是\(n\)个半平面，向量 \(\vec{c} \in \mathbb{R}^2\)（即二维实数域），\(m_1\)和\(m_2\)是解的额外限定约束条件。

输出：如果线性程序 \((H \cup {m_1, m_2}, \vec{c})\) 无解（infeasible），就声明无解并退出。否则，就点\(p\)，它是按字典序找到的、能使得目标函数\(f_{\vec{c}}(p)\)最大化的点。

步骤：

  • 令 \(v_0\) 是 \(C_0\) 角落上的点。

  • 令\(H\)是\(h_1, h_2, \dots, h_n\)这n个半平面

  • 从\(1\)到\(n\)，遍历\(i\)

    • 如果 \(h_{i-1}\)这个半平面的 optimal vertex \(v_{i-1}\)也在半平面 \(h_i\)上，那么\(h_i\)这个半平面的 optimal vertex \(v_i\)也同样是\(v_{i-1}\)。

    • 如果 \(h_{i-1}\) 这个半平面的 optimal vertex \(v_{i-1}\)也不在半平面 \(h_i\)上，那么找到\(\ell_i\)（半平面\(h_i\)的边界线）上的一个点\(p\)，它满足前\(i-1\)个线性条件约束（即\(H_{i-1}\)），并且它能够使得目标函数\(f_{\vec{c}}(p)\)最大化，此时，\(p\)就是要找的\(v_i\)。如果找不到这样的点\(p\)，就停止循环，报告这个线性程序无解，然后退出。

  • 循环到最后，报告\(v_n\)，这就是这个线性程序的解。

引理4.7 对一个有\(n\)个约束条件和\(2\)个变量的有界的线性程序，算法\(2DBoundedLP(H, \vec{c}, m_1, m_2)\)可以用\(O(n^2)\)的时间复杂度和线性的空间复杂度求得解。 证明： 根据前面的引理4.5及其结论，证明了可以找到这样的\(v_i\)是\(C_i\)的optimal vertex，因此可以得出引理4.7的正确性。而且如果某个\(C_i\)是空集，那么因为 \(C = C_n \subseteq C_i\) ，所以最终的\(C\)也就是空集。 易得空间复杂度是线性的。 同样根据引理4.6，对于每个\(i\)，计算一维线性程序的时间复杂度是\(O(i)\)，所以最终的时间复杂度是

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

虽然这样的算法是简单并且很好的，但还是太慢了。对于每次添加第\(i\)个约束时计算optimal vertex花费的诗句是\(O(i)\)，它不是紧确界（tight bound），而是渐近上限（即小于等于，\(\leqslant\)）。而根据引理4.5，只有当\(v_{i-1} \notin h_i\)时，花费的时间是\(\Theta(i)\)（渐近紧确界，即相等，\(=\)），而当\(v_{i-1} \in h_i\)时，只需要花费常数时间。但我们不能保证每次都有\(v_{i-1} \in h_i\)，书中此处举了一个例子，说明了每次第\(i\)个约束条件就会导致前一个的optimal vertex \(v_{i-1}\)不再满足这第\(i\)个约束条件，而且花费的时间就是\(\Theta(n^2)\)。接下来的4.3节说明了如何进行一定意义上的加速处理。 ==（举例的示意图在83页，页码是76）==

4.4 Randomized Linear Programming

在前面一节的结尾，作者举例说明了每次遍历到下一个半平面时，总是会导致optimal vertex改变，这样运行时间就会是最差的情况。但如果把遍历到的半平面的顺序做下改变，比如，从原先的\(H= \{h_1, h_2, \dots, h_n\}\)变为\(H = \{h_n, h_{n-1}, \dots, h_1\}\)，这样除了第一次添加半平面时需要计算一次optimal vertex，之后每次添加半平面，optimal vertex就不再发生变化，因此运行时间就会变成\(O(n)\)。 但很可惜的是，这样的顺序也许存在，但是我们不太好找出来，因为我们要在整个算法开始之前就要确定好添加半平面的顺序。 这就是一个引人入胜的现象了。我们不能保证\(H\)的顺序可以产生良好的运行时间，但是我们可以取\(H\)的一个随机顺序。虽然这样在运气不好的情况下仍然产生二次方的运行时间，但如果运气好，我们就能达到更快的运行时间。事实上，我们下面会证明在多数情况下，我们能取得较快的运行时间。 为了完整描述起见，把带有\(H\)随机顺序的算法罗列如下：

算法：\(2DRandomizedBoundedLP(H, \vec{c}, m_1, m_2)\) 输入：线性程序 \((H \cup {m_1, m_2}, \vec{c})\) ，\(H\)是\(n\)个半平面，向量 \(\vec{c} \in \mathbb{R}^2\)（即二维实数域），\(m_1\)和\(m_2\)是解的额外限定约束条件。 输出：如果线性程序 \((H \cup {m_1, m_2}, \vec{c})\) 无解（infeasible），就声明无解并退出。否则，就点\(p\)，它是按字典序找到的、能使得目标函数\(f_{\vec{c}}(p)\)最大化的点。 步骤：

  • 令 \(v_0\) 是 \(C_0\) 角落上的点。

  • 通过调用随机排序算法\(RandomPermutation(H[1\cdots n])\)，得到\(n\)个半平面\(H\)的一个随机排序序列\(h_1, h_2, \dots, h_n\)。

  • 从\(1\)到\(n\)，遍历\(i\)

    • 如果 \(h_{i-1}\)这个半平面的 optimal vertex \(v_{i-1}\)也在半平面 \(h_i\)上，那么\(h_i\)这个半平面的 optimal vertex \(v_i\)也同样是\(v_{i-1}\)。

    • 如果 \(h_{i-1}\) 这个半平面的 optimal vertex \(v_{i-1}\)也不在半平面 \(h_i\)上，那么找到\(\ell_i\)（半平面\(h_i\)的边界线）上的一个点\(p\)，它满足前\(i-1\)个线性条件约束（即\(H_{i-1}\)），并且它能够使得目标函数\(f_{\vec{c}}(p)\)最大化，此时，\(p\)就是要找的\(v_i\)。如果找不到这样的点\(p\)，就停止循环，报告这个线性程序无解，然后退出。

  • 循环到最后，报告\(v_n\)，这就是这个线性程序的解。

算法\(2DRandomizedBoundedLP\)，和算法\(2DBoundedLP\)，唯一的区别是，在遍历\(n\)个半平面之前，计算\(n\)个半平面的一个排序序列。 我们假设有一个随机数发生器\(Random(k)\)，它根据输入\(k\)产生一个范围在\([1,k]\)之间的一个随机数，并且花费的时间复杂度是常数时间。 随机排序的算法\(RandomPermutation\)如下。

算法：\(RandomPermutation(A)\) 输入：一个数组\(A[1 \cdots n]\) 输出：还是有同样元素的数组\(A[1 \cdots n]\)，但是是以一个随机顺序重新排序的 步骤：

  • 从\(n\)到\(2\)，遍历\(k\)

  • 对于每个当前的\(k\)，调用随机数发生器\(Random(k)\)，得到一个索引值\(rndindex\)

  • 交换元素\(A[k]\)和\(A[rndindex]\)。

这个新的线性程序的算法就是随机算法（randomized algorithm）。它的运行时间取决于算法中每次得到的随机序列。

那么这样的随机算法的运行时间（或者时间复杂度）是怎么样的呢？ 因为\(n\)个半平面的排列有\(n!\)种，所以运行时间就有\(n!\)种。而每种随机排列所产生的随机算法的运行时间是相互独立的（独立同分布），所以我们就需要分析它们的期望运行时间（expected running time），也就是这\(n!\)种运行时间的数学期望。而定理4.8说明了这样的期望运行时间是\(O(n)\)，而且我们对\(n\)个半平面的输入没有做任何的假设和限定。

引理4.8 二维平面空间上，有\(n\)个约束条件的线性规划问题，可以以\(O(n)\)的时间复杂度和最坏情况下线性的空间复杂的解决。

Lemma 4.8 The 2-dimensional linear programming problem with n constraints can be solved in \(O(n)\) randomized expected time using worst-case linear storage.

证明：我们之前已经观察到，需要的空间复杂度是线性的。 因为\(RANDOMPERMUTATION\)算法的时间复杂度是\(O(n)\)，所以剩下需要分析的就是添加半平面 \(h_1, h_2, \dots, h_n\)的运行时间。当optimal vertex不发生改变的时候，添加一个半平面需要的是常数时间（constant time）。但当optimal vertex发生改变的时候，我们需要解决一个1维的线性规划问题。现在来计算所有这些1维线性规划问题的时间界限。 假定 \(X_i\) 是一个随机变量，当 \(v_i \notin h_i\) 时取值为\(1\)，否则当 \(v_i \in h_i\) 时取值为 \(0\)，即

\[\begin{split} \begin{split} X_i = \left\{ \begin{matrix} {1 \space , v_i \notin h_i} \\ {0 \space , v_i \in h_i} \\ \end{matrix} \right.\end{split} \end{split}\]

回忆起对于有\(i\)个线性约束条件的1维线性规划问题，解决它的时间复杂度是\(O(i)\)。据此，有\(n\)个半平面（线性约束条件）的线性规划问题所需的时间复杂度就是，

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

使用数学期望的线性特性（linearity of expectation），随机变量之和的数学期望，等于随机变量数学期望的和（the expected value of a sum of random variables is the sum of the expected values of the random variables）。无论随机变量是否相关还是相互独立，这个线性特性都成立。所以，解决这个1维线性规划问题时间复杂度的数学期望就是，

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

\(E[X_i]\)是什么？它实际上就是\(v_{i-1} \notin h_i\)的概率。下面使用所谓的后向分析（backward analysis）技术来分析这个概率。 假设已经完成了算法步骤，计算出来了optimum vertex \(v_n\)。因为\(v_n\)是\(C_n\)上的一个vertex，所以它由至少两个半平面（half-plane）定义。此时我们后退一步，考察\(C_{n-1}\)。注意从\(C_n\)中移除\(h_n\)就得到了\(C_{n-1}\)。这种情况下，optimum vertex如何变化？（即\(C_n\)中对应的\(v_n\)，在从\(C_n\)中移除了\(h_n\)之后，对应的optimum vertex \(v_{n-1}\)是否发生变化）而\(v_{n-1}\)和\(v_n\)不同（即发生变化）的唯一情况是，\(v_n\)不是\(C_n\)中使得目标函数沿着\(\vec{c}\)方向上取得极值的顶点，而这只发生于此时移除的\(h_n\)是定义出\(v_n\)的（至少）两个半平面之一的时候。因为算法中添加半平面的顺序是随机的，所以\(h_n\)实际上是\(\{h_1, h_2, \dots, h_n\}\)中之一。因此，\(h_n\)是定义\(v_n\)的概率至多是\(2/n\)。这里说至多，原因是（一）\(v_n\)可能由多条半平面的边界直线的交点形成，因此这时候移除的\(h_n\)可能不会改变optimum vertex \(v_n\)（二）\(v_n\)还可能由\(m_1\)或\(m_2\)定义，它们不在前面提到的\(\{h_1, h_2, \dots, h_n\}\)之中。因为这两个原因，所以这个概率至多是\(2/n\)。

==（举例的示意图在85页，页码是78）==

为了计算\(E[X_i]\)的上界，我们固定前\(i\)个半平面构成的子集，这个子集构成了\(C_i\)。根据前面的分析，添加\(h_i\)时，使用后向分析。当添加\(h_i\)时需要重新计算一个新的optimal vertex的概率，就是我们把\(h_i\)从\(C_i\)中移除时optimal vertex发生改变的概率。而后者发生的概率，只有\(h_i\)是\(\{h_1, h_2, \dots, h_i\}\)中至多两个半平面之一的概率，而因为添加这些半平面的顺序是随机的，所以\(h_i\)是这两个特殊的半平面之一的概率就是\(2/i\)。我们推断出这个概率的条件，是基于这前\(i\)个半平面是\(H\)中的某个固定的子集的情况。但是因为这个推断出来的概率界限是对任何固定的子集适用的，所以我们就有 \(E[X_i] \leqslant 2/i\)。由此，我们得到了可以确定解决1维线性规划问题的总时间上界为

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

而且，我们也注意到这个算法其余部分的时间复杂度也是\(O(n)\)。 注意这里的数学期望只和算法的随机（顺序）选择有关。我们没有对可能的输入做平均。对任意一个有\(n\)个半平面的集合作为输入而言，算法运行时间的数学期望就是\(O(n)\)，并且没有“坏”的输入。

==（举例的示意图在86页，页码是79）==

4.5 Unbounded Linear Programs

在前面的一节里，我们通过手动添加两个额外的条件约束，避免处理无界的线性规划问题。但这种办法并不总是合适的。甚至有时候线性规划问题本身是有界的（bounded），但我们也不知道一个足够大的上界\(M\)。更进一步，在实际情况中，无界的线性规划问题是存在的，所以我们需要正确地解决它。

首先研究如何确定一个线性规划问题\((H, \vec{c})\)是无界的。正如我们前面提到的，一个无界的线性规划问题，意味着有一条射线 \(\rho\) 完全包含于可行区域（feasible region）\(C\)中，这样沿着射线 \(\rho\) 的方向，目标函数 \(f_{\vec{c}}\) 可以取得任意大的值。

如果我们以\(p\)表示射线的起点，以\(\vec{d}\)表示它的方向，那么就得到了射线 \(\rho\) 的参数化表达

\[ \rho = \{p + \lambda \cdot \vec{d} : \lambda \gt 0 \} \]

当且仅当 \(\vec{d} \cdot \vec{c} \gt 0\)时，目标函数 \(f_{\vec{c}}\) 取得任意大的值（为什么？下面是一个简单的证明）。另一方面，所有半平面（条件约束）\(H\)中的一个半平面\(h\)（\(h \in H\)），它的朝向有解区域（feasible region）的法向量\(\vec{\eta}(h)\)，有\(\vec{d} \cdot \vec{\eta} \geqslant 0\)。下面的推理4.9，说明了关于\(\vec{d}\)的两个条件，就足够判断一个线性规划问题是否有界。

这里首先简单地证明一下，为何当 \(\vec{d} \cdot \vec{c} \gt 0\)时，目标函数 \(f_{\vec{c}}\) （沿着射线 \(\rho\) 的方向）取得任意大的值。

首先 \(f_{\vec{c}}(x)\) 的定义如下，

\[ f_{\vec{c}}(x) = c_x \cdot p_x + c_y \cdot p_y \]

这里 \((p_x, p_y)\) 取射线 \(\rho\) 上的点。为了消除歧义，把前面射线参数的表达式改写为如下（即射线起点为 \(q\)），

\[\rho = \{q + \lambda \cdot \vec{d} : \lambda \gt 0 \} \]

这样射线上的任意一点的坐标可以写作 \(P = (p_x, p_y) = (q_x + \lambda d_x, q_y + \lambda d_y)\)。

那么 \(f_{\vec{c}}(x)\) 就可以写成

\[ f_{\vec{c}}(x) = c_x \cdot p_x + c_y \cdot p_y = c_x \cdot (q_x + \lambda d_x) + c_y \cdot (q_y + \lambda d_y) \]

化简得到

\[ f_{\vec{c}}(x) = c_x \cdot q_x + c_y \cdot q_y + \lambda \cdot (c_x d_x + c_y d_y) \]

而根据条件 \(\vec{d} \cdot \vec{c} \gt 0\)，以及 \(\lambda \gt 0\)，那么 \(f_{\vec{c}}(x)\) 可以进一步写成

\[ f_{\vec{c}}(x) > c_x \cdot q_x + c_y \cdot q_y \]

显然不等式右边是一个有限的数（因为\(c_x, c_y\)是固定的参数，\(q_x, q_y\)也是某个固定点的坐标），这就意味着此时 \(f_{\vec{c}}(x)\) 可以是一个无限大的值，而这正是通过沿着射线 \(\rho\) 上的点而取得的值。

推理4.9 一个线性规划问题 \((H, \vec{c})\)是无界的充要条件是，存在一个满足 \(\vec{d} \cdot \vec{c} \gt 0\) 的向量 \(\vec{d}\) ，使得\(H\)中的任意一个约束条件（半平面）\(h\) 朝向可行区域的法向量 \(\vec{\eta}(h)\) 有 \(\vec{d} \cdot \vec{\eta}(h) \geqslant 0\)，并且同时线性规划问题 \((H^{\prime}, \vec{c})\) 有解，这里 \(H^{\prime} = \{h \in H : \vec{\eta}(h) \cdot \vec{d} = 0\}\)。

Lemma 4.9 A linear program \((H, \vec{c})\) is unbounded if and only if there is a vector \(\vec{d}\) with \(\vec{d} \cdot \vec{c} \gt 0\) such that \(\vec{d} \cdot \vec{\eta}(h) \geqslant 0\) for all \(h \in H\) and the linear program \((H, \vec{c})\) is feasible, where \(H^{\prime} = \{h \in H : \vec{\eta}(h) \cdot \vec{d} = 0\}\).

注意，这里的\(\vec{\eta}(h)\)是半平面\(h\)的分界线的法向量，方向是朝向有可行域（feasible region）的一侧。所以，这里\(H^{\prime}\)实际上是半平面分界线和\(\vec{d}\)平行（共线）的那些半平面（约束条件）的集合，因为 \(\vec{d} \cdot \vec{\eta}(h) = 0\)，所以\(\vec{d}\)和半平面的法向量垂直，换言之就是和半平面的分界线平行（共线）。

证明：考虑线性规划问题 \((H, \vec{c})\) 和满足推理中设定条件的\(\vec{d}\)。因为线性规划问题 \((H^{\prime}, \vec{c})\) 有解，所以存在一个点 \(p_0 \in \textstyle \bigcap_{h \in H^{\prime}} h\) （即\(H^{\prime}\)中所有半平面的交集中的一个点）。考虑射线 \({\rho}_0 := \{p_0 + \lambda \cdot \vec{d} : \lambda \gt 0 \}\)。因为对于\(H^{\prime}\)中任意一个半平面\(h\)分界线的法向量 \(\vec{\eta}(h)\) 都有 \(\vec{d} \cdot \vec{\eta}(h) = 0\)（推理中的假定条件，即\(\vec{d}\)和半平面的分界线平行或共线），所以\(H^{\prime}\)中任意一个半平面 \(h\) 都完全包含射线 \({\rho}_0\)（因为\(\vec{d}\)和半平面的分界线平行或共线，又\(p_0\)在feasible region内，所以射线\(\rho_0\)就一定在每个半平面\(h\)中，这里\(h \in H^{\prime}\)）。进一步，因为有\(\vec{d} \cdot \vec{c} \gt 0\)，所以目标函数 \(f_{\vec{c}}\) 沿着 \({\rho}_0\) 的方向可以取得任意大的值。

对于任意一个属于\(H\)但不属于\(H^{\prime}\)的半平面\(h\)（即\(h \in H \backslash H^{\prime}\)），有 \(\vec{d} \cdot \vec{\eta}(h) > 0\)（因为根据推理给定的条件，\(\vec{d}\) 要么是和\(H\)中的任意一半平面分界线的法向量垂直，要么和它的夹角为锐角，即\(\vec{d} \cdot \vec{\eta} \geqslant 0\)，那么除去垂直的，剩下的自然都是夹角是锐角的，即\(\vec{d} \cdot \vec{\eta}(h) > 0\)）。

而这意味着，对于任意的一个\(h\)（\(h \in H \backslash H^{\prime}\)），存在一个\(\lambda_{h}\)，使得当所有的\(\lambda \geqslant \lambda_{h}\)时，射线\(p_0 + \lambda \cdot \vec{d}\) （的一部分）包含于半平面\(h\)（即有 \(p_0 + \lambda \cdot \vec{d} \in h\)）（为什么？还没想明白）。令 \(\lambda^{\prime} := \max_{h \in H \backslash H^{\prime}}{\lambda_{h}}\)（即\(\lambda^{\prime}\)是所有半平面\(h\)中对应的\(\lambda_h\)最大的一个），再令 \(p := p_0 + \lambda^{\prime} \cdot \vec{d}\)（\(\lambda^{\prime}\)就是刚刚的所有\(\lambda\)的最大值），此时，点 \(p\) 就一定包含于所有的半平面\(h\)（\(h \in H \backslash H^{\prime}\)）中，而这正是前面所得到的一条射线（的一部分）包含于所有半平面\(h\)的结论，即射线

\[ \rho = \{p + \lambda \cdot \vec{d}: \lambda \gt 0\} \]

是完全包含于每个半平面\(h\)中（\(h \in H\)），所以线性规划问题\((H, \vec{c})\)是无界的（unbounded）。

现在就可以像4.1节中叙述的一样，以解决一个1维线性规划问题，来测试2维线性规划问题 \((H, \vec{c})\) 是否无界。

首先旋转坐标系，使得 \(\vec{c}\) 的方向为垂直向上，即 \(\vec{c} = (0,1)\)。对任意方向的\(\vec{d} = (d_x, d_y)\)，当 \(\vec{d} \cdot \vec{c} \gt 0\) 时（夹角是锐角），它可以正规化为 \(\vec{d} = (d_x, 1)\)的形式（想象一下，把任意一个向量的起点移动到坐标原点，然后沿其方向延长，和直线 \(y=1\) 相交，就是点 \((d_x, 1)\)了，其实就是正规化），而这可以看做是 \(y=1\) 这条直线上坐标为 \(d_x\) 的点。对任一半平面（分界线）\(h\) 的法向量 \(\vec{\eta}(h) = ({\eta}_x, {\eta}_y)\)（这里的 \({\eta}_x, {\eta}_y\) 还是旋转之前坐标系里的坐标），下面的不等式

\[ \vec{d} \cdot \vec{\eta}(h) = d_x{\eta}_x + d_y{\eta}_y \geqslant 0 \]

可化为不等式 \(d_x{\eta}_x \geqslant -{\eta}_y\)（注意，这里的 \({\eta}_x, {\eta}_y\) 已经是旋转之后坐标系里的坐标了，也因此这里\(d_y=1\)了）。这样，我们就得到了一组 \(n\) 个线性不等式，也就是1维线性规划程序 \(\overline{H}\)（这里说线性规划程序这个术语不严谨，因为线性规划程序要包含约束和目标函数。不过这里我们只关注它是否可行，即feasibility，所以忽略目标函数，叙述起来比较方便）。

如果线性规划程序 \(\overline{H}\) 有一个解 \(d_{x}^{\ast}\)（即讨论这个固定不变的解\(d_{x}^{\ast}\)），我们把有严格解（solution is tight，即满足\(d_{x}^{\ast}{\eta}_x + {\eta}_y = 0\)）的那些半平面 \(h\) 的集合记作 \(H^{\prime}\)，当然有 \(H^{\prime} \subseteq H\)。我们仍然要验证 \(H^{\prime}\) 是有解的（因为前面说的是 \(\overline{H}\) 有解，不是 \(H^{\prime}\)有解）。虽然这时候我们面对的仍然是一个2为线性规划程序，但是它比较特殊，因为对任意的 \(h \in H^{\prime}\)，\(h\) 的法向量 \(\vec{\eta}(h)\) 是和 \(\vec{d} = (d_{x}^{\ast}, 1)\) 垂直的（正交），也就是说，此时 \(h\) 的分界线是和 \(\vec{d}\) 平行的。换句话说，就是 \(H^{\prime}\) 中的所有半平面的分界线是一系列平行分界线。这些平行的分界线和 \(x\) 轴相交，我们就会得到一个1维的线性规划程序 \(\overline{H^{\prime}}\)。（1）如果 \(\overline{H^{\prime}}\) 可行，那么就意味着 \(H^{\prime}\) 可行（因为 \(\overline{H^{\prime}}\) 是 \(H^{\prime}\) 在1维上的表示，所以\(\overline{H^{\prime}}\)可行就是\(H^{\prime}\)可行），那么（根据前面的推理4.9）原先的线性规划程序 \(H\) 就是无界的（unbounded），这时按照前面的推理，就可以以 \(O(n)\) 的时间复杂度，在可行域（feasible region）里画出一条射线 \(\rho\) （feasible ray）。（2）如果 \(\overline{H^{\prime}}\) 不可行，那么 \(H^{\prime}\) 不可行，因此 \(H\) 也不可行。

如果线性规划程序 \(\overline{H}\) 没有解，那么按照上面的推理，原先的线性规划程序 \((H, \vec{c})\) 就是有界的（bounded）。这种情况下能得到一些什么信息？回想1维线性规划程序的解决方案：\(\overline{H}\) 无解（不可行，即infeasible）的充要条件是，半直线（half-line）\(\overline{h_1}\) 的最大边界要比半直线（half-line）\(\overline{h_2}\) 的最小边界要更大（分析见下一段）。这时候两条半直线 \(\overline{h_1}\) 和 \(\overline{h_2}\) 没有交集。如果 \(h_1\) 和 \(h_2\) 分别是对应这两个约束的原始半平面（的分界线），那这就意味着 \((\{h_1, h_2\}, \vec{c})\) 是有界的。我们就能把 \(h_1\) 和 \(h_2\) 称作 certificates：它们证明了 \((H, \vec{c})\) 确实是有界的。

\(\overline{H}\) 无解（不可行，即infeasible）的充要条件是，半直线（half-line）\(\overline{h_1}\) 的最大边界要比半直线（half-line）\(\overline{h_2}\) 的最小边界要更大。

没搞懂什么意思？与之相关的内容在4.3节推理4.5下面，需要再研读思考。

（2023-02-04）想明白了。参考4.3解推理4.5下面的图，说的是1维的线性规划程序，就可以想象成位于 \(x\) 坐标轴上的区间。如果一个约束条件（对应这里的一个半平面 \(\overline{h_1}\)），得到的可行区域是半开区间 \(A = [a, +\infty)\)，而另一个约束条件（对应这里的另一个半平面 \(\overline{h_2}\)），得到的可行区域是半开区间 \(B = (-\infty, b]\)。那么 \(A\) 就是有左侧下限，而 \(B\) 就有右侧上限。又因为是一组约束条件（线性规划程序由一组多个约束条件组成），那么就可能有多个约束条件得到可行域是左闭右开的半区间，即有这样的多个\(A\)，那么这些左闭右开的半区间的左侧下限就有一个最大值。类似地，可能有多个约束条件得到可行域是左开右闭的半区间，即有这样的多个\(B\)，那么这些左开右闭的半区间的右侧上限就有一个最小值。这两个最值分别代表了一组约束条件中的两个条件，如果此时这个最大值比这个最小值还要大，这就意味着这两个约束条件的可行域（即1维坐标轴上的开区间）没有交集。

前面提到的 \(h_1, h_2\) 这样的两个certificates的用处是，在找到了它们之后，把它们当做算法2DRANDOMIZEDBOUNDEDLP中的两个约束\(m_1\)和\(m_2\)，而这意味着我们不需要再手动制造出约束来，使得解落于我们所允许的区间（区域）。

同时也要注意，有可能 \((\{h_1, h_2\}, \vec{c})\) 是有界的（bounded），但是没有按字典序的最小解。这种情况发生在，当1维线性规划程序由于有一个约束 \(h_1\) 的朝向其半平面内的法向量 \(\vec{\eta}(h_1)\)和目标函数的方向相反（即 \(\vec{\eta}(h_1) = -\vec{c} = (0, -1)\)），从而导致该1维线性规划程序不可行的时候。这种情况下，我们从剩下的约束条件（半平面）中找到一个满足法向量的\(x\)坐标满足 \(\eta_x(h_2) > 0\) 的半平面（约束条件）\(h_2\)。如果能找到这样的\(h_2\)，那么 \(h_1, h_2\) 就是certificates，能保证找到唯一的按字典序的最小解。如果找不到这样的\(h_2\)，要么线性规划程序不可行，要么就是找不到唯一的按字典序的最小解。要求解它，就要求解由一组满足\(\eta_x(h)=0\)的所有半平面\(h\)（即条件约束）组成的1维线性规划程序。而且如果它是可行的，我们就返回一条方向是\((-1,0)\)的射线\(\rho\)，这条射线\(\rho\)上的所有点都是可行的最优解（feasible optimal solution）。

现在可以给出2维线性规划程序的一般通用算法。

算法：\(2DRANDOMIZEDLP(H, \vec{c})\)

输入：一个线性规划程序\((H, \vec{c})\)，\(H\)是由\(n\)个半平面（条件约束）组成，并且\(\vec{c} \in \mathbb{R}^2\)（即目标函数的方向是2维实数域上的一个向量）。

输出：如果\((H, \vec{c})\)是无界的（unbounded），报告一条（包含于其可行域）的射线。如果\((H, \vec{c})\)是不可行的（infeasible），报告2个或3个certificate half-planes。否则就找到使得目标函数\(f_\vec{c}(p)\)最大化的并且按字典序的最小点\(p\)。

算法步骤：

  • 确定是否存在满足\(\vec{d}\cdot\vec{c}>0\)的这样一个向量\(\vec{d}\)，使得\(H\)中的所有\(h(h \in H)\)有\(\vec{d}\cdot\vec{\eta}(h) \geqslant 0\)。

  • 如果这样的\(\vec{d}\)存在

    • 计算线性规划程序\(H^{\prime}\)（即任意的\(h \in H'\)有\(\vec{d}\cdot\vec{\eta}(h) = 0\)），并确定其是否可行。

    • 如果\(H^{\prime}\)可行，报告一条射线\(\rho\)证明\((H, \vec{c})\)是无界的（unbounded），然后退出

    • 如果\(H^{\prime}\)不可行，报告\((H, \vec{c})\)不可行，然后退出

  • （如果这样的\(d\)不存在）找到两个半平面（约束条件）\(h_1, h_2 \in H\)是certificates，以证明\((H, \vec{c})\)是有界的（bounded），并且有一个唯一的按字典序的最小解。

  • 令\(v_2\)是这两个半平面\(h_1, h_2\)分界线\(\ell_1\)和\(\ell_2\)的交点

  • 令\(h_3, h_4, \dots, h_n\)是\(H\)中剩余半平面的一个随机排序序列

  • 令变量 \(i\)从\(3\)遍历到\(n\)

    • 如果上一次的optimal vertex \(v_{i-1}\)包含于当前的半平面\(h_i\)

      • 那么当前的optimal vertex \(v_i\)就是\(v_{i-1}\)，即有\(v_i = v_{i-1}\)

    • 如果上一次的optimal vertex \(v_{i-1}\)不在当前的半平面\(h_i\)之内

      • 那么找到\(\ell_i\)（半平面\(h_i\)的边界线）上的一个点\(p\)，它满足前\(i-1\)个线性条件约束（即\(H_{i-1}\)），并且它能够使得目标函数\(f_{\vec{c}}(p)\)最大化，此时，\(p\)就是要找的当前的optimal vertex \(v_i\)

      • 如果找不到这样的点\(p\)，令\(h_j, h_k\)是满足 \(h_j \bigcap h_k \bigcap \ell_i = \emptyset\)的certificates（这里\(j,k<0\)）。然后停止循环，报告这个线性程序无解，并且\(h_i, h_j, h_k\)是certificate，然后退出。

  • 变量\(i\)遍历至\(n\)结束，返回最终的optimal vertex \(v_n\)。

定理4.10 一个有\(n\)个约束的2维线性规划程序，可以以\(O(n)\)的随机平均时间复杂度解决，而且最坏情况下的空间复杂度是线性的。

Theorem 4.10 A 2-dimensional linear programming problem with \(n\) constraints can be solved in \(O(n)\) randomized expected time using worst-case linear storage.

4.6* Linear Programming in Higher Dimensions

前面小节所呈现的线性规划算法，可以推广到更一般的高维度上去。当维度不是太高的时候，和传统的算法相比，使用前面小节提到的算法更有利。

令\(H\)是\(d\)维实数域\(\mathbb{R}^d\)上\(n\)个闭合半空间的集合。给定向量\(\vec{c} = (c_1, c_2, \dots, c_d)\)，要求找到\(d\)维实数域\(\mathbb{R}^d\)上的一个点\(p=(p_1, p_2, \dots, p_n)\)，使得目标函数\(f_\vec{c}(p) := c_1p_1 + \cdots + c_dp_d\)取得最大值，同时要满足点\(p\)位于\(H\)中的任意一个半空间\(h\)中（\(h \in H\)）。当线性规划程序是有界的时候，为了确保解是唯一的，我们寻找的是按字典序的、使得目标函数\(f_\vec{c}(p)\)最大的最小点\(p\)。

和前面提到的二维平面上（解决方案）的版本一样，我们通过增量式地逐次增加约束条件，同时计算每次出现的optimal vertex。为此，我们需要确保每次的最优解（optimal solution）都是唯一的。像前面小节一样，我们首先确定线性规划程序是否无界。如果是有界的，求得一组\(d\)个certificates \(h_1, h_2, \dots, h_d \in H\)，确保解是有界的，并且有一个唯一的按字典序的最小解。我们暂时先关注该算法的主要部分，之后再讨论如何找到这\(d\)个certificates：\(h_1, h_2, \dots, h_d \in H\)。

通过检查线性规划程序是有界的，令\(h_1, h_2, \dots, h_d\)是\(d\)个certificate半空间（half-space），令\(h_{d+1}, h_{d+2}, \dots, h_n\)是\(H\)中剩余半空间的一个随机排序序列。进一步，令\(C_i\)是前\(i\)个半空间加入之后（形成一组约束条件）的可行域（可行空间），这里\(d \geqslant i \geqslant n\)：

\[ C_i := h_1 \cap h_2 \cap \cdots \cap h_i \]

令\(v_i\)表示可行域\(C_i\)的optimal vertex，即\(v_i\)使得目标函数\(f_\vec{c}\)（在前\(i\)个约束条件下）的值最大。推理4.5告诉了我们在二维平面上一种计算optimal vertex的简易方法：即当（增量式地）加入第\(i\)个半平面\(h_i\)（条件约束）的时候，当前的optimal vertex要么不变（和加入\(h_i\)之前的optimal vertex相同），要么新的optimal vertex就位于新加入的半平面\(h_i\)的分界线上。下面的推理4.11把这个结果推广到了更高维的情况下，并且它的证明也是推理4.5的一个直观的一般化推广。

推理4.11 令\(1 \leqslant i \leqslant n\)，令\(C_i\)是前\(i\)个约束条件的可行域，\(v_i\)是\(C_i\)的optimal vertex。这样就有

  • 如果\(v_{i-1} \in h_i\)，那么\(v_i = v_{i-1}\)

  • 如果\(v_{i-1} \notin h_i\)，那要么\(C_i = \emptyset\)，要么\(v_i \in g_i\)。这里\(g_i\)是\(h_i\)的分界面（hyperplane，超平面）。

如果我们用\(g_i\)表示半空间\(h_i\)的分界面，那么在\(g_i \cap C_{i-1}\)这个交集上找到的optimal vertex就是\(C_i\)的optimal vertex。

但如何找到\(g_i \cap C_{i-1}\)这个交集上找到的optimal vertex？在二维平面上这容易做到，因为所有的可能被限定到了一条直线上。首先看在三维空间上的情况。在三维空间里\(g_i\)是一个平面，并且\(g_i \cap C_{i-1}\)是一个二维平面上的凸多边形区域。在\(g_i \cap C_{i-1}\)上需要如何找到optimal vertex？答案是求解一个二维线性规划程序。定义在三维实数空间\(\mathbb{R}^3\)上的一个线性（目标）函数\(f_\vec{c}\)，在二维平面\(g_i\)上引出一个（对应的）线性函数，我们要做的就是在\(g_i \cap C_{i-1}\)上找到使得这个函数取得极值的点。如果目标函数的方向\(\vec{c}\)是和二维平面\(g_i\)垂直，那么\(g_i\)上所有的点都是同样的好（使目标函数取得同样大的极值），而按照我们的规则，我们需要找到的是按字典序的最小的点。通过正确选择目标函数可以达到这样的目的，比如，当平面\(g_i\)不是和\(x\)轴垂直的话，通过把向量\((-1, 0, 0)\)向平面\(g_i\)做投影而得到向量\(\vec{c}\)。

因此，在三维空间里，我们通过这样的方式找到\(g_i \cap C_{i-1}\)上的optimal vertex：计算所有\(i-1\)个半空间和\(g_i\)的交集，并且把下面四个向量向\(g_i\)做投影，选择第一个投影不为0的向量。

\[\begin{split} \vec{c}, \space \begin{pmatrix} -1 \\ 0 \\ 0 \end{pmatrix}, \space \begin{pmatrix} 0 \\ -1 \\ 0 \end{pmatrix}, \space \begin{pmatrix} 0 \\ 0 \\ -1 \end{pmatrix} \end{split}\]

而这会产生一个二维线性规划程序，我们可以用前面的2DRANDOMIZEDLP算法解决。

类似地，和上面解决三维空间的线性规划程序一样，\(g_i\)是\((d-1)\)维空间的一个超平面（hyperplane），我们需要在\(g_i \cap C_{i-1}\)上找到一个使得目标函数\(f_\vec{c}\)取得最大值的点。而这是一个\((d-1)\)维的线性规划程序，我们通过递归地调用解决\((d-1)\)维线性规划程序的算法版本来解决。这个递归算法的出口就是一个一维线性规划程序，它可以在线性时间里解决。

我们仍然需要确定这个线性规划程序是否有界，如果是有界的话，需要找到合适的certificates。我们首先验证推理4.9在任意维度空间也是适用的，实际上这个推理和它的证明不需要做修改（就可以推广到多维空间上）。这个推理也意味着，\(d\)维空间上线性规划程序有界的充要条件就是，某个\((d-1)\)维线性规划程序不可行（infeasible）。我们通过一个递归调用来求解这个\((d-1)\)维的线性规划程序。

如果这个\((d-1)\)维的线性规划程序是可行的（feasible），就能求得一个向量\(\vec{d}\)。此时这个\(d\)维线性规划程序要么沿着\(\vec{d}\)方向是无界的，要么就是不可行的（infeasible）。而这可以通过验证\((H^{\prime}, \vec{c})\)是否可行而确定，这里的\(H^{\prime}\)就是推理4.9中所定义的那样。而\(H^{\prime}\)中所有半空间的分界面（超平面）都是和\(\vec{d}\)平行的，所有就可以通过再解决一个\((d-1)\)维的线性规划程序来确定，这会再用到一个递归调用。

如果这个\((d-1)\)维的线性规划程序是不可行的（infeasible），那么它的解决办法就能找到\(k\)个certificate half-space \(h_1, h_2, \dots, h_k \in H (k < d)\)，从而证明\((H, \vec{c})\)是有界的。如果\(k < d\)，那么\((\{h_1, h_2, \dots, h_k\}, \vec{c})\)的optimal solution集合就是无界的。如果是这种情况，那么这些optimal solution就构成了一个\((d-k)\)维的子空间（subspace）。我们需要确定线性规划程序在这个\((d-k)\)维子空间的约束下，按照字典序是否是有界的。如果是无界的，我们就报告这个solution。如果是有界的，重复这个过程，直到求得一个有\(d\)个certificate的唯一解。

下面就是多维空间上的线性规划程序算法。这里仍然用\(g_i\)表示确定第\(i\)个超平面\(h_i\)的分界平面。

算法：\(RANDOMIZEDLP(H, \vec{c})\)

输入：一个线性规划程序\((H, \vec{c})\)，\(H\)是由\(n\)个半空间面（条件约束）组成，并且\(\vec{c} \in \mathbb{R}^d\)（即目标函数的方向是\(d\)维实数域上的一个向量）。

输出：如果\((H, \vec{c})\)是无界的（unbounded），报告一条（包含于其可行域）的射线。如果\((H, \vec{c})\)是不可行的（infeasible），报告至多\((d+1)\)个certificate half-planes。否则就找到使得目标函数\(f_\vec{c}(p)\)最大化的并且按字典序的最小点\(p\)。

算法步骤：

  • 确定是否存在满足\(\vec{d}\cdot\vec{c}>0\)的这样一个向量\(\vec{d}\)，使得\(H\)中的所有\(h(h \in H)\)有\(\vec{d}\cdot\vec{\eta}(h) \geqslant 0\)。

  • 如果这样的\(\vec{d}\)存在

    • 计算线性规划程序\(H^{\prime}\)（即任意的\(h \in H'\)有\(\vec{d}\cdot\vec{\eta}(h) = 0\)），并确定其是否可行。

    • 如果\(H^{\prime}\)可行，报告一条射线\(\rho\)证明\((H, \vec{c})\)是无界的（unbounded），然后退出

    • 如果\(H^{\prime}\)不可行，报告\((H, \vec{c})\)不可行，给出certificates，然后退出

  • （如果这样的\(d\)不存在）找到\(d\)个半空间（约束条件）\(h_1, h_2, \dots, h_d \in H\)是certificates，以证明\((H, \vec{c})\)是有界的（bounded）。

  • 令\(v_d\)是这\(d\)个半空间\(h_1, h_2, \dots, h_d\)分界线\(g_1\), \(g_2, \dots, g_d\)的交点。

  • 计算\(h_{d+1}, h_{d+2}, \dots, h_n\)是\(H\)中剩余半空间的一个随机排序序列

  • 令变量 \(i\)从\(d+1\)遍历到\(n\)

    • 如果上一次的optimal vertex \(v_{i-1}\)包含于当前的半空间\(h_i\)

      • 那么当前的optimal vertex \(v_i\)就是\(v_{i-1}\)，即有\(v_i = v_{i-1}\)

    • 如果上一次的optimal vertex \(v_{i-1}\)不在当前的半空间\(h_i\)之内

      • 那么找到\(g_i\)（确定半空间\(h_i\)的边界超平面）上的一个点\(p\)，它满足前\(i-1\)个线性条件约束（即\(H_{i-1} = (h_1, h_2, \dots, h_{i-1})\)），并且它能够使得目标函数\(f_{\vec{c}}(p)\)最大化，此时，\(p\)就是要找的当前的optimal vertex \(v_i\)

      • 如果找不到这样的点\(p\)，令\(H^\ast\)是\((d-1)\)维线性规划程序不可行的最多为\(d\)个的certificates，然后停止循环，报告这个线性程序无解，并且\(H^\ast \cup h_i\)是certificates，然后退出。

  • 变量\(i\)遍历至\(n\)结束，返回最终的optimal vertex \(v_n\)。

下面的定理说明了算法RANDOMIZEDLP的性能。尽管我们把\(d\)当做一个常量，这意味着运行时间的复杂度是\(O(n)\)，但查看运行时间如何依赖于\(d\)仍然是有帮助的。

定理4.12 对于每个给定的常量\(d\)，有\(n\)个约束条件的\(d\)维线性规划程序可以以\(O(n)\)的时间复杂度解决。

Theorem 4.12 For each fixed dimension \(d\), a \(d\)-dimensional linear programming problem with \(n\) constraints can be solved in \(O(n)\) expected time.

证明：

我们必须证明，存在一个约束\(C_d\)，这样算法才能花费最多\(C_dn\)的期望时间。接下来以\(d\)维空间做讨论。对二维空间而言，它的结果遵循定理4.10，所以这里假设\(d > 2\)。这里的推导过程，基本上和在二维空间中的推导过程一样。

首先，我们解决\((d-1)\)维空间下最多\(d\)个线性规划程序。根据推导的假设，这花费的时间是\(O(dn) + dC_{d-1}n\)。

算法花费\(O(d)\)的时间去计算\(v_d\)。测试\(v_{i-1}\)是否在\(h_i\)（\(v_{i-1} \in h_i\)）中将花费\(O(d)\)的时间。因此，如果不考虑算法中寻找使得目标函数\(f_{\vec{c}}(p)\)最大化点\(p\)的话，现在花费的时间就是\(O(dn)\)。

而算法中寻找使得目标函数\(f_{\vec{c}}(p)\)最大化点\(p\)，需要把目标函数\(f_{\vec{c}}\)对应的向量\(\vec{c}\)，以\(O(d)\)的时间向超平面\(g_i\)做投影，并且以\(O(di)\)的时间计算\(g_i\)和\(i-1\)个半空间的交集。进一步，我们要以\((d-1)\)的维度和\(i-1\)个半空间继续调用递归算法。

定义随机变量\(X_i\)，如果\(v_{i-1} \notin h_i\)，那么它的值就是\(1\)，否则就是\(0\)。因此，算法花费的总时间就是，

\[ O(dn) + dC_{d-1}n + \sum_{i=d+1}^{n}(O(di) + C_{d-1}(i-1)) \cdot E[X_i] \]

这里采用后向分析来确定\(E[X_i]\)的界。考虑已经添加了\(h_1, h_2, \dots, h_i\)的情况。可行域\(C_i\)上的optimal vertex是\(v_i\)，它是由半空间的\(d\)定义的（没想明白？？）。现在退回一步考虑。这个optimal vertex，只有当删除一个包含了\(v_i\)的半空间的时候，它才会发生改变。因为\(h_{d+1}, h_{d+2}, \dots, h_i\)是一个随机的排列，那么这种情况发生的概率就是\(d/(i-d)\)。

因此，我们得到了如下的算法运行期望时间的界（bound）。

\[ O(dn) + dC_{d-1}n + \sum_{i=d+1}^{n}(O(di) + C_{d-1}(i-1)) \frac{d}{i-d} \]

而这就是以\(C_dn\)为界的。这里\(C_d = O(C_{d-1}d)\)，所以对于和维度无关的常数\(c\)而已，\(C_d = O(c^dd!)\)。

所以当\(d\)作为一个常量出现时，可以说算法是按照线性时间运行的。但这仍然有误导性，原因是对应的常数系数\(C_d\)是一个关于\(d\)的函数，而且它随着\(d\)的增大而增大的非常快。所以，这种说法只有在维度比较低（\(d\)比较小）的情况下才更有意义。

4.7* Smallest Enclosing Discs

我们前面所采用的简单的随机方法，结果显示出人意料地强大。它不仅适用于线性规划问题，并且也能适用于其他一系列问题。本节就是这样的一个例子。

考虑这样一种情况，有一条机器手臂的基座固定在一个工作台上。这条机器手臂需要需要在不同的位置点上捡起物件，并且把它们放置在其他的位置点上。那么这条机器手臂的基座应该位于什么位置？这个位置可能是那些物件所在的位置点的中心附近的某个位置。更明确地说，这条机器手臂基座的一个好位置应该是包含了那些物件所在位置点的最小圆的中心点（圆心）。这个点能够使得机器手臂基座的位置点到其他任意一个物品位置点的所有距离中的最大值，取得最小的距离值。因此我们得到这样的一个问题：给定一个有\(n\)个位置点的集合\(P\)，每个点都位于工作台上，而且机器手臂必须够得着它，在这样的前提下，找到集合\(P\)的最小圆覆盖（最小包围圆，即smallest enclosing disc）。实际上这个最小圆覆盖是唯一的，下面的推理4.14(i)是一般化的描述。

和前面的小节一样，我们讲给出这个问题的一个随机的增量式的算法。首先，生成集合\(P\)中点的一个随机排序\(p_1, p_2, \dots, p_n\)，然后在每次加入一个点的同时，维护一个当前点集\(P_i\)的最小圆覆盖\(D_i\)。

在线性规划程序中，能帮助我们找到每次加入一个约束条件之后的optimal vertex的有利结论是，如果下一个半平面包含当前的optimal vertex，那么下一个optimal vertex就不会发生变化（即和当前的optimal vertex相同）。如果下一个半平面包含当前的optimal vertex，那么下一个optimal vertex就位于（下一个）半平面上。实际上，在最小圆覆盖问题中也是有类似的结论。

推理4.13 令变量\(i\)满足条件\(2 < i < n\)，令\(P_i\)是\(P\)中点的一个随机排列的前\(i\)个点集，\(D_i\)是\(P_i\)的最小圆覆盖。那么就有结论

  • 如果\(p_i \in D_i\)，那么\(D_i = D_{i-1}\)

  • 如果\(p_i \notin D_i\)，那么\(p_i\)就在\(D_i\)的边界（圆）上

下面会给出如何在随机的增量式的算法中使用它，而这个算法和线性规划问题算法相当类似。之后，我们会再证明这个推理。

==(图在第93页，页码是86)==

算法：\(MINIDISC(P)\)

输入：平面上\(n\)个点的集合\(P\)

输出：集合\(P\)的最小圆覆盖（smallest enclosing disc）

步骤：

  • 计算点集\(P\)的一个随机排列

  • 令\(D_2\)是集合\({p_1, p_2}\)的最小圆覆盖

  • 对变量\(i\)，从\(3\)遍历到\(n\)

    • 如果点\(p_i\)在第\(i-1\)个最小圆覆盖\(D_{i-1}\)上（\(p_i \in D_{i-1}\)），则第\(i\)个最小圆覆盖\(D_i\)就是\(D_{i-1}\)，即\(D_i = D_{i-1}\)

    • 如果点\(p_i\)不在第\(i-1\)个最小圆覆盖\(D_{i-1}\)上（\(p_i \notin D_{i-1}\)），调用函数\(MINIDISCWITHPOINT({p_1, p_2, \dots, p_{i-1}}, p_i)\)，并将返回值赋予\(D_i\)，即\(D_i = MINIDISCWITHPOINT({p_1, p_2, \dots, p_{i-1}}, p_i)\)。

  • 返回最终结果最小圆覆盖\(D_n\)

当点\(p_i\)不在第\(i-1\)个最小圆覆盖\(D_{i-1}\)上（\(p_i \notin D_{i-1}\)）时，如何解决是这个问题的关键所在。我们需要调用一个子函数来计算最小圆覆盖，条件是点\(p_i\)必须在圆的边界线上。那么如何实现这个子函数？令\(q := p_i\)，然后再次使用同样的框架，即以随机的顺序（依次）添加\(P_{i-1}\)中的每个点，并且维护\(P_{i-1} \cup \{q\}\)的最小圆覆盖，条件是\(q\)必须位于每个最小圆覆盖的边界线上。而每次添加一个点\(p_j\)之后的最小圆覆盖是否变化，会因为如下的结论而变的容易寻找：即，如果点\(p_j\)在已经包含在当前的最小圆覆盖里，那么添加了点\(p_j\)之后的最小圆覆盖不会变；如果点\(p_j\)在不包含在当前的最小圆覆盖里，那么添加了\(p_j\)之后的最小圆覆盖，必须满足的条件是，点\(p_j\)位于这个最小圆覆盖的边界线上。后者的情况中，相当于\(q\)和\(p_j\)都位于这个最小圆覆盖的边界线上。所以就有下面的子函数（子算法）。

算法：\(MINIDISCWITHPOINT(P,q)\)

输入：平面上\(n\)个点的集合\(P\)，有一个点\(q\)，并且点\(q\)位于最终找到的最小圆覆盖的边界线上。

输出：集合\(P\)的最小圆覆盖（smallest enclosing disc），并且点\(q\)位于最终找到的最小圆覆盖的边界线上。

步骤：

  • 计算集合\(P\)中点的一个随机排序序列\(p_1, p_2, \dots, p_n\)

  • 令\(D_1\)是\(q\)和\(p_1\)的最小圆覆盖，并且\(q\)位于这个最小圆覆盖的边界线上

  • 对变量\(j\)，从\(2\)遍历到\(n\)

    • 如果点\(p_j\)位于当前最小圆覆盖\(D_{j-1}\)里，那么加入点\(p_j\)之后的最小圆覆盖\(D_j\)就是\(D_{j-1}\)，即\(D_j = D_{j-1}\)。

    • 如果点\(p_j\)位于当前最小圆覆盖\(D_{j-1}\)外，那么就调用子函数\(MINIDISCWITH2POINTS({p_1, p_2, \dots, p_{j-1}}, p_j, q)\)

  • 返回最终的最小圆覆盖\(D_n\)

现在遗留问题是，子函数\(MINIDISCWITH2POINTS\)如何实现？也就是，给定一个点的集合，以及两个额外的点\(q_1\)和\(q_2\)，如何找到它们的最小圆覆盖，并且使得点\(q_1\)和\(q_2\)位于这个圆的边界线上？很简单，我们再一次使用同样的办法，即按照一种点的集合的随机排序，逐次添加点，并且动态维护一个optimal disc。当点\(p_k\)被加入时，如果它就位于当前的最小圆覆盖内部，那么什么也不需要做（也就是说加入了\(p_k\)之后的新的最小圆覆盖和之前的相同，不需要改变）；如果它位于当前的最小圆覆盖外部，那么它必须落在新的最小圆覆盖的边界线上。后者的情况就是，点\(q_1\)，\(q_2\)和\(p_k\)同时位于所求的最小圆覆盖的边界线上，而这样的圆有且仅有一个。下面的步骤就是这个子函数的实现。

算法：\(MINIDISCWITH2POINTS(P, q_1, q_2)\)

输入：平面上\(n\)个点的集合\(P\)，有两个点\(q_1\)和\(q_2\)，并且这两个点位于所求的\(P\)的最小圆覆盖的边界线上。

输出：集合\(P\)的最小圆覆盖（smallest enclosing disc）\(D_n\)，并且点\(q_1\)和\(q_2\)位于最终找到的最小圆覆盖的边界线上。

步骤：

  • 令\(D_0\)是\(q_1\)和\(q_2\)的最小圆覆盖，并且\(q_1\)和\(q_2\)都位于\(D_0\)的边界线上

  • 令变量\(k\)从\(1\)遍历到\(n\)

    • 如果点\(p_k\)位于\(D_{k-1}\)内部，那么\(D_k = D_{k-1}\)

    • 如果点\(p_k\)位于\(D_{k-1}\)外部，那么\(q_1\)，\(q_2\)和\(p_k\)这三个点就位于所寻找的\(D_k\)的边界线上（求这三点组成的任意两条线段的垂直平分线的交点，就是\(D_k\)的圆心）

  • 返回最终集合\(P\)的最小圆覆盖\(D_n\)

至此，求解一个点集\(P\)的最小圆覆盖的完整算法就实现了。但在分析它（的时间和空间复杂度）之前，需要先证明在算法中用到的一些结论是正确的。比如，如果我们添加的一个点位于当前optimal disc的外部，那么这个点就必须落在新的最小圆覆盖的边界线上。

注意，下面的推理4.14在原书上的第三点结论似乎不正确。如果\(p\)位于原先最小圆覆盖的外部，那么新的最小圆覆盖会经过点\(p\)，而此时很可能\(R\)中的点，除了\(p\)，其它点都已经位于这个最小圆覆盖的内部了。这也就是说，此时的\(P\)已经包含了\(R\)中除了\(p\)的其他所有点，因此\(P \cap R \neq \emptyset\)，这就和定理的前提已经矛盾了。因此我对其第三点描述做了修改，以正视听。

推理4.14（已更正书中错误） 令\(P\)是平面上一个点的集合，令\(R\)是\(P\)的一个子集（可能为空），即\(R \subseteq P\)。令点\(p\)是集合\(P\)中的一个点。那么就有以下结论成立：

  • (i) 如果有一个圆包含\(K = P \backslash R\)中所有的点（即\(P\)除去\(R\)中所有的点的集合），并且\(R\)中所有的点都位于这个圆的边界线上，那么这样的最小的圆是唯一的，我们把它记作\(MD(P) = md(P, R)\)

  • (ii) 如果点\(p\)位于除去它的\(P\)和\(R\)的上述最小圆的内部（即\(p \in MD(P\backslash{\{p\}}) = md(P\backslash{\{p\}, R})\)），那么\(MD(P) = md(P, R)\)就是\(MD(P) = md(P\backslash{\{p\}, R})\)，即\(md(P, R) = md(P\backslash{\{p\}, R})\)

  • (iii) 如果点\(p\)位于除去它的\(P\)和\(R\)的上述最小圆的外部（即\(p \notin MD(P\backslash{\{p\}}) = md(P\backslash{\{p\}, R})\)），那么\(MD(P) = md(P, R')\)就是\(MD(P\backslash{\{p\}}) = md(P\backslash{\{p\}}, \{R \backslash \Delta\} \cup \{p\})\)，其中\(\Delta\)是\(R\)的一个子集，即\(\Delta \subseteq R\)，\(R' = \{R \backslash \Delta\} \cup \{p\}\)。也就是说，包含了点\(p\)的\(P\)的最小圆覆盖，\(p\)就一定在它的边界线上，而之前\(R\)中的点，可能有一部分或者全部都位于这个新的最小圆覆盖上了。

==(图在第95页，页码是88)==

证明：

(i) 采用反证法。假设存在两个半径相同但圆心不同的圆都是所求的最小圆覆盖，记圆心分别为\(x_0\)和\(x_1\)，两个圆分别为\(D_0\)和\(D_1\)。根据条件，\(P\)中所有的点都要被包含在这样的最小圆覆盖中，那么就必须有\(P\)中所有的点都位于\(D_0\)和\(D_1\)的交集中，即\(P \subseteq \{D_0 \cap D_1\}\)。现在按照如下方式定义一系列连续的圆的集合\(\{D(\lambda) | 0 \leqslant \lambda \leqslant 1 \}\)。令\(z\)是\(D_0\)和\(D_1\)的边界（分别记作\(\partial{D_0}\)和\(\partial{D_1}\)）的交点。圆\(D(\lambda)\)的圆心为\(x(\lambda) := (1-\lambda)x_0 + \lambda{x_1}\)（即圆\(D(\lambda)\)的圆心位于\(x_0\)和\(x_1\)的连线上），圆\(D(\lambda)\)的半径为\(r(\lambda) := d(x(\lambda),z)\)（即圆心和\(\partial{D_0}\)和\(\partial{D_1}\)的交点之间的距离）。显然，对任意的\(\lambda\)（\(0 \leqslant \lambda \leqslant 1\)），都有\(D_0 \cap D_1 \subset D(\lambda)\)，特别地，当\(\lambda = 1/2\)时，有\(D_0 \cap D_1 \subset D(1/2)\)。更进一步，圆\(D(1/2)\)的边界\(\partial{D(1/2)}\)也经过\(\partial{D_0}\)和\(\partial{D_1}\)的所有交点。又根据条件\(R \subset \partial{D_0} \cap \partial{D_1}\)（因为条件是\(R\)位于最小圆覆盖的边界上，那么根据这里的假设，\(R\)就必须既在\(\partial{D_0}\)上，又在\(\partial{D_1}\)上，那么\(R\)就必须在它们两个的交集上，即交点上），这就意味着\(R\)也必须落在\(\partial{D(1/2)}\)上（因为\(\partial{D(1/2)}\)经过\(\partial{D_0}\)和\(\partial{D_1}\)的所有交点）。换句话说，\(D(1/2)\)也是一个包含\(P\)，并且使得\(R\)都位于其边界上的圆覆盖，而且\(D(1/2)\)的半径是严格地小于\(D_0\)和\(D_1\)的半径，这就出现了一个半径更小的且满足条件要求的最小圆覆盖，这就和开始的假设相互矛盾，因此，这样的最小圆覆盖（即\(md(P, R)\)）就是唯一的。

关于此结论的另一个证明方法比较直观：# 最小圆覆盖问题 算法步骤与证明+代码模板

这个直观方法是，仍然假设有两个存在两个半径相同但圆心不同的圆都是所求的最小圆覆盖，记圆心分别为\(x_0\)和\(x_1\)，两个圆分别为\(D_0\)和\(D_1\)。根据条件，\(P\)中所有的点都要被包含在这样的最小圆覆盖中，那么就必须有\(P\)中所有的点都位于\(D_0\)和\(D_1\)的交集中，即\(P \subseteq \{D_0 \cap D_1\}\)。然后以两个圆的交点连线的弦长作为直径，做一个新的圆，这个圆必然也包含\(P\)，而这个圆的半径又严格小于\(D_0\)和\(D_1\)，那么按照定义它才是一个最小圆覆盖，这就说明假设和前提矛盾，故由此得证。

(ii) 令\(D := md(P\backslash{\{p\}, R})\)，那么如果\(p \in D\)（点\(p\)在\(D\)中），那么显然\(D\)就包含\(P\)中所有的点，并且\(R\)中所有的点都位于其边界上。那么有没有可能存在一个比\(D\)更小的圆\(D'\)，使得它包含\(P\)中所有的点，并且\(R\)中所有的点都位于其边界上？答案是没有。原因是如果存在半径更小的圆\(D'\)满足\(D' := md(P, R)\)，那么这个圆\(D'\)实际上也就满足\(D' := md(P\backslash{\{p\}}, R)\)，因为\(P\)包含\(P\backslash{\{p\}}\)，而这就产生了另外一个最小圆覆盖，和\(D := md(P\backslash{\{p\}, R})\)矛盾。所以此时也有\(D := md(P, R)\)。

(iii) 这条结论因为书中描述的关系，我认为书中的描述证明并不正确，故按照自己的叙述证明如下。假设\(MD(P\backslash{\{p\}}) = md(P\backslash{\{p\}}, R)\)的圆心为\(x\)，半径为\(r\)（即包含点集\(P\backslash{\{p\}}\)的最小圆覆盖圆心为\(x\)，半径为\(r\)），做\(x\)和点\(p\)之间的连线\(\overline{xp}\)，取\(y\)为连线\(\overline{xp}\)上任意一点，即\(y = y(\lambda) = (1-\lambda)x + \lambda{p}\)，其中\((0 \leqslant \lambda \leqslant 1)\)。以\(y(\lambda)\)为圆心，\(r + \overline{xy(\lambda)}\)为半径做圆，显然当圆心是\(y(0)\)，半径是\(r + \overline{xy(0)}\)时，就是之前的最小圆，\(p\)位于这个圆的外面。而当圆心是\(y(1)\)，半径是\(r + \overline{xy(1)}\)时，这个圆一定包含点\(p\)（因为圆心就是\(p\)）。由于连续性，中间必然存在一个值\({\lambda}^{\ast}\)，使得以\(y({\lambda}^{\ast})\)为圆心，\(r + \overline{xy({\lambda}^{\ast})}\)为半径的圆经过点\(p\)（即点\(p\)就位于该圆的边界线上）。而这个圆必然包含\(P\backslash{\{p\}}\)中所有的点，因为该圆的半径是\(r + \overline{xy({\lambda}^{\ast})}\)，大于\(r\)。因此当加入点\(p\)到\(P\backslash{\{p\}}\)中时，如果\(p\)位于原先最小圆覆盖的外面，那么新的最小圆覆盖就一定经过\(p\)，即\(MD(P) = md(P, R')\)就是\(MD(P\backslash{\{p\}}) = md(P\backslash{\{p\}}, \{R \backslash \Delta\} \cup \{p\})\)。

定理4.15 可以用\(O(n)\)的时间复杂度和最坏情况下线性的空间复杂度，求解平面上有\(n\)个点的最小圆覆盖。

Theorem 4.15 The smallest enclosing disc for a set of \(n\) points in the plane can be computed in \(O(n)\) expected time using worst-case linear storage.

证明：

算法\(MINIDISCWITH2POINTS\)的时间复杂度是\(O(n)\)，因为循环中的每个步骤都花费的是常数时间，并且是线性空间复杂度。而算法\(MINIDISCWITHPOINT\)和\(MINIDISC\)也都花费线性空间复杂度，所以剩下的就是如何确定它们花费的时间复杂度。

对于算法\(MINIDISCWITHPOINT\)，如果不考虑调用\(MINIDISCWITH2POINTS\)，那么时间复杂度就是\(O(n)\)。那么调用\(MINIDISCWITH2POINTS\)的概率是多少？我们依然使用后验概率来分析。假如固定一个点集\(\{p_1, p_2, \dots, p_i\}\)，令\(D_i\)是这个点集\(\{p_1, p_2, \dots, p_i\}\)的最小圆覆盖，并且有额外的一个点\(q\)位于其边界线上。如果我们从点集\(\{p_1, p_2, \dots, p_i\}\)中移除一个点，那么最小圆覆盖如何变化？答案是，如果移除的点是位于最小圆覆盖边界上的三个点之一时，最小圆覆盖就会发生变化。

为什么是三个点，不是四个点或更多或更少？原因是，三个点就可以确定一个圆的圆心，使得这个圆就能同时经过这三个点，所以不是更少。而如果是四个或更多，那么移除这些点中的任意一个，最小圆覆盖不会发生变化，因为当前的最小圆覆盖依然经过这些点。

因为边界线上已经有一个点是\(q\)，那么从点集\(\{p_1, p_2, \dots, p_i\}\)中移除这个点是剩下的两个点的概率就是\(2/i\)，所以算法\(MINIDISCWITHPOINT\)的期望运行时间就是

\[ O(n) + \sum_{i=2}^{n} O(i) \cdot \frac{2}{i} = O(n) \]

而使用同样的办法分析，得到算法\(MINIDISC\)的时间复杂度也是\(O(n)\)。

实际上有不同的方法来提高算法\(MINIDISC\)的性能。首先，不必在算法\(MINIDISCWITHPOINT\)中每次都使用一个新的随机排序序列。可以在算法\(MINIDISC\)的开头，计算出一个随机排序序列，然后把这个随机排序序列当做参数传递给\(MINIDISCWITHPOINT\)。进一步，可以用一个单一的算法\(MINIDISCWITHPOINTS(P, R)\)来取代三个不同的算法的结合，就像推理4.14中叙述的那样。

4.8 Notes and Comments

本章通过引出一个铸造问题而学习了一种对应的算法难题。而制造过程中的其他问题也引出了其他各类算法难题。在过去的几十年里，计算几何对这些问题进行了研究。著有书籍的有Dutta等人。Janardan & Woo和Bose & Toussaint进行了相应的调研。

计算半平面相交的问题是一个由来已久又被充分研究的问题。在第11章中会讲到，这个问题相当于在平面上计算凸包的点。在计算几何领域里，这个问题有很长的历史，Preparata & Shamos列出了许多的解决方案。在第一章里，讲述了更多关于在二维平面上计算凸包的内容。

在二维平面上或三维空间里，计算半平面或半空间的可以在时间复杂度\(O(nlogn)\)内完成。但当维数变高之后，这就变成了一个需要算力的问题。原因是，一般情况下，（低维）多面体之间产生的交（集）的面（faces）可能高达\(\Theta(n^{\lfloor{d/2}\rfloor})\)。所以，如果目标仅仅是只寻找一个可行点（feasible point），那么显式地计算（多面体的）交（集）就不是很有吸引力的办法了。

线性规划是数值分析和组合分析中的基本问题，而调研这一领域的内容就超过了本章的范围，我们在本章也只提到了单纯形算法（simplex algorithm）和它的变种，以及Khachiyan和Karmarkar给出的多项式时间复杂度的解法。关于线性规划，更多的内容可以参考Chvatal和Schrijver的书籍。

Megiddo是第一个把线性规划当做是计算几何中问题来对待的，他说明了检查半空间的交集是否为空，要比计算交集要简单很多。对于线性规划，他给出了第一个确定性的算法，时间复杂度是\(O(C_dn)\)，\(C_d\)是一个只和维度有关的因子。他的这个算法对于任意给定的维度，都是\(n\)的线性复杂度。实际上，在他算法中\(C_d = 2^{2^d}\)，后来有改进到\(C_d = 3^{d^2}\)。最近出现了一些更简单实际的随机算法。有许多随机算法的时间复杂度是亚指数增长的（subexponential），但依旧不是关于维度的多项式时间。寻找一个严格的多项式时间的算法（组合多项式复杂度的），仍旧是线性规划领域中一个主要的待解决的开放课题。

本章给出的二维或更高维度的简单随机增量式算法，源自Seidel。和本章中讲述的不同的是，他在处理无界的线性规划程序时，用符号代替了本章中的因子\(M\)。这也许更有效和优雅，而且可以显示无界的\(d\)维线性规划程序和\(d-1\)维的可行线性规划程序之间的关系。在Seidel的版本中，因子\(C_d\)是\(O(d!)\)。

计算最小圆覆盖的一般化算法源自Welzl。他同时也说明了如何在高维空间里找到一个点集的最小覆盖球（smallest enclosing ball），或者最小覆盖椭圆（ellipse），或者最小覆盖椭圆体（ellipsoid）。

进一步，Sharir和Welzl总结了相关的技术，并提出了\(LP-type\)问题，这些问题可以使用【189】和【354】中提到的算法解决。一般来说，这种技术适用于这样的优化问题：当一个新的约束加入时，解决方案不变，或者新的约束部分定义了解决方案而使得算法的维度有所降低。同时，也显示\(LP-type\)问题的某些特性，引出了所谓的\(Helly-type\)定理。

随机化通常是一种使得算法能够变简单和有效的技术。后面的章节里会遇到更多。（使用了这种技术后）我们付出的代价是运行时间就是一个期望值，就像我们观察到的一样，某些情况下，算法的运行时间会更长。由于这个原因，有些人认为随机算法是不可靠的，并且不应该被使用。比如在核电站中或重症监护站（intensive care station）中使用的计算机。

另一方面，确定性的算法只有在理论上是完美的。在实际当中，任何稍微复杂一点的程序都可能包含bug，哪怕我们忽略这一点，也会有硬件故障或软件错误的风险，比如核心内存中单个比特在\(\alpha\)射线的影响下发生翻转错误。而由于随机化算法通常更简单，并且代码更少，那么产生这种不幸情况的可能就更少。因此，在要求的时间里，随机化算法无法计算出正确答案的概率，通常不比确定性算法失败的可能性更大。此外，通常我们允许一个较大的期望运行时间，以此来降低随机化算法的运行时间超过运行时间的概率。

4.9 Exercises

4.10 References

4.10.1 Latex 数学字体示例

\(\mathnormal{ABCDEabcde1234}\)

\(\mathrm{ABCDEabcde1234}\)

\(\mathit{ABCDEabcde1234}\)

\(\mathbf{ABCDEabcde1234}\)

\(\mathsf{ABCDEabcde1234}\)

\(\mathtt{ABCDEabcde1234}\)

\(\mathcal{ABCDEabcde1234}\)

\(\mathscr{ABCDEabcde1234}\)

\(\mathfrak{ABCDEabcde1234}\)

\(\mathbb{ABCDEabcde1234}\)

4.10.2 渐近界

\(\Theta\) ：渐近紧确界，表示既为上界，也为下界（tight），即相等（\(=\)）、准确的复杂度。

\(O\) ：渐近上界，表示渐近上界（tightness unknown），即小于等于（\(\le\)），近似复杂度。

\(\Omega\) ：渐近下界，表示渐近上界（tightness unknown），即大于等于（\(\ge\)），近似复杂度。

\(o\) ：非渐近紧确上界，表示上界（not tight），即明确小于（\(\lt\)），准确计算得出。

\(\omega\) ：非渐近紧确上界，表示下界（not tight），即明确大于（\(\gt\)），准确计算得出。

时间复杂度的五个记号

算法导论-渐近记号

4.10.3 Latex Mathematical Symbols

A PDF Link of Latex Mathematical Symbols

4.10.4 线性规划

线性规划简介（一） - 知乎

线性规划简介（一） - 知乎

5 Orthogonal Range Searching Querying a Database

本章引言部分，以个人数据管理为例子，说明了数据库（database）和几何之间存在某种意义上的联系。在所举例子中，对个人数据中两项类别（field）数据的联合查询，相当于在二维平面上对代表个人数据的点集的查询（query）。进一步，如果查询的是三项类别的数据，那么就相当于在三维空间里对代表个人数据的点集的查询。

查询数据记录（record）中某些数据域（field）处于某些设定的值的范围内，这样的查询可以转换为对\(d\)维空间上，沿着\(d\)个坐标轴平行的超空间（box）里，查询其中对应的数据点集。这样的查询，在计算几何中，叫做区域范围查询，或者正交范围查询。

A query asking to report all records whose fields lie between specified values then transforms to a query asking for all points inside a \(d\)-dimensional axis-parallel box. Such a query is called a rectangular range query, or an orthogonal range query.

5.1 1-Dimensional Range Searching

在讲解二维或高维矩形区域查询之前，先介绍了一维矩形区域查询。即给定实数轴上的一系列实数，然后查询位于一维矩形区域（即区间）\([x:x^{\prime}]\)中的点。

这里用到的数据结构是平衡二叉搜索树（Balanced Binary Search Tree，BBST），记作\(\mathcal{T}\)。记\(P := \{p_1, p_2, \dots, p_n\}\)是实数轴上要查询的点集。这个二叉树的叶子节点就是这些待查询的点，而二叉树的中间节点，都是用来指导如何搜索的分叉节点（但实际上也都属于待查询的点中的一部分）。把存储在分叉节点\(\nu\)中的分割值（splitting value）记作\(x_{\nu}\)。那么根据二叉树的性质，节点\(\nu\)的左子树中的所有节点的值，都是小于或等于\(x_{\nu}\)，而右子树中所有节点的值，都是严格地大于\(x_{\nu}\)。

报告出区间\([x:x^{\prime}]\)的中点的步骤如下。我们在这个平衡二叉搜索树\(\mathcal{T}\)中搜索\(x\)和\(x'\)，当搜索结束时，记它们对应的叶子节点分别是\(\mu\)和\({\mu}'\)（比如从根节点出发，开始搜索\(x\)，沿着分叉一直搜索到叶子节点，把这个叶子节点记作\(\mu\)）。那么区间\([x:x^{\prime}]\)的中点就是\(\mu\)和\({\mu}'\)之间的那些叶子节点（上的值），而且还可能包括\(\mu\)和\({\mu}'\)这两个叶子节点（上的值）。

这里以一个具体的平衡二叉搜索树为例，说明了查询一个区间\([18:77]\)的情况。

从图中可以看到的是，\(\mu\)和\({\mu}'\)之间的节点如何找到？就是从通向\(\mu\)和\({\mu}'\)的搜索路径中的某个节点开始的（左或右）子树的叶子节点。为了找到要从（通向\(\mu\)和\({\mu}'\)的搜索路径中的）那个节点开始去报告其左或右子树，首先要找到通向\(\mu\)和\({\mu}'\)的搜索路径是从哪个节点开始分叉的。把这个分叉的节点记作\({\nu}_{split}\)，记一个节点\(\nu\)的左孩子和右孩子分别为\(lc(\nu)\)和\(rc(\nu)\)。

算法：\(FindSplitNode(\mathcal{T}, x, x')\)

输入：一棵（平衡二叉搜索）树，两个值\(x\)和\(x'\)代表搜索区间，且有\(x \leqslant x'\)。

输出：通向\(x\)和\(x'\)的搜索路径上分叉的节点\(\nu\)，或者是两条路径终止的叶子节点。

步骤：

  • 令\(\nu\)是树的根节点（root）

  • （循环）如果\(\nu\)不是叶子节点，并且其节点的值\(x_{\nu}\)严格小于区间下限\(x\)（\(x_{\nu} \lt x\)），或大于等于区间上限\(x\)（\(x_{\nu} \geqslant x'\)），就做如下操作

    • 如果当前的节点的值\(x_{\nu}\)大于等于区间上限\(x'\)，即\(x_{\nu} \geqslant x'\)

      • 把左孩子\(lc(\nu)\)赋予\(\nu\)（即向左走）

    • 如果当前的节点的值\(x_{\nu}\)严格小于区间下限\(x\)，即\(x_{\nu} \lt x\)

      • 把右孩子\(rc(\nu)\)赋予\(\nu\)（即向左走）

  • 返回循环结束时的节点\(\nu\)

从这个分割节点\({\nu}_{split}\)开始，如果沿着通向\(x\)（搜索下限）的路径，就报告这个路径上每个节点的右子树的所有叶子节点（因为它们的值都是位于搜索区间里的）；如果沿着通向\(x'\)（搜索上限）的路径，就报告这个路径上每个节点的左子树的所有叶子节点（因为它们的值也都是位于搜索区间里的）。最后，检查这两条路径结束的叶子节点的值，查看它们的值是否位于搜索区间\([x:x']\)。

下面的\(1DRangeQuery\)算法利用了子函数\(ReportSubtree\)。这个子函数从一个给定的节点开始，遍历整个树，并报告其所有叶子节点。因为任何二叉树的中间节点都比其叶子节点少，所以这个子函数需要的时间和其报告的节点的数量成线性关系。

算法：\(1DRangeQuery(\mathcal{T}, [x: x'])\)

输入：一棵（平衡二叉搜索）树，两个值\(x\)和\(x'\)代表搜索区间\([x: x']\)，且有\(x \leqslant x'\)。

输出：\(\mathcal{T}\)中所有位于搜索区间\([x: x']\)的点。

步骤：

  • 调用函数\(FindSplitNode(\mathcal{T}, x, x')\)找到分叉节点\({\nu}_{split}\)

  • 如果节点\({\nu}_{split}\)已经是叶子节点

    • 检查该节点所存储的值是否位于搜索区间并被报告出来

  • 否则

    • 沿着通向\(x\)的搜索路径，报告路径上每个点的右子树的所有叶子节点

      • 令\(\nu\)是分叉节点\({\nu}_{split}\)的左孩子

      • （循环）如果\(\nu\)不是叶子节点

        • 如果当前节点的值\(x_{\nu}\)大于或等于搜索下限，即\(x_{\nu} \geqslant x\)

          • 调用子函数\(ReportSubtree\)来报告其右子树，\(ReportSubtree(rc(\nu))\)

          • 然后令\(\nu\)为左孩子（即向左走）

        • 如果当前节点的值\(x_{\nu}\)小于搜索下限\(x\)，即\(x_{\nu} < x\)

          • 令为右孩子（即向右走）

      • 循环结束时，\(\nu\)是叶子节点，检查其值是否位于搜索区间并报告

    • 沿着通向\(x'\)的搜索路径，报告路径上每个点的左子树的所有叶子节点

      • 令\(\nu\)是分叉节点\({\nu}_{split}\)的右孩子

      • （循环）如果\(\nu\)不是叶子节点

        • 如果当前节点的值\(x_{\nu}\)小于或等于搜索上限，即\(x_{\nu} \leqslant x'\)

          • 调用子函数\(ReportSubtree\)来报告其左子树，\(ReportSubtree(lc(\nu))\)

          • 然后令\(\nu\)为右孩子（即向右走）

        • 如果当前节点的值\(x_{\nu}\)大于搜索上限\(x'\)，即\(x_{\nu} > x'\)

          • 令为左孩子（即向左走）

      • 循环结束时，\(\nu\)是叶子节点，检查其值是否位于搜索区间并报告

时间复杂度：

引理5.1 算法\(1DRangeQuery\)报告了位于搜索区间中的每一个点，不多也不少。

证明：

关于证明算法中报告的任何一个点\(p\)，都位于搜索区间中。如果\(p\)是在搜索\(x\)或\(x'\)的路径终点的叶子节点上被报告出来时，\(p\)本身就会经过检查，确保它一定位于搜索区间\([x:x']\)当中。如果不是，那么\(p\)一定是在子函数\(ReportSubtree\)中被报告出来的。

假设\(p\)是在通往搜索\(x\)的路径上在节点\(\nu\)处，通过调用\(ReportSubtree(rc(\nu))\)报告出来的（因为此时一定是通过调用某个节点的右孩子）。首先\(\nu\)和它的右孩子\(rc(\nu)\)都位于分叉节点\({\nu}_{split}\)的左子树上，而\(p\)又在\(rc(\nu)\)为根的树上，所以\(p\)一定小于\(x_{{\nu}_{split}}\)，即\(p \leqslant x_{{\nu}_{split}}\)。而搜索\(x'\)的路径是在\({\nu}_{split}\)的右子树上，所以就有\(p < x'\)。而搜索\(x\)的路径在\(\nu\)处向左，\(p\)又是在\(\nu\)的右子树上，所以\(p > x\)。所以最终就有\(p \in [x:x']\)。

如果\(p\)是通往搜索\(x'\)的路径上通过\(ReportSubtree\)调用而被报告出来的，分析和上面类似，结论也一致。

而关于位于搜索区间中的每一个点都被报告出来的证明，书上的分析始终没有看懂（2023-02-18），需要之后再看。

关于性能。因为是平衡二叉搜索树，所以建立它的空间复杂度是\(O(n)\)，时间复杂度是\(O(n\log{n})\)。在最坏情况下，给定集合中的所有点都可能在搜索区域中，此时查询时间复杂度就是\(\Theta(n)\)，而且如果要报告所有的点，那么这个时间复杂度实际上是避免不了的。所以这里的分析不仅仅考虑\(n\)（点集中点的个数），还考虑报告了的点的个数\(k\)。这就是第二章提到过的output-sensitive。

据此，回想前面提到\(ReportSubtree\)花费的时间是和报告的点个数线性相关的，所以它的时间复杂度是\(O(k)\)。而通向\(x\)或\(x'\)的搜索路径，其长度（深度）是\(O(\log{n})\)，访问每个节点花费常数时间\(O(1)\)，所以访问通向\(x\)或\(x'\)的搜索路径上的所有节点时间是\(O(\log{n})\)。最终得到该算法的时间复杂度就是\(O(\log{n} + k)\)。

下面的定理总结了一维区域搜索的特点。

定理5.2 \(P\)是一维空间上的点集，它可以用平衡二叉搜索树来存储，而建立该树的空间复杂度和时间复杂度分别是\(O(n)\)和\(O(n\log{n})\)。报告给定搜索区间中的点的查询时间复杂度是\(O(k + \log{n})\)，这里\(k\)是报告出来的点的个数。

5.2 Kd-Trees

本节讲解的是二维矩形区域搜索问题。记\(P\)是平面上\(n\)个点的集合。这里做了个假设，任何两个点的\(x\)坐标或\(y\)坐标都是不同的。虽然实际中这样的假设是不可能的，但后面5.5节可以耍个把戏来解决这个问题。

二维平面的矩形区域搜索，就是查询点集\(P\)位于矩形区域\([x:x'] \times [y:y']\)中的点。一个点\(p := (p_x, p_y)\)位于这个区域的充要条件是，

\[ p_x \in [x:x'] \space \space \mathbf{and} \space \space p_y \in [y:y'] \]

可以说二维平面上的矩形区域搜索，就就包含了两个一维矩形区域的子搜索，分别是对点的\(x\)坐标和\(y\)坐标的搜索。

一维空间（就是实数轴）搜索利用的是一个平衡二叉搜索树。它的构建是，首先把（一维）点分为大致相同数量的两部分，一部分小于等于分割点的值，另一部分大于分割点的值。这个分割点就是树的根，而这两部分就通过递归的方式存储在它的左右两棵子树里。

在二维空间中，每个点有\(x\)和\(y\)坐标。因此可以对点集，先按\(x\)坐标做分割，然后再按\(y\)坐标分割，之后再按照\(x\)做峰，依次类推，直到将整个点集分割完毕。

更具体地，首先用一条垂直线把点集\(P\)一分为二，每部分点的个数基本相等。这条分割线就是树的根。位于分割线左边或分割线上的点，属于左边点的集合\(P_{left}\)，它是存在根的左子树里。位于分割线右边的点，属于右边点的集合\(P_{right}\)，它是存在根的右子树里。

然后，把根的左子树里的点\(P_{left}\)，用一条水平直线再分为两个大致相同数量的子集。这个水平直线就是根的左孩子。位于这条水平直线下面或水平直线上的点，存在这个左孩子的左子树里，而位于水平直线上方的点，存在这个左孩子的右子树里。对于跟的右子树里的点\(P_{right}\)，按照上面同样的方式分割并存储在树里。而到根的孙子节点时，就再用垂直线划分并继续构建树，依次类推。

所以，当树的节点所在的深度是奇数时，使用了垂直线划分；当树的节点所在的深度是偶数时，使用了水平线划分。最终形成的这样一棵二叉树，就叫做\(\mathbf{KD-tree}\)。这个名字最早的意义是\(k\)-dimensional tree，但现在这个原始的意义已经在演化过程中丢失了，所以上面提到的这个二叉树现在叫做2-dimensional kd-tree，而不是2d-tree了。

下面的递归算法\(BuildKDTree\)用来构建一棵二维KD-tree。它有两个参数，点集\(P\)和一个整数。集合\(P\)就是用来构建KD-tree的点的集合，而整数是子树的根的深度，它会递归地调用这个算法来继续构建，所以它很重要，因为它决定了子树里的点是用垂直线还是水平线分割的。在最开始，这个整数是0（即调用这个算法的地方），表明点集\(P\)对应的KD-tree的根节点的深度是0。

算法：\(BuildKDTree(P, depth)\)

输入：点集\(P\)，以及当前的深度\(depth\)。

输出：存储点集\(P\)的KD-tree的根节点

步骤：

  • 如果\(P\)中只有一个点

    • 直接返回存储这个点的树的节点

  • 如果深度\(depth\)是偶数

    • 按照当前点集\(P\)中点的\(x\)坐标，用一条垂直线\(\ell\)把点集\(P\)分为两部分。记\(P_1\)是位于垂直线\(\ell\)左边或垂直线\(\ell\)上的点的集合，记\(P_2\)是位于垂直线\(\ell\)右边的点的集合。

  • 如果深度\(depth\)是奇数

    • 按照当前点集\(P\)中点的\(y\)坐标，用一条水平线\(\ell\)把点集\(P\)分为两部分。记\(P_1\)是位于水平线\(\ell\)下方或水平线\(\ell\)上的点的集合，记\(P_2\)是位于垂直线\(\ell\)上方的点的集合。

  • 递归调用函数本身\(BuildKDTree(P_1, depth + 1)\)，得到的节点\({\nu}_{left}\)。

  • 递归调用函数本身\(BuildKDTree(P_2, depth + 1)\)，得到的节点\({\nu}_{right}\)。

  • 创建一个节点\(\nu\)用来存储当前的垂直线或水平线\(\ell\)，然后连接前面得到的两个节点，\({\nu}_{left}\)是\(\nu\)的左孩子，\({\nu}_{right}\)是\(\nu\)的右孩子。

  • 返回节点\(\nu\)。

时间复杂度：\(O(n\log{n})\)

空间复杂度：\(O(n)\)

这里的约定是，如果一个点位于垂直或水平分割线上，就把它归于垂直分割线的左边，或水平分割线的下方。垂直或水平分割线实际上就是\(x\)或\(y\)坐标的中间值（median）。所以\(n\)个数的中间数是就是第\(\lceil n/2 \rceil\)小的数，即向下取整。

在建造KD-tree的过程中，最耗时的步骤是用来寻找分割线的递归调用。它根据depth的值是偶数还是奇数，从而决定分割线应该是垂直还是水平，进而决定寻找的是\(x\)坐标的中间数（median）还是\(y\)坐标的中间数。寻找这样中间数的线性时间算法比较复杂，而把点集按照\(x\)坐标和\(y\)坐标分别排序，然后把这两个列表当做两个参数传递，就是更好的办法。这样不仅容易查找\(x\)或\(y\)坐标的中间数，还能以线性时间，基于给定的两个点的列表再构建排序列表，继而给递归调用使用。

所以构建树的时间复杂度是

\[\begin{split} T(n) = \left\{ \begin{array}{lcl} O(1), & \mathrm{if} & n = 1 \\ O(n) + 2 T(\lceil n/2 \rceil), & \mathrm{if} & n > 1 \\ \end{array}\right. \end{split}\]

所以时间复杂度最终是\(O(n\log{n})\)，它也包含了把点按\(x\)和\(y\)坐标预排序所花费的时间。

对于空间复杂度，KD-tree的每个叶子节点存储的是独一无二的点，所以共有\(n\)个叶子。而KD-tree也是二叉树，所以中间节点和叶子节点也花费同样的存储空间。所以最终的空间复杂度就是\(O(n)\)。

引理5.3 对\(n\)个点的点集构建一个KD-tree，时间复杂度是\(O(n\log{n})\)，空间复杂度是\(O(n)\)。

Lemma 5.3 A kd-tree for a set of n points uses \(O(n)\) storage and can be constructed in \(O(n\log{n})\) time.

KD-tree中，某一个节点，代表了平面区域上的由一条或几条分割线限定的一块区域。如，根的左孩子的左孩子，这个节点代表了在第一条垂直分割线左边的半个平面里，再由一条水平分割线划分的这半个平面的下方区域。所以可以把树中的任意一个节点\(\nu\)都认为是平面上的一个矩形区域，当然这个矩形区域可能是由分割线（splitting line）限定的有限区域，也可能是某几边没有分割线限定的半开放区域。

书中举例说明了一个节点代表的（开放或封闭的）矩形区域。

==(图在第108页，页码是102)==

把KD-tree中一个节点对应的区域记作\(region(\nu)\)，那么平面上的一个点位于以\(\nu\)为根的子树上的充要条件就是这个点位于区域\(region(\nu)\)中。因此，当待查询的区域和\(region(\nu)\)相交时，就要查找以\(\nu\)为根的子树。因此，查询算法就是，遍历KD-tree，但只访问那些和待查询区域相交的节点。当一个节点对应的矩形区域完全包含于待查询区域时，就报告以这个节点为根的子树的所有叶节点；（而不是完全包含时）就一直检查是否和待查询区域相交，直到叶子节点，然后判断对应的点是否在待查询区域当中，再决定是否报告。

书中这里举例，说明了在一个KD-tree中，那个节点对应的矩形区域被待查询区域完全包含，从而报告其所有叶子节点。而另一些节点直到遍历到叶子节点时（遍历路径上的点对应的矩形区域肯定都和待查询区域相交），才判断点是否在待查询区域当中。

==(图在第109页，页码是103)==

下面给出的查询算法，以KD-tree的根节点，和待查询区域\(R\)作为输入参数，使用了\(ReportSubtree(\nu)\)子函数报告其以\(\nu\)为根的树的所有叶子节点。\(lc(\nu)\)和\(rc(\nu)\)分别表示节点\(\nu\)的左孩子和右孩子。

注意，算法中的待查询区域可以是任意区域，不必一定是矩形区域。

算法：\(SearchKDTree(\nu, R)\)

输入：一个KD-tree的根节点\(\nu\)，以及待查询的区域\(R\)

输出：\(\nu\)节点下面所有位于待查询区域\(R\)当中的叶子节点

步骤：

  • 检查节点\(\nu\)是叶子节点

    • 如果节点\(\nu\)存储的点位于待查询的区域\(R\)当中，就报告它

  • 检查节点\(\nu\)的左孩子节点所代表的（封闭或开放的）区域\(region(lc(\nu))\)是否完全被待查询区域\(R\)所包含

    • 如果完全被包含，就调用子函数\(ReportSubtree(lc(\nu))\)

    • 如果不是完全被包含但是和\(R\)相交，就递归调用\(SearchKDTree(lc(\nu), R)\)

  • 检查节点\(\nu\)的右孩子节点所代表的（封闭或开放的）区域\(region(rc(\nu))\)是否完全被待查询区域\(R\)所包含

    • 如果完全被包含，就调用子函数\(ReportSubtree(rc(\nu))\)

    • 如果不是完全被包含但是和\(R\)相交，就递归调用\(SearchKDTree(rc(\nu), R)\)

算法中的主要的计算来自于检查节点\(\nu\)代表的区域和待查询区域是否相交。可以在算法开始前对每个节点计算并存储，但更好的办法是维护一个区域，然后根据当前的水平或垂直分割线来更新它。

这里以一条垂直分割线举例。在算法中，垂直分割线代表的是偶数深度。那么此时节点\(\nu\)的左孩子所对应的区域就可以通过下面的办法计算

\[ region(lc(\nu)) = region(\nu) \cap \ell(\nu)^{left} \]

这里\(\ell(\nu)\)是节点\(\nu\)存储的分割线，而\(\ell(\nu)^{left}\)表示的是分割线左边的半平面（包括分割线）。

==(图在第109页，页码是103)==

引理5.4 使用一个和坐标轴平行的矩形区域，对一个存储\(n\)个点的KD-tree进行查询，时间复杂度是\(O(\sqrt{n} + k)\)，这里\(k\)是报告的点的个数。

Lemma 5.4 A query with an axis-parallel rectangle in a kd-tree storing \(n\) points can be performed in \(O(\sqrt{n}+k)\) time, where k is the number of reported points.

证明：

首先，遍历一棵子树，并且报告其叶子节点，这样的时间复杂度是和其报告的叶子节点数成线性关系，即\(O(k)\)，这里\(k\)是报告的点的个数。

其次需要确定的是由算法遍历访问到的、和待查询区域\(R\)相交的节点（就是和节点对应的区域相交）的个数。这些节点对应的区域和待查询区域\(R\)相交，但不是被待查询区域\(R\)完全包含。分析的办法就是，考察一条垂直线和树中节点（对应的区域）的相交关系。考察垂直线就是确定待查询区域\(R\)的左右边界相交的区域个数（即节点个数），类似地，考察水平线就能确定待查询区域\(R\)的上下边界相交的区域个数（亦即节点个数）。

记\(\ell\)是任意的一条垂直线，\(\mathcal{T}\)是一棵KD-tree，\(\ell(root(\mathcal{T}))\)是KD-tree的根节点存储的垂直分割线。那么\(\ell\)要么就和\(\ell(root(\mathcal{T}))\)左边的那些区域相交，要么就和\(\ell(root(\mathcal{T}))\)右边的那些区域相交，而不可能同时相交。这似乎意味着这条垂直线和树中相交的区域个数（即节点个数）是\(Q(n) = 1 + Q(n/2)\)。但这个递归式子不成立，原因是根的孩子节点对应的分割线是水平的。比如，\(\ell\)和根的左孩子对应的区域\(region(lc(root(\mathcal{T})))\)相交，那么\(\ell\)就和\(lc(root(\mathcal{T}))\)的下面两个孩子都相交了，因此递归式子不成立了。

为了解决这个问题，就需要再向下看一层，即考察根的孙子节点。根的孙子节点因为深度是偶数（这里就是2），所以存储的分割线也是垂直线，记作\(\ell_{gc}\)，那么\(\ell\)肯定就和\(\ell_{gc}\)的左边或右边的区域相交（同理不可能和左右的区域同时相交），所以相交的区域就是至多\(\lceil\lceil n/2 \rceil / 2 \rceil\)，亦即\(n/4\)。所以递归下去，就得到如下的递归式子。

\[\begin{split} Q(n) = \left\{ \begin{array}{lcl} O(1), & \mathrm{if} & n = 1 \\ 2 + 2 Q(n/4), & \mathrm{if} & n > 1 \\ \end{array}\right. \end{split}\]

这个递归式子最终得到\(Q(n) = O(\sqrt{n})\)，这意味着一条垂直线（即待查询区域的左或右边界线）和KD-tree中\(O(\sqrt{n})\)个节点（对应的区域）相交。同理，一条水平线（即待查询区域的上或下边界线）和KD-tree中\(O(\sqrt{n})\)个节点（对应的区域）相交。那么最终就得到，待查询区域和KD-tree中\(O(\sqrt{n})\)个节点（对应的区域）相交。因此也得到花费在和相交的节点上的时间就是\(O(\sqrt{n})\)。

至此，得到总的时间复杂度就是\(O(\sqrt{n} + k)\)。

前面的分析实际上是悲观的分析，因为是按照垂直线或水平线和相交的区域个数计算的。但实际情况中，待查询区域可能比较小，因此边界也比较短，所以实际上相交的区域个数就比较少。

定理5.5 二维平面上有\(n\)个点的集合\(P\)，可以用\(O(n\log{n})\)的时间复杂度和\(O(n)\)的空间复杂度，建立对应于\(P\)的KD-tree。用平面上一个矩形做查询，时间复杂度是\(O(\sqrt{n} + k)\)，这里\(k\)是报告的点的个数，即位于矩形区域内的点的个数。

Theorem 5.5 A kd-tree for a set P of \(n\) points in the plane uses \(O(n)\) storage and can be built in \(O(n\log{n})\) time. A rectangular range query on the kd-tree takes \(O(\sqrt{n} + k)\) time, where \(k\) is the number of reported points.

同样地，KD-tree可以用在三维或更高维空间上，而且算法和二维平面上的十分相似。在\(d\)维空间中\(\mathbb{R}^d\)，有\(d\)个维度，任意一个点有\(d\)个坐标，即\(x = (x_1, x_2, \dots, x_d)\)。我们首先把点集\(P\)中所有的点，沿着第一个维度的坐标轴\(x_1\)-axis的一个垂直（超）平面（hyperplane），大致分为两个数量相等的集合。这第一个超平面就是整个KD-tree的根节点。然后对这两个点集，分别用一个和第二个维度的坐标轴\(x_2\)-axis垂直的一个超平面再分别进行划分。这两个超平面就是根的孩子节点（此时深度是1）。依次类推，直到深度是\(d-1\)的时候，空间是用和\(x_d\)坐标轴垂直的超平面划分的。再往下到\(d\)深度，就又开始用和第一个维度坐标轴垂直的超平面进行划分。这种递归直到最后只剩一个点才结束。

可以看到，\(d\)维的KD-tree也是一个二叉树，而且叶子节点就是对应的点集\(P\)中的\(n\)个点。所以，构建这样的\(d\)维的KD-tree，时间复杂度和空间复杂度也都分别是\(O(n\log{n})\)和\(O(n)\)。而因为\(d\)维的关系，查询一个\(d\)维子空间里点的个数的时间复杂度是\(O(n^{1-1/d}+k)\)，分析和前面二维空间上的分析十分类似。

5.3 Range Trees

上一节提到的KD-tree，在最终需要报告的点较少的情况下，查询的时间复杂度较高（\(O(\sqrt{n} + k)\)），本节介绍的是range tree，查询时间复杂度是\(O(\log^2{n} + k)\)，有更好的查询性能。当然代价是空间复杂度增加了，从KD-tree的\(O(n)\)增加到了\(O(n\log{n})\)。

二维的区域搜索本质上是由两个一维搜索组成的，KD-tree也是由此启发而得出。对于range tree，也是类似的思路，但有所不同。

记\(P\)是\(n\)个点的集合，将在它上做矩形区域查询。查询矩形区域是\([x:x'] \times [y:y']\)。首先查询\(x\)坐标位于区间\([x:x']\)的点，然后再考虑\(y\)坐标。如果只考虑\(x\)坐标的查询，那么它就是一维区域查询，5.1节介绍了相关的算法。简单地说，就是用基于\(x\)坐标的点的二叉搜索树。大体步骤是：从根向下，找到通往\(x\)和\(x'\)的路径的分叉节点\({\nu}_{split}\)，之后左边的路径通往\(x\)，右边通往\(x'\)。报告左边路径上每个节点\(\nu\)的右子树的所有叶子节点，也报告右边路径上每个节点\(\nu\)的左子树的所有叶子节点。最后检查两条路径的叶子节点\(\mu\)和\({\mu}'\)是否也位于查询区间里。

==（这里显示查询的示意图在112页，页码是106）==

把每个节点\(\nu\)的叶子节点所存储的点的集合，叫做\(P\)的正则子集，canonical subset，记作\(P(\nu)\)。根节点的正则子集就是\(P\)集合本身。前面提到的\(x\)坐标位于查询区间的点的集合，是一系列不相交并集（disjoint union），也就是一些根节点是\(\nu\)的子树的正则子集的集合。我们感兴趣的是其中点的\(y\)坐标位于区间\([y:y']\)的点，而这正是令一个一维区域查询。如果我们有一个基于\(y\)坐标的点的二叉搜索树，我们就能完成对\(y\)坐标的查询。这就引出下面的数据结构，它用来在平面上对有\(n\)个点的集合做矩形区域查询。

  • 主树\(\mathcal{T}\)是一个平衡二叉搜索树，是基于\(P\)中点的\(x\)坐标

  • 主树\(\mathcal{T}\)中的每个中间节点\(\nu\)，它的正则子集\(P(\nu)\)都有一个对应的基于其集合中点的\(y\)坐标的平衡二叉搜索树\(\mathcal{T}_{assoc}(\nu)\)，而节点\(\nu\)有一个指向\(\mathcal{T}_{assoc}(\nu)\)的指针，这个\(\mathcal{T}_{assoc}(\nu)\)结构就叫做\(\nu\)的联合结构（associated structure）。

上述数据结构就是区域树（range tree，也叫范围树）。存有指向联合结构的节点的数据结构，叫做多层次数据结构（multi-level data structures）。主树\(\mathcal{T}\)叫做第一层树，而联合结构是第二层树。这种多层次数据结构在第10章和第16章分别有示例。下面的图是这种多层次数据结构的一个示意图。

==（多层次数据结构的示意图在113页，页码是107）==

下面是构建区域树的递归算法，它的输入是以\(x\)坐标排序的\(n\)个点的集合\(P\)（\(P := \{p_1, p_2, \dots, p_n\}\)），返回的是\(P\)对应的二维区域树\(\mathcal{T}\)。（这里仍然假设每两个点的\(x\)或\(y\)坐标不同，在5.5节则会解除这个限制）。

算法：\(Build2DRangeTree(P)\)

输入：平面上有\(n\)个点的集合\(P\)。

输出：二维区域树\(\mathcal{T}\)的根节点

步骤：

  • 创建联合数据结构：以\(P\)中每个点的\(y\)坐标为值构建二叉搜索树\(\mathcal{T}_{assoc}(\nu)\)。该二叉搜索树\(\mathcal{T}_{assoc}(\nu)\)的每个叶子节点存储的不仅仅是每个点对应的\(y\)坐标，而且还要有这个点本身。

  • 如果\(P\)中只有一个点

    • 创建一个叶子节点\(\nu\)，它存储的就是这一个点，然后令\(\mathcal{T}_{assoc}\)为\(\nu\)的联合数据结构。

  • 否则（如果\(P\)中不止一个点）

    • 把\(P\)按照点的\(x\)坐标分为两个子集，一个子集\(P_{left}\)中每个点的\(x\)坐标小于或等于\(P\)中所有点的\(x\)坐标的中间值\(x_{mid}\)，另一个子集\(P_{right}\)中每个点的\(x\)坐标严格大于\(x_{mid}\)。

    • 调用该递归函数本身\(Build2DRangeTree(P_{left})\)，返回这个子集对应的二维区域树的根节点\({\nu}_{left}\)。

    • 调用该递归函数本身\(Build2DRangeTree(P_{right})\)，返回这个子集对应的二维区域树的根节点\({\nu}_{right}\)。

    • 创建节点\(\nu\)，存储\(x_{mid}\)（以及它对应的点本身），然后分别连接\({\nu}_{left}\)和\({\nu}_{right}\)，使得它们分别是\(\nu\)的左孩子和有孩子，然后令\(\mathcal{T}_{assoc}\)为\(\nu\)的联合数据结构。

  • 返回根节点\(\nu\)

需要注意的是，每个叶子节点不仅仅存储了点的\(y\)坐标，而且也存储了这个点本身，因为最终要报告的是点本身。

引理5.6 平面上一个有\(n\)个点的集合，对应的区域树的空间复杂度是\(O(n\log{n})\)。

Lemma 5.6 A range tree on a set of \(n\) points in the plane requires \(O(n\log{n})\) storage.

证明：

只有通往一个叶子节点\(p\)（\(p\)也是集合\(P\)中的点）路径上的每个节点，才会在对应的联合结构里面存储点\(p\)。所以树\(\mathcal{T}\)中每层上的所有节点中，只会有一个节点对应的联合结构中存储了点\(p\)。而由于一维区域树使用的是线性空间，所以树树\(\mathcal{T}\)中任意深度上一层所有节点，空间复杂度都是\(O(n)\)。整个树的深度是\(O(\log{n})\)，所以二者相乘是\(O(n\log{n})\)。因此整个树\(\mathcal{T}\)的空间复杂度就是\(O(n\log{n})\)。

==（总共有\(O(n\log{n})\)个联合结构存储了\(p\)的示意图在114页，页码是108）==

上述算法的描述，并不能得到构建树的最优时间复杂度是\(O(n\log{n})\)（但实际可以达到）。为了证明这一点，要小心研究。如果以未排序的\(n\)个点的键值来构建二叉搜索树，时间复杂度就已经是\(O(n\log{n})\)了，那么算法中仅每次递归调用中构建二叉搜索树\(\mathcal{T}_{assoc}\)就占用了\(O(n\log{n})\)的时间。但如果\(P\)中的点已经是按照其\(y\)坐标预先排序的，那么按照自底向上的办法构建就只会花费线性时间。所以构建过程中，维护两个列表，分别是基于\(x\)坐标和\(y\)左边的\(P\)中点的排序。这样的话，花费在主树\(\mathcal{T}\)中的每个节点上的时间，就是线性正比于其正则子集的规模。这说明总的构造时间复杂度就和空间复杂度相同，即\(O(n\log{n})\)。由于预排序也花费\(O(n\log{n})\)的时间，所以最终总的构建时间就是\(O(n\log{n})\)。

对于查询算法，首先会选择到\(O(\log{n})\)个正则子集（原因就是沿着通往查询下限\(x\)或上限\(x'\)的路径上有\(O(n\log{n})\)个节点），它们的\(x\)坐标位于\([x:x']\)区间中。然后，再在这些子集中，报告\(y\)坐标位于\([y:y']\)区间中的点。之后，再在每个存储了正则子集的联合结构（也是二叉搜索树）上，应用之前的一维查询算法\(1DRangeQuery\)。唯一不同的是，子函数\(ReportSubtree\)被替换为\(1DRangeQuery\)。

算法：\(2DRangeQuery(\mathcal{T}, [x:x']\times[y:y'])\)

输入：一个二维区域树\(\mathcal{T}\)，一个查询范围\([x:x']\times[y:y']\)。

输出：二维区域树\(\mathcal{T}\)中所有位于查询范围中的点

步骤：

  • 使用（前面提到的）函数\(FindSplitNode(\mathcal{T}, x, x')\)找到分叉节点\({\nu}_{split}\)。

  • 如果分叉节点\({\nu}_{split}\)是叶子节点

    • 检查该节点\({\nu}_{split}\)中存储的对应点是否应该被报告出来

  • 否则

    • （沿着通往\(x\)坐标查询下限的路径上）

    • 令\(\nu\)是分叉节点\({\nu}_{split}\)的左孩子

    • 循环：如果检查\(\nu\)不是叶子节点

      • 如果该节点\(\nu\)对应的\(x\)坐标值\(x_{\nu}\)大于或等于查询下限\(x\)

        • 调用一维查询函数\(1DRangeQuery(\mathcal{T}_{assoc}(rc(\nu)), [y:y'])\)，查询的值是\(y\)坐标

        • 然后令\(\nu\)向左孩子移动：\(\nu \longleftarrow lc(\nu)\)

      • 否则（即该节点\(\nu\)对应的\(x\)坐标值\(x_{\nu}\)严格小于查询下限\(x\)）

        • 然后令\(\nu\)向右孩子移动：\(\nu \longleftarrow rc(\nu)\)

    • （此时\(\nu\)是叶子节点了）检查\(\nu\)对应的点是否应该被报告

    • （沿着通往\(x'\)坐标查询上限的路径上）

    • 令\(\nu\)是分叉节点\({\nu}_{split}\)的右孩子

    • 循环：如果检查\(\nu\)不是叶子节点

      • 如果该节点\(\nu\)对应的\(x\)坐标值\(x_{\nu}\)小于或等于查询下限\(x\)

        • 调用一维查询函数\(1DRangeQuery(\mathcal{T}_{assoc}(rc(\nu)), [y:y'])\)，查询的值是\(y\)坐标

        • 然后令\(\nu\)向右孩子移动：\(\nu \longleftarrow rc(\nu)\)

      • 否则（即该节点\(\nu\)对应的\(x\)坐标值\(x_{\nu}\)严格大于查询上限\(x'\)）

        • 然后令\(\nu\)向左孩子移动：\(\nu \longleftarrow lc(\nu)\)

    • （此时\(\nu\)是叶子节点了）检查\(\nu\)对应的点是否应该被报告

引理5.7 对存有\(n\)个点的二维区域树，按照和坐标轴平行的矩形区域查询，时间复杂度是\(O(\log^2n + k)\)，这里\(k\)是最终报告的点的个数。

A query with an axis-parallel rectangle in a range tree storing \(n\) points takes \(O(\log^2n + k)\) time, where \(k\) is the number of reported points.

证明：

树中每层花费常数时间判断搜索路径的走向，然后有可能调用\(1DRangeQuery\)。根据引理5.2，这个单独的递归调用时间复杂度就是\(O(\log{n} + k_{\nu})\)，\(k_{\nu}\)是对应路径上某个节点的正则子集中所报告的点的个数。所以最终花费的时间是搜索路径上每个节点\(\nu\)花费的时间的总和：

\[ \sum_{\nu}O(\log{n} + k_{\nu}) \]

这里显然有\(\sum_{\nu}k_{\nu}\)就是\(k\)（因为就报告的点的总个数，就等于在那些节点\(\nu\)中报告的点的个数之和）。而通向\(x\)或\(x'\)的路径深度就是\(O(\log{n})\)，所以\(\sum_{\nu}O(\log{n})\)就是\(O(\log{n}) \times O(\log{n})\)，即\(O(\log^2n)\)。

所以最终的查询时间复杂度就是\(O(\log^2n + k)\)。

综上所述，有如下定理。

定理5.8 记\(P\)是二维平面上有\(n\)个点的集合，它对应的区域树，可以用\(O(n\log{n})\)的时间复杂度和\(O(n\log{n})\)的空间复杂度构建出来。以一个矩形区域查询该区域树，时间复杂度是\(O(\log^2n + k)\)，这里\(k\)是最终报告的点的个数。

Theorem 5.8 Let \(P\) be a set of \(n\) points in the plane. A range tree for \(P\) uses \(O(n\log{n})\) storage and can be constructed in \(O(n\log{n})\) time. By querying this range tree one can report the points in \(P\) that lie in a rectangular query range in \(O(\log2 n+k)\) time, where \(k\) is the number of reported points.

实际上在5.6节中，通过一种叫做分散层叠（fractional cascading）的技术，可以把查询时间复杂度降低到\(O(\log{n} + k)\)。

5.4 Higher-Dimensional Range Trees

从二维区域树可以很直观地推广到高维区域树。

令\(P\)是\(d\)维空间上\(n\)个点的集合，即任意一个点的坐标是\(X = (x_1, x_2, \dots, x_d)\)。首先按照每个点的第一个坐标\(x_1\)来构建一棵平衡二叉搜索树。这棵第一次层次的树（即主树）上，任意一个节点\(\nu\)的正则子集\(P(\nu)\)，是由以\(\nu\)为根节点的子树的所有叶子节点构成的。再在这样的每一个节点\(\nu\)上，我们继续构建对应的联合数据结构\(\mathcal{T}_{assoc}(\nu)\)（也就是第二层次的树）。这个\(\mathcal{T}_{assoc}(\nu)\)是一个\((d-1)\)维的区域树，是根据\(\nu\)对应子树的所有叶子节点的后\((d-1)\)个坐标确定的。类似地，这样的一棵区域树，也是通过递归的方式建立起来的，即这是一棵基于第二个坐标\(x_2\)建立的平衡二叉搜索树，每个节点又有指向一棵\((d-2)\)维的区域树，（这棵\((d-2)\)维的区域树）又是由后面\((d-2)\)个坐标确定的。以此类推，直到按照最后一个坐标\(x_n\)建立一个一维的区域树（即平衡二叉搜索树）。

查询算法也是和二维情况下类似。首先在第一层次（first-level）的树上确定\(O(\log{n})\)个节点，这些节点的正则子集包含的点的第一个坐标\(x_1\)都位于\(x_1\)坐标的查询范围内。然后，在这些正则子集上执行对第二层次（second-level）的树区域范围查询。同样地，在每个第二层次（second-level）的树上，确定了\(O(\log{n})\)个节点（亦即正则子集），这就意味着在第二个层次上总共有\(O(\log^2{n})\)个正则子集，它们包含的点的第一个和第二个坐标（\(x_1, x_2\)）都是落于查询范围内的。以此类推，直到最后的一维的区域树（正则子集），它们的叶节点上的点不仅最后一个坐标\(x_d\)落于查询范围内，而且也是最终要报告的点。

定理5.9 记\(P\)是\(d\)维空间上有\(n\)个点的集合，这里\(d \geqslant 2\)。可以以\(O(n\log^{d−1} n)\)的空间复杂度，并以\(O(n\log^{d−1} n)\)的时间复杂度建立一棵对应于\(P\)的区域树。查询\(P\)中位于一个矩形区域的点的时间复杂度是\(O(log^d n+k)\)，这里\(k\)是最后报告的点的个数。

Theorem 5.9 Let \(P\) be a set of \(n\) points in \(d\)-dimensional space, where \(d \geqslant 2\). A range tree for \(P\) uses \(O(n\log^{d−1} n)\) storage and it can be constructed in \(O(n\log^{d−1} n)\) time. One can report the points in \(P\) that lie in a rectangular query range in \(O(log^d n+k)\) time, where \(k\) is the number of reported points.

证明：

记\(T_d(n)\)是在\(d\)维空间上构建\(n\)个点对应的区域树的时间复杂度。根据定理5.8有\(T_2(n) = O(n\log{n})\)。而构建一棵区域树分为两个部分，一个是构建一个平衡二叉搜索树，时间复杂度是\(O(n\log{n})\)，二是构建联合数据结构。在第一层次（first-level）的树的任意深度上，\(P\)中的任意一个点，最多只会在联合数据结构里面存储一次。在某个深度上，构建所有节点对应的联合数据结构的时间复杂度是\(O(T_{d-1}(n))\)，这也是构造对应于根节点的联合数据结构所需的时间复杂度，而且时间复杂度至少是线性的（因为有可能更慢）。因此，总的构建时间是

\[ T_d(n) = O(n\log{n}) + O(\log{n}) \cdot T_{d-1}(n) \]

由于\(T_2(n) = O(n\log{n})\)，所以最终有

\[ T_d(n) = O(n\log^{d−1} n) \]

使用和上面类似的分析，可以得到空间复杂度也是\(O(n\log^{d−1} n)\)。

记\(Q_d(n)\)是在\(d\)维空间上，查询有\(n\)个点的一棵区域树的时间复杂度（不包括报告点的时间）。而查询一棵\(d\)维空间上的区域树的步骤，分为：一、搜索第一层次的树，这花费的时间是\(O(\log{n})\)；二、对对数级的\((d-1)\)维的树进行搜索，因此有

\[ Q_d(n) = O(n\log{n}) + O(\log{n}) \cdot Q_{d-1}(n) \]

而\(Q_2(n) = O(\log^2n)\)，因此有\(Q_d(n) = O(\log^dn)\)，再加上报告点的时间，即\(O(k)\)，那么最终的时间就是\(O(log^d n+k)\)。

而同样地，和二维空间上的情况类似，也可以使用一个对数因子来实现性能的提高。

5.5 General Sets of Points

前面讨论的情况里，都假设每两个点的\(x\)和\(y\)坐标不相同（二维情况下，高维情况类似），但这不现实。为了解决这个问题，引入了一个叫做合成数（composite-number）的概念。这是因为前面的讨论里，并没有假设坐标都是实数（即可以是其他类型），我们只需要它们都是有序的。

所以，使用所谓合成数空间（composite-number space）中的元素来替换每个坐标（坐标原先是实数）。这种合成数元素实际上是一对实数，可以用\((a|b)\)表示，这里\(a\)和\(b\)都是实数。我们使用字典序来定义合成数空间的顺序，对任意两个合成数\((a|b)\)和\((a'|b')\)，如果有\((a|b) < (a'|b')\)，那么它的充要条件是

\[ (a|b) < (a'|b') \Leftrightarrow a < a', or \space (b < b' \space if \space a = a') \]

这样就能允许\(P\)中的任意两个点可以有相同的\(x\)坐标或\(y\)坐标。把每个点的坐标首先做替换，\(p := (p_x, p_y)\)替换为\(\hat{p} := ((p_x | p_y), (p_y | p_x))\)，这样任意两个点的\(x\)或\(y\)坐标即不会相同了。（前提是没有两个点的\(x\)和\(y\)坐标同时相等，即没有重复的点）。这样可以想象一下，如果原先有两个点的\(x\)坐标相同，但它们的\(y\)坐标不同，那么对应的前面提到的替换后的合成数的\(x\)坐标也是不同的。

这样，替换之后就得到一个新的点集\(\hat{P}\)，该点集中任意两个点的\(x\)坐标不同，任意两个点的\(y\)坐标不同。然后将它们排序，再构建对应于\(\hat{P}\)的KD-tree和二维空间上的区域树。

现在假设要查询的矩形范围是\([x:x'] \times [y:y']\)，但我们已经把\(P\)映射到了\(\hat{P}\)上了，因此查询的实际上是合成数空间对应的区域树（或KD-tree）。所以也要把查询范围转换为合成数查询范围。即，

\[ \hat{R} := [(x|-\infty):(x'|+\infty)] \times [(y|-\infty):(y'|+\infty)] \]

剩下需要证明的就是，在\(\hat{R}\)中最终报告出来的点，是否就对应于应该在\(R\)中报告出来的点。下面的引理证明了这一点。

引理5.10 记\(p\)是矩形\(R\)中的一个点，那么\(p\)位于\(R\)中的充要条件就是\(\hat{p}\)位于\(\hat{R}\)，即

\[ p \in R \Leftrightarrow \hat{p} \in \hat{R} \]

证明：

令\(p := (p_x, p_y)\)，\(R := [x,x'] \times [y:y']\)。根据定义，\(p\)位于\(R\)中的充要条件就是\(x \leqslant p_x \leqslant x'\)，并且\(y \leqslant p_y \leqslant y'\)。而又根据合成数的定义以及合成数比较的定义，充要条件可以写作\((x|-\infty) \leqslant (p_x|p_y) \leqslant (x'|+\infty)\)，并且\((y|-\infty) \leqslant (p_y|p_x) \leqslant (y'|+\infty)\)，而这正是\(\hat{p} \in \hat{R}\)的定义，因此命题得证。

实际上，在实际操作过程中，完全没有必要存储一个新的合成数，而是直接用对应的合成数的定义去比较、去排序，而存储的点仍然是原先的点，从它上就可以立即得到对应的合成数。

5.6* Fractional Cascading

前面第5.3小节说明了对一个平面上的区域树进行查询的时间复杂度是\(O(\log^2n + k)\)，本节使用名为分散层叠（fractional cascading）的技术，把查询时间复杂度降低到\(O(\log{n} + k)\)。

（接下来一段是回顾如何查询一个区域树，具体见5.3小节）平面上对应的点集\(P\)的区域树是一个二层次的数据结构。主树是基于点集\(P\)的\(x\)坐标的一棵平衡二叉搜索树，树中每个节点\(\nu\)又有一个联合数据结构\(\mathcal{T}_{assoc}(\nu)\)，它又是一个基于点集\(P\)的一个正则子集的\(y\)坐标的一棵平衡二叉搜索树。查询矩形区域\(R := [x,x'] \times [y:y']\)，就是先选定主树中\(O(\log{n})\)个节点，它们的正则子集中的点的\(x\)坐标位于区间\([x,x']\)，然后再依次查询这\(O(\log{n})\)个节点的\(y\)坐标，使得其位于区间\([y:y']\)。查询每个节点的联合数据结构\(\mathcal{T}_{assoc}(\nu)\)的时间复杂度是\(O(\log{n} + k_{\nu})\)，而总共要查询的联合数据结构的个数（即前面选定的主数中的节点个数）是\(O(\log{n})\)个，所以（二者相乘）总时间复杂度就是\(O(\log^2n + k)\)。

如果能把查询每个联合数据结构的时间复杂度降低到\(O(1 + k_{\nu})\)，那么总的时间复杂度就能降低到\(O(\log{n} + k)\)。问题的关键在于，在一维空间（即一维坐标轴）同样的范围内反复查询，我们利用一个的查询结构来加速其他的查询。

下面通过一个例子，讲述了分散层叠的思路。记\(S_1\)和\(S_2\)是两个集合，集合中的每个元素都是含有一个键（key）的对象（object）。这两个集合的元素分别存储在数组\(A_1\)和\(A_2\)中，并且是有序的。假设我们想要报告\(S_1\)和\(S_2\)中key位于区间\([y:y']\)中的元素对象，就按照下面的办法。首先在\(A_1\)中通过二分查找找到大于等于\(y\)的最小元素，然后从它开始，向右遍历数组，遍历过程中，报告元素，直到遇到第一个大于\(y'\)的元素才停止。使用同样的办法报告\(S_2\)中的元素。如果报告元素的总个数是\(k\)，那么整个的查询时间，就是\(O(k)\)，加上两次二分查找的时间（即在\(A_1\)和\(A_2\)中分别搜索\(y\)）。然而，如果\(S_2\)的key是\(S_1\)的子集的话，我们就可以采用以下的办法，避免第二次（即在\(A_2\)上）的二分查找：即在\(A_1\)中的每个元素上增加一项（一个成员，亦即指针）。如果\(A_1[i]\)存储的元素对象的key是\(y_i\)，那么就在该元素上存储一个指针，这个指针指向\(A_2\)中key大于等于\(y_i\)的最小元素。如果不存在这样的最小元素，那么这个指针就指向空。

==（此处关于\(A_1\)中元素指向\(A_2\)中元素指针的示意图在119页，页码是113）==

如何利用这种数据结构来查询\(S_1\)和\(S_2\)中key位于区间\([y:y']\)中的元素对象？查询\(S_1\)的办法仍然和之前相同，先用二分查找法搜索\(S_1\)对应的\(A_1\)，找到大于等于\(y\)最小元素，然后向右遍历并报告，直到遇到第一个大于\(y'\)的元素停止。不同的是查询\(S_2\)。假设搜索\(A_1\)时，找到大于等于\(y\)最小元素是\(A_1[i]\)，那么因为\(S_2\)是\(S_1\)的子集，所以\(A_1[i]\)指向的\(A_2\)中的元素对应的key值也一定是大于等于\(y\)，然后就从\(A_2\)的那个元素开始向右遍历，并依次报告，直到遇到第一个大于\(y'\)的元素停止。由此就避免了在数组\(A_2\)上的二分查找，并且\(S_2\)中的查询时间就是\(O(1 + k)\)，这里\(k\)是最终报告的元素个数。

根据前面的图，下面通过查询区间\([20:65]\)，举例说明了查询过程。原文在119页（页码是113），此处省略。

下面回到区域树的讨论。这里观察的关键点在于，正则子集\(P(lc(\nu))\)和\(P(rc(\nu))\)都是\(P(\nu)\)的子集，因此根据前面的讨论，就有可能加速查询过程。这里细节更加复杂一点，原因是我们要快速搜索的\(P(\nu)\)的子集是两个，而不是类似前面讨论的一个（即\(A_2\)是\(A_1\)的子集）。

令\(\mathcal{T}\)是平面上有\(n\)个点的集合\(P\)的区域树。\(P(\nu)\)的正则子集存储在对应的联合数据结构中。但和更早之前讨论不同的是，这里我们使用一个排序数组\(A(\nu)\)当做一个联合数据结构，而不是一个二叉搜索树（如5.3节讨论）。这个排序数组\(A(\nu)\)是以每个点的\(y\)坐标大小排序的。进一步，\(A(\nu)\)中的每个元素包含两个指针，一个指向\(A(lc(\nu))\)，另一个指向\(A(rc(\nu))\)。也就是说，这两个指针，一个指向节点\(\nu\)的左孩子的正则子集对应的联合数据结构（这里是一个排序数组），另一个指向节点\(\nu\)的右孩子的正则子集对应的联合数据结构（也是一个排序数组）。

假设\(A(\nu)[i]\)存储一个点\(p\)，那么就有一个指针从\(A(\nu)[i]\)指向\(A(lc(\nu))\)中的一个元素，这个元素存储的点\(p'\)的\(y\)坐标是\(A(lc(\nu))\)中大于等于点\(p\)的\(y\)坐标\(p_y\)的最小值。正如前面所述，因为\(P(lc(\nu))\)是\(P(\nu)\)的子集，所以如果\(p\)的\(y\)坐标是\(P(\nu)\)中大于等于搜索下限\(y\)的最小点，那么同样就有\(p'\)的\(y\)坐标是\(P(lc(\nu))\)中大于等于搜索下限\(y\)的最小点。而类似上面所述，从\(A(\nu)[i]\)指向\(A(lc(\nu))\)中的一个元素的指针，其指针指向的元素一样也是\(P(rc(\nu))\)中大于等于搜索下限\(y\)的最小点。

这种修改版本的区域树叫做层次化区域树。

==（一个层次化区域树示例的示意图在120页，页码是114，两个图分别是层次化区域树，以及对应的联合数据结构，这里是数组）==

下面阐述如何通过搜索一个层次化的区域树来找到对应区域\([x:x'] \times [y:y']\)中的值。

和之前的讨论一样，在主树\(\mathcal{T}\)在搜索\(x\)坐标的上下限\([x:x']\)，以此确定位于区间\([x:x']\)中的\(O(\log{n})\)个节点对应的正则子集（每个点的\(x\)坐标位于搜索区间\([x:x']\)中）。正如前面小节所述，假设搜索上限\(x\)和下限\(x'\)的路径分叉于节点\({\nu}_{split}\)，而这些\(O(\log{n})\)个节点，是位于以下两种路径的节点的左或右孩子，

  • 从\({\nu}_{split}\)向左通向搜索下限\(x\)的路径上，每个节点的右孩子的正则子集

  • 从\({\nu}_{split}\)向右通向搜索上限\(x'\)的路径上，每个节点的左孩子的正则子集

在\({\nu}_{split}\)节点处，搜索对应的联合数据结构（数组）\(A({\nu}_{split})\)，找到大于等于\(y\)的最小节点，这可以在\(O(\log{n})\)时间里完成。当在主树\(\mathcal{T}\)中搜索上下限\(x\)和\(x'\)时，同步记录在联合数据结构中找到的大于等于\(y\)的最小节点。借助对应数组中的指针，这些节点可以以常数时间来维护和更新。

现在假设\(\nu\)是这些\(O(\log{n})\)个节点中的一个，我们要报告的是\(\nu\)对应的数组\(A(\nu)\)中点的\(y\)坐标位于搜索区间\([y:y']\)中的点。因为\(\nu\)是这些\(O(\log{n})\)个节点中的一个，那么它的父节点\(parent(\nu)\)就在通向\(x\)或\(x'\)的路径上，而前面也提到，\(parent(\nu)\)就有指向它的左孩子和右孩子对应的正则子集中大于等于\(y\)的最小节点。在数组\(A(\nu)\)中找到这样的点只需要花费常数时间，然后从它开始向右遍历，并依次报告，直到对应节点的\(y\)坐标大于搜索上限\(y'\)。因此，报告\(A(\nu)\)中节点的\(y\)坐标位于区间\([y:y']\)的时间复杂度就是\(O(1 + k_{\nu})\)，这里\(k_{\nu}\)是在节点\(\nu\)（对应的\(A(\nu)\)）上最后报告的点的个数。所以，最终总的查询时间复杂度就是\(O(\log{n} + k)\)。

对高维区域采样分散层叠计算，同样可以把查询时间复杂度降低一个对数因子。回忆要实现一个\(d\)维区域树的搜索，首先选定第\(d\)个坐标位于搜索区间的节点，以及它们对应的\(O(\log{n})\)个正则子集，然后在每个正则子集上，解决一个\((d-1)\)维的区域搜索。类似地，按照递归的办法，最终得到一个二维区域搜索，就可以使用前面讨论的办法解决，也因此有下面的定理。

定理5.11 \(P\)是\(d\)维空间上有\(n\)个点的集合（\(d \geqslant 2\)）。可以用 \(O(n\log^{d−1}n)\) 的时间复杂度和 \(O(n\log^{d−1}n)\) 的空间复杂度，建立一个\(P\)对应的层次化的区域树。在这个层次化的区域树中查询位于一个矩形区域（超平面）中的点，花费的时间复杂度是\(O(\log^{d−1} + k)\)，这里\(k\)是最终位于该搜索区域中的点的个数。

Theorem 5.11 Let \(P\) be a set of \(n\) points in \(d\)-dimensional space, with \(d \geqslant 2\). A layered range tree for \(P\) uses \(O(n\log^{d−1}n)\) storage and it can be constructed in \(O(n\log^{d−1}n)\) time. With this range tree one can report the points in \(P\) that lie in a rectangular query range in \(O(\log^{d−1} + k)\) time, where \(k\) is the number of reported points.

5.7 Notes and Comments

二十世纪七十年代（1970s)，计算几何发展早期，正交区域搜索是该领域最重要的研究课题之一，有许多人致力于该领域。这产生了许多研究成果，下面讨论其中一些。

在第14章的meshing这部分，讨论了正交区域搜索最早的数据结构，quadtree。不幸的是，在最坏情况下quadtree的表现比较差。1975年Bentley首先提出了KD-tree，它是对quadtree的优化和提升。Samet的书中详细介绍了quadtree，KD-tree以及它们的应用。一些年以后，一些研究者各自独立发现了range tree（区域树）。Lueker和Willard使用分散层叠技术把区域树的查询时间复杂度提升到了\(O(\log{n} + k)\)。而分散层叠技术不仅仅用于区域树，还用于其他对同一个建进行多次查询的场景。Chazelle和Guibas讨论了这种技术的广泛应用。除此之外，分散层叠技术还用于dynamic setting。Chazelle描述了，二维区域搜索最有效的数据结构是分层区域树（layered range tree）的一个变种，他在保持查询时间复杂度是\(O(\log{n} + k)\)的情况下，把空间复杂度降低到了\(O(n\log{n}/\log{\log{n}})\)，并且证明了是最优解。如果搜索区域有一侧是没有边界的，比如\([x:x'] \times [y:+\infty]\)，那么利用第10章讨论的priority search tree，就能达到线性空间的复杂度以及\(O(\log{n})\)的时间复杂度。Chazelle也研究了高维空间里正交区域搜索的最好结果：他给出了用多重对数时间的（polylogarithmic）复杂度和\(O(n(\log{n}/\log{\log{n}})^{d-1})\)的空间复杂度实现对\(d\)维空间的搜索，而且证明了是最优解。当然，空间和时间的平衡调整，也是可能的。

只有基于某些模型的计算情况下，lower-bound的结果才是合理的，这就运行在特定的一些情况下来优化结果。Overmars描述了当点位于\(U \times U\)这样的格子中时，一种更加高效的区域搜索数据结构。这种数据结构根据所允许的预处理时间，最终得到查询时间复杂度是\(O(\log{\log{U}} + k)\)或\(O(\sqrt{U} + k)\)。这些结果使用了Willard早前描述的数据结构。和一般普通的情况相比较，一旦对象的坐标点被限定位于格点上时，在计算几何中许多问题就可以得到更好的时间上下限，例如nearest neighbor searching problem，point location和line segment intersection。

在数据库中，区域搜索（range query）是多重查询中被认为最基本的三种类型之一。其他两种分别是完全匹配查询（exact match queries）和部分匹配查询partial match queries。完全匹配查询是这样一种查询：一个有特定属性（坐标）的对象（点）是否位于数据库中？或诸如此类的查询（such and such）。对于完全匹配查询，显而易见的数据结构就是对坐标使用字典序比较的平衡二叉树，利用它可以得到\(O(\log{n})\)的查询时间复杂度。如果维数增加（即属性个数增加），那么用来表示查询性能的表达式不仅仅使用\(n\)，还要使用维数\(d\)。比如，一个用来做完全匹配查询的二叉树，查询时间复杂度就是\(O(d\log{n})\)，因为它要花费\(O(d)\)的时间来比较两个点，而这个时间复杂度可以很容易地降低到\(O(d + \log{n})\)，而且它是最优的。部分匹配查询只指定了坐标中的一个或几个维度上的值，然后要求查询所有点中和所指定的值相同的点。比如，在平面上，部分匹配查询只指定\(x\)坐标或\(y\)坐标的值。用几何语言描述的话，平面上的部分匹配查询就是查询在位于一条水平或垂直直线上的点。利用\(d\)维的KD-tree，指定\(s (s \lt d)\)个坐标值，部分匹配查询的时间复杂度就是\(O(n^{1-s/d} + k)\)，这里\(k\)是最终报告的点的个数。

在许多应用领域当中，我们的输入数据不是点的集合，而是某些对象（比如多边形）的集合。如果要查询一个对象（多边形）是否完全包含在搜索区域\([x:x'] \times [y:y']\)中，那么就可以把原始查询转化为查询更高维度下的点的查询操作（见Exercise 5.13）。有时候，也要查询一个多边形是否是部分位于搜索区域里，这种问题就是第10章要讨论的截窗问题（windowing problem）。

区域搜索的变种包括当搜索区域的类型是其他的样子的时候，比如圆或三角形。这变种中的许多类型可以划分树（partition trees）解决，它在第16章讨论。

5.8 Exercises

（Omitted）

5.9 References

参考文章：范围搜索 (Range Query)

邓俊辉主页

6 Point Location Knowing Where You Are

引言部分，通过经纬度的说明，引出了点的位置查询（point location query）的概念：给定一个地图以及一个点的坐标，在地图中找到包含该点的区域。这里的地图，就是平面上的细分。