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

1 //叉积,可以用来判断方向和求面积
2 
3 doublecross(Point a,Point b,Point c){
4 
5     return (c.x-a.x)*(b.y-a.y) - (b.x-a.x)*(c.y-a.y);
6 
7 }
8 
9  
10 
11  
12 
13 //求多边形的面积
14 
15 double S(Point p[],intn){
16 
17     double ans = 0;
18 
19     p[n] = p[0];
20 
21     for(int i=1;i<n;i++)
22 
23        ans += cross(p[0],p[i],p[i+1]);
24 
25     if(ans < 0) ans = -ans;
26 
27     return ans / 2.0;
28 
29 }

1 //求多边形的重心
2 
3 Point grabity(Point p[],intn){
4 
5 Point G;
6 
7     double sum_area=0;
8 
9     for(int i=2;i<n;i++){
10 
11         double area = cross(p[0],p[i-1],p[i]);
12 
13         sum_area+=area;
14 
15         G.x+=(p[0].x+p[i-1].x+p[i].x)*area;
16 
17         G.y+=(p[0].y+p[i-1].y+p[i].y)*area;
18 
19 }
20 
21     G.x=G.x/3/sum_area,G.y=G.y/3/sum_area;
22 
23     returnG;
24 
25 }

1 /**
2 
3 求二维凸包Andrew算法，将所有的点按x小到大(x相等，y小到大)排序
4 
5 删去重复的点，得到一个序列p1,p2...,然后把p1,p2放入凸包中，从p3
6 
7 开始当新点再前进方向左边时(可以用叉积判断方向)继续，否则，依次
8 
9 删除最近加入凸包的点，直到新点再左边。
10 
11 **/
12 
13  
14 
15 int ConvexHull(Point *p,int n,Point *stack){
16 
17     sort(p,p+n);
18 
19     n=unique(p,p+n)-p;
20 
21     int m=0;
22 
23     for(int i=0;i<n;i++){//如果不希望凸包的边上有输入的点则把两个等号去掉
24 
25         while(m>1&&cross(stack[m-2],p[i],stack[m-1])<=0) m--;
26 
27         stack[m++]=p[i];
28 
29 }
30 
31     int k=m;
32 
33     for(int i=n-2;i>=0;i--){
34 
35         while(m>k&&cross(stack[m-2],p[i],stack[m-1])<=0)m--;
36 
37         stack[m++]=p[i];
38 
39 }
40 
41     if(n>1) m--;
42 
43     returnm;
44 
45 }

//判断符号，提高精度

int dcmp(doublex){

    if(fabs(x)<eps) return 0;

    else return x < 0 ? -1 : 1;

}

structPoint{

    doublex,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, doublep){

    return Vector(A.x*p, A.y*p);

}

Vector operator / (Vector A, doublep){

    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)
}

doubleLength(Vector A){

    return sqrt(Dot(A, A));    //向量的长度
}

doubleAngle(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,doublea){

    return Vector (A.x*cos(a)-A.y*cos(a),A.x*sin(a)+A.y*cos(a));

}

 

double trans(doubleang){

    return ang/180*acos(-1.0);

}

//旋转卡壳求凸包的直径，平面距离最远的点对的距离

double rotatint_calipers(Point *p,intn){

    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])));

    }

    returnans;

}

1 double rotating_calipers(Point *p,intn){
2 
3     int k = 1;
4 
5     double ans = 0x7FFFFFFF;
6 
7     p[n] = p[0];
8 
9     for(int i=0;i<n;i++){
10 
11         while(fabs(cross(p[i],p[i+1],p[k])) < fabs(cross(p[i],p[i+1],p[k+1])))
12 
13              k = (k+1) %n;
14 
15         double tmp = fabs(cross(p[i],p[i+1],p[k]));
16 
17         double d   = dist(p[i],p[i+1]);
18 
19         ans = min(ans,tmp/d);
20 
21 }
22 
23     returnans;
24 
25 }

1 //求线段的中垂线  
2 
3 inline Line getMidLine(const Point &a, const Point &b) {  
4 
5     Point mid = (a +b);  
6 
7     mid.x/=2.0;  
8 
9     mid.y/=2.0;  
10 
11     Point tp = b-a;  
12 
13     return Line(mid, mid+Point(-tp.y, tp.x));  
14 
15 } 

structLine  

{  

    Point s,e;  

    Line(){}  

    Line(Point _s,Point _e)  

    {  

        s =_s;  

        e =_e;  

    }  

    bool operator ==(Line v)  

    {  

        return (s == v.s)&&(e ==v.e);  

    }  

    voidinput()  

    {  

        s.input();  

        e.input();  

    }  

    //两线段相交判断  

    //2 规范相交  

    //1 非规范相交  

    //0 不相交  

    intsegcrossseg(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 不相交  

    intlinecrossseg(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);  

    }  

};