对半平面交算法正确性解释的探索

Piwry

2020-06-12 21:06:58

Algo. & Theory

话说现存的教材对半平面交算法正确性的解释都十分随意...于是无法接受用一个自己都搞不懂的算法的我就捣鼓了这个问题好几天..

本文尝试从半平面交的表面性质入手，并尝试用一个能保证正确性且直观的方式 （为什么说方式而不是构造..? 因为我实在找不到一个好且严谨的构造方法QAQ） 理解算法。

同时文中提出了一种可行的在线做法 （肯定有大佬早就想到类似的了...） 并用同种方式解释了网上一种现较冷门的、常数稍小一些的离线做法。

在都是废话的前言也顺便说明下变量的命名约定：

typedef pair<int, int> pad;/*pad awa*/

const double pi =acos(-1), eps =1e-6;

/*点或者说向量*/
struct vect{
    double x, y;
    vect(){}
    vect(double xx, double yy):x(xx), y(yy){}
    vect operator + (vect v){ return vect(x+v.x, y+v.y); }
    vect operator - (vect v){ return vect(x-v.x, y-v.y); }
    vect operator * (double mu){ return vect(mu*x, mu*y); }
    double operator / (vect v){ return x*v.y-y*v.x; }/*叉积*/
};

/*直线*/
struct line{
    vect u, v;
    double angle;
    line(){}
    line(vect uu, vect vv):u(uu), v(vv){ angle =atan2(vv.y-uu.y, vv.x-uu.x); }
};

line hull[];/*储存 凸包 / 环（定义见下） / 答案*/

line ls[];/*代表加入直线序列*/

int l, r;/*左右指针，左闭右开*/

以及一些函数的声明和实现：

/*(a > b)*/
inline bool gtr(double a, double b){ return (a-b > eps); }

/*(a == b)*/
inline bool eq(double a, double b){ return (a-b < eps && a-b > -eps); }

/*点是否在有向直线右侧，是返回 true，否（在左侧或直线上）返回 false*/
inline bool onright(line f, vect w){ return (gtr((w-f.u)/(f.v-f.u), 0)); }

/*求两条直线交点（无交情况未定义）*/
vect getIntersection(line f, line g){
    double w =((g.u-f.u)/(f.u-f.v))/((f.u-f.v)/(g.u-g.v));
    return g.u+(g.u-g.v)*w;
}

/*用于排序的比较函数*/
int cmp(line A, line B){
    if(eq(A.angle, B.angle)) return onright(B, A.u);/*有向直线最左的会排在最后面，并被保留*/
    else return (gtr(B.angle, A.angle));
}

最后是说明上的一些约定：

  1. 文中的半平面交默认是有向直线左侧交，包括边界。

  2. 实现中所用的计较都是直接使用 atan2 的返回值，并未做其他处理（因此值域为 (-\pi, \pi]）、

  3. 这里的增量法实现都是按直线极角从小到大加入的。

  4. 话说下面以及上面提到的 "增量法"、"标准双端队列做法"，一种典型示范就是刘汝佳书中有关半平面交章节的示例代码。或者关于这种做法的实现细节及说明，可以先看看之前的 计算几何日报（来自 \text{wjyyy} 大佬awa）

Part 1. 引出及一些思考

1.1 存储方式

观察上图中半平面交的形状，有一个很显然的思路就是通过保存交的多边形形状（例如保存顶点）来表示交（其实也可以知道这个多边形总是凸包）；

而对于无穷大的交，我们可以给所有交在 "无穷远"（相对于数据） 处围上一个边界 （话说这方法是不是被提过很多次了...） ，将这种情况也转化为多边形。

无交则只需判断存储的点集是否为空即可。

1.2 一种简易的在线维护思路

现在考虑给已知半平面交添加一条有向直线的情况。

显然新交不包含新加入直线右侧的所有点（在直线上的点保不保留没大区别）。

同时，还发现这些点还满足：

  1. 它们一定是在凸包边界上 "连续" 的一段
  2. 我们将凸包所有边按极角排序，第一个极角大于和第一个极角小于插入直线的极角的两条边的交点一定是需要弹出的点（即使在最坏情况下（只弹出一个点））。

（如下图；对更多点情况可平移直线并由凸包性质发现）

实现时，我们先按 性质二 尝试找到一个需要弹出的点。如果这个点不在加入直线的右侧，则这条直线就不会对交产生影响；否则只需逆/顺时针依次判断即可（这些点是相邻的），并对于一个方向在第一个不满足要求（在加入直线右侧）的点停下。

（至于点在当前直线上...其实这时候不特判去除也不会影响贡献。但在下一部分的做法中，不判这种情况可能会在某些特殊情况下出错）

同时，为了保证精度，我们仍保存加入直线的信息（通常用两个点表示直线...加入后并不重求点为交点），实际的凸包点在判断时求出。而直线删除的条件是它的两个交点都被弹出（一个交点被弹出的话保留，由于相邻的直线（边）被更新了，那个交点下次求出来也是不同的，就相当于弹出了）；同时注意对可能将所有直线弹出的加入操作（即无交）需要特判。（或者干脆存点...这段可以无视（

为了能支持二分查找而快速找到上述两条直线，保存的直线需按极角排序，考虑到还可能要在序列中间插入点，我们可能需要用平衡树来维护。

（由于在竞赛没什么实用价值代码就暂时咕咕了 QAQ （话说半平面交也真的常考吗...））

这种在线做法也可以很方便地改成离线。不过单纯照搬的话其实有很多部分是可以简化一下的。（其实这里应该早就有人看出这种做法和标准增量法很相似了...）

Part 2. 对常用做法的一种解释

2.1 一种较直观的 模型（应该？） 解释

我们先将存下的直线想象成一条由一段段线段组成的 "链"：

或者具体地来讲，是：

（其实也不太严谨QAQ，看看下图理解下吧）

再定义 "环"：

如果这样一个环没有我们不需要的直线，或者说它正好形成了一个凸包，那么这个环就可以直接作为答案了。

可以发现这其实就是 [情况1]。

2.2 对几种环的处理

现在我们考虑下怎么排除 [情况2], [情况3], [情况4] 的无用线，将它们转化成 [情况1]：

首先可以看出对于 [情况2], [情况3]，只要判断线段的交点在不在 队尾、队头 元素的右侧即可。

就像这样：

/*情况2*/
while(l < r-1 && onright(hull[r-1], getIntersection(hull[l], hull[l+1]))) ++l;
/*情况3*/
while(l < r-1 && onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;

而 [情况4] 就比较棘手了。我们先观察队首队尾元素的交点（包括和它们相邻元素的交点）排列情况：

（红点是头尾元素和其相邻元素的交点，E 则是头尾元素的交点（可能在延长线上））

很容易想到共有四种情况。其中三种情况（[A, B], [C, D], [D, B]）我们仍能通过刚才的方法判断。唯独当交点取 [A, C] 时的情况，不论是队首还是队尾元素都判不出来。

但可以发现，如果交点取 [A, C]，此时首尾元素实际上就是有交的（不是延长线上）。这时保留的元素实际上就是这样：

（借用了上面的图一下...）

因此这种情况就可以停止循环了。此时剩下的元素正好就形成了一个凸包。

代码就像这样：

while(l < r-1){
    if(onright(hull[r-1], getIntersection(hull[l], hull[l+1]))) ++l;
    else if(onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;
    else break;/*已经没有更新了*/
}

这种处理方式对 [情况2], [情况3] 也适用。

2.3 离线算法的实现

这里的实现实际上和标准双端队列做法很相似：我们维护一个栈，按极角顺序加入直线；如果栈顶两元素交点在当前直线右侧就弹出栈顶（因为是求链，这里没有对队头的判断）。退栈至少保留一个元素（避免越界问题）。

这其中，弹出是为了维护链（当前直线和栈顶的交点用来拓展链是不满足要求的，因为极角跨度太大；否则相反。这可通过加入直线的极角单调性得出，也可见下图）；而弹出栈顶却不舍弃加入直线则是为了去除无用直线。

从半平面交的意义理解，可以知道这样的处理是 “最优的”（"外围" 的直线对最终的交没有贡献）。

其实这部分的逻辑就像在下图中找出最中间那一圈链：

参考代码：

sort(ls, ls+totl, cmp);
int l =0, r =0;
for(int i =0; i < totl; ++i){
    while(i < totl-1 && eq(ls[i].angle, ls[i+1].angle)) ++i;/*去重*/
    /*这里后面会解释（交点在当前直线上的情况）*/
    while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;hull[r++] =ls[i];
}

而获得链之后，我们只需向上一节那样处理即可（总的代码可见后文）。

2.4 对无交的判断

在上一节的做法中，我们在现存的直线不大于一条时就不判断而是直接压入当前直线了。

如果是有交还好说，但如果实际上没交的话要怎么判断呢？（最后总会保留不少于两条直线）

先很容易想到一种特殊的无交情况（这种情况留到下面可能会比较麻烦..可能会导致求交点时出错）：有两根方向相反且平行的直线。

首先可以考虑在加入边时弹出元素的 while 循环判断：

while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))){
    --r;
    if(r-l > 1 && eq(ls[i].angle-hull[r-1].angle, pi)) return pad(0, 0);/*返回空*/
}
if(r-l > 1 && eq(ls[i].angle-hull[r-1].angle, pi)) return pad(0, 0);/*外面的判断避免开始几个元素不进入 while 循环*/
hull[r++] =ls[i];

但我们其实可以发现，如果有这样的一组直线，那么在加入其中第二根时它一定会趋于将栈中的元素弹完。

实际上它们的位置关系是这样的：

因为凸包的单调性，其他的交点一定要比这组直线先加入的那条线要 "远"，因此一定满足弹出的条件。

还有人可能想到通过延长边来构造点使其到达当前线的左侧：

不过这种情况下，延长的那条线极角一定要比当前线要大(可以看出这是向直线正方向延长并找到的一个点，如果极角小于当前线那么正方向是在一直远离当前线的)，是比当前线后加入的；因而这种情况自然也就不存在了。

因此，在这过程中也会将和它平行的直线弹出，除非那条直线处于栈底。

所以实际上我们只需紧贴着在入栈的前面判一次就行了：

（一般的无交可以用下面的通用方法判断。这里只是为了避免在求栈顶交点时带入平行线）

while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;
if(r-l > 1 && eq(ls[i].angle-hull[r-1].angle, pi)) return pad(0, 0);
hull[r++] =ls[i];

接下来考虑对更一般的无交情况的判断。

首先在每次开始前，我们先将 4 条限制边界的直线也放入序列。这样我们就可以保证有交的情况中最后栈中至少剩下三条线（这也方便了我们处理无穷交的情况）

而对于无交的情况，实际上最后求出来是这样的一个条链：

（除了中间两个元素（\overrightarrow{AB}, \overrightarrow{BC}），其他元素的情况可能因具体实现方式有差异）

（且从半平面交的意义考虑，我们去除的都是对最终的交没有贡献的直线，因此在无交时总不会忽略这种情况。）

即（相邻两个元素极角差超过 \pi 导致的）整个链不成环。

同时这里还可以发现用 2.2 的方式处理这种链后，队列中仅会剩下中间的两个元素（\overrightarrow{AB}, \overrightarrow{BC}）。

因此只要最后判断若栈中只剩下两个元素就返回无交即可。

 

同时，我们注意交退化成直线或点的情况，可以想到在弹栈时如果在栈顶交点在当前直线上时也弹出，就会把这种情况判成无交。因此对于交点恰好在加入直线上的情况，我们不再弹出。

2.5 对常见做法的重解释

现在再回来看一般的双端队列做法。

添加直线时的实现和上面讲到的几乎完全一致。

同时因为我们在添加直线时还顺便判了队头，最终得出的链应该为 [情况3]（或直接求得答案）。因而最后我们只需简单地对队尾弹出即可。

示例代码：

/*这里求得的交包括边界，主要取决于 onright() 函数*/
inline pad getHPI(line ls[], int totl, line hull[]){
    sort(ls, ls+totl, cmp);
    int l =0, r =0;
    for(int i =0; i < totl; ++i){
        while(i < totl-1 && eq(ls[i].angle, ls[i+1].angle)) ++i;/*去重*/
        while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;
        while(r-l > 1 && onright(ls[i], getIntersection(hull[l], hull[l+1]))) ++l;
        if(r-l > 0/*首次循环可能为空*/ && eq(ls[i].angle-hull[r-1].angle, pi)) return pad(0, 0);/*判平行的无交*/
        hull[r++] =ls[i];
    }
    while(r-l > 2 && onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;/*仅弹出队尾即可*/
    if(r-l < 3) return pad(0, 0);
    else return pad(l, r);/*返回答案所在区间*/
}

（其实也可以想想会使先判队头做法出问题的情况转化成链是什么样子，不过本质也是特判这里就不展开了qaq）

2.6 另一种冷门的离线做法

这种做法最初是我在（由洛谷模板一篇题解导向的）一本通里发现的。

原书的半平面交章节仅仅是整理了已有论文的解释。但其实论文中的 "排序增量算法" 仅仅讲了思路而没有实现（其实那篇论文（2006 《朱泽园--半平面交的新算法及其实用价值》）中提到的算法实现起来也挺坑...可以看 这里 的下面第二个回复），且下文的样例程序也和上面完全没有关系...

具体地讲，书中的算法是将标准双端队列做法在加入直线时对队头的判断完全去掉，并在加入完直线后用 2.2 对 [情况4] 的处理方式处理。

可以想到，如果在加入直线时不判队头，最后求得的环可能是 [情况2], [情况3], [情况4] 的任意一种。因而只需用 2.2 中讲到的对 [情况4] 的处理方式就可以得到答案了。

示例代码：

/*这里求得的交包括边界，主要取决于 onright() 函数*/
inline pad getHPI(line ls[], int totl, line hull[]){
    sort(ls, ls+totl, cmp);
    int l =0, r =0;
    for(int i =0; i < totl; ++i){
        while(i < totl-1 && eq(ls[i].angle, ls[i+1].angle)) ++i;/*去重*/
        while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;
        if(r > 0/*首次循环可能为空*/ && eq(ls[i].angle-hull[r-1].angle, pi)) return pad(0, 0);/*判平行的无交*/
        hull[r++] =ls[i];
    }
    while(r-l > 1){
        if(onright(hull[r-1], getIntersection(hull[l], hull[l+1]))) ++l;
        else if(onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;
        else break;/*已经没有更新了*/
    }
    if(r-l < 3) return pad(0, 0);
    else return pad(l, r);/*返回答案所在区间*/
}

这种写法的常数应该也和原来的写法差不多...要说有什么优势的话..或许更好记？（不用背 while 位置？）

尾声

从开始写这篇文章前，我曾试图找到一个能准确描述半平面交的模型，同时还要易于理解，但最终我思考了许久也没得出结果，因此只好转而使用仅便于理解的描述方式（就是文中的那种）。

说实话这样的结果让我很沮丧，不过即使这样可能也会比完全不对正确性做解释要稍微好一点（也是为了记录我对半平面交的思考），于是我就硬着头皮把这篇文章写完了。

仅仅是形象却难以用数学语言定义的说法导致这篇文章中许多用到数学语言的地方令我感到很别扭，在修订时我时不时有想删去这些段落的想法；不过想了想删去反而可能导致讲不清 “这到底是个什么东西”，于是最后还是保留了XD（

如果还有时间可能还会对文章可能有的表意不明的地方再做修订，不过相信读者自己手玩一下还是能搞懂描述的核心意思的qwq。

（说实话我写这篇文章的目的就是为了让更多人在敲（半平面交）代码时能明白每一行代码的作用，如果不写某行代码可能导致的后果；而不是仅仅背住部分特殊数据，背记着模板敲出代码（不过也因人而异啦XD））