计算几何模板合集

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
const double Pi=acos(-1.0);
const int N=3e5+10;
inline int dcmp(double x){return (x<-eps)?-1:(x>eps?1:0);}
inline double Abs(double x){return x*dcmp(x);}

namespace CG{
struct pt{
double x,y;pt(double _x=0,double _y=0){x=_x;y=_y;}
};
inline bool cmpx(const pt &a,const pt &b){return (a.x!=b.x)?a.x<b.x:a.y<b.y;}

typedef pt vec;
inline double Len(const vec &a){return sqrt(a.x*a.x+a.y*a.y);}//模长
inline double angle(const vec &a){return atan2(a.y,a.x);}
inline vec operator +(const vec &a,const vec &b){return vec(a.x+b.x,a.y+b.y);}
inline vec operator -(const vec &a,const vec &b){return vec(a.x-b.x,a.y-b.y);}
inline vec operator *(const vec &a,double b){return vec(a.x*b,a.y*b);}
inline vec operator /(const vec &a,double b){return vec(a.x/b,a.y/b);}

inline double operator *(const vec &a,const vec &b){return a.x*b.x+a.y*b.y;}//点积
inline double operator ^(const vec &a,const vec &b){return a.x*b.y-a.y*b.x;}//叉积

inline vec rotate(vec a,double theta){//将点或向量a绕原点(向量是顶点)逆时针旋转theta的弧度
double x=a.x*cos(theta)-a.y*sin(theta);
double y=a.x*sin(theta)+a.y*cos(theta);
return pt(x,y);
}
inline vec rotate_90(vec a){return pt(a.y,-a.x);}
inline pt rotate_P(pt a,pt b,double theta){return rotate(a-b,theta)+b;}//将点a绕点b逆时针旋转theta

//命名技巧：点P(point)，线段S(segment)，射线R(ray)，直线L(line)
struct line{
pt s,t;
line(pt _s=pt(0,0),pt _t=pt(0,0)){s=_s;t=_t;}
};
inline double maxx(const line &L){return max(L.s.x,L.t.x);}
inline double maxy(const line &L){return max(L.s.y,L.t.y);}
inline double minx(const line &L){return min(L.s.x,L.t.x);}
inline double miny(const line &L){return min(L.s.y,L.t.y);}
inline double ang(const line &L){return angle(L.t-L.s);}

inline bool operator <(const line &a,const line &b){
double a1=angle(a.t-a.s),a2=angle(b.t-b.s);
if(dcmp(a2-a1)!=0) return dcmp(a2-a1)>0;
return dcmp((b.t-a.s)^(a.t-a.s))>0;
}
inline bool operator ==(pt a,pt b){return (!dcmp(a.x-b.x))&&(!dcmp(a.y-b.y));}
inline double dis_PP(pt a,pt b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}//点a与点b间的距离
inline bool judge_PL(pt a,line b){return !dcmp((a-b.s)^(b.t-b.s));}//判断点是否在直线上
inline bool judge_PS(pt a,line b){return (!dcmp((a-b.s)^(b.t-b.s)))&&(dcmp((a-b.s)*(a-b.t))<=0);}//判断点是否在线段上

inline pt Footprint(pt a,line b){//点A关于直线ST的垂足
pt x=a-b.s,y=a-b.t,z=b.t-b.s;
double s1=x*z,s2=-1.0*(y*z);//分别求出AS，AT关于ST的投影
return b.s+z*(s1/(s1+s2));
}
inline pt mirror(pt a,line b){return a+(Footprint(a,b)-a)*2.0;}//点a关于直线b的对称点
inline double dis_PL(pt a,line b){return Abs((a-b.s)^(a-b.t))/Len(b.t-b.s); }//点与直线的距离:面积除以底边长
inline double dis_PS(pt a,line b){//点与线段的距离
pt x=a-b.s,y=a-b.t,z=b.t-b.s;
if(dcmp(x*z)<0) return Len(x);//距离左端点最近
if(dcmp(y*z)>0) return Len(y);//距离右端点最近
return dis_PL(a,b);
}
inline pt point_PS(pt a,line b){//点a在线段b上距离最近的点
pt x=a-b.s,y=a-b.t,z=b.t-b.s;
if(dcmp(x*z)<0) return b.s;//距离左端点最近
if(dcmp(y*z)>0) return b.t;//距离右端点最近
return Footprint(a,b);
}

inline pt cross_LL(line a,line b){//直线的交点
pt x=a.t-a.s,y=b.t-b.s,z=a.s-b.s;
return a.s+x*((y^z)/(x^y));
}
inline bool judge_cross_SL(line a,line b){//判断线段a与直线b是否相交
if(!dcmp((a.t-a.s)^(b.t-b.s))) return false;
return judge_PS(cross_LL(a,b),a);//看交点是否在线段上即可
}
inline bool judge_cross_SS(line a,line b){//判断两线段是否相交
if(maxx(a)<minx(b)||maxy(a)<miny(b)) return false;
if(maxx(b)<minx(a)||maxy(b)<miny(a)) return false;
double s1=(a.t-a.s)^(b.s-a.s),s2=(a.t-a.s)^(b.t-a.s);
double s3=(b.t-b.s)^(a.s-b.s),s4=(b.t-b.s)^(a.t-b.s);
return dcmp(s1)*dcmp(s2)<=0&&dcmp(s3)*dcmp(s4)<=0;//a的端点在b的两侧且b的端点在a的两侧
}
}
using namespace CG;

namespace Polygon{//多边形
struct polygon{
vector<pt> pts;
inline pt& operator [](int x){return pts[x];}
pt cur;
}
inline void clear(){pts.clear();}
inline int nxt(int x){return x<pts.size()-1?x+1:0;}
inline int include(pt p){//点在多边形上：0:在多边形外,1:在多边形内，2：在多边形的边上
int cnt=0;
for(int i=0;i<pts.size();++i){
pt s=pts[i],t=pts[nxt(i)];
line L=line(s,t);
if(judge_PS(p,L)) return 2;
if(dcmp(p.y-miny(L))>=0&&dcmp(maxy(L)-p.y)>0){
double nx=s.x+((p.y-s.y)/(t.y-s.y)*(t.x-s.x));
if(dcmp(nx-p.x)>0) cnt++;
}
}
return cnt&1;
}

inline double area(){//面积
double ans=0;
for(int i=1;i<pts.size()-1;++i) ans+=(pts[i]-pts[0])^(pts[nxt(i)]-pts[0]);
return Abs(ans)/2;
}
inline double peri(){//周长
double ans=0;
for(int i=0;i<pts.size();++i) ans+=dis_PP(pts[i],pts[nxt(i)]);
return ans;
}
inline bool Left(pt x,pt l,pt r){return dcmp((l-x)^(r-x))<=0; }//xl是否在xr左侧
inline void rever(){reverse(pts.begin(),pts.end());}//转顺时针为逆时针
inline int include_bs(pt p){//二分法判断点是否在多边形内
int n=pts.size();
if(!Left(pts[0],p,pts[1])) return 0;
if(!Left(pts[0],pts[n-1],p)) return 0;
if(judge_PS(p,line(pts[0],pts[1]))) return 2;
if(judge_PS(p,line(pts[0],pts[n-1]))) return 2;
int l=1,r=n-2,ans=1;
while(l<=r){
int mid=(l+r)>>1;
if(!Left(pts[0],pts[mid],p)) l=mid+1,ans=mid;
else r=mid-1;
}
if(!Left(pts[ans],p,pts[nxt(ans)])) return 0;
if(judge_PS(p,line(pts[ans],pts[nxt(ans)]))) return 2;
return 1;
}
inline void insert(pt p){pts.push_back(p);}
};
inline bool disjoint(polygon A,polygon B){//多边形AB是否相离
for(int i=0;i<A.pts.size();++i){
line x=line(A[i],A[A.nxt(i)]);
for(int j=0;j<B.pts.size();++j){
line y=line(B[j],B[B.nxt(j)]);
if(judge_cross_SS(x,y)) return false;
}
}
if(A.include_bs(B[0])) return false;
if(B.include_bs(A[0])) return false;
return true;
}
}
using namespace Polygon;

namespace Hull{// 凸包、旋转卡壳、半平面交、闵可夫斯基和
inline void Andrew(polygon &p){//Andrew算法求凸包
vector<pt> q(p.pts.size()<<1);
sort(p.pts.begin(),p.pts.end(),cmpx);
int top=-1,n=p.pts.size();
q[++top]=p[0];
for(int i=1;i<n;++i){
while(top&&dcmp((q[top-1]-q[top])^(p[i]-q[top]))>=0) top--;
q[++top]=p[i];
}
int now=top;
for(int i=n-2;i>=0;--i){
while(top>now&&dcmp((q[top-1]-q[top])^(p[i]-q[top]))>=0) --top;
q[++top]=p[i];
}
if(n>1) --top;
p.clear();
for(int i=0;i<=top;++i) p.insert(q[i]);
}
inline double S(const pt &x,const pt &y,const pt &z){return Abs((y-x)^(z-x));}
inline double diam(polygon &p){//旋转卡壳求凸包直径
int n=p.pts.size();double ans=0;
for(int i=0,j=1;i<n;++i){
for(;;j=p.nxt(j))
if(dcmp(S(p[j],p[i],p[p.nxt(i)])-S(p[p.nxt(j)],p[i],p[p.nxt(i)]))>=0) break;
ans=max(ans,dis_PP(p[j],p[i]));ans=max(ans,dis_PP(p[j],p[p.nxt(i)]));
}
return ans;
}
inline polygon SI(vector<line> q){//半平面交算法 S&I
int n=q.size();
sort(q.begin(),q.end());
vector<line> li(n+1),L(n+1);vector<pt> p(n+1);
int len=0;
for(int i=0;i<n;++i){
if(i>0&&!dcmp(ang(q[i])-ang(q[i-1]))) continue;
L[++len]=q[i];
}
int l=1,r=0;
for(int i=1;i<=len;++i){
while(l<r&&dcmp((L[i].t-p[r])^(L[i].s-p[r]))>0) --r;
while(l<r&&dcmp((L[i].t-p[l+1])^(L[i].s-p[l+1]))>0) ++l;
li[++r]=L[i];if(r>l) p[r]=cross_LL(li[r],li[r-1]);
}
while(l<r&&dcmp((li[l].t-p[r])^(li[l].s-p[r]))>0) --r;
while(l<r&&dcmp((li[r].t-p[l+1])^(li[r].s-p[l+1]))>0) ++l;

p[r+1]=cross_LL(li[r],li[l]),++r;
polygon P;
for(int j=l+1;j<=r;++j) P.insert(p[j]);
return P;
}
inline polygon merge(polygon A,polygon B){//闵可夫斯基和
int n=A.pts.size(),m=B.pts.size(),tot1=0,tot2=0;
vector<pt> p(n+1),q(m+1);
for(int i=1;i<n;++i) p[++tot1]=A[i]-A[i-1];p[++tot1]=A[0]-A[n-1];
for(int i=1;i<m;++i) q[++tot2]=B[i]-B[i-1];q[++tot2]=B[0]-B[m-1];
int i=1,j=1,tot=0;pt last=pt(0,0);
polygon C;C.pts.resize(n+m+1);C[0]=last=A[0]+B[0];
while(i<=n&&j<=m){
if(dcmp(p[i]^q[j])>=0) C[++tot]=last+p[i++],last=C[tot];
else C[++tot]=last+q[j++],last=C[tot];
}
while(i<=n) C[++tot]=last+p[i++],last=C[tot];
while(j<=m) C[++tot]=last+q[j++],last=C[tot];
Andrew(C);
return C;
}
}
using Hull::Andrew;
using Hull::diam;
using Hull::merge;

namespace Circle{//圆：三点确定一圆、最小圆覆盖
struct circle{
pt o;double r;
circle(pt _o=pt(0,0),double _r=0){o=_o;r=_r;}
circle(pt A,pt B,pt C){//三点确定一圆
pt D=(A+B)/2,E=(B+C)/2;
o=cross_LL(line(D,D+rotate_90(B-A)),line(E,E+rotate_90(C-B)));
r=dis_PP(o,A);
}
inline bool include(pt p){return dcmp(r-dis_PP(p,o))>=0;}//点在圆内
inline void print(int x){
printf("%.*lf\n",x,r);
printf("%.*lf %.*lf",x,o.x,x,o.y);
}
};
inline circle mincov(const vector<pt> &p){//最小圆覆盖
int n=p.size();
circle c=circle(0,0);
for(int i=0;i<n;++i){
if(!c.include(p[i])){
c=circle(p[i],0);
for(int j=0;j<i;++j){
if(!c.include(p[j])){
c=circle((p[i]+p[j])/2.0,dis_PP(p[i],p[j])/2.0);
for(int k=0;k<j;++k)
if(!c.include(p[k])) c=circle(p[i],p[j],p[k]);
}
}
}
}
return c;
}
}
using namespace Circle;

posted @ 2021-04-12 19:56  cjTQX  阅读(88)  评论(0编辑  收藏  举报