计算几何（基础）

点积和叉积

如果涉及的运算都是整数，可以将返回值改为int。

叉积可用于求平行四边形面积和三角形面积。

1
2
3
4
5
6
7
8
9
10
11
12
13
double cro(point a,point b)
{
    return a.x*b.y-a.y*b.x;
}
double dot(point a,point b)
{
    return a.x*b.x+a.y*b.y;
}
int sgn(double x)
{
    if(fabs(x)<eps)return 0;
    else return x<0?-1:1;
}

向量旋转

设向量为$(x, y)$，绕起点逆时针旋转，设旋转角度为$\theta$，那么旋转后的向量为：
$x^{‘} = xcos \theta - ysin \theta$
$y^{‘} = xsin \theta + ycos \theta$

以点代向量：

1
2
3
4
point Roate(point A, double rad)
{
    return point(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}

求单位法向量

1
2
3
4
point Normal(point A)
{
    return point(-A.y/length(A), A.x/length(A));
}

点到直线的投影

对于给定的三个点p1，p2，p。发现点p到直线p1p2的投影点x：

alt text

使用点积和叉积可以轻松推出来这个表达式，但是我忘了咋推的了（

1
2
3
4
5
point PointToLine(point p,point p1,point p2)
{
    double k = dot(p2-p1,p-p1)/length2(p2-p1);
    return p1+(p2-p1)*k;
}

粘贴一个完整代码，内含point类的运算符重载：

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
#include <bits/stdc++.h>
using namespace std;
struct point
{
    double x,y;
    point():x(0),y(0) {}
    point(double x,double y):x(x),y(y) {}
    point(const point &p):x(p.x),y(p.y) {}
};
double length2(point a)
{
    return a.x*a.x + a.y*a.y;
}
double dot(point a,point b)
{
    return a.x*b.x + a.y*b.y;
}
double cro(point a,point b)
{
    return a.x*b.y - a.y*b.x;
}
point operator+(point a,point b)
{
    return point(a.x+b.x,a.y+b.y);
}
point operator-(point a,point b)
{
    return point(a.x-b.x,a.y-b.y);
}
point operator*(point a,double b)
{
    return point(a.x*b,a.y*b);
}
point PointToLine(point p,point p1,point p2)
{
    double k = dot(p2-p1,p-p1)/length2(p2-p1);
    return p1+(p2-p1)*k;
}
istream &operator>>(istream &in,point &p)
{
    in>>p.x>>p.y;
    return in;
}
ostream &operator<<(ostream &out,const point p)
{
    printf("%.10lf %.10lf",p.x,p.y);
    return out;
}
int main()
{
    point p1,p2,p;
    cin>>p1>>p2;
    int q;
    cin>>q;
    while(q--)
    {
        cin>>p;
        cout<<PointToLine(p,p1,p2)<<endl;
    }
    return 0;
}

点关于直线的对称点

求对称点，首先要求出点到直线的投影点，然后用中点坐标公式一推就可以了。

1
2
3
4
5
point SymmetricPoint(point p, point p1, point p2)
{
    point q = PointToLine(p, p1, p2);
    return point(2*q.x-p.x,2*q.y-p.y);
}

判断两个向量的位置关系

判断向量$p_0p_2$与向量$p_0p_1$的位置关系：

（为了美观，将point类的重载函数、点积叉积函数和main函数删去了）

若$\mathbf{A} \times \mathbf{B} > 0 $则B在A的逆时针方向。

若$\mathbf{A} \times \mathbf{B} < 0 $则B在A的顺时针方向。

若叉积等于0，则说明三点共线，使用点积判断同方向还是反方向。

这个方法可以判断点在直线的左侧还是右侧，或者在直线上。

点是p2，线段是p0p1

alt text

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
#define eps 1e-8
typedef struct point
{
    double x;
    double y;
    point(double xx=0.0,double yy=0.0):x(xx),y(yy) {}
};
enum {COUNTER_CLOCKWISE,CLOCKWISE,ONLINE_BACK,ONLINE_FRONT,ON_SEGMENT};
double length(point a);

int LinePointStatus(point p0,point p1,point p2)
{
    if(fabs(cro(p1-p0,p2-p0))<eps)
    {
        if(dot(p1-p0,p2-p0)<0)
        {
            return ONLINE_BACK;
        }
        else
        {
            if(p2==p0||p2==p1||length(p1-p0)>length(p2-p0))
                return ON_SEGMENT;
            else return ONLINE_FRONT;
        }
    }
    else if(dot(p1-p0,p2-p0)>0)
    {
        return COUNTER_CLOCKWISE;
    }
    else
    {
        return CLOCKWISE;
    }
}

判断两条直线关系

点p1和p2构成一条直线，点p3和p4构成另一条直线：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
int LineRelation(point p1,point p2,point p3,point p4)
{
    point v1(p2-p1), v2(p4-p3);
    if(sgn(dot(v1,v2))==0)
    {
        // 两直线正交
        return 1;
    }
    if(sgn(cro(v1,v2))==0)
    {
        // if(LinePointStatus(p1,p2,p3)==ON_SEGMENT||
        //     LinePointStatus(p1,p2,p3)==ONLINE_BACK||
        //     LinePointStatus(p1,p2,p3)==ONLINE_FRONT) return 0;  //两直线重合
        return 2; //两直线平行
    }
    return 0;   //普通的相交
}

两条直线的交点

依然是点积叉积推出来的

1
2
3
4
5
6
7
8
9
10
point operator/(point a,double b)
{
    return point(a.x/b,a.y/b);
}
point CrossPoint(point p1,point p2,point p3,point p4)
{
    double s1 = cro(p2-p1,p3-p1);
    double s2 = cro(p2-p1,p4-p1);
    return point(p3.x*s2-p4.x*s1,p3.y*s2-p4.y*s1)/(s2-s1);
}

两条线段是否相交

先求出交点，再判断是否在这两个线段上。线段有重合部分有的时候也被视为相交了（但我个人不认同），看情况而定吧。

这题NEUOJ上需要开long double才能过。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
bool PointOnLine(point p,point p1,point p2)
{
    int status=LinePointStatus(p1,p2,p);
    return status==ON_SEGMENT||status==ONLINE_BACK||status==ONLINE_FRONT;
}
bool SegmentCross(point p1,point p2,point p3,point p4)
{
    if(LineRelation(p1,p2,p3,p4)==2)
    {
        // 线段平行的时候特殊判断
        int st1=LinePointStatus(p1,p2,p3),st2=LinePointStatus(p1,p2,p4);
        if(st1==ON_SEGMENT||st2==ON_SEGMENT)return true;
        if(st1==ONLINE_BACK&&st2==ONLINE_FRONT)return true;
        if(st2==ONLINE_BACK&&st1==ONLINE_FRONT)return true;
        return false;
    }
    point p(CrossPoint(p1,p2,p3,p4));
    return LinePointStatus(p1,p2,p)==ON_SEGMENT&&LinePointStatus(p3,p4,p)==ON_SEGMENT;
}

几个距离

点到直线距离

使用叉积计算出平行四边形面积，再除以底就是了。

1
2
3
4
double pointLineDistance(point p,point p1,point p2)
{
    return fabs(cro(p-p1,p2-p1))/length(p1-p2);
}

点到线段距离

如果垂足在线段内部，那么就是点到直线距离；反之则是点到某一端点的距离。

1
2
3
4
5
double pointSegmentDistance(point p,point p1,point p2)
{
    if(LinePointStatus(p1,p2,PointToLine(p,p1,p2))==ON_SEGMENT)return pointLineDistance(p,p1,p2);
    return min(length(p-p1),length(p-p2));
}

线段与线段的距离

就是4个点到线段的距离，取最小

1
2
3
4
5
6
double SegmentsDistance(point p1,point p2,point p3,point p4)
{
    if(SegmentCross(p1,p2,p3,p4))return 0;
    return min(min(pointSegmentDistance(p1,p3,p4),pointSegmentDistance(p2,p3,p4)),
        min(pointSegmentDistance(p3,p1,p2),pointSegmentDistance(p4,p1,p2)));
}

多边形的面积

随便找一个点，然后与n多边形的点两两组成n个三角形，面积就为三角形的叉积和。注意求叉积的时候不能取绝对值，因为有负面积。适用于所有多边形包括凹多边形。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
double PolygonArea(point*p,int n)
{
    double ans=0;
    for(int i=0;i<n;i++)ans+=cro(p[i],p[(i+1)%n]);
    return ans/2;
}
const int N=1e5+1;
point p[N];

int main()
{
    int n;
    cin>>n;
    for(int i=0;i<n;i++)cin>>p[i];
    printf("%.10lf\n",PolygonArea(p,n));
    return 0;
}

点与多边形的关系

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int PointInPolygon(point t,point*p,int n)
{
    for(int i=0;i<n;i++)if(t==p[i])return 1;
    for(int i=0;i<n;i++)
    {
        if(LinePointStatus(p[i],p[(i+1)%n],t)==ON_SEGMENT)return 1;
    }
    int cnt=0;
    for(int i=0;i<n;i++)
    {
        int j = (i+1)%n;
        int c = sgn(cro(t-p[j],p[i]-p[j]));
        int u = sgn(p[i].y-t.y);
        int v = sgn(p[j].y-t.y);
        if(c>0&&u<0&&v>=0)cnt++;
        if(c<0&&u>=0&&v<0)cnt--;
    }
    if(cnt!=0)
    {
        // 点在内部
        return 2;
    }
    return 0;
}

判断多边形是不是凸多边形

直接每个角的叉积大于等于0就是凸多边形。注意叉积的两个向量有顺序，下面的代码假设多边形以逆时针顺序给出：

1
2
3
4
5
6
7
8
bool isConvex(point*p,int n)
{
    for(int i=0;i<n;i++)
    {
        if(cro(p[(i+1)%n]-p[i],p[(i-1+n)%n]-p[i])<0)return false;
    }
    return true;
}

凸包

给定一些点，求能把所有点包含在内的面积最小的多边形。相当于一个橡皮筋把整个多边形包住的最小收缩。

使用Andrew算法求凸包前，要将点进行排序：按照x从小到大的顺序，x相同则按照y从小到大。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void AndrewConvex(point* num,point* ans,int n,int&m)
{
    m=0;
    for(int i=0;i<n;i++)
    {
        while(m>1&&sgn(cro(ans[m-1]-ans[m-2],num[i]-ans[m-1]))<0)m--;
        ans[m++]=num[i];
    }
    int j=m;
    for(int i=n-2;i>=0;i--)
    {
        while(m>j&&sgn(cro(ans[m-1]-ans[m-2],num[i]-ans[m-1]))<0)m--;
        ans[m++]=num[i];
    }
    if(n>1)m--;
}

LuoguP2742 圈奶牛：

求面积最小的多边形，能把所有奶牛圈起来，输出这个多边形的周长。

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
#include <bits/stdc++.h>
using namespace std;
struct point
{
    double x,y;
    point():x(0),y(0) {}
    point(double x,double y):x(x),y(y) {}
};
#define eps 1e-8
int sgn(double x)
{
    if(fabs(x)<0)return 0;
    return x<0?-1:1;
}
double cro(point a,point b)
{
    return a.x*b.y-a.y*b.x;
}
point operator-(point a,point b)
{
    return point(a.x-b.x,a.y-b.y);
}
void AndrewConvex(point* num,point* ans,int n,int&m)
{
    m=0;
    for(int i=0;i<n;i++)
    {
        while(m>1&&sgn(cro(ans[m-1]-ans[m-2],num[i]-ans[m-1]))<0)m--;
        ans[m++]=num[i];
    }
    int j=m;
    for(int i=n-2;i>=0;i--)
    {
        while(m>j&&sgn(cro(ans[m-1]-ans[m-2],num[i]-ans[m-1]))<0)m--;
        ans[m++]=num[i];
    }
    if(n>1)m--;
}
double length(point a,point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool operator<(const point a,const point b)
{
    return a.x<b.x||(fabs(a.x-b.x)<eps&&a.y<b.y);
}
const int N=1e6+1;
point p1[N],p2[N];

int main()
{
    int n;
    cin>>n;
    for(int i=0;i<n;i++)cin>>p1[i].x>>p1[i].y;
    int m=0;
    sort(p1,p1+n);
    AndrewConvex(p1,p2,n,m);
    double ans=0;
    for(int i=0;i<m;i++)ans+=length(p2[i],p2[(i+1)%m]);
    printf("%.2lf",ans);
}

旋转卡壳

用于求凸包直径，即多边形的两点距离最大值。

注意得先求凸包，再求旋转卡壳。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
double trangleSquare(point a, point b, point c)
{
    // 三角形面积的2倍，其实应该就是平行四边形
    return fabs(cro(b - a, c - a));
}
double ConvexDiameter(point* p, int n)
{
    if (n == 2)return length(p[0]-p[1]);
    int cur = 0;
    double ans = 0.0;
    for (int i = 0; i < n; i++)
    {
        while (trangleSquare(p[i], p[(i + 1) % n], p[cur]) <= trangleSquare(p[i], p[(i + 1) % n], p[(cur + 1) % n]))
            cur = (cur + 1) % n;
        ans = max(ans, max(length(p[i]-p[cur]), length(p[i + 1]-p[cur])));
    }
    return ans;
}

（注：洛谷上这个代码TLE了2个样例，但很多题解的思路都是这个，我怀疑是卡常了，因为有的旋转卡壳题这个代码能过）

牛客多校2024-8-13 H题

给出两个多边形，一个固定一个随便变换（平移、旋转、翻转），但必须两个多边形至少一个交点，求B的轨迹所形成的周长。

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
#include <bits/stdc++.h>
using namespace std;
struct point
{
    double x,y;
    point():x(0),y(0) {}
    point(double x,double y):x(x),y(y) {}
};
#define eps 1e-8
int sgn(double x)
{
    if(fabs(x)<0)return 0;
    return x<0?-1:1;
}
double cro(point a,point b)
{
    return a.x*b.y-a.y*b.x;
}
point operator-(point a,point b)
{
    return point(a.x-b.x,a.y-b.y);
}
void AndrewConvex(point* num,point* ans,int n,int&m)
{
    m=0;
    for(int i=0;i<n;i++)
    {
        while(m>1&&sgn(cro(ans[m-1]-ans[m-2],num[i]-ans[m-1]))<0)m--;
        ans[m++]=num[i];
    }
    int j=m;
    for(int i=n-2;i>=0;i--)
    {
        while(m>j&&sgn(cro(ans[m-1]-ans[m-2],num[i]-ans[m-1]))<0)m--;
        ans[m++]=num[i];
    }
    if(n>1)m--;
}
double length(point a)
{
    return sqrt((a.x*a.x)+(a.y*a.y));
}
bool operator<(const point a,const point b)
{
    return a.x<b.x||(fabs(a.x-b.x)<eps&&a.y<b.y);
}
const int N=1e6+5;
point A[N], B[N], C[N];
double trangleSquare(point a, point b, point c)
{
    return fabs(cro(b - a, c - a));
}
double ConvexDiameter(point* p, int n)
{
    if (n == 2)return length(p[0]-p[1]);
    int cur = 0;
    double ans = 0.0;
    for (int i = 0; i < n; i++)
    {
        //while (Line2Point(p[i], p[(i + 1) % n], p[cur]) <= Line2Point(p[i], p[(i + 1) % n], p[(cur + 1) % n]))
        while (trangleSquare(p[i], p[(i + 1) % n], p[cur]) <= trangleSquare(p[i], p[(i + 1) % n], p[(cur + 1) % n]))
            cur = (cur + 1) % n;
        ans = max(ans, max(length(p[i]-p[cur]), length(p[i + 1]-p[cur])));
    }
    return ans;
}
istream& operator>>(istream& in, point& p)
{
    in >> p.x >> p.y;
    return in;
}

void solve()
{
    #define pi  3.1415926535897932384626433
    int n;
    cin >> n;
    for (int i = 0; i < n; i++)cin >> A[i];
    double ans = 0;
    for (int i = 0; i < n; i++)ans += length(A[i]-A[(i + 1) % n]);
    int m;
    cin >> m;
    for (int i = 0; i < m; i++)cin >> B[i];
    int k = 0;
    sort(B, B + m);
    AndrewConvex(B, C, m, k);
    ans += 2 * pi * ConvexDiameter(C, k);
    printf("%.15lf\n", ans);
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int t = 1;
    cin>>t;
    while (t--)solve();
    cout.flush();
    return 0;
}

最近点对

一堆点，求距离最近的两个点的距离。其实就是分治法。

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
double distance(point a,point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool PointCMPy(const point a,const point b)
{
    return a.y<b.y;
}
double PointsMinDistance(point*p,point*pt,int l,int r)
{
    double ans=inf;
    if(l==r)return ans;
    if(l+1==r)return distance(p[l],p[r]);
    int mid=(l+r)/2;
    ans=min(PointsMinDistance(p,pt,l,mid),PointsMinDistance(p,pt,mid+1,r));
    int k=0;
    for(int i=l;i<=r;i++)
    {
        if(fabs(p[mid].x-p[i].x)<=ans)pt[k++]=p[i];
    }
    sort(pt,pt+k,PointCMPy);
    for(int i=0;i<k;i++)
    {
        for(int j=i+1;j<k;j++)
        {
            if(pt[j].y-pt[i].y>=ans)break;
            ans=min(ans,distance(pt[i],pt[j]));
        }
    }
    return ans;
}

多边形内的晶格点

给定一个多边形，您的任务是计算多边形内部及其边界上的格点数量。 格点是坐标为整数的点。

皮克定理：指一个计算点阵中顶点在格点上的多边形面积公式，该公式可以表示为S=a+b÷2－1，其中a表示多边形内部的点数，b表示多边形落在格点边界上的点数，S表示多边形的面积。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
void CountSquareLattice(point*p,int n,long long&inside,long long&boundary)
{
    boundary=abs(CountBoundaryLattice(p,n));
    inside=abs(AllSquare(p,n)+1-boundary/2);
}
long long CountBoundaryLattice(point*p,int n)
{
    long long ans=n;
    for(int i=0;i<n;i++)ans+=CountLineLattice(p[i],p[(i+1)%n]);
    return ans;
}
long long CountLineLattice(point p1,point p2)
{
    if(p1==p2)return 0;
    if(p1.x==p2.x)return abs(p1.y-p2.y)-1;
    if(p1.y==p2.y)return abs(p1.x-p2.x)-1;
    if(abs(p1.x-p2.x)%abs(p1.y-p2.y)==0)return abs(p1.y-p2.y)-1;
    if(abs(p1.y-p2.y)%abs(p1.x-p2.x)==0)return abs(p1.x-p2.x)-1;
    return gcd(abs(p1.x-p2.x),abs(p1.y-p2.y))-1;
}
long long gcd(long long a,long long b)
{
    return b?gcd(b,a%b):a;
}

直接调用CountSquareLattice函数就可以获得多边形内部的格点数和边缘上的格点数。