【算法模板】计算几何：点类与线类

计算几何（Computational Geometry）是计算机科学中的一个分支，专注于研究几何问题的算法和数据结构。它主要涉及几何对象（如点、线段、多边形等）的表示、操作和分析，旨在高效地解决几何问题。

基本结构和类型定义

#include <bits/stdc++.h>
using namespace std;

// 使用别名定义浮点数类型为Real
using Real = double;
// 定义精度常量，用于浮点数比较
const Real EPS = 1e-8;
1
2
3
4
5
6
7

解释：

• #include <bits/stdc++.h>: 引入所有标准库，以便简化代码书写。
• using Real = double;: 定义浮点数类型的别名，方便后续修改。
• const Real EPS = 1e-8;: 定义用于浮点数比较的精度常量，避免因精度问题导致的错误比较。

基本概念

点和向量的定义

首先定义表示点和向量的结构体Point。点和向量在几何中是基本的表示形式，通常用于描述位置和方向。

// 定义点和向量的结构体
struct Point {
    Real x, y;  // 点的x和y坐标

    // 构造函数，默认初始化为原点(0, 0)
    Point(Real x = 0, Real y = 0) : x(x), y(y) {}

    // 向量的加法
    Point operator+(const Point &p) const { return {x + p.x, y + p.y}; }

    // 向量的减法
    Point operator-(const Point &p) const { return {x - p.x, y - p.y}; }

    // 向量的乘法 (数乘)
    Point operator*(Real c) const { return {x * c, y * c}; }

    // 向量的除法 (数除)
    Point operator/(Real c) const { return {x / c, y / c}; }

    // 判断两个点是否相等，考虑浮点数精度
    bool operator==(const Point &p) const { return fabs(x - p.x) < EPS && fabs(y - p.y) < EPS; }

    // 比较运算符，用于点的排序
    bool operator<(const Point &p) const { return x < p.x - EPS || (fabs(x - p.x) < EPS && y < p.y - EPS); }

    // 输出点的坐标，格式为(x, y)
    friend ostream& operator<<(ostream &os, const Point &p) { return os << "(" << p.x << ", " << p.y << ")"; }
};
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

解释：

• Point结构体用于表示二维空间中的一个点或一个向量。它包含两个成员变量x和y，分别表示点的横坐标和纵坐标。
• 我们为该结构体定义了一些基本的运算符重载，如加法、减法、乘法和除法，以方便进行向量运算。
• 运算符==和<用于点的比较，考虑到浮点数运算的精度问题，使用了EPS作为误差容忍度。
• 重载了<<运算符，以便可以直接输出点的坐标。

向量的基本操作

向量是点的扩展，可以用来进行更多的几何操作，例如点积、叉积、向量旋转等。

// 计算向量的点积（内积），a · b = a.x * b.x + a.y * b.y
Real dot(const Point &a, const Point &b) { 
    return a.x * b.x + a.y * b.y; 
}

// 计算向量的叉积，a × b = a.x * b.y - a.y * b.x，结果为标量
Real cross(const Point &a, const Point &b) { 
    return a.x * b.y - a.y * b.x; 
}

// 计算向量的模长（欧几里得距离）
Real length(const Point &v) { 
    return sqrt(dot(v, v)); 
}

// 计算两点之间的距离
Real distance(const Point &a, const Point &b) { 
    return length(a - b); 
}

// 判断两向量是否平行
bool areParallel(const Point &a, const Point &b) { 
    return fabs(cross(a, b)) < EPS; 
}

// 判断两向量是否垂直
bool arePerpendicular(const Point &a, const Point &b) { 
    return fabs(dot(a, b)) < EPS; 
}

// 计算两向量的夹角
Real angleBetween(const Point &a, const Point &b) {
    return acos(dot(a, b) / (length(a) * length(b)));
}

// 向量旋转任意角度（逆时针）
Point rotate(const Point &v, Real theta) {
    return {v.x * cos(theta) - v.y * sin(theta), v.x * sin(theta) + v.y * cos(theta)};
}

// 点绕点旋转任意角度
Point rotateAround(const Point &p, const Point &origin, Real theta) {
    return origin + rotate(p - origin, theta);
}
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

解释：

• dot函数计算两个向量的点积，它是标量，反映了两个向量的相似程度。
• cross函数计算两个向量的叉积，它的结果也是标量，反映了两个向量之间的面积（方向）。
• length函数计算向量的模长，实际上是计算点到原点的距离。
• distance函数计算两点之间的距离，这是通过计算它们之间的向量的模长得到的。
• areParallel和arePerpendicular分别判断两向量是否平行或垂直。
• angleBetween计算两个向量之间的夹角，结果是弧度值。
• rotate函数用于将向量逆时针旋转指定的角度theta。
• rotateAround则是让点p绕点origin旋转指定角度theta。

线与线段

在几何中，线段和直线都是由两点确定的，因此我们使用两个Point来表示一条直线或线段。

struct Line {
    Point a, b;  // 直线的两个端点
    
    // 构造函数
    Line(Point a = {}, Point b = {}) : a(a), b(b) {}
    
    // 输出直线的端点，格式为<点a, 点b>
    friend ostream& operator<<(ostream &os, const Line &l) { return os << "<" << l.a << ", " << l.b << ">"; }
};
1
2
3
4
5
6
7
8
9

解释：

• Line结构体表示由两个点a和b确定的一条直线或线段。

线的几何操作

以下函数实现了直线与线段的各种几何操作，包括点到直线的距离、投影，判断相对位置关系，以及计算直线的交点等。

// 计算点p在直线上的投影
Point projection(const Line &l, const Point &p) {
    Point v = l.b - l.a;
    return l.a + v * (dot(p - l.a, v) / dot(v, v));
}

// 计算点p到直线的距离
Real distanceToPoint(const Line &l, const Point &p) {
    return fabs(cross(l.b - l.a, p - l.a)) / length(l.b - l.a);
}

// 判断点p是否在线段ab上
bool onSegment(const Line &l, const Point &p) {
    return fabs(cross(l.b - l.a, p - l.a)) < EPS && dot(p - l.a, p - l.b) < EPS;
}

// 判断点p相对于直线的左右位置 (-1: 左侧, 0: 共线, 1: 右侧)
int pointPosition(const Line &l, const Point &p) {
    Real crossProduct = cross(l.b - l.a, p - l.a);
    if (fabs(crossProduct) < EPS) return 0;  // 共线
    return crossProduct > 0 ? -1 : 1;  // 左侧为负值，右侧为正值
}

// 判断两直线是否平行
bool isParallel(const Line &l1, const Line &l2) {
    return areParallel(l1.b - l1.a, l2.b - l2.a);
}

// 判断两直线是否垂直
bool isPerpendicular(const Line &l1, const Line &l2) {
    return arePerpendicular(l1.b - l1.a, l2.b - l2.a);
}

// 计算两直线的夹角
Real angleBetween(const Line &l1, const Line &l2) {
    return angleBetween(l1.b - l1.a, l2.b - l2.a);
}

// 判断两线段是否相交
bool intersects(const Line &l1, const Line &l2) {
    int d1 = pointPosition(l1, l2.a);
    int d2 = pointPosition(l1, l2.b);
    int d3 = pointPosition(l2, l1.a);
    int d4 = pointPosition(l2, l1.b);
    return d1 * d2 <= 0 && d3 * d4 <= 0;
}

// 计算两直线的交点
Point intersection(const Line &l1, const Line &l2) {
    Real cross1 = cross(l1.b - l1.a, l2.b - l2.a);
    Real cross2 = cross(l2.a - l1.a, l2.b - l2.a);
    if (fabs(cross1) < EPS) throw runtime_error("Lines are parallel");  // 平行无交点
    return l1.a + (l1.b - l1.a) * (cross2 / cross1);
}
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

解释：

• projection函数计算点在直线上的投影，即点在直线上最近的点。
• distanceToPoint计算点到直线的距离。
• onSegment判断点是否在线段上，通过检查点是否共线且在线段的两端点之间。
• pointPosition判断点相对于直线的左右位置，返回值表示左侧、共线或右侧。
• isParallel和isPerpendicular判断两直线是否平行或垂直。
• angleBetween计算两直线之间的夹角。
• intersects判断两线段是否相交。
• intersection计算两直线的交点，如果两直线平行则抛出异常。

完整代码

#include <bits/stdc++.h>
using namespace std;

using Real = double;
const Real EPS = 1e-8; 

struct Point {
    Real x, y;

    Point(Real x = 0, Real y = 0) : x(x), y(y) {}

    Point operator+(const Point &p) const { return {x + p.x, y + p.y}; }
    Point operator-(const Point &p) const { return {x - p.x, y - p.y}; }
    Point operator*(Real c) const { return {x * c, y * c}; }
    Point operator/(Real c) const { return {x / c, y / c}; }

    bool operator==(const Point &p) const { return fabs(x - p.x) < EPS && fabs(y - p.y) < EPS; }
    bool operator<(const Point &p) const { return x < p.x - EPS || (fabs(x - p.x) < EPS && y < p.y - EPS); }

    friend ostream& operator<<(ostream &os, const Point &p) { return os << "(" << p.x << ", " << p.y << ")"; }
};

Real dot(const Point &a, const Point &b) { return a.x * b.x + a.y * b.y; }
Real cross(const Point &a, const Point &b) { return a.x * b.y - a.y * b.x; }
Real length(const Point &v) { return sqrt(dot(v, v)); }
Real distance(const Point &a, const Point &b) { return length(a - b); }

bool areParallel(const Point &a, const Point &b) { return fabs(cross(a, b)) < EPS; }
bool arePerpendicular(const Point &a, const Point &b) { return fabs(dot(a, b)) < EPS; }
Real angleBetween(const Point &a, const Point &b) { return acos(dot(a, b) / (length(a) * length(b))); }
Point rotate(const Point &v, Real theta) { return {v.x * cos(theta) - v.y * sin(theta), v.x * sin(theta) + v.y * cos(theta)}; }
Point rotateAround(const Point &p, const Point &origin, Real theta) { return origin + rotate(p - origin, theta); }

struct Line {
    Point a, b;

    Line(Point a = {}, Point b = {}) : a(a), b(b) {}
    
    friend ostream& operator<<(ostream &os, const Line &l) { return os << "<" << l.a << ", " << l.b << ">"; }
};

Point projection(const Line &l, const Point &p) {
    Point v = l.b - l.a;
    return l.a + v * (dot(p - l.a, v) / dot(v, v));
}

Real distanceToPoint(const Line &l, const Point &p) {
    return fabs(cross(l.b - l.a, p - l.a)) / length(l.b - l.a);
}

bool onSegment(const Line &l, const Point &p) {
    return fabs(cross(l.b - l.a, p - l.a)) < EPS && dot(p - l.a, p - l.b) < EPS;
}

int pointPosition(const Line &l, const Point &p) {
    Real crossProduct = cross(l.b - l.a, p - l.a);
    if (fabs(crossProduct) < EPS) return 0;
    return crossProduct > 0 ? -1 : 1;
}

bool isParallel(const Line &l1, const Line &l2) {
    return areParallel(l1.b - l1.a, l2.b - l2.a);
}

bool isPerpendicular(const Line &l1, const Line &l2) {
    return arePerpendicular(l1.b - l1.a, l2.b - l2.a);
}

Real angleBetween(const Line &l1, const Line &l2) {
    return angleBetween(l1.b - l1.a, l2.b - l2.a);
}

bool intersects(const Line &l1, const Line &l2) {
    int d1 = pointPosition(l1, l2.a);
    int d2 = pointPosition(l1, l2.b);
    int d3 = pointPosition(l2, l1.a);
    int d4 = pointPosition(l2, l1.b);
    return d1 * d2 <= 0 && d3 * d4 <= 0;
}

Point intersection(const Line &l1, const Line &l2) {
    Real cross1 = cross(l1.b - l1.a, l2.b - l2.a);
    Real cross2 = cross(l2.a - l1.a, l2.b - l2.a);
    if (fabs(cross1) < EPS) throw runtime_error("Lines are parallel");
    return l1.a + (l1.b - l1.a) * (cross2 / cross1);
}
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