计算几何模板

Posted on | In

点、线段

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
const double eps=1e-8;
inline int sgn(double x){
    if(fabs(x)<eps) return 0;
    else if(x<0) return -1;
    else return 1;
}

struct Point{
    double x,y;
    Point(){
        x=0;
        y=0;
    }
    Point(double a,double b):x(a),y(b){}
    inline Point operator-(const Point& b)const{
        return Point(x-b.x,y-b.y);
    }
    inline Point operator+(const Point& b)const{
        return Point(x+b.x,y+b.y);
    }
    inline double dot(const Point& b)const{
        return x*b.x+y*b.y;
    }
    inline double cross(const Point&b,const Point& c)const{
        return (b.x-x)*(c.y-y)-(c.x-x)*(b.y-y);
    }
    double operator^(Point a){return x*a.y-y*a.x;}
    double operator*(Point a){return x*a.x+y*a.y;}
};

double dist(const Point& a, const Point& b){
    return sqrt((b-a).dot(b-a));
}



//得线段交点
Point LineCross(const Point &a, const Point &b, const Point &c, const Point &d) {
     double u = a.cross(b, c), v = b.cross(a, d);
     return Point((c.x * v + d.x * u) / (u + v), (c.y * v + d.y * u) / (u + v));
}

//判断线段和直线相交 poj3304
//s-e是直线
bool check(const Point& s,const Point& e){
    if(sgn(dist(s,e))==0) return false;	//注意特殊情况，线段退化为点时会造成共线的误判
    for(int i=0;i<n;i++){
        if(sgn(s.cross(lines[i][0],e))*sgn(s.cross(lines[i][1],e))>0)
           return false;
    }
    return true;
}

凸包

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
Point stk[maxn];
bool operator<(const Point& a,const Point& b){
    return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
int Andrew(){	//逆时针扫描
    sort(ps,ps+n);
    int top=0;
    for(int i=0;i<n;i++){
        while(top>1&&sgn((stk[top-1]-stk[top-2])^(ps[i]-stk[top-1]))<=0) top--;
        stk[top++]=ps[i];
    }
    int k=top;
    for(int i=n-2;i>=0;i--){
        while(top>k&&sgn((stk[top-1]-stk[top-2])^(ps[i]-stk[top-1]))<=0) top--;
        stk[top++]=ps[i];
    }
    return top;	//凸包点数
}//结果点集逆时针存放在stk中

半平面交

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
//poj3335求多边形有无核

struct Line {
    Point s, e;
    Line() {}
    Line(Point a, Point b) {s = a, e = b;}
};
double getAngle(Line a){
    return atan2(a.e.y-a.s.y,a.e.x-a.s.x);
}
bool operator<(Line a,Line b){
    double A=getAngle(a);
    double B=getAngle(b);
    if(fabs(A-B)<eps) return ((a.e-a.s)^(b.e-a.s))>0;  //把靠右的排在前面
    else return A<B;
}
Point getIntersectPoint(Line a,Line b){     //用线段b的两点坐标组合成交点坐标
    double u=a.s.cross(a.e,b.e);
    double v=a.e.cross(a.s,b.s);
    return Point((v*b.e.x+u*b.s.x)/(u+v),(v*b.e.y+u*b.s.y)/(u+v));
}


bool onRight(Line a,Line b,Line c){ //判断bc交点是否在a的右边
    Point o=getIntersectPoint(b,c);
    if(((a.e-a.s)^(o-a.s))<0)   //这里不能取等，也许是因为能留下一条边就尽量留？这样tail-head的合法判定更易满足
        return true;
    return false;
}

Line lines[maxn];
Line que[maxn];
Point ps[maxn];
int n;

bool HalfPlaneIntersect(){
    sort(lines,lines+n);
    int head=0,tail=0,cnt=0;
    for(int i=0;i<n-1;i++){
        if(fabs(getAngle(lines[i])-getAngle(lines[i+1]))<eps){
            continue;
        }
        lines[cnt++]=lines[i];
    }
    lines[cnt++]=lines[n-1];

    for(int i=0;i<cnt;i++){
        while(tail-head>1&&onRight(lines[i],que[tail-1],que[tail-2]))   //用新边优化队尾
            tail--;
        while(tail-head>1&&onRight(lines[i],que[head],que[head+1])) //用新边优化队首
            head++;
        que[tail++]=lines[i];
    }
    while(tail-head>1&&onRight(que[head],que[tail-1],que[tail-2]))
        tail--;
    while(tail-head>1&&onRight(que[tail-1],que[head],que[head+1]))
        head++;
    if(tail-head<3)
        return false;
    return true;
}

bool judge() {  //判断点的顺序是逆时针
  double ans = 0;
  for(int i=1; i<n-1; i++) {
    ans += ((ps[i]-ps[0])^(ps[i+1]-ps[0]));
  }
  return ans>0;
}
int t;
int main(){
    cin>>t;
    while(t--){
        cin>>n;
        for(int i=0;i<n;i++)
            cin>>ps[i].x>>ps[i].y;
        if(judge()){
            for(int i=0;i<n;i++)
                lines[i]=Line(ps[i],ps[(i+1)%n]);
        }else{
            for(int i=0;i<n;i++){
                lines[i]=Line(ps[(i+1)%n],ps[i]);
            }
        }
        if(HalfPlaneIntersect()){
            printf("YES\n");
        }else
            printf("NO\n");
    }
    return 0;
}

旋转卡壳

求多边形最长直径

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
//注意它的思想，实际情况中可能求的不是直径
double rotating_caliper(Point* ps){
    double res=0;
    if(n==2){
        res=dist(ps[0],ps[1]);
    }else{
        int j=2;
        ps[n]=ps[0];	//省的再对i取模……
        for(int i=0;i<n;i++){
            while(ps[i].cross(ps[i+1],ps[j])<ps[i].cross(ps[i+1],ps[j+1]))
                j=(j+1)%n;
            res=max(res,max(dist(ps[i],ps[j]),dist(ps[i+1],ps[j])));
        }
    }
    return res;
}