计算几何

· · 个人记录

1.向量

向量的定义

在二维及以上维度里,既有大小又有方向的量叫做向量

向量可以用有向线段表示

如图,从A到B的有向线段表示了一个向量,记作 \overrightarrow{AB}\vec{a} ,读作向量 AB 或向量 a

注意: 向量可以用有向线段表示,但向量并不等同于有向线段,所有方向相同且长度相等的有向线段都可以表示同一个向量。

正因如此,我们可以将有向线段的起点平移至平面直角坐标系的原点,并用其终点的坐标来表示向量。

如图,\overrightarrow{OA} 可以表示为 (2,3)

因为坐标系中的点坐标可以表示向量,所以我们可以将坐标系中的点 A 看作是向量 \overrightarrow{OA},从而就可以让坐标系中的点直接参与向量运算。

向量的模

向量的长度叫做向量的模,记作 | \overrightarrow{AB} |

如上图,| \overrightarrow{OA} | = \sqrt{13}

零向量

长度为 0 的向量叫做零向量,记作 0

零向量的起点和终点重合,故零向量的方向不确定,也就是说零向量的方向是任意的,但零向量的大小是确定的。

相反向量

若两个向量大小相等,方向相反,则称这两个向量互为相反向量,记作 \vec{a}- \vec{a}

特别的,我们规定,零向量的相反向量为它本身。

表示成坐标形式就是 \vec{a} = (x,y) , -\vec{a} = (-x,-y)

共线向量

方向相同或相反的非零向量叫做共线向量(或平行向量),记作 \vec{a} /\mskip-4mu/ \vec{b}

由于零向量的方向不确定,我们规定,零向量与任意向量平行,即 \forall \, \vec{a},0 /\mskip-4mu/ \vec{a}

2.向量的运算

以下默认 \vec{a} = (x_1,y_1) , \vec{b} = (x_2,y_2)

向量的加法

向量既有大小,又有方向,其加减运算不能简单的将大小相加。

向量的加法满足平行四边形定则,即以参与运算的两个向量作为邻边作平行四边形,这两个向量的和即为该平行四边形的对角线,如下图:

因为平行四边形对边平行且相等,所以我们可以将向量的加法理解为将 \vec{b} 的起点平移至 \vec{a} 的终点,则由 \vec{a} 的起点指向 \vec{b} 的终点的有向线段表示的向量就是 \vec{a} + \vec{b} ,如下图:

表示成坐标形式就是

\vec{a} + \vec{b} = (x_1 + x_2 , y_1 + y_2)

向量的加法满足交换律和结合律,即:

\vec{a} + \vec{b} = \vec{b} + \vec{a} (\vec{a} + \vec{b}) + \vec{c} = \vec{a} + (\vec{b} + \vec{c})

加法代码:

Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}

向量的减法

与实数运算中减去一个数等于加上这个数的相反数类似,减去一个向量等于加上这个向量的相反向量,即 \vec{a} - \vec{b} = \vec{a} + (-\vec{b})

由此我们可以将减法转化成加法。

从图中我们可以看出,\overrightarrow{OC} = \overrightarrow{AB}

所以 \vec{b} - \vec{a} 的几何意义就是从 \vec{a} 的终点指向 \vec{b} 的终点的向量。

根据这个几何意义,\overrightarrow{AB} 可以这样计算:

\overrightarrow{AB} = B - A

减法代码:

Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}

向量的数乘

向量的数乘即向量与一个实数相乘,其实质是向量的缩放。

实数 k 与向量 \vec{a} 的数乘记作 k\vec{a} ,结果是一个向量,其模为 | k | \cdot | \vec{a} |

k>0 时, k\vec{a}\vec{a} 同向
k<0 时, k\vec{a}\vec{a} 反向
k=0 时, k\vec{a} = 0

表示成坐标形式就是

k\vec{a} = (k x_1,k y_1)

数乘满足以下运算律:

\lambda,\mu \in \mathbb{R}

结合律:\lambda(\mu \vec{a}) = (\lambda \mu)\vec{a} = \mu(\lambda \vec{a})

向量对于数的分配律(第一分配律):(\lambda + \mu)\vec{a} = \lambda \vec{a} + \mu \vec{a}

数对于向量的分配律(第二分配律):\lambda(\vec{a} + \vec{b}) = \lambda \vec{a} + \lambda \vec{b}

数乘代码:

Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}

向量的数量积

又称点积、内积,记作 \vec{a} \cdot \vec{b} ,结果是一个数量(没有方向)

定义: \vec{a} \cdot \vec{b} = | \vec{a} | | \vec{b} | \cos\theta\theta\vec{a}\vec{b} 的夹角,且 0 \le \theta \le \pi

用坐标表示为: \vec{a} \cdot \vec{b} = x_1 x_2 + y_1 y_2

坐标表示的推导:

如图,A(x_1,y_1) , B(x_2,y_2)

根据余弦定理,| \vec{a} | ^2 + | \vec{b} | ^2 - 2 | \vec{a} | | \vec{b} | \cos\theta = |AB| ^2

{x_1}^2 + {y_1}^2 + {x_2}^2 + {y_2}^2 - 2 | \vec{a} | | \vec{b} | \cos\theta = {x_1}^2 - 2 x_1 x_2 + {x_2}^2 + {y_1}^2 - 2 y_1 y_2 + {y_2}^2

\therefore | \vec{a} | | \vec{b} | \cos\theta = x_1 x_2 + y_1 y_2 \because \vec{a} \cdot \vec{b} = | \vec{a} | | \vec{b} | \cos\theta \therefore \vec{a} \cdot \vec{b} = x_1 x_2 + y_1 y_2

点积满足以下运算律:

\lambda \in \mathbb{R}

交换律:\vec{a} \cdot \vec{b} = \vec{b} \cdot \vec{a}

关于数乘的结合律: (\lambda \vec{a}) \cdot \vec{b} = \lambda(\vec{a} \cdot \vec{b}) = \vec{a} \cdot (\lambda \vec{b})

分配律:\vec{a} \cdot (\vec{b} + \vec{c}) = \vec{a} \cdot \vec{b} + \vec{a} \cdot \vec{c}

点积可用于判定两向量是否垂直: \vec{a} \perp \vec{b} \Leftrightarrow \vec{a} \cdot \vec{b} = 0

点积代码:

double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}

向量的向量积

又称叉积、外积,记作 \vec{a} \times \vec{b} ,结果为向量

叉积的模长: | \vec{a} \times \vec{b} | = | \vec{a} | | \vec{b} | \sin\theta (即两向量作为邻边围成的平行四边形的面积)

叉积的方向垂直于参与运算的两向量所在的平面,且方向与参与运算的两向量的相对位置有关。具体来说:

如图,对于 \vec{a} \times \vec{b} ,若 \vec{b}\vec{a} 的逆时针方向(即以 \vec{a} 所在的射线作为始边,以 \vec{b} 所在的射线作为终边的角 \theta 满足 0 < \theta < \pi),则叉积的方向垂直于平面向读者方向。

若若 \vec{b}\vec{a} 的顺时针方向(即以 \vec{a} 所在的射线作为始边,以 \vec{b} 所在的射线作为终边的角 \theta 满足 -\pi < \theta < 0)反之。

\vec{a} /\mskip-4mu/ \vec{b} ,则 \vec{a} \times \vec{b} = 0

向量的叉积同样可以用坐标表示:\vec{a} \times \vec{b} = x_1 y_2 - x_2 y_1

用坐标求出的值若为正数,则 \vec{b}\vec{a} 的逆时针方向

若为负数,则 \vec{b}\vec{a} 的顺时针方向

若为 0,则 \vec{a} /\mskip-4mu/ \vec{b}

利用这个性质我们可以判断两个向量的位置关系。

向量的叉积满足以下运算律:

\lambda \in \mathbb{R}

\vec{a} \times \vec{b} = (-\vec{b}) \times \vec{a} = \vec{b} \times (-\vec{a}) = (-\vec{a}) \times (-\vec{b})

关于数乘的结合律:(\lambda \vec{a}) \times \vec{b} = \lambda(\vec{a} \times \vec{b}) = \vec{a} \times (\lambda \vec{b})

左分配律:\vec{a} \times (\vec{b} + \vec{c}) = \vec{a} \times \vec{b} + \vec{a} \times \vec{c}

右分配律:(\vec{b} + \vec{c}) \times \vec{a} = \vec{b} \times \vec{a} + \vec{c} \times \vec{a}

叉积代码:

double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}

3.求直线交点

如图,求直线 ABCD 的交点

S_{\text{平行四边形ACDE}} = |\overrightarrow{AC} \times \overrightarrow{CD}| , S_{\text{平行四边形ABGE}} = |\overrightarrow{AB} \times \overrightarrow{CD}| \because \triangle AFH \backsim \triangle ABI \therefore \frac{AF}{AB} = \frac{FH}{BI} = \frac{S_{\text{平行四边形ACDE}}}{S_{\text{平行四边形ABGE}}} = \frac{|\overrightarrow{AC} \times \overrightarrow{CD}|}{|\overrightarrow{AB} \times \overrightarrow{CD}|} \therefore F = A + \overrightarrow{AF} = A + \frac{AF}{AB} \overrightarrow{AB} = A + \frac{|\overrightarrow{AC} \times \overrightarrow{CD}|}{|\overrightarrow{AB} \times \overrightarrow{CD}|} \overrightarrow{AB} = A + \frac{|\overrightarrow{AC} \times \overrightarrow{CD}|}{|\overrightarrow{AB} \times \overrightarrow{CD}|} (B-A)

代码:

Vec getcrp(Line A,Line B) //crp=cross point 交点
{
    double k=cross(B.a-A.a,B.b-B.a)/cross(A.b-A.a,B.b-B.a); //A.a和A.b表示直线A的两个端点,直线B同理
    return A.a+(A.b-A.a)*k;
}

4.求多边形面积

如图,对于平面上的一个多边形,我们可以从某个点向其他所有顶点连线,将其分割成若干个三角形,三角形面积的和即为多边形面积。

那么三角形面积怎么求呢?用向量。

如图,从某一点向其他所有点连向量,因为向量的叉积可以求两向量作为邻边围成的平行四边形的面积,而三角形的面积是平行四边形的面积的一半,这样我们就可以用叉积求出多边形的面积。

那么可能有人会问了:如果多边形不是凸多边形怎么办呢?

确实,如果多边形非凸,我们会算上一些多余部分。

但是我们可以按顺序依次算叉积,即计算 \overrightarrow{AB} \times \overrightarrow{AC} , \overrightarrow{AC} \times \overrightarrow{AD} , \overrightarrow{AD} \times \overrightarrow{AE} …

因为用坐标计算的叉积,若第二个向量在第一个向量的顺时针方向,算出来的结果是负数,那么加上这个负数就相当于减去的多余部分。

如图,当我们加上 \overrightarrow{AC} \times \overrightarrow{AD} 后,就相当于减去了红色部分,也就是减去了多余部分。

因此,只要我们依次计算叉积,就可以保证求出来的面积就是多边形的面积,也就可以保证这样求面积的方法是普遍适用的。(当然,很少有题目会让你求非凸多边形的面积)

5.判断点与线、线与线的位置关系

直线的左半平面与右半平面

在计算几何中,直线由有顺序的两个点确定,如下图为直线AB

那么我们定义,直线AB的左半平面为以射线AB为始边,射线AC为终边且 0 < \angle BAC < \pi 的所有射线AC的集合,即下图中阴影部分。

相似的,直线AB的右半平面的定义为以射线AB为始边,射线AC为终边且 -\pi < \angle BAC < 0 的所有射线AC的集合,即上图中白色部分。

特别的,直线AB上的点不属于任何一个半平面。

点与线的位置关系(向量与向量的位置关系)

判断点与线的位置关系,其实就是判断点在直线的左半平面,直线上还是右半平面。

如图,判断点C与直线AB的位置关系。

前面在讲叉积时说过,用坐标算出来的叉积结果可以判断参与运算的两个向量之间的位置关系,那么我们可以用同样的办法判断点与直线的位置关系。

计算 \overrightarrow{AB} \times \overrightarrow{AC}

若为正数,则 \overrightarrow{AC}\overrightarrow{AB} 的逆时针方向,即C在直线AB的左半平面。

若为负数,则 \overrightarrow{AC}\overrightarrow{AB} 的顺时针方向,即C在直线AB的右半平面。

若为 0,则 \overrightarrow{AC} /\mskip-4mu/ \overrightarrow{AB},即C在直线AB上。

线与线的位置关系

将直线上的两个点表示成向量,用向量的运算判断。

平行

\vec{a} \times \vec{b} = 0 \Leftrightarrow a /\mskip-4mu/ b

相交

两条直线不是平行就是相交

垂直

\vec{a} \cdot \vec{b} = 0 \Leftrightarrow a \perp b

6.二维凸包

定义:能包含平面上给定的所有点的最小凸多边形。

性质:凸包上的边按逆时针顺序排列,第 i+1 条边一定在第 i 条边的逆时针方向

求法:Graham扫描法

首先选取y坐标最小的点(若有多个y坐标相等的点则选取其中x坐标最小的点),将该点作为极点建立极坐标系,将平面上所有点按极角从小到大排序。

由于凸包第 i+1 条边一定在第 i 条边的逆时针方向,我们可以用单调栈维护凸包。

怎么判断新加入的点是否满足性质呢?

我们可以把加入栈的点之间依次连接向量,将栈顶的点与新点连接向量,然后判断两个向量之间的位置关系。

如图,A,B 为已入栈的点,C 为新点。

那么我们判断 \overrightarrow{AB}\overrightarrow{BC} 的位置关系

\overrightarrow{BC}\overrightarrow{AB} 的顺时针方向或与 \overrightarrow{AB} 共线,则边AB不满足条件,将B弹出栈。

否则边AB满足条件,停止循环。

最后在将栈顶和栈底的点连边,就可以构成一个完整的凸包。

代码(Luogu P2742):

#include <bits/stdc++.h>
#define maxn 1000005
using namespace std;
struct Vec{double x,y;}v[maxn];
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double angle(Vec a){return atan2(a.y,a.x);}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec O;//极点
bool cmp(Vec a,Vec b)
{
    double k1=angle(a-O);
    double k2=angle(b-O);
    if (k1!=k2) return k1<k2;
    if (a.x!=b.x) return a.x<b.x;
    return a.y<b.y; 
}
int n;
int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++)
    {
        scanf("%lf%lf",&v[i].x,&v[i].y);
        if (i==1) O=v[i];
        if (v[i].y<O.y||(v[i].y==O.y&&v[i].x<O.x)) O=v[i];
    }
    sort(v+1,v+1+n,cmp);
    Vec s[maxn];
    int top=0;
    for (int i=1;i<=n;i++)
    {
        while(top>1&&cross(s[top]-s[top-1],v[i]-s[top])<=0) top--;
        //top>1是因为若栈中只剩一个点是无法判断的,此时新点必须入栈才能进一步判断之后的点
        s[++top]=v[i];
    }
    double ans=0.0;
    for (int i=1;i<=top;i++)
    {
        int tar=i%top+1;//逆时针方向的下一个点
        ans+=dist(s[tar]-s[i]);
    }
    printf("%.2lf",ans);
    return 0;
} 

7.判断点是否在凸包内

如图,判断点P是否在凸包内。

将凸包上的边按逆时针顺序连接向量,并将凸包所有节点与P连接向量。

P在凸包内,则有 \overrightarrow{P_iP_{i+1}} \times \overrightarrow{P_iP} \ge 0 ,即P在每一条边的逆时针方向或边上。

否则,P不在凸包内。

如图,判断点 A 在凸包内,凸包中的点按逆时针顺序排列,编号为 0n-1

首先将点 A 在阴影部分的情况判掉。

然后,令 l=1 , r=n-1

接下来,我们判断点 A\overrightarrow{P_0P_{mid}} 的左侧还是右侧,若在左侧,令 l=mid ,否则 r=mid

最后,P_lP_r 会变成相邻的两个点,如图:

此时点 A 在射线 OP_lOP_r 之间,此时只需要判断点 A 是否在 \overrightarrow{P_lP_r} 的左侧或 \overrightarrow{P_lP_r} 上,是则点在凸包内,否则点在凸包外。

Code:

//d为需判断的点,E为凸包顶点坐标
//1为在凸包内,0反之
if (cross(E[N-1]-E[0],d-E[0])>0){puts("0");continue;}
if (cross(E[1]-E[0],d-E[0])<0){puts("0");continue;}
int l=1,r=N-1;
while(l+1<r)
{
    int mid=(l+r)>>1;
    if (cross(E[mid]-E[0],d-E[0])>0) l=mid;
    else r=mid;
}
if (cross(E[r]-E[l],d-E[l])>=0) puts("1");
else puts("0");

8.闵可夫斯基和

定义:两个多边形A,B的闵可夫斯基和定义为 \{a+b \ | \ a \in A,b \in B\}

其中 a,b 指终点坐标对应的点在多边形内的向量 \overrightarrow{OP}

此处我们只讨论凸包的闵可夫斯基和。

其实就是将其中一个多边形的某个顶点放在另一个多边形上绕一圈,如图:

红色凸包与蓝色凸包的闵可夫斯基和即为图中整个多边形。

由图中可看出,闵可夫斯基和对应的边就是将两个凸包的边按极角排序后一次收尾相接,根据这个性质我们可以求出闵可夫斯基和的形状。

那么闵可夫斯基和的位置怎么求呢?

因为闵可夫斯基和的定义是向量相加,所以闵可夫斯基和y坐标最小的点即为两个凸包y坐标最小的点的和(具体看代码)

最后记得将相邻的且极角相等的边合并,不然会有奇怪的错误

例题:P4557 [JSOI2018] 战争

给出两个凸包A,B,求给出的向量p是否满足 \forall a \in A\text{且}b \in B , a \ne b+p

我们考虑将这个式子移一下项,可得:p \ne a-b

即:p \ne a+(-b)

我们发现,原条件等价于p \notin \{a+(-b) \ | \ a \in A,b \in B\}

回看一下闵可夫斯基和的定义,我们发现,就是将B的坐标都乘以 -1 后与 A 求闵可夫斯基和,然后判断p是否在闵可夫斯基和内。

Code:

#include <bits/stdc++.h>
#define maxn 100005
using namespace std;
const double pi=acos(-1);//圆周率 
struct Vec{double x,y;};//向量结构体 
bool operator == (Vec a,Vec b){return a.x==b.x&&a.y==b.y;}
bool operator != (Vec a,Vec b){return !(a==b);}
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double angle(Vec a){return atan2(a.y,a.x);}//极角 
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}//叉积 
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec O;
bool cmp(Vec a,Vec b)//求凸包时的排序 
{
    double k1=angle(a-O);
    double k2=angle(b-O);
    if (k1!=k2) return k1<k2;
    if (a.x!=b.x) return a.x<b.x;
    return a.y<b.y;
}
void hull(Vec *v,int n,Vec *e,int &num)//求凸包 
{
    sort(v+1,v+1+n,cmp);
    static Vec st[maxn];
    int top=0;
    for (int i=1;i<=n;i++)
    {
        while(top>1&&cross(st[top]-st[top-1],v[i]-st[top])<=0) top--;
        st[++top]=v[i];
    }
    for (int i=1;i<=top;i++)
    {
        int tar=(i%top)+1;
        e[++num]=st[tar]-st[i];//记录凸包的边 
    }
}
int n,m,q;
Vec v1[maxn],v2[maxn];
Vec e[maxn<<1];//存储两个凸包的边,记得开两倍 
Vec E[maxn<<1];
int num,N;
bool cmp2(Vec a,Vec b)
{
    double k1=angle(a);
    double k2=angle(b);
    //因为atan2函数求出的极角是-π到π之间的,而排序时要按照极角在0到2π之间的值排序,所以将小于0的极角加上2π 
    if (k1<0) k1+=pi*2;
    if (k2<0) k2+=pi*2;
    return k1<k2;
}
int main()
{
    scanf("%d%d%d",&n,&m,&q);
    for (int i=1;i<=n;i++)
    {
        scanf("%lf%lf",&v1[i].x,&v1[i].y);
        if (i==1) O=v1[i];
        if (v1[i].y<O.y||(v1[i].y==O.y&&v1[i].x<O.x)) O=v1[i];
    }
    hull(v1,n,e,num);
    for (int i=1;i<=m;i++)
    {
        scanf("%lf%lf",&v2[i].x,&v2[i].y);
        v2[i].x*=-1;v2[i].y*=-1;
        if (i==1) O=v2[i];
        if (v2[i].y<O.y||(v2[i].y==O.y&&v2[i].x<O.x)) O=v2[i];
    }
    hull(v2,m,e,num);
    sort(e+1,e+1+num,cmp2);//按极角排序 
    for (int i=1;i<=num;i++)
    {
        if (i==1||angle(e[i])!=angle(e[i-1])) E[++N]=e[i];
        else E[N]=E[N]+e[i];//合并极角相等的边 
    }
    E[0]=v1[1]+v2[1];//闵可夫斯基和的y坐标最小的点 
    for (int i=1;i<=N;i++) E[i]=E[i-1]+E[i];//求出各顶点的坐标 
    for (int i=1;i<=q;i++)
    {
        Vec d;
        scanf("%lf%lf",&d.x,&d.y);
        if (cross(E[N-1]-E[0],d-E[0])>0){puts("0");continue;}
        if (cross(E[1]-E[0],d-E[0])<0){puts("0");continue;}
        int l=1,r=N-1;
        while(l+1<r)
        {
            int mid=(l+r)>>1;
            if (cross(E[mid]-E[0],d-E[0])>0) l=mid;
            else r=mid;
        }
        if (cross(E[r]-E[l],d-E[l])>=0) puts("1");
        else puts("0");
    }
    return 0;
}

练习:CF87E Mogohu-Rea Idol

题解link

9.旋转卡壳

一种重要算法,可用于求凸包直径,也就是求凸包上任意两点间距离的最大值。

做法:先求出凸包,然后逆时针枚举凸包上的边,对于每条边找到凸包上距离该边最远的点,设为P_k,用该边两个端点与P_k的距离更新答案。

由于逆时针旋转后,新的P_k^\prime一定在原来的P_k的逆时针方向,所以我们不用每次都重新枚举一遍点,直接从上次的P_k开始逆时针枚举即可。

如何判断P_kP_{k+1}哪个更优?

如图,AB 为当前边。

要比较P_kP_{k+1}哪个更优,实质就是比较 P_kEP_{k+1}D 哪个更长。

因为 \triangle ABP_k\triangle ABP_{k+1} 同底,我们可以用比较 S_{\triangle ABP_k}S_{\triangle ABP_{k+1}} 来比较 P_kEP_{k+1}D

面积怎么求?用叉积算出平行四边形面积在除以二即可(其实不除也行,因为比较面积的两倍跟比较面积是等价的)

代码(Luogu P1452):

#include <bits/stdc++.h>
#define maxn 1000005
using namespace std;
struct Vec{double x,y;}v[maxn];
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double angle(Vec a){return atan2(a.y,a.x);}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec O;
bool cmp(Vec a,Vec b)
{
    double k1=angle(a-O);
    double k2=angle(b-O);
    if (k1!=k2) return k1<k2;
    if (a.x!=b.x) return a.x<b.x;
    return a.y<b.y; 
}
int n;
int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++)
    {
        scanf("%lf%lf",&v[i].x,&v[i].y);
        if (i==1) O=v[i];
        if (v[i].y<O.y||(v[i].y==O.y&&v[i].x<O.x)) O=v[i];
    }
    sort(v+1,v+1+n,cmp);
    Vec s[maxn];
    int top=0;
    for (int i=1;i<=n;i++)
    {
        while(top>1&&cross(s[top]-s[top-1],v[i]-s[top])<=0) top--;
        s[++top]=v[i];
    }
    double ans=0.0;
    int k=2;
    for (int i=1;i<=top;i++)
    {
        int tar=i%top+1;
        while(cross(s[tar]-s[i],s[k%top+1]-s[i])>cross(s[tar]-s[i],s[k]-s[i])) k=k%top+1;
        ans=max(ans,sqrdist(s[k]-s[i]));
        ans=max(ans,sqrdist(s[k]-s[tar]));
    } 
    printf("%.0lf",ans);
    return 0;
} 

10.半平面交

求若干个半平面的交集,这里的有效半平面定义为直线的左半平面。

求法:

首先将所有的直线按极角排序,极角相同的只保留最左侧的那条直线。

显然最后的交是一个凸多边形,我们使用deque(双端队列)来维护。

考虑新加入一个半平面(一条直线)

(注:图中直线AB(蓝色)为deque中倒数第二条直线,直线CD(绿色)为deque尾端的直线,直线EF(红色)为新加入的直线,以下除特殊说明外同理)

如图,若deque中最后两条直线的交点在新加入的直线的右侧,显然deque尾端的直线对于半平面的交已经没有了影响,可以弹出。

如图,若deque中最后两条直线的交点在新加入的直线的左侧,那么此时deque尾端的直线显然都还对半平面的交有影响,不必弹出,此时直接将新的直线入队。

这样看起来,我们用单调栈貌似也能实现目的,为什么我们要用deque呢?

(注:图中直线AB(蓝色)为deque头端的直线,直线CD(绿色)为deque众第二条直线,直线EF(红色)为新加入的直线)

如图,deque中的前两条直线的交点也有可能在新加入的直线的右侧,此时显然deque头端的直线对半平面的交没有了影响,所以我们要弹出队头直线。

(注:此处图示中直线的含义与情况三图示相同)

如图,若deque中的前两条直线的交点在新加入的直线的左侧,同情况二,直接将新的直线入队。

由于deque中前两条直线的交点可能在队尾直线的右侧,此时队头直线是多余的。

并且最后两条直线的交点也有可能在队头直线的右侧,此时队尾直线是多余的。

所以,我们还要将队头和队尾多余的直线弹出。

例题:Luogu P4196

这题说是让我们求凸多边形的交,但实际上,每个凸多边形可以看成若干个半平面的交,多边形的每条边所在的直线代表着一个半平面,这样我们就将原问题转化为了半平面交。

代码(注意,若要使用本代码解决其他半平面交问题,请确保构造的直线的有效半平面为左半平面):

#include <bits/stdc++.h>
#define maxn 1005
using namespace std;
struct Vec{double x,y;}v[maxn];
struct Line{Vec a,b;double k;}L[maxn],l[maxn]; 
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double angle(Vec a){return atan2(a.y,a.x);}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec getcrp(Line A,Line B)
{
    double k=cross(B.a-A.a,B.b-B.a)/cross(A.b-A.a,B.b-B.a);
    return A.a+(A.b-A.a)*k;
}
bool cmp(Line A,Line B)
{
    if (A.k!=B.k) return A.k<B.k;
    return cross(B.a-A.a,B.b-A.a)>0;//若斜率相等,就将在最左侧的直线排到最前面 
}
bool check(Line l,Vec p){return cross(l.b-l.a,p-l.a)<0;}//判断点是否在直线的右侧 
int n,N,num;
int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++)
    {
        int m;
        scanf("%d",&m);
        for (int j=1;j<=m;j++) scanf("%lf%lf",&v[j].x,&v[j].y);
        for (int j=1;j<=m;j++)
        {
            int tar=j%m+1;
            N++;
            L[N].a=v[j];L[N].b=v[tar];L[N].k=angle(v[tar]-v[j]);
        }
    }
    sort(L+1,L+1+N,cmp);
    for (int i=1;i<=N;i++)
        if (i==1||L[i].k!=L[i-1].k) l[++num]=L[i];
    Line q[maxn*2];
    int h=1,t=0;
    Vec crp[maxn];
    for (int i=1;i<=num;i++)
    {
        while(h<t&&check(l[i],crp[t])) t--;//情况一 
        while(h<t&&check(l[i],crp[h+1])) h++;//情况三 
        q[++t]=l[i];
        if (h<t) crp[t]=getcrp(q[t-1],q[t]);
    }
    while(h<t&&check(q[t],crp[h+1])) h++;//弹出多余直线 
    while(h<t&&check(q[h],crp[t])) t--;//弹出多余直线 
    crp[h]=getcrp(q[h],q[t]);
    double ans=0.0;
    for (int i=h+1;i<t;i++) ans+=cross(crp[i]-crp[h],crp[i+1]-crp[h]);//求面积 
    printf("%.3lf",ans/2);
    return 0;
}