Templates for XCPC - 计算几何基础/凸包/半平面交

点/向量

struct Point{
    double x,y;
    Point operator+(const Point& b)const{return {x+b.x,y+b.y};}
    Point operator-(const Point& b)const{return {x-b.x,y-b.y};}
    Point operator*(double p)const{return {x*p,y*p};}
};
typedef Point Vector;

直线

struct Line{
    Point p;
    Vector v;
    double ang;
    Line(){}
    Line(Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);}
};

叉积

double cross(Point a1,Point a2,Point b1,Point b2){
    return (a2.x-a1.x)*(b2.y-b1.y)-(b2.x-b1.x)*(a2.y-a1.y);
}

判断点是否在线段上

bool isOnSeg(Point p,Point A,Point B){
    if(dcmp(cross(p,A,p,B)))return 0;
    return dcmp(min(A.x,B.x)-p.x)<=0&&
           dcmp(p.x-max(A.x,B.x))<=0&&
           dcmp(min(A.y,B.y)-p.y)<=0&&
           dcmp(p.y-max(A.y,B.y))<=0;
}

判断点是否在多边形里

int isInPolygon(Point p,Point* s,int m){
    int wn=0;
    for(int i=1;i<=m;i++){
        if(isOnSeg(p,s[i],s[i+1]))return -1;
        int k=dcmp(cross(s[i],s[i+1],s[i],p));
        int d1=dcmp(s[i].y-p.y);
        int d2=dcmp(s[i+1].y-p.y);
        if(k>0&&d1<=0&&d2>0)wn++;
        if(k<0&&d2<=0&&d1>0)wn--;
    }
    return wn!=0;
}

判断线段是否有公共点

bool isSegsIntersect(Point A,Point B,Point C,Point D){
    int d1=dcmp(cross(A,B,A,C));
    int d2=dcmp(cross(A,B,A,D));
    int d3=dcmp(cross(C,D,C,A));
    int d4=dcmp(cross(C,D,C,B));
    if(d1*d2<0&&d3*d4<0)return 1;
    if(d1==0&&isOnSeg(C,A,B))return 1;
    if(d2==0&&isOnSeg(D,A,B))return 1;
    if(d3==0&&isOnSeg(A,C,D))return 1;
    if(d4==0&&isOnSeg(B,C,D))return 1;
    return 0;
}

Andrew 算法求凸包

int Andrew(int n,Point *p,Point *s){
    int m=0;
    sort(p+1,p+n+1,[&](Point a,Point b){
        return a.x<b.x||(a.x==b.x&&a.y<b.y);
    });
    for(int i=1;i<=n;i++){
        while(m>=2&&cross(s[m-1],s[m],s[m],p[i])<=0)m--;
        s[++m]=p[i];
    }
    int k=m;
    for(int i=n-1;i>=1;i--){
        while(m>k&&cross(s[m-1],s[m],s[m],p[i])<=0)m--;
        s[++m]=p[i];
    }
    if(m>1)m--;
    s[m+1]=s[1];
    return m;
}

半平面交（S&I 算法）

int dcmp(double x){
    if(fabs(x)<eps)return 0;
    return x<0?-1:1;
}
double cross(Vector a,Vector b){
    return a.x*b.y-a.y*b.x;
}
struct Line{
    Point p;
    Vector v;
    double ang;
    Line(){}
    Line(Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);}
};
bool cmp(const Line& a,const Line& b){
    if(dcmp(a.ang-b.ang)==0)return dcmp(cross(a.v,b.p-a.p))<0;
    return a.ang<b.ang;
}
bool onRight(Line L,Point p){
    return dcmp(cross(L.v,p-L.p))<=0;
}
Point getIntersection(Line a,Line b){
    Vector u=a.p-b.p;
    double t=cross(b.v,u)/cross(a.v,b.v);
    return a.p+a.v*t;
}
int n;
Line L[N],q[N];
Point p[N];
int main(){
    ios::sync_with_stdio(false);cin.tie(0);
    if(!(cin>>n))return 0;
    int cnt=0;
    for(int i=1;i<=n;i++){
        int m;cin>>m;
        vector<Point> poly(m);
        for(int j=0;j<m;j++)cin>>poly[j].x>>poly[j].y;
        for(int j=0;j<m;j++)L[++cnt]=Line(poly[j],poly[(j+1)%m]-poly[j]);
    }
    sort(L+1,L+cnt+1,cmp);
    int tot=0;
    //极角去重
    for(int i=1;i<=cnt;i++){
        if(i>1&&dcmp(L[i].ang-L[i-1].ang)==0)continue;
        L[++tot]=L[i];
    }
    int head=1,tail=0;
    for(int i=1;i<=tot;i++){
        //新加入的线把队尾的交点切掉了
        while(head<tail&&onRight(L[i],getIntersection(q[tail-1],q[tail])))tail--;
        //新加入的线把队首的交点切掉了
        while(head<tail&&onRight(L[i],getIntersection(q[head],q[head+1])))head++;
        q[++tail]=L[i];
        //判断平行反向导致闭合空间为空
        if(head<tail&&dcmp(cross(q[tail-1].v,q[tail].v))<=0){
            cout<<fixed<<setprecision(3)<<0.000<<"\n";
            return 0;
        }
    }
    //首尾互相更新去冗余
    while(head<tail&&onRight(q[head],getIntersection(q[tail-1],q[tail])))tail--;
    while(head<tail&&onRight(q[tail],getIntersection(q[head],q[head+1])))head++;
    //不足以构成多边形
    if(tail-head+1<3||dcmp(cross(q[tail].v,q[head].v))<=0){
        cout<<fixed<<setprecision(3)<<0.000<<"\n";
        return 0;
    }
    int m=0;
    for(int i=head;i<tail;i++)p[++m]=getIntersection(q[i],q[i+1]);
    p[++m]=getIntersection(q[tail],q[head]);//收尾相接的交点
    //计算多边形面积
    double ans=0;
    for(int i=1;i<=m;i++)ans+=cross(p[i],p[i%m+1]);
    ans=fabs(ans)/2.0;
    cout<<fixed<<setprecision(3)<<ans<<"\n";
    return 0;
}

旋转卡壳

求凸包直径

if(m==2){
    cout<<dist(s[1],s[2])<<'\n';
    return 0;
}
s[m+1]=s[1];
int ans=0,cur=3;
for(int i=1;i<=m;i++){
    while(sqr(s[i],s[i+1],s[(cur)%m+1])>=sqr(s[i],s[i+1],s[cur])){
        cur=cur%m+1;
    }
    ans=max(ans,max(dist(s[i],s[cur]),dist(s[i+1],s[cur])));
}
cout<<ans<<endl;

求凸包最小外接矩形

double ans=1e18;
int l,r,t;
int I,L,R,T;//边、最右、最左、最远点
int main(){
    cin>>n;
    for(int i=1;i<=n;i++){
        cin>>p[i].x>>p[i].y;  
    }
    Andrew();
    //最右、最左、最远点
    l=2,r=2,t=2;
    for(int i=1;i<=m;i++){
        // t：距离边 (i,i+1) 最远的点（控制矩形高度）
        while(sqr(s[i],s[i+1],s[t%m+1])>=sqr(s[i],s[i+1],s[t]))t=t%m+1;
        // r：沿方向 s[i]->s[i+1] 投影最靠前的点（控制矩形右侧边界）
        while(dot(s[i],s[i+1],s[i+1],s[r%m+1])>=dot(s[i],s[i+1],s[i+1],s[r]))r=r%m+1;
        if(i==1)l=r;// 防止劣化
        // l：沿方向 s[i+1]->s[i]（反向）投影最靠前的点（控制矩形左侧边界）
        while(dot(s[i+1],s[i],s[i],s[l%m+1])>=dot(s[i+1],s[i],s[i],s[l]))l=l%m+1;
        // 矩形面积 = 宽度 * 高度
        // 宽度 = dist(i,i+1) + 右边投影长度 + 左边投影长度
        // 高度 = 三角形面积*2 / 底边长 = sqr(i,i+1,t) / dist(i,i+1)
        double res=
        (
            dist(s[i],s[i+1])                      // 底边长度
            +dot(s[i],s[i+1],s[i+1],s[r])/dist(s[i],s[i+1])   // 右侧超出部分在底边上的投影
            +dot(s[i+1],s[i],s[i],s[l])/dist(s[i],s[i+1])     // 左侧超出部分在底边上的投影
        )*sqr(s[i],s[i+1],s[t])/dist(s[i],s[i+1]);  // 乘以高度（顶点t到底边的距离）
        if(res<ans){
            ans=res;
            I=i;L=l;R=r;T=t;      // 记录最优解时的四个指针
        }
    }
    printf("%.5f\n",ans);
    Point W[4]={};
    double len=dist(s[I],s[I+1]);               // 底边长度
    Vector dir=(s[I+1]-s[I])*(1.0/len);         // 底边方向的单位向量
    W[0]=s[I]+dir*dot(dir,s[R]-s[I]);           // 右下顶点：R 在底边上的投影
    W[3]=s[I]+dir*dot(dir,s[L]-s[I]);           // 左下顶点：L 在底边上的投影
    Vector norm={-dir.y,dir.x};                  // 底边的法向量（逆时针旋转90°）
    double len2=sqr(s[I],s[I+1],s[T])/len;      // 矩形高度（T 到底边的距离）
    W[1]=W[0]+norm*len2;                         // 右上顶点
    W[2]=W[3]+norm*len2;                         // 左上顶点
    for(int i=0;i<4;i++){
        printf("%.5f %.5f\n",W[i].x,W[i].y);                    
    }
    return 0;
}