计算几何总结

基础的东西

一些定义

typedef double db;
struct poi{
    db x,y;
    poi(db X=0,db Y=0) {x=X,y=Y;}
    db len() { return sqrt(x*x+y*y); }
};
poi operator + (poi a,poi b) {return poi(a.x+b.x,a.y+b.y);}
poi operator - (poi a,poi b) {return poi(a.x-b.x,a.y-b.y);}
poi operator * (poi a,db p) {return poi(a.x*p,a.y*p);}
poi operator / (poi a,db p) {return poi(a.x/p,a.y/p);}

判断正负0

//判断正负0
int sign(db x)
{
    if(fabs(x)<eps) return 0;
    if(x>0) return 1;
    return -1;
}

 向量旋转

poi rotate(db c){return (poi){x*cos(c)-y*sin(c),x*sin(c)+y*cos(c)};}
    //逆时针旋转-pi~pi

判断线段和直线是否相交 2984. 线段

//跨立实验：判断一条线段的两端是否在另一条线段的两侧（两个端点与另一线段的叉积乘积为负）。需要正反判断两侧。
int panxdj(poi a,poi b,poi c,poi d)//ab线段，cd直线
{
    if(sign(c.x-d.x)==0&&sign(c.y-d.y)==0) return 0;//直线不能是一个点
    if(sign(cross(a-c,d-c)*cross(b-c,d-c))>0) return 0;
    return 1;
}

判断点在直线上

int cxd(poi a,poi b,poi p)
//p-a x p-b ==0 && p-a * p-b <=0
{
    if(sign(cross(p-a,p-b))==0&&sign(dot(p-a,p-b))<=0) return 1;
    return 0;
}

射线法判断点在多边形内

db dbarea()
{
    db res=0;
    for(int i=1;i<=n;i++)
    res+=cross(a[i],a[i+1]);
    return res;
}

int cdb(poi p)
{
    int cnt=0;
    for(int i=1;i<=n;i++)
    {
        if(cxd(a[i],a[i+1],p)) return 1;
        int d1=sign(a[i].y-p.y),d2=sign(a[i+1].y-p.y);
        int det=sign(cross(a[i]-p,a[i+1]-p));//cross
        if((det>=0&&d1<0&&d2>=0)||(det<=0&&d1>=0&&d2<0)) cnt++;
    }
    return cnt&1;
}


for(int i=1;i<=n;i++)
scanf("%lf %lf",&a[i].x,&a[i].y);
a[n+1]=a[1];
if(sign(dbarea())<0)
{
    reverse(a+1,a+n+1);
    a[n+1]=a[1];
}

判断两直线关系 

先判重合，再判平行，再求交点

db cross(poi a,poi b)
{
    return a.x*b.y-b.x*a.y;
}
poi zhixianjd(poi a,poi b,poi c,poi d)
{
    poi p;
    db s1=cross(c-a,d-a),s2=cross(d-b,c-b);//s1 ac*ad s2 bd*bc
    p=s1/(s1+s2)*(b-a)+a;
    return p;
}
int pandianxian(poi a,poi b,poi p)
{
    if(sign(cross(p-a,p-b))==0) return 1;
    return 0;
}
if(pandianxian(a,b,c)&&pandianxian(a,b,d))
{
    cout<<"LINE"<<endl;
    continue;
}
if(sign(cross(b-a,d-c))==0)
{
    cout<<"NONE"<<endl;
    continue;
}
poi p=zhixianjd(a,b,c,d);

 求凸包

 

poi dian[N];
int cmp(poi l,poi r)
{
    int det=sing(cross(l-dian[1],r-dian[1]));
    if(det) return det>0;
    return (l-dian[1]).len()<(r-dian[1]).len();
}
int n,m;
poi p[N];
void graham()
{
    int k=1;
    for(int i=2;i<=n;i++)
    if(dian[i].x<dian[k].x||(sing(dian[i].x-dian[k].x)==0&&dian[i].y<dian[k].y)) k=i;
    if(k^1) swap(dian[1],dian[k]);
    sort(dian+2,dian+n+1,cmp);
    //for(int i=1;i<=n;i++) cout<<dian[i].x<<" "<<dian[i].y<<endl;
    m=1;
    p[1]=dian[1];
    for(int i=2;i<=n;i++)
    {
        while(m-1&&sing(cross(p[m]-dian[i],p[m-1]-p[m]))>=0) m--;
        p[++m]=dian[i];
    }
}

 

半平面交

 

int l=1,r=1;
int cmp(qw l,qw r)
{
    db t1=atan2(l.b.y-l.a.y,l.b.x-l.a.x),t2=atan2(r.b.y-r.a.y,r.b.x-r.a.x);
    if(sing(t1-t2)!=0) return t1<t2;
    return sing(cross(l.b-l.a,r.a-l.a))<0;//r.a
}
db ang(qw xian)
{
    return atan2(xian.b.y-xian.a.y,xian.b.x-xian.a.x);
}
void bpingj()
{
    sort(line+1,line+n+1,cmp);
    q[1]=line[1];
    for(int i=2;i<=n;i++)
    {
        poi a=line[i].a,b=line[i].b;
        if(sing(ang(line[i])-ang(line[i-1]))==0) continue;
        while(l<r&&sing(cross(b-a,jd[r]-a))<=0) r--;
        while(l<r&&sing(cross(b-a,jd[l+1]-a))<=0) l++;
        q[++r]=line[i];
        if(l<r) jd[r]=zhixianjd(q[r-1].a,q[r-1].b,a,b);
    }
    while(l<r&&sing(cross(q[l].b-q[l].a,jd[r]-q[l].a))<=0) r--;
    jd[r+1]=zhixianjd(q[l].a,q[l].b,q[r].a,q[r].b);
    r++;
}

最小圆覆盖

cir get_cir(poi i,poi j,poi k)
{
    poi p1,p2;
    p1=(i+j)/2+(j-i).rotate(pi/2);
    p2=(i+k)/2+(k-i).rotate(pi/2);
    poi jd=zhixianjd((i+j)/2,p1,(i+k)/2,p2);
    cir c;
    c.a=jd;
    c.r=(jd-i).len();//i
    return c;
}
    
random_shuffle(dian+1,dian+n+1);
cir c=(cir){dian[1],0};
for(int i=2;i<=n;i++)
if(sing((dian[i]-c.a).len()-c.r)==1)
{
    c=(cir){dian[i],0};
    for(int j=1;j<i;j++)
    if(sing((dian[j]-c.a).len()-c.r)==1)
    {
        c=(cir){(dian[i]+dian[j])/2,(dian[i]-dian[j]).len()/2};
        for(int k=1;k<j;k++)
        if(sing((dian[k]-c.a).len()-c.r)==1)
        {
             c=get_cir(dian[i],dian[j],dian[k]);
        }
    }
}

 

三维凸包

#include<bits/stdc++.h>
using namespace std;
const int N=405;
typedef double db;
const db eps=1e-12;
struct poi{
    db x,y,z;
    poi(db X=0,db Y=0,db Z=0){x=X,y=Y,z=Z;}
    db len(){return sqrt(x*x+y*y+z*z);}
    db rand_eps(){return ((db)rand()/RAND_MAX-0.5)*eps;}
    void shake() { x+=rand_eps(),y+=rand_eps(),z+=rand_eps(); }
}dian[N];
poi operator + (poi a,poi b){return (poi){a.x+b.x,a.y+b.y,a.z+b.z};}
poi operator - (poi a,poi b){return (poi){a.x-b.x,a.y-b.y,a.z-b.z};}
db dot(poi a,poi b)
{
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
poi cross(poi a,poi b)
{
    return (poi){a.y*b.z-b.y*a.z,a.z*b.x-b.z*a.x,a.x*b.y-b.x*a.y};
}
struct plane{
    int v[3];
    poi a,b,c;
    poi normal()//求法向量
    {
        return cross(c-a,b-a);
    }
    db area()//求面积
    {
        return normal().len()/2;
    }
    int above(poi p)
    {
        return dot(p-a,normal())>=0;
    }
}pl[N],lb[N];
bool g[N][N];
int n,m;
void get_tubao3()
{
    pl[++m]=(plane){1,2,3,dian[1],dian[2],dian[3]};
    pl[++m]=(plane){3,2,1,dian[3],dian[2],dian[1]};
    for(int i=4;i<=n;i++)
    {
        int cnt=0;
        for(int j=1;j<=m;j++)
        {
            int t=pl[j].above(dian[i]);
            if(!t) lb[++cnt]=pl[j];
            for(int k=0;k<3;k++)
            g[pl[j].v[k]][pl[j].v[(k+1)%3]]=t;
        }
        for(int j=1;j<=m;j++)
        for(int k=0;k<3;k++)//不要写xxx+xxx==1 这样的话面就重复了，一开始加两个一样的面是为了让4一定有面做底，另一个面和4一起组成另外三个面
        if(g[pl[j].v[k]][pl[j].v[(k+1)%3]]&&!g[pl[j].v[(k+1)%3]][pl[j].v[k]])
            lb[++cnt]={pl[j].v[k],pl[j].v[(k+1)%3],i,dian[pl[j].v[k]],dian[pl[j].v[(k+1)%3]],dian[i]};
        m=cnt;
        for(int j=1;j<=cnt;j++) pl[j]=lb[j];
    }
}
int main()
{
    cin>>n;
    for(int i=1;i<=n;i++)
    scanf("%lf %lf %lf",&dian[i].x,&dian[i].y,&dian[i].z);
    for(int i=1;i<=n;i++) dian[i].shake();
    get_tubao3();
    db res=0;
    for(int i=1;i<=m;i++) res+=pl[i].area();
    printf("%.6lf",res);
    return 0;
}

 

旋转卡壳

 

    db res=0;
    for(int i=1,j=3;i<=m;i++)
    {
        poi a=p[i],b=p[i%m+1];
        while(sing(fabs(cross(a-b,a-p[j%m+1]))-fabs(cross(a-b,a-p[j])))>0) j=j%m+1;
        res=max(res,(p[j]-a).len());
        res=max(res,(p[j]-b).len());
    }

 

扫描线

区间合并

 

        ll num=0;
        int st=la[1].x,ed=la[1].y;
        for(int j=2;j<=cnt;j++)
        if(la[j].x<=ed) ed=max(ed,la[j].y);
        else res+=1ll*(ed-st)*(r-l+1),st=la[j].x,ed=la[j].y;
        num+=1ll*(ed-st)*(r-l+1);

 

自适应辛普森积分

db spi(db a,db b)
{
    return (f(a)+4*f((a+b)/2)+f(b))*(b-a)/6;
}
db xps(db a,db b)
{
    db s=spi(a,b),l=spi(a,(a+b)/2),r=spi((a+b)/2,b);
    if(fabs(s-l-r)<=eps) return s;
    return xps(a,(a+b)/2)+xps((a+b)/2,b);
}

f是自己定义的函数