二维几何基本运算
在开始写算法模板前,要定义好相关的结构体
定义点与向量
二维向量运算:
#define IN inline
struct Point{
double x,y;
int id=0;
//排序为Andrew凸包算法做准备
IN bool operator<(const Point &t)const {//按照x排序
if(x!=t.x) return x<t.x;
return y<t.y;//x相同按照y排序
}
IN Point operator+(const Point &t)const {return {x+t.x,y+t.y};}
IN Point operator-(const Point &t)const {return {x-t.x,y-t.y};}
IN Point operator*(const double &t)const {return {x*t,y*t};}//数乘
IN double operator*(const Point &t)const {return x*t.y-y*t.x;} //叉积
IN double operator^(const Point &t)const {return x*t.x+y*t.y;}
friend ostream & operator<<(ostream &out,const Point &t);
};
ostream & operator<<(ostream &out,const Point &t){out<<t.x<<' '<<t.y;return out;}
IN double sqr(double x){return x*x;}
IN double dis(Point a){return sqr(a.x)+sqr(a.y);}//距离平方
$a\times b$ 表示的几何意义是$|a|\cdot|b|\cdot sin(\theta) $,其中 $\theta $ 为 $a,b$ 的夹角,定义为从向量 $a$ 逆时针旋转到向量 $b$ 的角度
叉乘在后续的算法模板中,大量应用于判断点与向量的位置关系以及两个向量之间的位置关系
三维向量运算
#define IN inline
struct Point{
double x,y;
int id=0;
//排序为Andrew凸包算法做准备
IN bool operator<(const Point &t)const {//按照x排序
if(x!=t.x) return x<t.x;
return y<t.y;//x相同按照y排序
}
IN Point operator+(const Point &t)const {return {x+t.x,y+t.y};}
IN Point operator-(const Point &t)const {return {x-t.x,y-t.y};}
IN Point operator*(const Point &t)const{return {y*t.z-z*t.y,z*t.x-x*t.z,x*t.y-y*t.x};}
IN Point operator*(const double &t)const{return{x*t,y*t,z*t};}
IN double operator^(const Point &t)const{return x*t.x+y*t.y+z*t.z;}
friend ostream & operator<<(ostream &out,const Point &t);
};
ostream & operator<<(ostream &out,const Point &t){out<<t.x<<' '<<t.y<<' '<<t.z;return out;}
IN double sqr(double x){return x*x;}
IN double dis(Point a){return sqr(a.x)+sqr(a.y);}//距离平方
三维向量的叉乘求的是与两个向量垂直的法向量
定义点向式直线
//u为直线上一点,v为方向向量
//在模板里面的方向向量v的水平指向均指向x轴正向
struct Line{
Point u,v;
double a=atan2(v.y,v.x);
IN bool operator<(const Line t)const {//按照极角排序
if(a==t.a)return (t.u-u)*v>0;
//若构建下半面的半平面,极角相同,从左到右排序,改为>0
//若构造上半面的半平面,极角相同,从右到左排序,改为<0
return a<t.a;
}
IN bool operator !=(const Line t)const {return a!=t.a;}//判断极角是否相同
};
IN Point Rotation(Point v){return {-v.y,v.x};}//逆时针旋转90度
IN Point Inter(Line l1,Line l2){//求交点
//此处没有避免平行和0向量的情况,需要在代码solve部分解决
double t=(l2.v*(l2.u-l1.u))/(l2.v*l1.v);
return l1.u+l1.v*t;
}
计算n边形的面积
已知 $n$ 边形所有顶点的坐标
用 Point d[n]
表示这些点,则可以通过高斯面积公式
计算其面积为:
IN double S(Point d[]){
double s=0;
for(int i=1;i<n;++i){
s+=d[i]*d[i+1];//注意内部的加法是不带绝对值的
}
s+=d[n]*d[1];
return fabs(s)/2;//确保面积为正数
}
注意:叉乘不满足交换律,所以千万不能写反!
计算平行面卡壳最短距离
用 $O(n^2)$ 维护所有的向量,之后求出两两向量之间构成的法向量(叉乘求出)
枚举所有的点到法向量上的最长和最短投影,两个投影相减即为答案
复杂度:$O(n^5)$
凸包相关
凸包构造
Andrew
算法:
按 $先x,后y$ 的优先顺序排序
从前往后扫一次,求下半凸包
再从后往前扫一次,求上半凸包
当前点在直线下方时弹栈,否则加入
//n记录点数,top记录栈大小,栈尾有重复的1号点
//a[]表示输入的点集,stk[]表示从[1,top]顺时针构造凸包的点集
int top,n;
Point a[maxn],stk[maxn];
void Andrew(){
sort(a+1,a+1+n);
for(int i=1;i<=n;++i){
while(top>1 && ((a[i]-stk[top-1])*(stk[top]-stk[top-1]))<=0)--top;
stk[++top]=a[i];
}
int tmp=top;
for(int i=n;i>=1;--i){
while(top>tmp && ((a[i]-stk[top-1])*(stk[top]-stk[top-1]))<=0)--top;
stk[++top]=a[i];
}
--top;//去掉重复的1号点
}
$a\times b =|a|*|b|*sin(\theta)$ 其中 $\theta$ 表示 向量 $a$ 逆时针旋转到达 $b$ 时转过的角度
凸包是顺时针的,所以
(stk[top]-stk[top-1])
的向量指向 $x$ 轴负方向当这个角度是负数时代表点在向量左侧
旋转卡壳
求二维凸包直径
一般来说,要先用Andrew
算法求凸包
求出凸包之后,用比较面积的方式求出离初始直线最远的一点
之后发现直线每旋转一个点,最远点也会随之旋转一个点
扫描时间为 $O(1)$,扫描 $n$ 次,所以时间复杂度为 $O(n)$
//返回距离的平方
//top,stk[]与Andrew算法中相同
double get_long(){//旋转卡壳求直径
double ans=0;
for(int i=1,j=2;i<=top;++i){
while(fabs((stk[i]-stk[i+1])*(stk[i]-stk[j]))<fabs((stk[i]-stk[i+1])*(stk[i]-stk[j%top+1])))
j=j%top+1;//寻找最远点
ans=max(ans,max(dis(stk[i]-st[j]),dis(stk[i+1]-stk[j])));
}
return ans;
}
求K远点对距离
由于一次旋转卡壳算法可以求出最远点对
所以在 $K$ 较小的时候可以作 $K$ 次旋转卡壳
用大小为 $K$ 的小根堆维护答案
每次求出当前最远的点对 $(u,v)$ 并把剩余所有点到这两点的距离加入堆中,之后删除这两点
重复操作 $K$ 次后,堆顶即为答案
复杂度稳定 $O(nklogn+nk)$
//输出为距离的平方,所以定义的都是整数
//n记录点数,top记录栈大小,K记录要求的第k远
//a[]记录输入,stk[]记录每次求出的凸包顶点
//vis[]记录是否删除
//q为小根堆
typedef long long ll;
ll top,n,k;
Poll a[maxn],stk[maxn];
ll vis[maxn];
priority_queue<ll,vector<ll>,greater<ll>>q;
void Andrew(){
for(ll i=1;i<=n;++i){
if(vis[a[i].id])continue;//去掉的点不进循环
while(top>1 && ((a[i]-stk[top-1])*(stk[top]-stk[top-1]))<=0)--top;
stk[++top]=a[i];
}
ll tmp=top;
for(ll i=n;i>=1;--i){
if(vis[a[i].id])continue;//去掉的点不进循环
while(top>tmp && ((a[i]-stk[top-1])*(stk[top]-stk[top-1]))<=0)--top;
stk[++top]=a[i];
}
--top;//去掉重复的1号点
}
void insert(ll x){//插入大小为k的小根堆
q.push(x);
if(q.size()>k)q.pop();
}
void get_long(){
ll ans=0,r=0,l=0;
for(ll i=1,j=2;i<=top;++i){
while(fabs((stk[j]-stk[i])*(stk[i+1]-stk[i]))<fabs((stk[j%top+1]-stk[i])*(stk[i+1]-stk[i])))
j=j%top+1;
if(ans<max(dis(stk[j]-stk[i+1]),dis(stk[j]-stk[i]))){//开始获取最长距离的两点
ans=max(dis(stk[j]-stk[i+1]),dis(stk[j]-stk[i]));
r=stk[j].id;
if(dis(stk[j]-stk[i+1])>dis(stk[j]-stk[i]))l=stk[i+1].id;
else l=stk[i].id;
}
}
vis[l]=1,vis[r]=1;
for(ll i=1;i<=n;++i){
if(vis[i])continue;
insert(dis(a[i]-a[r]));
insert(dis(a[i]-a[l]));
}
insert(ans);
}
void solve(){
cin>>n>>k;
for(ll i=1;i<=n;++i){
cin>>a[i].x>>a[i].y;
}
sort(a+1,a+1+n);
for(ll i=1;i<=n;++i)a[i].id=i;
for(ll i=1;i<=k;++i){
top=0;
Andrew();
get_long();
}
cout<<q.top()<<'\n';
}
求最小矩形覆盖
旋转卡壳求出最高点
利用往右的点一定是点积由正到负,往左的点的点积一定由负到正的特点求出左右边界点
利用直线相交求出交点与面积,或者直接求面积
//top,stk[]与Andrew算法中意义相同
//l,r表示最小矩形的左右边界上的一点
double ans=1e18;
void cal(int i,int l,int r,int j){//面积,求交点参考1.0.2
double t1=fabs((stk[i]-stk[i+1])*(stk[i]-stk[j]));
double t2=fabs((stk[i]-stk[i+1])^(stk[i]-stk[r]))+fabs((stk[i]-stk[i+1])^(stk[i]-stk[l]));
double t3=fabs((stk[i]-stk[i+1])^(stk[i]-stk[i+1]));
ans=min(ans,t1*t2/t3);
}
void get_s(){
for(int i=1,j=2,r=2,l=2;i<=top;++i){
while(fabs((stk[i]-stk[i+1])*(stk[i]-stk[j]))<fabs((stk[i]-stk[i+1])*(stk[i]-stk[j%top+1])))j=j%top+1;
l=max(l,j);
while(((stk[i+1]-stk[i])^(stk[r%top+1]-stk[r]))>=0)r=r%top+1;
while(((stk[i+1]-stk[i])^(stk[l%top+1]-stk[l]))<=0)l=l%top+1;
cal(i,l,r,j);//计算面积或者交点
}
}
求最大内接四边形
旋转卡壳可以求凸包内最大内接四边形
最大内接四边形的对角线一定是对踵点
通过对角线确定最大内接四边形
在维护对角线的时候,利用最大内接四边形的特点,用类似旋转卡壳的思想维护另一条对角线
时间复杂度:$O(nlogn+n)$ 排序复杂度加旋转卡壳复杂度
//u,v为一条对角线,l,r为另一条对角线
double ans=0;
void get_s_in(){
for(int u=1,v=2,r=1,l=2;u<=top;++u){
while((stk[u%top+1]-stk[u])*(stk[v%top+1]-stk[v])<0){//找对角线
if(l==v)l=l%top+1;//保证l在v的后面
v=v%top+1;
}
while(l%top+1!=u&&fabs((stk[u]-stk[v])*(stk[l]-stk[v]))<fabs((stk[u]-stk[v])*(stk[l%top+1]-stk[v])))
l=l%top+1;
while(r%top+1!=v&&fabs((stk[u]-stk[v])*(stk[r]-stk[v]))<fabs((stk[u]-stk[v])*(stk[r%top+1]-stk[v])))
r=r%top+1;
ans=max(ans,fabs((stk[u]-stk[v])*(stk[r]-stk[v]))+fabs((stk[u]-stk[v])*(stk[l]-stk[v])));
if(r==u)r=u%top+1;
}
}
对踵点:对于凸包上的一点 $P$ ,过 $P$ 作一条不经过凸包内部的直线 $l$ ,在凸包上能够找到另一个点 $P’$ 使得过 $P’$ 作的 $l$ 发平行线同样不经过凸包内部, $P,P’$ 互为对踵点。
旋转卡壳算法的本质其实就是利用了对踵点的特点,即 $P$ 顺时针移动位置时, $P’$ 必然也会顺时针移动
半平面交
半平面交求面积
在点向直线的基础上,选择构造下半平面
先对极角进行排序
在保证极角不同的情况下,对于当前预插入向量:
弹出队尾,若上一个交点在向量右侧则弹出
弹出队首,若第一个交点在向量右侧则弹出
之后插入当前向量,求出新的交点
所有向量插入完毕后,删除多余的元素
//输出半平面交的面积
//e[]存输入的向量
//q[]为单调队列
//(l,r]维护单调队列的有效范围,左开右闭,便于更新交点
//ans[]存合法点集
Line e[maxn];
solve(){
sort(e+1,e+1+tot);//排序
int l = 0, r = 0;
for (int i = 1; i <= tot; ++i)
if (e[i] != e[i - 1]) {
// 注意要先检查队尾
while (r - l > 1 && e[i].v * (e[i].u - ans[r]) >0) // 如果上一个交点在向量右侧则弹出队尾
--r;
while (r - l > 1 && e[i].v * (e[i].u - ans[l + 2]) >0) // 如果第一个交点在向量右侧则弹出队首
++l;
q[++r] = e[i];
if (r - l > 1) ans[r] = inter(q[r], q[r - 1]); // 求新交点
}
while (r - l > 1 && q[l+1].v * (q[l + 1].u - ans[r]) > 0) // 注意删除多余元素
--r;
if(r<=l+2)cout<<"0.000\n";//剩余的点只有1个或者没有的情况,直接输出0
else{
ans[r + 1] = inter(q[l + 1], q[r]); // 再求出新的交点
++r;
double sum=0;
for(int i=l+3;i<=r;++i)sum+=ans[i-1]*ans[i];//逆时针三角剖分求出面积
sum+=ans[r]*ans[l+2];
cout<<fixed<<setprecision(3);
cout<<fabs(sum/2);
}
}
半平面交确定可见直线
从上方往下看,作上半面的半平面交
从下往上看,作下半面半平面交
之后按照极角排序
对于一个预插入的直线
- 若栈中元素小于两个直接插入
- 若当前直线与栈顶
top
直线的交点在上上一个top-1
直线的下方则弹出栈顶元素 - 否则插入当前直线
//从大到小输出可见直线的序号 //e[]记录输入直线 //stk[]维护动态栈 //ans[]记录最终答案 int stk[maxn] void solve(){ sort(e+1,e+1+n); int top=0; for(int i=1;i<=n;++i){ if(i==n||e[i]!=e[i+1]){//排除相同的直线 //此处的inter写内部,防止top为0时出现与0向量相交的情况 while(top>1&&(inter(e[i],e[stk[top]])-e[stk[top-1]].u)*e[stk[top-1]].v>=0)top--; stk[++top]=i; ans[top]=e[i].id; } } for(int i=1;i<=top;++i)ans[i]=e[stk[i]].id; sort(ans+1,ans+1+top); for(int i=1;i<=top;++i) cout<<ans[i]<<' '; }
平面最近点对
用分治法可以在 $O(nlogn)$ 的时间内求出平面最近点对
首先对所有的点按照 $x$ 的坐标从小到大排序
之后取所有点的中位数 $mid$,并不断计算点 $mid$ 两边区域中,各点对的最小距离
假设求出的最小距离为 $d$
而由于从两边区域各取一点构成的点对距离可能小于 $d$ ,所以要进一步讨论
对 $y$ 坐标进行归并排序
要讨论横坐标的范围为 $[mid.x-d,mid.x+d]$
而纵坐标的讨论范围为 $[mid.y-d/2,mid.y+d/2]$
此时发现,在这个区间内合法的点最多只能有 $6$ 个,暴力枚举即可
//a[]:输入的点集,tmp[]:用于合并的临时数组,b[]:提取出的合法枚举区间内的点集
Point a[maxn],tmp[maxn],b[maxn];
void Merge(ll l,ll r){
ll mid=(l+r)>>1,i=l,j=mid+1;
for(ll k=l;k<=r;++k){
if(i<=mid&&(j>r||a[i].y<a[j].y))tmp[k]=a[i++];
else tmp[k]=a[j++];
}
for(ll i=l;i<=r;++i)a[i]=tmp[i];
}
ll get_close(ll l, ll r) {
if (r - l<=0) return inf;//到自己的距离定义为无穷
if(r-l==1){if(a[l].y>a[r].y)swap(a[l],a[r]);return dis(a[l]-a[r]);}//区间内只有两点,返回距离并排序
// 分解
ll mid = (l + r) >> 1;
ll t=a[mid].x,cnt=0;//t:中位数的x坐标 cnt:满足条件的点的个数
ll d=min(get_close(l, mid),get_close(mid+1, r));//获取区间最小距离
// 合并
Merge(l,r);//对y归并排序
for(ll i=l;i<=r;++i){
if((a[i].x-t)*(a[i].x-t)<d){//横向距离比d小,
b[++cnt]=a[i];
}
}
for(ll i=1;i<=cnt;++i){
for(ll j=i+1;j<=cnt&&(b[j].y-b[i].y)*(b[j].y-b[i].y)<d;++j)
d=min(d,dis(b[j]-b[i]));
}
return d;
}
void solve(){
ll n;cin>>n;
for(ll i=1;i<=n;++i)cin>>a[i].x>>a[i].y;
sort(a+1,a+1+n);
cout<<get_close(1,n)<<'\n';
}
随机增量法
随机增量法通常与随机化结合,追求均摊时间复杂度最优化
类似于数学归纳法,本质是将一个问题转化为规模刚好小一层的子问题,解决后加入当前对象
$$
T(n)=T(n-1)+g(n)
$$
最小圆覆盖
由于三点确定一个圆
设前 $i-1$ 个点构成的圆半径为 $r$
对于第 $i$ 个点,若在圆内或者圆上,则属于合法点,不作任何处理
若在圆外,则需要构造新的圆且点 $i$ 必然在圆上
以第 $i$ 个点为圆心,半径为 $0$ 开始枚举 $i$ 之前的一点 $j$
确定 $j$ 后,以 $i,j$ 中点作为圆心,$i,j$ 距离的一半作为半径作出新的圆并开始枚举 $j$ 前的一点 $k$
若 $k$ 在圆内同样不处理
若在圆外则以 $i,j,k$ 三点作为新的圆
看上去是 $O(n^3)$ 的复杂度,实际上随机化数据后,均摊复杂度为 $O(n)$
//n:输入点数
//a[]:输入的点
int n;
Point a[maxn];
Point GetCircle(Point a,Point b,Point c){//中垂线交点求圆心
Point m1=(a+b)*0.5;
Point m2=(c+b)*0.5;
Line l1={m1,Rotation(a-b)};
Line l2={m2,Rotation(c-b)};
return Inter(l1,l2);
}
void GetMinCircle(){//求最小圆
Point o=a[1];//圆心
double r=0;
for(int i=1;i<=n;++i){
if(dis(o-a[i])<=r)continue;//在圆内或者圆上
o=a[i],r=0;
for(int j=1;j<i;++j){
if(dis(o-a[j])<=r)continue;
o=(a[i]+a[j])*0.5;//中d
r=dis(a[j]-a[i])/4;
for(int k=1;k<j;++k){
if(dis(o-a[k])<=r)continue;
o=GetCircle(a[i],a[j],a[k]);
r=dis(o-a[k]);
}
}
}
cout<<fixed<<setprecision(5);
cout<<o<<' ';
cout<<sqrt(r)<<'\n';
}
void solve(){
cin>>n;
for(ll i=1;i<=n;++i)cin>>a[i].x>>a[i].y;
random_shuffle(a+1,a+1+n);//打乱a数组
GetMinCircle();
}