计算几何初步

# 计算几何初步 ## 基础知识 ### 1.点 没什么可说的 用$(x,y)$表示 ### 2.向量及运算 就是矢量 可以表示一条有向线段 模长 即向量的长度 向量的运算 (1)加/减 直接相加减 $(x_1\pm x_2,y_1\pm y_2)$ (2)数乘 直接乘 $(x*a,y*a)$ (3)点乘 （答案为标量）可以求投影的长度 $|a|*|b|*cos\theta=x_1*x_2+y_1*y_2$ (4)叉积 （答案为标量）也叫叉乘 求的是两个向量间围成的有向面积 $x_1*y_2-x_2*y_1$  如图即为向量$A$和向量$B$的叉积 注意叉积是没有交换律的 $A*B$的正负表示$B$在$A$的哪个方向上 顺时针$<0$ 逆时针$>0$ 如图即为逆时针的情况 所以$B*A=-A*B<0$ 如果要求围成的三角形的面积 则需要将叉积除以$2$ ```cpp struct Vec{ double x,y; Vec(){x=0;y=0;} friend Vec operator-(Vec a,Vec b){a.x-=b.x,a.y-=b.y;return a;} friend Vec operator+(Vec a,Vec b){a.x+=b.x,a.y+=b.y;return a;} friend double operator*(Vec a,Vec b){return a.x*b.y-b.x*a.y;} friend double operator^(Vec a,Vec b){return a.x*b.x+a.y*b.y;} } ``` ### 3.线段 和向量不同的是 线段要保存起点和终点 （可以理解为向量的起点默认为原点） 常用的表示方法是通过起点$+$向量表示法 即起点$+$这个向量 $=$终点 ### 4.常见问题 #### (1)判断线段和线段的位置关系 通过叉积判断即可 #### (2)点$A$是否在线段$BC$上 判断$|AB|+|AC|==|BC|$即可 #### (3)凸多边形的面积 就等于相邻的点与原点构成的向量的叉积之和再除以$2$  如图 注意叉积是有向的 则通过加加减减就能形成凸多边形的面积 #### (3)求两个线段的交点 可以通过解析几何来求 不过这里有一种用向量$+$相似的简便方法  如图 我们要求$AB$和$CD$的交点$E$ 只需要找到$t=\frac{AE}{AB}$ 再把$AB*t$就行了 构造向量$AC$ 然后求$AC$与$AF$的叉积$S_1$ 求$AB$和$CD$的叉积$S_2$ 则$t=\frac{S_1}{S_2}$ （这里的线段可以用起点$+$向量表示法） ```cpp inline Point get_Intersect_Point(Line a,Line b){ Point c=a.x-b.x; double t=(b.y*c)/(a.y*b.y); a.y.x*=t;a.y.y*=t; return a.x+a.y; } ``` #### (4)将向量/点排序 这里的排序是指按照顺时针或者逆时针排序 这种排序方法又称极角排序 要先选择一个极点 然后直接定义$<$为叉积的大小 （这里的向量为相对极点的向量） 如果要按照逆时针排序 则返回叉积$>0$ 否则返回$<0$ 也可以直接比较 $atan2(y,x)$的大小 按照从小到大排序为逆时针（三、四、一、二象限） 否则为顺时针 （注意这里选的极点为$(0,0)$） #### (5)判断点/向量在不在凸包内部 首先把凸包弄成到原点的向量形式然后排序 然后找这个向量（点）$q$逆时针方向上第一个凸包向量$v$ 判断$q$是否在和$v$下一个向量$v+1$所构成三角形的内部即可 ```cpp int in(Vec q){ if((q-bs)*(D[1]-bs)>0||(q-bs)*(D[tot]-bs)<=0)return 0; int pos=lower_bound(D+1,D+tot+1,q,cmp)-D-1; return ((q-bs)-(D[pos]-bs))*(D[pos%tot+1]-D[pos])<=0; } ``` #### (6)旋转点 可以随机打乱数据以构造出随机数据的效果 假设旋转$\theta$ 那么 $$ x'=x\cos\theta-y\sin \theta $$ $$ y'=x\sin \theta+y\cos \theta $$ ## 常用算法 ### 1.凸包 给定二维平面上的点集，凸包就是将最外层的点连接起来构成的凸多边形，它能包含点集中所有的点。  可以想象成一根橡皮筋围成的点 用$Graham$扫描法可以求一个凸包 我们首先要选择一个一定在凸包上的点 可以选择$y$坐标最小的 $y$相同时$x$最小的点 然后 以这个点作为极点 把其他点极角排序 然后维护一个队列 每次加入点时 通过叉积比较最后三个点形成的两个向量是否按照顺时针关系 如果满足就加进去 否则不断$r--$ 最后得到的就是一个凸包了 [P2742 【模板】二维凸包 / [USACO5.1\]圈奶牛Fencing the Cows](https://www.luogu.org/problem/P2742) 求凸包周长 每两个点之间算一下$dis$就行了 ```cpp #include<bits/stdc++.h> using namespace std; const int maxn=10005; struct Vec{ double x,y; Vec(){x=0;y=0;} friend Vec operator-(Vec a,Vec b){a.x-=b.x,a.y-=b.y;return a;} friend Vec operator+(Vec a,Vec b){a.x+=b.x,a.y+=b.y;return a;} friend double operator*(Vec a,Vec b){return a.x*b.y-b.x*a.y;} }O,p[maxn],s[maxn]; inline double dis(Vec a,Vec b){return sqrt((double)(b.y-a.y)*(b.y-a.y)*1.0+(double)(b.x-a.x)*(b.x-a.x)*1.0);} inline bool cmp(Vec a,Vec b){ Vec c=a-p[1],d=b-p[1];//按极点进行极角排序 return c*d==0?dis(O,c)<dis(O,d):c*d>0; } inline bool check(Vec a,Vec b,Vec c){ Vec d=b-a,e=c-b;//三个点必须保证连成的两个向量是逆时针方向即叉积>=0 return d*e<=0; } int n,top; double ans,mid; int main() { scanf("%d",&n); for(int i=1;i<=n;i++) { scanf("%lf%lf",&p[i].x,&p[i].y); if(p[i].y<p[1].y)swap(p[1],p[i]);//找极点 } sort(p+2,p+1+n,cmp);//按极点进行极角排序 s[++top]=p[1];//极点一定在凸包上 for(int i=2;i<=n;i++) { while(top>1&&check(s[top-1],s[top],p[i]))top--; s[++top]=p[i]; } s[++top]=p[1];//最后一个点和第一个点求距离 for(int i=1;i<top;i++)ans+=dis(s[i],s[i+1]); printf("%.2lf",ans); } ``` 还可以动态维护凸包 就是用$set$维护凸包中的点 每次插入时分别向左向右找到第一个满足凸包条件的点 用叉积判断就行了 [P2521 [HAOI2011]防线修建](https://www.luogu.org/problem/P2521) 动态维护凸包的周长 把删除点倒过来改成插入点就行了 写的比较早了 码风有点奇怪 ```cpp #include<bits/stdc++.h> using namespace std; int n,x,y,m,q; double now; struct Q{ int op,num; double ans; }opt[300000+19]; struct P{ int x,y; }p[200000]; bool vis[200000]; inline bool operator <(P a,P b) { if(a.x==b.x)return a.y<b.y; return a.x<b.x; } inline P operator -(P a,P b) { P t;t.x=a.x-b.x;t.y=a.y-b.y; return t; } inline double operator *(P a,P b) { return (a.x*b.y)-(b.x*a.y); } inline double dis(P a,P b) { return sqrt((double)(a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } set<P>s; inline void insert(int a,int b) { P x;x.x=a;x.y=b; set<P>::iterator r=s.lower_bound(x),l=r,t; l--; //叉积不满足交换律 (*l-x)*(*r-x)表示后一个向量在前一个向量的哪个方向 >0逆时针 <0顺时针 if((*l-x)*(*r-x)<0)return;//在凸包内 （l-x表示从x指向l的向量 r-x表示从x指向r的向量） now-=dis(*l,*r); while(r!=s.end()) { t=r;r++; if((*r-x)*(*t-x)>0)break; now-=dis(*t,*r); s.erase(t); } while(l!=s.begin()) { t=l;l--; if((*t-x)*(*l-x)>0)break; now-=dis(*t,*l); s.erase(t); } s.insert(x); l=r=t=s.find(x); l--;r++; now+=dis(*l,x)+dis(x,*r); } int main() { scanf("%d%d%d",&n,&x,&y); scanf("%d",&m); for(int i=1;i<=m;i++)scanf("%d%d",&p[i].x,&p[i].y); scanf("%d",&q); for(int i=1;i<=q;i++) { scanf("%d",&opt[i].op); if(opt[i].op==1) { scanf("%d",&opt[i].num); vis[opt[i].num]=1; } } P cap,S,T;cap.x=x;cap.y=y;S.x=0;S.y=0;T.x=n;T.y=0; now+=dis(S,cap);now+=dis(cap,T); s.insert(cap);s.insert(S);s.insert(T); for(int i=1;i<=m;i++) if(!vis[i])insert(p[i].x,p[i].y); for(int i=q;i>=1;i--) { if(opt[i].op==2){ opt[i].ans=now; } else insert(p[opt[i].num].x,p[opt[i].num].y); } for(int i=1;i<=q;i++) if(opt[i].op==2)printf("%.2lf\n",opt[i].ans); } ``` ### 2.旋转卡壳 如果要求凸包的直径 我们就需要借助这个算法了 引入概念 (1) 凸包的切线：就是过凸包上一个点且与凸包只有一个交点的直线 或者是凸包的一条边所在直线 (2) 凸包的直径：凸包上最远的两个点的距离 (3) 对踵点：两条平行切线对应的点对 可能是点--点 点--边 边--边 但实际情况下 点--边是最容易写的 所以只考虑这种情况 我们每次按顺序枚举一条边 然后找到对踵点（一定是离它最远的点）算一下边上两个点中的一个点到对踵点的距离 取$max$就行了 找对踵点可以通过叉积来找 因为底相同 高越大 面积越大 只需要找到最大的叉积对应的那个点就行了  但是这样枚举是$O(n^2)$的 考虑性质 对踵点和选择的边同方向变化 所有只需要维护一个指针作为对踵点 每次移动边时用叉积判断一下这个点要移动多少个就行了  [P1452 Beauty Contest](https://www.luogu.org/problem/P1452) 板子题 求凸包直径 ```cpp #include<bits/stdc++.h> typedef long long ll; using namespace std; const int maxn=52000; struct Vec{ ll x,y; Vec(){x=0;y=0;} friend Vec operator-(Vec a,Vec b){a.x-=b.x,a.y-=b.y;return a;} friend Vec operator+(Vec a,Vec b){a.x+=b.x,a.y+=b.y;return a;} friend double operator*(Vec a,Vec b){return a.x*b.y-b.x*a.y;} }O,p[maxn],s[maxn]; inline ll dis(Vec a,Vec b){return (b.y-a.y)*(b.y-a.y)+(b.x-a.x)*(b.x-a.x);} inline bool cmp(Vec a,Vec b){ Vec c=a-p[1],d=b-p[1]; return c*d==0?dis(O,c)<dis(O,d):c*d>0; } inline ll S(Vec a,Vec b,Vec c){ Vec d=b-a,e=c-b; return d*e; } int n,top; void get_convex() { sort(p+2,p+1+n,cmp); s[++top]=p[1]; for(int i=2;i<=n;i++) { while(top>1&&S(s[top-1],s[top],p[i])<=0)top--; s[++top]=p[i]; } } ll ans; void work() { int j=2;s[top+1]=s[1];//补成循环 if(top==2)ans=dis(s[1],s[2]); for(int i=1;i<=top;i++) { while(S(s[i],s[i+1],s[j])<S(s[i],s[i+1],s[j+1])){//找对踵点 j++;if(j>top)j-=top;//在凸包上旋转 } ans=max(ans,dis(s[i],s[j])); } } int main() { scanf("%d",&n); for(int i=1;i<=n;i++) { scanf("%lld%lld",&p[i].x,&p[i].y); if(p[i].y==p[1].y?p[i].x<p[1].x:p[i].y<p[1].y)swap(p[1],p[i]); } get_convex(); work(); printf("%lld",ans); } ``` 旋转卡壳还可以求一个凸包的最小外接矩形的面积 这时 我们不仅要像求直径那样维护最高点 同时还要维护一个最靠左边的点和最靠右的点 这里可以用点积根据点在当前底边上的投影来判断 越靠左投影越小  [P3187 [HNOI2007\]最小矩形覆盖](https://www.luogu.org/problem/P3187) 模板题 不过需要输出四个点的坐标 需要通过投影、单位向量什么的算一下 细节有些毒瘤而且还卡精度 然后就是注意最左边的点一开始要等于凸包上最后一个点 ```cpp #include<bits/stdc++.h> using namespace std; const int maxn=1e5+9; const double eps=1e-5; struct Vec{ double x,y; Vec(){x=0;y=0;} friend Vec operator-(Vec a,Vec b){a.x-=b.x,a.y-=b.y;return a;} friend Vec operator+(Vec a,Vec b){a.x+=b.x,a.y+=b.y;return a;} friend double operator*(Vec a,Vec b){return a.x*b.y-b.x*a.y;} friend double operator^(Vec a,Vec b){return a.x*b.x+a.y*b.y;} }O,p[maxn],s[maxn],v[maxn]; inline double dis(Vec a,Vec b){return sqrt((b.y-a.y)*(b.y-a.y)+(b.x-a.x)*(b.x-a.x));} inline bool cmp(Vec a,Vec b){ Vec c=a-p[1],d=b-p[1]; return c*d==0?dis(O,c)<dis(O,d):c*d>0; } inline double S(Vec a,Vec b,Vec c){//面积 Vec d=b-a,e=c-a; return d*e; } inline double L(Vec a,Vec b,Vec c){//投影 Vec d=b-a,e=c-a; if(!dis(O,d))return 0.0; return (d^e)/dis(O,d); } int n,top; void get_convex() { sort(p+2,p+1+n,cmp); s[++top]=p[1]; for(int i=2;i<=n;i++) { while(top>1&&S(s[top-1],s[top],p[i])<=0)top--; s[++top]=p[i]; } } double ans=1e9; void work() { int j=2,q=2,P=top;//j为最高点 q为右边 P为左边 矩形的底边与i-i+1边重合 投影即为边长 s[top+1]=s[1]; for(int i=1;i<=top;i++) { while(S(s[i],s[i+1],s[j])<S(s[i],s[i+1],s[j+1]))j=j%top+1; while(L(s[i],s[i+1],s[q])<L(s[i],s[i+1],s[q+1]))q=q%top+1; while(L(s[i],s[i+1],s[P])>L(s[i],s[i+1],s[P+1]))P=P%top+1; double ss=S(s[i],s[i+1],s[j]); double w=L(s[i],s[i+1],s[q])-L(s[i],s[i+1],s[P]); double a=dis(s[i],s[i+1]); if(w*ss/a+eps<ans) { ans=w*ss/a; Vec I,tmp,u=s[i+1]-s[i]; double mo=dis(O,u); I.x=u.x/mo;I.y=u.y/mo; double l=L(s[i],s[i+1],s[P]); I.x=I.x*l;I.y=I.y*l; v[1]=p[i]+I; I.x=u.x/mo;I.y=u.y/mo; l=L(s[i],s[i+1],s[q]); I.x=I.x*l;I.y=I.y*l; v[2]=p[i]+I; I.x=u.x/mo;I.y=u.y/mo; l=L(s[i],s[i+1],s[j]); I.x=I.x*l;I.y=I.y*l; tmp=p[i]+I;tmp=s[j]-tmp; v[3]=v[1]+tmp; v[4]=v[2]+tmp; } } } inline bool cmp2(Vec a,Vec b){ Vec c=a-v[1],d=b-v[1]; return c*d==0?dis(O,c)<dis(O,d):c*d>0; } int main() { scanf("%d",&n); for(int i=1;i<=n;i++) { scanf("%lf%lf",&p[i].x,&p[i].y); if(p[i].y==p[1].y?p[i].x<p[1].x:p[i].y<p[1].y)swap(p[1],p[i]); } get_convex(); work(); for(int i=1;i<=4;i++) if(v[i].y==v[1].y?v[i].x<v[1].x:v[i].y<v[1].y)swap(v[1],v[i]); sort(v+1,v+1+4,cmp2); printf("%.5lf\n",round(ans)); for(int i=1;i<=4;i++) printf("%.5lf %.5lf\n",fabs(round(v[i].x))<eps?0.00000:round(v[i].x),fabs(round(v[i].y))<eps?0.00000:round(v[i].y)); } ``` ### 3.半平面交 半平面 就是平面的一半 即一条直线的左边或者右边 半平面交就是一堆半平面的交集 如图 现在有一些直线  这就是它们的半平面交（所有直线的左边）  形象上看 半平面交是多边形的内核 在这个区域内 我们可以看到多边形中任何一个点 有$S\&I$算法 $O(nlogn))$可以求半平面交 先把所有线段极角排序 然后维护一个双端队列$q$和队列中相邻线段的交点$p$ 每次插入时比较最后两条线段的交点是否在当前线段的右边 如果在 则队尾不合法 需要删去 同理 比较队头两条线段的交点的位置关系删除队头 然后插入当前线段即可共线向量直接特判 保留较靠左边的一个 注意最后需要特判队头和队尾是否能接上 即将队头插入队尾时是否合法 因为有这种情况 无法闭合  判断交点时即这种情况  用双端队列为这种情况 这时队头队尾都要删除  [P2283 [HNOI2003]多边形](https://www.luogu.org/problem/P2283) 模板题 不过数据好像有些奇怪 [[CQOI2006]凸多边形](https://www.luogu.org/problem/P4196) 就是求所有多边形上的线段的半平面交的面积 最后求一下交点围成的凸多边形的面积即可 ```cpp #include<bits/stdc++.h> using namespace std; const int maxn=2000; const double eps=1e-7; struct Point{ double x,y; friend Point operator+(Point a,Point b){a.x+=b.x;a.y+=b.y;return a;} friend Point operator-(Point a,Point b){a.x-=b.x;a.y-=b.y;return a;} friend double operator*(Point a,Point b){return a.x*b.y-b.x*a.y;} Point(double xx=0,double yy=0){x=xx;y=yy;} }a[maxn],p[maxn]; struct Line{ Point x,y;//起点 终点 friend bool operator<(Line a,Line b){return atan2(a.y.y,a.y.x)<atan2(b.y.y,b.y.x);}//以原点为极点 }v[maxn],q[maxn]; int n,l,r,top; inline Point get_Intersect_Point(Line a,Line b){ Point c=a.x-b.x; double t=(b.y*c)/(a.y*b.y); a.y.x*=t;a.y.y*=t; return a.x+a.y; } bool right(Line a,Point o){ return a.y*(o-a.x)<=eps; } double Ans; void work() { l=r=0; sort(v+1,v+1+top); q[l=r=1]=v[1]; for(int i=2;i<=top;i++) { while(l<r&&right(v[i],p[r-1]))r--; while(l<r&&right(v[i],p[l]))l++; q[++r]=v[i]; if(fabs(q[r].y*q[r-1].y)<eps){//共线 r--; if(!right(q[r],v[i].x))q[r]=v[i]; } if(l<r)p[r-1]=get_Intersect_Point(q[r-1],q[r]);//求交点 } while(l<r&&right(q[l],p[r-1]))r--; if(r-l>1)p[r]=get_Intersect_Point(q[l],q[r]); double ans=0; for(int i=l;i<r;++i)ans+=p[i]*p[i+1]; ans+=p[r]*p[l]; ans=ans/2; printf("%.3lf\n",ans); } int main() { int t; scanf("%d",&t); while(t--) { scanf("%d",&n); for(int i=1;i<=n;i++)scanf("%lf%lf",&a[i].x,&a[i].y); for(int i=1;i<=n;i++)v[++top]=Line{a[i],a[i%n+1]-a[i]};//起点 方向 } work(); } ``` ### 4.闵可夫斯基和 定义:两个图形$A,B$的闵可夫斯基和$C=\{a+b|a\in A,b\in B\}$ 也就是两个图形边界上所有点两两相加 可以看作是一个图形绕着另一个图形的每个点转一圈然后最外圈的点形成的新图形  比较暴力的做法就是把所有边拆开 排个序然后首尾相连即可 起点是两个图形的起点相加 可以用归并排序 具体来说就是把两个凸包分别用差分向量表示 然后比较两个向量 如果$B$在$A$的顺时针方向则把$B$加进去 ```cpp void Minkovski(){ for(int i=1;i<n;++i)s1[i]=A[i+1]-A[i];s1[n]=A[1]-A[n]; for(int i=1;i<m;++i)s2[i]=B[i+1]-B[i];s2[m]=B[1]-B[m]; D[tot=1]=A[1]+B[1]; int p1=1,p2=1; while(p1<=n&&p2<=m)++tot,D[tot]=D[tot-1]+(s1[p1]*s2[p2]>=0?s1[p1++]:s2[p2++]); while(p1<=n)++tot,D[tot]=D[tot-1]+s1[p1++]; while(p2<=m)++tot,D[tot]=D[tot-1]+s2[p2++]; } ``` [P4557 [JSOI2018]战争](https://www.luogu.com.cn/problem/P4557) 应该算是模板题了 考虑$v$是可行向量集合 则满足$B+w=A$ 因此 $w=A-B$ 求出$A$和$-B$的和 （$B$凸包上每个点取相反数） 然后每次判断给出的向量在不在这个凸包内就行了 ```cpp #include<bits/stdc++.h> #define int long long using namespace std; const int maxn=100005; struct Vec{ int x,y; Vec(int x=0,int y=0):x(x),y(y){} friend Vec operator-(Vec a,Vec b){a.x-=b.x,a.y-=b.y;return a;} friend Vec operator+(Vec a,Vec b){a.x+=b.x,a.y+=b.y;return a;} friend int operator*(Vec a,Vec b){return a.x*b.y-b.x*a.y;} }O,bs,A[maxn],B[maxn],s[maxn],s1[maxn],s2[maxn],D[maxn]; int top; inline int dis(Vec a,Vec b){return 1ll*(b.y-a.y)*(b.y-a.y)+1ll*(b.x-a.x)*(b.x-a.x);} inline bool cmp1(Vec A,Vec B) {return A.y<B.y||(A.y==B.y&&A.x<B.x);} inline bool cmp(Vec a,Vec b){ Vec c=a-bs,d=b-bs;//按极点进行极角排序 return c*d==0?dis(O,c)<dis(O,d):c*d>0; } inline bool check(Vec a,Vec b,Vec c){ Vec d=b-a,e=c-b;//三个点必须保证连成的两个向量是逆时针方向即叉积>=0 return d*e<=0; } void convex(Vec *C,int &n){ sort(C+1,C+1+n,cmp1); bs=C[1];top=0; sort(C+2,C+1+n,cmp); s[++top]=C[1]; for(int i=2;i<=n;++i){ while(top>1&&check(s[top-1],s[top],C[i]))--top; s[++top]=C[i]; } for(int i=1;i<=top;i++)C[i]=s[i]; n=top;C[n+1]=C[1]; } int n,m,q,tot; void Minkovski(){ for(int i=1;i<n;++i)s1[i]=A[i+1]-A[i];s1[n]=A[1]-A[n]; for(int i=1;i<m;++i)s2[i]=B[i+1]-B[i];s2[m]=B[1]-B[m]; D[tot=1]=A[1]+B[1]; int p1=1,p2=1; while(p1<=n&&p2<=m)++tot,D[tot]=D[tot-1]+(s1[p1]*s2[p2]>=0?s1[p1++]:s2[p2++]); while(p1<=n)++tot,D[tot]=D[tot-1]+s1[p1++]; while(p2<=m)++tot,D[tot]=D[tot-1]+s2[p2++]; } int in(Vec qly){ if((qly-bs)*(D[1]-bs)>0||(qly-bs)*(D[tot]-bs)<=0)return 0; int pos=lower_bound(D+1,D+tot+1,qly,cmp)-D-1; return ((qly-bs)-(D[pos]-bs))*(D[pos%tot+1]-D[pos])<=0; } signed main(){ O=Vec(0,0); scanf("%lld%lld%lld",&n,&m,&q); for(int i=1;i<=n;++i){ int x,y;scanf("%lld%lld",&x,&y); A[i]=Vec(x,y); } convex(A,n); for(int i=1;i<=m;++i){ int x,y;scanf("%lld%lld",&x,&y); B[i]=Vec(-x,-y); } convex(B,m); Minkovski(); convex(D,tot); bs=D[1]; while(q--){ int x,y;scanf("%lld%lld",&x,&y); Vec qly=Vec(x,y); cout<<in(qly)<<endl; } } ```