计算几何

本文最后更新于 2024年9月12日 晚上

计算几何

1. 求凸包

Graham扫描法

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
//以左下点为极点求凸包周长
 
#include<bits/stdc++.h>
using namespace std;
 
struct point{
    double x,y;
    point friend operator -(point a,point b)
    {return {a.x-b.x,a.y-b.y};}
}p[105],s[105];//p为所有点的集合,s为凸包点集
double dis(point a,point b)
{
    point c=a-b;
    return sqrt(c.x*c.x+c.y*c.y);
}
double cross(point a,point b)
{
    return a.x*b.y-a.y*b.x;
}
//按照极角(polar angle)从小到大排序(以 p1为极点)
//极角相同的点按照到的距离从小到大排序。
int cmp(point a,point b)
{
    double x=cross(a-p[1],b-p[1]);
 
    if(x>0) return 1;
    if(x==0&&dis(a,p[1])<dis(b,p[1])) return 1;
    return 0;
}
double multi(point p1,point p2,point p3)//计算向量叉积
{
    return cross(p2-p1,p3-p1);
}
int main()
{
    int N;
    while(scanf("%d",&N),N)
    {
        for(int i=1;i<=N;i++) cin>>p[i].x>>p[i].y;
 
        if(N==1)
        {
            printf("0.00\n");
            continue;
        }
        else if(N==2)
        {
            printf("%.2lf\n",dis(p[1],p[2]));
            continue;
        }
 
        int k=1;
        for(int i=2;i<=N;i++)
        if(p[i].y<p[k].y||(p[i].y==p[k].y&&p[i].x<p[k].x))k=i;
        swap(p[1],p[k]);//找到处于最左下方的点,赋给p[1]
 
        sort(p+2,p+1+N,cmp);//对点按极角排序,p[1]为极点
 
        s[1]=p[1];
        s[2]=p[2];
        int cnt=2;
        //while循环把发现不是凸包顶点的点移除出去,因为当逆时针遍历凸包时，我们应该在每个顶点向左转
        //因此当while循环发现在一个顶点处没有向左转时，就把该顶点移除出去。
        for(int i=3;i<=N;i++)//求凸包
        {
            while(cnt>=2&&multi(s[cnt-1],s[cnt],p[i])<=0) cnt--;
            s[++cnt]=p[i];
        }
        double sum=0;
        for(int i=1;i<cnt;i++)//求凸包周长
            sum+=dis(s[i],s[i+1]);
        printf("%.2f\n",sum+dis(s[1],s[cnt]));
        /*
        double area=0;
        for(int i=2;i<=cnt-1;i++)//求凸包面积，这里因为就是凸多边形所以可以用某个顶点作为顶点分割三角形
            area+=fabs(multi(s[1],s[i],s[i+1]));
        printf("%d\n",int(area/100));
        */
        /*标准求多边形面积,别漏了最后一个和起点的叉积
        for(int i=1;i<=cnt;i++)//求凸包面积
            area+=cross(s[i],s[i==cnt?1:i+1]);
        if(area<0)
            area=-area;
        */
    }
    return 0;
}

分治法：求解上凸包和下图包合并

  1. 按横坐标从小到大（再总坐标从小到大）排序，首先排序后的起终点必然在凸包里
  2. 连接p1pn，这条直线将点集分成了上下两个部分，则在这两个部分中分别求得上、下凸包
  3. 对于上凸包，找到上部分离直线最远的点pmax,连接p1max,pmaxpn，则直线p1pmax右将上部分分成了两个凸包部分
  4. 重复，求下凸包类似
  5. 判断点位于直线的左侧还是右侧可以用叉积
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
//求能围住所有的点且距离所有点的距离>=L的围墙的长度
 
#include<iostream>
#include<cmath>
#include<cstdio>
#include<algorithm>
using namespace std;
 
struct Point
{
    double x, y;

    Point operator-(Point & p)
    {
        Point t;
        t.x = x - p.x;
        t.y = y - p.y;
        return t;
    }

    double cross(Point p)
    {
        return x*p.y - p.x*y;
    }

    double dist(Point & p)
    {
        return sqrt((x - p.x)*(x - p.x) + (y - p.y)*(y - p.y));
    }
};
 
bool cmp(Point p1, Point p2)//先横坐标从小到大,再纵坐标从小到大
{
    if (p1.x != p2.x)
        return p1.x < p2.x;
    return p1.y < p2.y;
}
 
Point point[1005];
int convex[1005];
int N, L;
 
//返回极点个数+1
int getConvexHull()
{
    sort(point, point + N, cmp);
    int temp;
    int total = 0;
    for (int i = 0; i < N; i++)//求解下凸包,<=关系则如果存在共线的点只取端点
    {
        while (total > 1 && (point[convex[total - 1]] - point[convex[total - 2]]).cross(point[i] - point[convex[total - 1]]) <= 0)
            total--;
        convex[total++] = i;
    }
    //求解上凸包
    temp = total;
    for (int i = N - 2; i >= 0; i--)
    {
        while (total > temp && (point[convex[total - 1]] - point[convex[total - 2]]).cross(point[i] - point[convex[total - 1]]) <= 0)
            total--;
        convex[total++] = i;
    }
    return total;//返回组成凸包的点的个数，实际上多了一个,是起点(开始与最后都是起点),所以组成凸包的点个数是total-1
}
 
int main()
{
    scanf("%d%d", &N, &L);
    for (int i = 0; i < N; i++)
        scanf("%lf%lf", &point[i].x, &point[i].y);

    int num = getConvexHull();
    double ans = 0.0;
    for (int i = 0; i < num - 1; i++)
        ans += point[convex[i]].dist(point[convex[i + 1]]);
    ans += 3.14159265358927 * 2 * L;
    printf("%.0f\n", ans);
    return 0;
}


2. 判断线段关系

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
#include<bits/stdc++.h>
using namespace std;

struct Point
{
    double x, y;
};

//算叉乘
double cross(Point A, Point B, Point C)
{
    return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x);
}


bool is_intersect(Point A, Point B, Point C, Point D)
{
    if (min(A.x, B.x) <= max(C.x, D.x) &&
        min(C.x, D.x) <= max(A.x, B.x) &&
        min(A.y, B.y) <= max(C.y, D.y) &&
        min(C.y, D.y) <= max(A.y, B.y) &&
        cross(A, B, C) * cross(A, B, D) <= 0 &&
        cross(C, D, A) * cross(C, D, B) <= 0)
        return true;

    return false;
}
bool is_para(Point A, Point B, Point C, Point D){
    return (B.y-A.y)*(D.x-C.x)==(B.x-A.x)*(D.y-C.y);
}

3. 旋转卡壳求凸包直径

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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstdio>

using namespace std;

const int MAX = 200005;
const double eps = 1e-7;

struct Point {
    double x, y;

    Point operator+(const Point &b) const {
        return {x + b.x, y + b.y};
    }

    Point operator-(const Point &b) const {
        return {x - b.x, y - b.y};
    }

    double operator^(const Point &b) const {
        return x * b.y - y * b.x;
    }

    bool operator<(const Point &b) const {
        if (x != b.x)
            return x < b.x;
        return y < b.y;
    }
};

Point p[MAX];
Point s[MAX];
int top;

void selMin(int n);
int cmp(Point a, Point b);
bool equal(double a, double b);
double dis(Point a, Point b);
void graham(int n);
double s_sqr(Point a, Point b, Point c);
double diameter();

int main() {
    int n;
    cin >> n;
    for (int i = 0; i < n; i++)
        cin >> p[i].x >> p[i].y;
    selMin(n);
    sort(p + 1, p + n, cmp);
    graham(n);
    printf("%.6f", sqrt(diameter())) ;
    return 0;
}

void selMin(int n) {
    Point Min = p[0];
    int IDMin = 0;
    for (int i = 0; i < n; i++)
        if (p[i] < Min) {
            Min = p[i];
            IDMin = i;
        }
    swap(p[0], p[IDMin]);
}

int cmp(Point a, Point b) {
    double x = (a - p[0]) ^ (b - p[0]);
    if (x > 0)
        return 1;
    if (equal(x, 0) && (dis(a, p[0]) < dis(b, p[0])))
        return 1;
    return 0;
}

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

void graham(int n) {
    top = 1;
    s[0] = p[0];
    s[1] = p[1];
    for (int i = 2; i < n; i++) {
        while (top > 1 && ((p[i] - s[top]) ^ (s[top - 1] - s[top])) <= 0)
            top--;
        s[++top] = p[i];
    }
}

double s_sqr(Point a, Point b, Point c) {
    return fabs((a - b) ^ (c - b));
}

double diameter() {
    double diam = 0;
    int j = 2;
    s[++top] = s[0];
    if (top < 3)
        return dis(s[0], s[1]);
    for (int i = 0; i < top - 1; i++) {
        while (s_sqr(s[i], s[i + 1], s[j]) <
               s_sqr(s[i], s[i + 1], s[(j + 1) % top]))
            j = (j + 1) % top;
        diam = max(diam, max(dis(s[i], s[j]), dis(s[i + 1], s[j])));
    }
    return diam;
}
bool equal(double a, double b){
    return fabs(a - b) < eps;
}