计算几何的基本代码模板

代码来自刘汝佳的训练指南

计算π值

const double PI=acos(-1);

角度化弧度

double torad(double deg)
{
    return deg/180*PI;
}

坐标点和运算

struct Point
{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y){}
};
bool operator < (const Point &a,const Point &b)//排序用，按照x的坐标从小到大，如果x相同，那么按照y从小到大
{
    return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
Vector operator +(Vector A,Vector B){return Vector(A.x+B.x,A.y+B.y);}//坐标点相加
Vector operator -(Vector A,Vector B){return Vector(A.x-B.x,A.y-B.y);}//相减

三态函数，用来减少精度问题

int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    else
        return x<0?-1:1;
}

点积

double Dot(Vector A,Vector B) {return A.x*B.x+A.y*B.y;}

叉积

#define Vector Point
double Cross(Vector A,Vector B)
{
    return A.x*B.y-A.y*B.x;
}

andrew算法 求凸包

int ConvexHull(Point *p,int n,Point *ch)//p是所有点，n所有点的个数，ch里面记录形成凸包的点，返回凸包点的个数
{
    sort(p,p+n);
    int m=0;
    for(int i=0;i<n;i++)
    {
        while(m>1&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)
            m--;
        ch[m++]=p[i];
    }
    int k=m;
    for(int i=n-2;i>=0;i--)
    {
        while(m>k&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)
            m--;
        ch[m++]=p[i];
    }
    if(n>1)
        m--;
    return m;
}

某点绕着原点逆时针旋转

#define Vector Point
Vector Rotate(Vector A,double rad)//旋转
{
    return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

计算多边形有向面积

double PolyonArea(Point* p,int n)//多边形面积
{
    double area=0;
    for(int i=1;i<n-1;i++)
        area+=Cross(p[i]-p[0],p[i+1]-p[0]);
    return area/2;
}

判断点是否在线段上面

bool OnSegment(Point p,Point a1,Point a2)
{
    return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}

判断点是否在多边形内部

int isPointInPolygon(Point p,Point *ch,int n)
{
    int wn=0;
    for(int i=0;i<n;i++)
    {
        if(OnSegment(p,ch[i],ch[(i+1)%n]))return -1;
        int k=dcmp(Cross(ch[(i+1)%n]-ch[i],p-ch[i]));
        int d1=dcmp(ch[i].y-p.y);
        int d2=dcmp(ch[(i+1)%n].y-p.y);
        if(k>0&&d1<=0&&d2>0) wn++;
        if(k<0&&d2<=0&&d1>0) wn--;
    }
    if(wn!=0)
        return 1;
    return 0;
}

线段是否相交（不考虑端点的情况）

bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2)
{
    double c1 = Cross(a2-a1,b1-a1), c2 = Cross(a2-a1,b2-a1),
    c3 = Cross(b2-b1,a1-b1), c4 = Cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}

两点距离的平方

double Sq_dist(Point a,Point b)
{
    return (b.x-a.x)*(b.x-a.x)+(b.y-a.y)*(b.y-a.y);
}

旋转卡壳模板

double rotating_calipers(Point *ch,int n)//ch为顺序凸包上的点,n为点的个数
{
    int q=1;
    double ans=0;
    ch[n]=ch[0];
    for(int p=0;p<n;p++)
    {
        while(Cross(Vector(ch[q+1]-ch[p+1]),Vector(ch[p]-ch[p+1]))>Cross(Vector(ch[q]-ch[p+1]),Vector(ch[p]-ch[p+1])))
            q=(q+1)%n;
        ans=max(ans,max(Sq_dist(ch[p],ch[q]),Sq_dist(ch[p+1],ch[q+1])));
    }
    return ans;//返回凸包直径的平方
}

//由旋转卡壳的原理可以知道，
//以凸包上某条边为底，顺时针或逆时针以凸包上的点为顶点，得到的三角形的面积是一个有且仅有一个最大值的函数
//由于三角形面积先增加，再减少
//while循环的目的就是找到当前以当前顶点q+1的三角形面积，比前以一个顶点q形成一个三角形面积小的时候顶点
//则该店顶点q为距离边{p,p+1}最远的顶点

to be continue~