计算几何-优快云博客

本文链接：https://blog.youkuaiyun.com/wu6789/article/details/115666169

大纲对计算几何的要求：

矢量及其运算

点线面之间的位置判断

常见图形的面积计算

二维凸包的求法及其应用

半平面交

问题分解

复杂的算法都是由许多简单的算法组合而成的

计算几何最基本的算法：
求直线的斜率
求2条直线的交点
判断2条线段是否相交
求叉积

计算几何经典算法
求凸包
求最近点对
判断点是否在多边形内等等

几何题的类型

纯粹的计算求解题
比如求直线的斜率时，直线的斜率为无穷大，求2条直线的交点时，2直线平行，等等。
解这一类题除了需要有扎实的解析几何的基础，还要全面地看待问题，仔细地分析题目中的特殊情况。

存在性问题
这一类问题可以用计算的方法来直接求解，如果求得了可行解，则说明是存在的，否则就是不存在的，但是模型的效率同模型的抽象化程度有关，模型的抽象化程度越高，它的效率也就越高，几何模型的的抽象化程度是非常低的，而且存在性问题一般在一个测试点上有好几组测试数据，几何模型的效率显然是远远不能满足要求的，这就需要对几何模型进行一定的变换，转换成高效率的模型。

求几何中的最佳值问题
这类问题是几何题中比较难的问题，一般没有什么非常有效的算法能够求得最佳解，最常用的是用近似算法去逼近最佳解，近似算法的优劣也完全取决于得出的解与最优解的近似程度。

存储方法

一般为了避免分类讨论和减少精度误差，通常不使用解析几何的方法

点：(x,y) 坐标

直线/线段/射线：两个点，带方向

圆：圆心和半径

多边形：按顺时针/逆时针存储所有的顶点

直线：存直线上的两个点，线段：存两个端点，射线：存端点和射线上的一个点

向量的表示

向量：把起点平移到坐标原点，记录下终点的坐标，存储同点一样，用一个数对表示，(x,y)

向量的缩放：(x,y) -> (kx,ky)

把要研究的图形放在平面直角坐标系或极坐标系下，这样解决问题就会方便很多。

在 n 维空间下，矢量经常被表达为 n 个数的元组 a=(a1,a2,a3...,an)。

在高等代数中，n 维矢量一般表示n×1 的矩阵。

在二维空间下则以 (x,y) 一对整数表示。

精度问题

实数运算会导致精度问题，浮点误差

常用方法：设定一个eps，例如 10^-8

判断a b相等，fabs(a-b)<eps

不相等>eps

a<b：a+eps < b

a<=b： a<b+eps

单位向量

矢量 $\frac{a}{\left | a \right |}$ 是与 a 同向的单位矢量，即模长是 1 的矢量。

a 同向，但长度是 l 的矢量，为 $l \frac{a}{\left | a \right |}$ 。

与 a 共线但方向相反，长度是 l 的矢量为 $-l \frac{a}{\left | a \right |}$ 。

矢量的垂直

垂直的矢量有 $\theta = \frac{\pi}{2}$ ，即 $cos \theta = 0$ .

所以 a 和 b 垂直定义为 $a \cdot b = 0$

点积

点积 $a \cdot b = \left | a \right| \left | b \right| cos \theta$ ，a的模长乘b的模长再成ab夹角的cos值，其中 θ 为 a 和 b 之夹角。

a 与 b 的点积的值实际上就是 a 的模乘 b 在 a 上的投影的模，但是若其投影与 a 方向相反则为负。

两个 n 维矢量的点积是一个标量，有
$(a_1,a_2,a_3,...,a_n)\cdot(b_1,b_2,b_3,...,b_n)=a_1b_1+a_2b_2+a_3b_3+...+a_nb_n$

向量的模（即长度） $\left|a \right| = \sqrt{a \cdot a} = \sqrt{a^2}$

点积满足交换律。

点积满足分配律 $a \cdot (b+c)=a \cdot b + a \cdot c$

两个点(x1,y1) (x2,y2)，点积 $(x1,y1) \cdot (x2,y2)=(x_1 i + y_1 j) \cdot (x_2 i + y_2 j)=(x_1 i x_2 i + x_2 i y_1 j + x_1 i y_2 j + y_1 j y_2 j)=x_1 x_2 i i+y_1 y_2 j j=x_1 x_2 + y_1 y_2$

（i,j分别是x y 轴的单位向量，根据定义 ii=1 ，ij=0）

应用

两个向量垂直，等价于他们的点积为0

两个向量夹角<90,等价于他们的点积>0

c落在以AB为x轴的坐标系的y轴的左侧还是右侧

叉积

叉积的结果是向量，而不是标量

$a \times b = \left |a\right|\left | b \right|sin \theta$

几何意义：以两个向量为邻边构成的平行四边形的有向面积（底乘高），当b在a的左侧为正，右侧为负

c落在以AB为x轴的坐标系的x轴的上侧还是下侧

通常判断位置关系只使用叉积就行，比如求凸包

有时候需要判断在哪个象限，需要用到点积和叉积，比如极角排序

$a \times b = -b \times a$

叉积满足分配律 $a \times (b+c)=a \times b + a \times c$

设 $a=(x_1,y_1) , b=(x_2,y_2)$ ，

$(x_1,y_1)\times(x_2,y_2)=x_1 y_2 -x_2 y_1$

应用

两个向量平行，等价于他们的叉积为0

求三角形面积，fabs(两个边向量的叉积/2)

点到直线的距离，连接顶点以及直线上的两点，组成三角形，求出面积和直线上两点间的距离，会得到高，即点到直线的距离

投影

已知线段AB，和一点C，求C落到AB上的投影位置O的坐标

$=A+unit(B-A)*\frac{\overrightarrow{AC}\cdot \overrightarrow{AB}}{\left| AB \right|}$

反射

点相对于线段的对称点

已知线段AB，和一点C，求C相对于AB的对称点

先求出C的垂足O，C‘=2O-C

计算几何就是搭积木的过程

两条直线交点

求AB和CD的交点

取AB上的一点O，O=A+x(B-A)

$\overrightarrow{OD}\times \overrightarrow{OC}$ 是关于x的一次函数

O=A+x(B-A)

利用叉积计算出两个三角形ACD、BCD的面积比，求出x，进而得到O

套路

平移旋转法（OI）

向量(x,y)仰角为theta，逆时针旋转alpha角度，得到的新向量

point turn(db k1){
	return (point){
		x*cos(k1)-y*sin(k1),x*sin(k1)+y*cos(k1)
	};
}

用途，求凸包，有一条线段两个端点x坐标相同，这时候sin值是不存在的，程序会出问题

解决方法：所有点绕原点旋转alpha角度，alpha是一个随机数

点线面的表示

(x[i],y[i])代表第i个点的xy坐标值，写代码不方便，中括号太多

解决方法：

面向对象

点类（加减乘除，点积，叉积，旋转，距离）

直线类（交点，判平行）

多边形类

圆类

写代码会比较方便

单元操作

点到直线的投影
直线/圆和直线/圆的交点
点/圆到圆的切线
点和多边形位置关系

点和多边形的位置关系

判断点在多边形内部还是外部

射线法

存在问题：射线与多边形的端点相交，或与多边形的某条边重叠

解决方法1：所有点旋转一个随机度数，有精度误差

解决方法2：射线的角度随机

解决方法3

把多边形每条线段都标记为左闭右开

旋转向量

a=(x,y)，把a旋转 $\theta$ 弧度后，得到的新向量 $(xcos\theta-ysin\theta,xsin\theta+y cos \theta)$ ，证明：利用三角函数和角公式

直线的交点

AB和CD交于点O

先使用叉积求出ABC ABD的面积，得到OC和OD的比值，然后就可以求出O点的坐标了

多边形面积

已知一个简单多边形，顶点按逆时针顺序依次为 p1 p2 p3 ... pn

其面积为 $\frac{1}{2}(p_n \times p_1 + \sum_{i=1}^{n-1}p_i \times p_{i+1} )$

证明：考虑平面上的一个点，我们以它为起点引一条反向延长线过原点的射线。如果该点在多边形外，这条射线会经过偶数条多边形的边，最终这个点的面积不会被算入答案中:如果在多边形内则会经过奇数条边，相抵消后恰好计算了一次它的面积。

圆的交点

用余弦定理算出 $\theta$ ，然后把D点绕圆心A旋转即可求出交点C的坐标。

两圆的公切线

BE=两圆的半径差，用AB和BE算出∠BAE，AB缩到小圆半径，并旋转得到C，BA缩小到大圆半径，并旋转得到D

内公切线也类似

凸包

凸多边形是指所有内角大小都在 [0,pi]范围内的简单多边形

凸包定义为：对于给定点集合X ，所有包含X的凸集的交集。

凸包就是能包围给定点的最小的凸多边形。

凸包用最小的周长围住了给定的所有点

求凸包

Graham扫描法:把点集按极角排序后用栈求出凸包。（极角序）

struct Node
{
       int x,y;
}p[MAX],S[MAX];//p储存节点的位置，S是凸包的栈 
inline bool cmp(Node a,Node b)//比较函数，对点的极角进行排序 
{
       double A=atan2((a.y-p[1].y),(a.x-p[1].x));
       double B=atan2((b.y-p[1].y),(b.x-p[1].x));
       if(A!=B)return A<B;
       else    return a.x<b.x; //这里注意一下，如果极角相同，优先放x坐标更小的点 
}
long long Cross(Node a,Node b,Node c)//计算叉积 
{
       return 1LL*(b.x-a.x)*(c.y-a.y)-1LL*(b.y-a.y)*(c.x-a.x);
}
void Get()//求出凸包 
{
       p[0]=(Node){INF,INF};int k;
       for(int i=1;i<=n;++i)//找到最靠近左下的点 
              if(p[0].y>p[i].y||(p[0].y==p[i].y&&p[i].x<p[0].x))
               {p[0]=p[i];k=i;}
       swap(p[k],p[1]);   
       sort(&p[2],&p[n+1],cmp);//对于剩余点按照极角进行排序 
       S[0]=p[1],S[1]=p[2];top=1;//提前在栈中放入节点 
       for(int i=3;i<=n;)//枚举其他节点 
       {
              if(top&&Cross(S[top-1],p[i],S[top])>=0)
                        top--;//如果当前栈顶不是凸包上的节点则弹出 
              else  S[++top]=p[i++];//加入凸包的栈中 
       }
       //底下这个玩意用来输出凸包上点的坐标 
       //for(int i=0;i<=top;++i)
       //    printf("(%d,%d)\n",S[i].x,S[i].y);
}

https://blog.youkuaiyun.com/qq_30974369/article/details/76405546

Andrew算法

是Graham扫描法的变种。
把所有点按横坐标为第一关键字、纵坐标为第二关键字升序排序。（水平序）

显然排序后最小的元素和最大的元素一定在凸包上。而且因为是凸多边形，我们如果从一个点出发逆时针走，轨迹总是“左拐”的，一旦出现右拐，就说明这一段不在凸包上。因此我们可以用一个单调栈来维护上下凸壳。

因为从左向右看，上下凸壳所旋转的方向不同，为了让单调栈起作用，我们首先升序枚举求出下凸壳，然后降序求出上凸壳

找出最左侧和最右侧的点，从左侧一直枚举到最右侧的点，按顺序枚举，如果该点在当前凸包前进方向的左侧则把它加入栈，否则不断退栈直到满足前面的条件。枚举完我们就得到了点集的下凸包。
倒序枚举每个点，用同样的方法求出，上凸包。时间复杂度的瓶颈在排序上。

注意:一定要按照纵坐标为第二关键字排序，否则多于三个点横坐标相等的情况可能会挂。

旋转卡壳

由于凸多边形的性质，固定一条边，其他端点到这条边的距离都是单峰函数，可以用三分求出峰值

旋转卡壳利用凸多边形的性质简化运算

例求包含凸包的面积最小的矩形

方法一：显然矩形一定有一条边在凸包的边上，求出峰值和距离最远的两个垂足，得到矩形，多个矩形比较

方法二：用一条直线在凸包的边上按一定方向滚动，利用单调性确定峰值（旋转卡壳）

while i+1 比 i 优秀 { i修改为i+1}

例给定凸包，和一个方向，求凸包在这个方向上的两个切点

上半凸包下半凸包各一个切点

对于上半凸包，找到一个点，凸包上的一条边斜率大于给定方向，下一条边斜率小于给定方向，二分

下半凸包类似

https://www.luogu.com.cn/problem/P1452

合并凸包

如何快速的合并两个凸包

一边归并一边用前面的方法维护当前的凸包

与凸包相关的询问

查询一个点是否在凸包内
在上下凸包上分别按照水平序二分出在哪两个点之间，判一下是否在该线段内侧。

一条斜率为k的直线从正无穷远处向下平移，问碰到的第一个点?

在凸包上二分斜率

动态凸包

只有插入，离线且对答案贡献独立: CDQ分治。
例题: NOI2007 货币兑换(洛谷P4027)

只有插入，在线或对答案贡献不独立:用平衡树维护。
例题: CodeForces 70D Professor's task

要求插入和删除，离线:把每个点视为在时间[l,r]内有效，然后外面套一个线段树分治。
例题: BZ0J4311 向量

要求插入和删除，在线:平衡树套可持久化平衡树维护，详见陈立杰的《可持久化数据结构研究》。

闵可夫斯基和

上定义两个点集A,B的闵可夫斯基和为 A+B= {a+b|a∈A,b∈B}。

结论1 jie闵可夫斯基和一定是凸集

凸集的定义， $if \ \ \alpha, \beta \in A,then \ \ x \alpha + (1-x) \beta \in A$ ，两个点在集合里，两点间的任意点都在集合里

$if \ \ \alpha_1, \alpha_2 \in A, \ \ \beta_1, \beta_2 \in B,then \ \ \alpha_1+ \beta_1,\alpha_2+\beta_2 \in A+B$

$(\alpha_1+\beta_1)x+(1-x)(\alpha_2 + \beta_2)=(\alpha_1x+(1-x)\alpha_2)+(\beta_1 x+(1-x)\beta_2)$

结论2

从形态上看，一条出现在任意凸包上的边，一定会出现在闵可夫斯基和的凸包上

计算闵可夫斯基和：

形态：A凸包，边 $\overrightarrow{A_i A_{i+1}}$ ，B凸包，边 $\overrightarrow{B_i B_{i+1}}$ ，按斜率排序 O(n)

位置：通过A和B的x坐标最小值，y坐标最小值确定闵可夫斯基和的xy坐标最小值

例 JSOI2018 P4557战争

一个点集的领地为它的凸包(包括边界)

给出两个点集A,B以及q组询问

每组询问给出一个向量，问把A沿着这个向量平移后，领地是否会和B有交

alpha是一个偏移向量，A偏移以后的集合A+alhpa

需要判断A+alpha和B是否相交

考虑哪些alpha会使得集合A+alpha和集合B相交

在A中存在点x，在B中存在点y，满足y=x+alhpa，即 alpha=y-x

所有alpha的集合C， $C=\{y-x|y \in B,x \in A \}$

可以看出是闵可夫斯基和的形式，对A中的所有点取负，和B做闵可夫斯基和

判断偏移点是否在凸包内集合

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
const ll N=1e5+10;
struct Node
{
	ll x,y;
	Node operator - (Node A) {return (Node){x-A.x,y-A.y};}
	Node operator + (Node A) {return (Node){x+A.x,y+A.y};}
	ll operator * (Node A) const {return x*A.y-y*A.x;}
	ll len() const {return x*x+y*y;}
}A[N],C1[N],C2[N],s1[N],s2[N],bs;
ll cmp1(const Node&A,const Node&B) {return A.y<B.y||(A.y==B.y&&A.x<B.x);}
ll cmp2(const Node&A,const Node&B) {return A*B>0||(A*B==0&&A.len()<B.len());}
ll n,m,sta[N],top,q,tot;
void Convex(Node *A,ll &n)
{
	sort(A+1,A+n+1,cmp1);
	bs=A[1];sta[top=1]=1;
	for(ll i=1;i<=n;i++) A[i]=A[i]-bs;
	sort(A+2,A+n+1,cmp2);
	for(ll i=2;i<=n;sta[++top]=i,i++)
		while(top>=2&&(A[i]-A[sta[top-1]])*(A[sta[top]]-A[sta[top-1]])>=0) top--;
	for(ll i=1;i<=top;i++) A[i]=A[sta[i]]+bs;
	n=top;A[n+1]=A[1];
}
void Minkowski()
{
	for(ll i=1;i<n;i++) s1[i]=C1[i+1]-C1[i];s1[n]=C1[1]-C1[n];
	for(ll i=1;i<m;i++) s2[i]=C2[i+1]-C2[i];s2[m]=C2[1]-C2[m];
	A[tot=1]=C1[1]+C2[1];
	ll p1=1,p2=1;
	while(p1<=n&&p2<=m) ++tot,A[tot]=A[tot-1]+(s1[p1]*s2[p2]>=0?s1[p1++]:s2[p2++]);
	while(p1<=n) ++tot,A[tot]=A[tot-1]+s1[p1++];
	while(p2<=m) ++tot,A[tot]=A[tot-1]+s2[p2++];
}
ll in(Node a)
{
	if(a*A[1]>0||A[tot]*a>0) return 0;
	ll ps=lower_bound(A+1,A+tot+1,a,cmp2)-A-1;
	return (a-A[ps])*(A[ps%tot+1]-A[ps])<=0;
}
int main()
{
	cin>>n>>m>>q;
	for(ll i=1;i<=n;i++)
		scanf("%lld%lld",&C1[i].x,&C1[i].y);
	Convex(C1,n);
	for(ll i=1;i<=m;i++)
	{
		scanf("%lld%lld",&C2[i].x,&C2[i].y);
		C2[i].x=-C2[i].x;C2[i].y=-C2[i].y;
	}
	Convex(C2,m);
	Minkowski();
	Convex(A,tot);
	bs=A[1];for(ll i=tot;i>=1;i--) A[i]=A[i]-A[1];
	while(q--)
	{
		scanf("%lld%lld",&A[0].x,&A[0].y);
		printf("%lld\n",in(A[0]-bs));
	}
	return 0;
}

例 BZOJ2564

给出两个点集A,B，求它们的闵可夫斯基和的凸包的面积。|A|,|B|< 10^5，坐标的绝对值不超过10^8。
先分别求出点集A和点集B的凸包。记ai bi分别表示A或B的凸包上的第i个点。
首先显然a1+b1在A + B的凸包上。假设我们已知ai+bj在A+ B的凸包上，可以发现凸包上的下一个点就是
a(i+1) + bj,ai+b(j+1)中更凸的一个。因此拿两个指针用类似归并排序的方式扫一遍即可。