计算几何常用的函数/方法

时间:2022-01-12 06:16:13

(一)求多边形的面积(用叉积计算)

代码如下:

//叉积,可以用来判断方向和求面积
double cross(Point a,Point b,Point c){
    return (c.x-a.x)*(b.y-a.y) - (b.x-a.x)*(c.y-a.y);
}


//求多边形的面积
double S(Point p[],int n){
    double ans = 0;
    p[n] = p[0];
    for(int i=1;i<n;i++)
       ans += cross(p[0],p[i],p[i+1]);
    if(ans < 0) ans = -ans;
    return ans / 2.0;
}

(二)求多边形的重心

代码如下:

//求多边形的重心
Point grabity(Point p[],int n){
    Point G;
    double sum_area=0;
    for(int i=2;i<n;i++){
        double area = cross(p[0],p[i-1],p[i]);
        sum_area+=area;
        G.x+=(p[0].x+p[i-1].x+p[i].x)*area;
        G.y+=(p[0].y+p[i-1].y+p[i].y)*area;
    }
    G.x=G.x/3/sum_area,G.y=G.y/3/sum_area;
    return G;
}

(三)andrew算法求凸包

代码如下:

/**
求二维凸包Andrew算法,将所有的点按x小到大(x相等,y小到大)排序
删去重复的点,得到一个序列p1,p2...,然后把p1,p2放入凸包中,从p3
开始当新点再前进方向左边时(可以用叉积判断方向)继续,否则,依次
删除最近加入凸包的点,直到新点再左边。
**/

int ConvexHull(Point *p,int n,Point *stack){
    sort(p,p+n);
    n=unique(p,p+n)-p;
    int m=0;
    for(int i=0;i<n;i++){//如果不希望凸包的边上有输入的点则把两个等号去掉
        while(m>1&&cross(stack[m-2],p[i],stack[m-1])<=0) m--;
        stack[m++]=p[i];
    }
    int k=m;
    for(int i=n-2;i>=0;i--){
        while(m>k&&cross(stack[m-2],p[i],stack[m-1])<=0)m--;
        stack[m++]=p[i];
    }
    if(n>1) m--;
    return m;
}


(四)比较函数提高精度:

//判断符号,提高精度
int dcmp(double x){
    if(fabs(x)<eps) return 0;
    else return x < 0 ? -1 : 1;
}

(五)向量/以及常见运算重载

struct Point{
    double x,y;
    Point():x(0),y(0){}
    Point(double _x,double _y):x(_x),y(_y){}
    bool operator <(const struct Point &tmp) const{
        if(x==tmp.x) return y<tmp.y;
        return x<tmp.x;
    }
};

typedef Point Vector;
Vector operator + (Vector A, Vector B){
    return Vector(A.x+B.x, A.y+B.y);
}
Vector operator - (Point A, Point B){
    return Vector(A.x-B.x, A.y-B.y);
}
Vector operator * (Vector A, double p){
    return Vector(A.x*p, A.y*p);
}
Vector operator / (Vector A, double p){
    return Vector(A.x/p, A.y/p);
}
bool operator == (Vector A,Vector B){
    return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0;
}
double Dot(Vector A, Vector B){//向量相乘
    return A.x*B.x + A.y*B.y;  //a*b*cos(a,b)
}
double Length(Vector A){
    return sqrt(Dot(A, A));    //向量的长度
}
double Angle(Vector A, Vector B){
    return acos(Dot(A, B) / Length(A) / Length(B));    //向量的角度
}
double Cross(Vector A, Vector B){//叉积
    return A.x*B.y - A.y*B.x;
}
/**
向量(x,y) 绕起点逆时针旋转a度。
x' = x*cosa - y*sina
y' = x*sina + y*cosa
**/
Vector Rotate(Vector A,double a){
    return Vector (A.x*cos(a)-A.y*cos(a),A.x*sin(a)+A.y*cos(a));
}

double trans(double ang){
    return ang/180*acos(-1.0);
}

(六)旋转卡壳求凸包的直径,平面最远的点对

代码如下:

//旋转卡壳求凸包的直径,平面距离最远的点对的距离
double rotatint_calipers(Point *p,int n){
    int k=1;
    int ans = 0;
    p[n]=p[0];
    for(int i=0;i<n;i++){
        while(fabs(Cross(p[i+1],p[k],p[i]))<fabs(Cross(p[i+1],p[k+1],p[i])))
            k=(k+1)%n;
        ans = max(ans,max(dis(p[i],p[k]),dis(p[i+1],p[k])));
    }
    return ans;
}

(七)旋转卡壳求凸包的宽度,即找一组距离最近的平行线似的凸包的点在两根线的内侧

代码如下:

double rotating_calipers(Point *p,int n){
    int k = 1;
    double ans = 0x7FFFFFFF;
    p[n] = p[0];
    for(int i=0;i<n;i++){
        while(fabs(cross(p[i],p[i+1],p[k])) < fabs(cross(p[i],p[i+1],p[k+1])))
             k = (k+1) % n;
        double tmp = fabs(cross(p[i],p[i+1],p[k]));
        double d   = dist(p[i],p[i+1]);
        ans = min(ans,tmp/d);
    }
    return ans;
}

(八)求线段的中垂线

//求线段的中垂线  
inline Line getMidLine(const Point &a, const Point &b) {  
    Point mid = (a + b);  
    mid.x/=2.0;  
    mid.y/=2.0;  
    Point tp = b-a;  
    return Line(mid, mid+Point(-tp.y, tp.x));  
} 

(九)直线相关

struct Line  
{  
    Point s,e;  
    Line(){}  
    Line(Point _s,Point _e)  
    {  
        s = _s;  
        e = _e;  
    }  
    bool operator ==(Line v)  
    {  
        return (s == v.s)&&(e == v.e);  
    }  
    void input()  
    {  
        s.input();  
        e.input();  
    }  
    //两线段相交判断  
    //2 规范相交  
    //1 非规范相交  
    //0 不相交  
    int segcrossseg(Line v)  
    {  
        int d1 = sgn((e-s)^(v.s-s));  
        int d2 = sgn((e-s)^(v.e-s));  
        int d3 = sgn((v.e-v.s)^(s-v.s));  
        int d4 = sgn((v.e-v.s)^(e-v.s));  
        if( (d1^d2)==-2 && (d3^d4)==-2 )return 2;  
        return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) ||  
            (d2==0 && sgn((v.e-s)*(v.e-e))<=0) ||  
            (d3==0 && sgn((s-v.s)*(s-v.e))<=0) ||  
            (d4==0 && sgn((e-v.s)*(e-v.e))<=0);  
    }  
    //直线和线段相交判断  
    //-*this line   -v seg  
    //2 规范相交  
    //1 非规范相交  
    //0 不相交  
    int linecrossseg(Line v)  
    {  
        int d1 = sgn((e-s)^(v.s-s));  
        int d2 = sgn((e-s)^(v.e-s));  
        if((d1^d2)==-2) return 2;  
        return (d1==0||d2==0);  
    }  
};