ACM
February 21, 2021

计算几何

计算几何

基础

高精度圆周率

1
const double pi = acos(-1.0);

偏差值

1
const double eps = 1e-8;

sgn

1
2
3
4
5
int sgn(double x)	// 判断x是否等于0
{
	if (fabs(x) < eps) return 0;
	else return x < 0 ? -1 : 1;
}

dcmp

1
2
3
4
5
int dcmp(double x, double y)	// 比较两个浮点数，0为相等，-1为小于，1为大于
{
	if (fabs(x - y) < eps) return 0;
	else return x < y ? -1 : 1;
}

点和向量

点

1
2
3
4
5
6
struct Point
{
	double x, y;
	Point() {}
	Point(double x, double y) : x(x), y(y) {}
};

两点间距离

1. 

1
2
3
4
double Dist(Point A, Point B)
{
	return hypot(A.x - B.x, A.y - B.y);
}

2. 

1
2
3
4
double Dist(Point A, Point B)
{
	return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

向量

1
typedef Point Vector;

向量的运算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
struct Point
{
	double x, y;
	Point() {}
	Point(double x, double y) : x(x), y(y) {}

	Point operator+(Point B) { return Point(x + B.x, y + B.y); }	// 加

	Point operator-(Point B) { return Point(x - B.x, y - B.y); }	// 减

	Point operator*(double k) { return Point(x * k, y * k); }	// 乘

	Point operator/(double k) { return Point(x / k, y / k); }	// 除

	bool operator==(Point B) { return sgn(x - B.x) == 0 && sgn(y - B.y) == 0; }	// 等于

};

点积

定义

A · B = |A||B|cosθ

几何定义

θ: 表示向量 b 与向量 a 的夹角

向量 b 在向量 a 上的投影长度乘以向量 a 的模长

公式

A · B = A.x * B.x + A.y * B.y

代码

1
2
3
4
double Dot(Vector A, Vector B)
{
	return A.x * B.x + A.y * B.y;
}

应用

  1. 判断向量 A 与向量 B 的夹角是钝角还是锐角

    • 若 Dot(A, B) > 0: 向量 A 与向量 B 的夹角为锐角
    • 若 Dot(A, B) < 0: 向量 A 与向量 B 的夹角为顿角
    • 若 Dot(A, B) = 0: 向量 A 与向量 B 的夹角为直角

  2. 求向量 A 的模长

    1
    2
    3
    4
    double Len(Vector A)
    {
    	return sqrt(Dot(A, A));
    }

  3. 求向量 A 与向量 B 的夹角大小

    1
    2
    3
    4
    double Angle(Vector A, Vector B)
    {
    	return acos(Dot(A, B) / Len(A) / Len(B));
    }

叉积

定义

A × B = |A||B|sinθ

几何意义

θ: 向量 A 旋转到向量 B 所经过的夹角

|A × B| 在数值上等于由向量和向量构成的平行四边形的面积

代码

1
2
3
4
double Cross(Vector A, Vector B)
{
	return A.x * B.y - A.y * B.x;
}

应用

  1. 判断向量与向量的方向关系

    • 若 A × B > 0，B 在 A 的逆时针方向
    • 若 A × B < 0，B 在 A 的顺时针方向
    • 若 A × B = 0，B 与 A 共线，可能是同方向，也可能是反方向

  2. 计算两向量构成的平行四边形的有向面积

    3 个点 A, B, C 以 A 为公共点，得到两个向量 B - A, C - A 它们构成的平行四边形的面积如下

    1
    2
    3
    4
    double Area2(Point A, Point B, Point C)
    {
    	return Cross(B - A, C - A);
    }

    同理 Area2(A, B, C) / 2 就是计算以 A, B, C 三点构成三角形的面积

  3. 向量旋转

    向量 A 逆时针旋转的角度为 rad

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

    有时需要求单位法向量，即逆时针旋转 90°，然后取单位值

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

  4. 用叉积检查两个向量是否平行或重合

    1
    2
    3
    4
    bool Parallel(Vector A, Vector B)
    {
    	return sgn(Cross(A, B)) == 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
32
33
34
35
struct Line	// 直线
{
	Point p1, p2;	// 线上的两个点
	Line() {}
	// 直接用两个点来构造直线
	Line(Point p1, Point p2) : p1(p1), p2(p2) {}
	// 根据一个点和倾斜角angle确定直线，0≤angle≤pi
	Line(Point p, double angle)
	{
		p1 = p;
		if (sgn(angle - pi / 2) == 0)
			p2 = (p1 + Point(0, 1));
		else
			p2 = (p1 + Point(1, tan(angle)));
	}
	// ax + by + c = 0
	Line(double a, double b, double c)
	{
		if (sgn(a) == 0)
		{
			p1 = Point(0, -c / b);
			p2 = Point(1, -c / b);
		}
		else if (sgn(b) == 0)
		{
			p1 = Point(-c / a, 0);
			p2 = Point(-c / a, 1);
		}
		else
		{
			p1 = Point(0, -c / b);
			p2 = Point(1, (-c - a) / b);
		}
	}
};

线段的表示

1
typedef Line Segment;

点和直线的位置关系

用直线 v 上的两点 p1 和 p2 与点 p 构成两个向量，用叉积的正负判断方向，得到相对的位置关系点

1
2
3
4
5
6
7
8
9
int Point_line_relation(Point p, Line v)
{
	int c = sgn(Cross(p - v.p1, v.p2 - v.p1));
	if (c < 0) // 1：p在v的左边
		return 1;
	if (c > 0) // 2：p在v的右边
		return 2;
	return 0; // 0：p在v上
}

点和线段的位置关系

判断点 p 是否在线段 v 上，先用叉积判断是否共线，然后用点积看 p 和 v 的两个端点产生的角是否为钝角

1
2
3
4
bool Point_on_seg(Point p, Line v)	// 0为点不在线段v上；1为点在线段v上
{
	return sgn(Cross(p - v.p1, v.p2 - v.p1)) == 0 && sgn(Dot(p - v.p1, p - v.p2)) <= 0;
}

点到直线的距离

已知点 p 和直线 v(p1, p2)，求 p 到 v 的距离。首先用叉积求 p, p1, p2 构成的平行四边形的面积，然后用面积除以平行四边形的底边长，也就是线段 (p1, p2) 的长度，就得到了平行四边形的高，即点 p 到直线 v 的距离

1
2
3
4
double Dis_point_line(Point p, Line v)
{
	return fabs(Cross(p - v.p1, v.p2 - v.p1)) / Dist(v.p1, v.p2);
}

点在直线上的投影

1
2
3
4
5
Point Point_line_proj(Point p, Line v)	// 点p在直线v上的投影
{
	double k = Dot(v.p2 - v.p1, p - v.p1) / Len2(v.p2 - v.p1);
	return v.p1 + (v.p2 - v.p1) * k;
}

点关于直线的对称点

1
2
3
4
5
Point Point_line_symmetry(Point p, Line v)	// 点p关于直线v的对称点
{
	Point q = Point_line_proj(p, v);
	return Point(2 * q.x - p.x, 2 * q.y - p.y);
}

点到线段的距离

对于点 p 到线段 v(p1, p2) 的距离，在以下 3 个距离中取最小值：从 p 出发对线段 v 做垂线，如果交点在 v 上，这个距离就是最小值；p 到 p1 的距离，p 到 p2 的距离

1
2
3
4
5
6
double Dis_point_seg(Point p, Segment v)
{
	if (sgn(Dot(p - v.p1, v.p2 - v.p1)) < 0 || sgn(Dot(p - v.p2, v.p1 - v.p2)) < 0)
		return min(Dist(p, v.p1), Dist(p, v.p2));
	return Dis_point_line(p, v);
}

两条直线的位置关系

1
2
3
4
5
6
7
8
9
10
11
int Line_relation(Line v1, Line v2)
{
	if (sgn(Cross(v1.p2 - v1.p1, v2.p2 - v2.p1)) == 0)
	{
		if (Point_line_relation(v1.p1, v2) == 0)
			return 1; // 1：重合
		else
			return 0; // 0：平行
	}
	return 2; // 2：相交
}

求两条直线的交点

1
2
3
4
5
6
Point Cross_point(Point a, Point b, Point c, Point d)	// Line: ab, Line: cd
{
	double s1 = Cross(b - a, c - a);
	double s2 = Cross(b - a, d - a);
	return Point(c.x * s2 - d.x * s1, c.y * s2 - d.y * s1) / (s2 - s1);
}

判断两个线段是否规范相交

这里利用叉积有正负的特点。如果一条线段的两端在另一条线段的两侧，那么两个端点与另一线段产生的两个叉积的正负相反，也就是说两个叉积相乘为负。如果两条线段互相满足这一点，那么就是规范相交的。

规范相交：交点在线段内部

非规范相交：交点在某条线段的端点

1
2
3
4
5
6
7
8
9
10
bool Cross_segment(Point a, Point b, Point c, Point d)
{
	// 规范相交
	double c1 = Cross(b - a, c - a), c2 = Cross(b - a, d - a);
	double d1 = Cross(d - c, a - c), d2 = Cross(d - c, b - c);
	return sgn(c1) * sgn(c2) < 0 && sgn(d1) * sgn(d2) < 0; // 1: 相交，0: 不相交

	// 非规范相交
	return max(a.x, b.x) >= min(c.x, d.x) && max(c.x, d.x) >= min(a.x, b.x) && max(a.y, b.y) >= min(c.y, d.y) && max(c.y, d.y) >= min(a.y, b.y) && sgn(Cross(b - a, c - a)) * sgn(Cross(b - a, d - a)) <= 0 && sgn(Cross(d - c, a - c)) * sgn(Cross(d - c, b - c)) <= 0;
}

多边形

代码

1
Point p[N];

应用

  1. 判断点是否在多边形内部

    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
    int Point_in_polygon(Point pt, Point* p, int n)	// 点pt，多边形Point* p
    {
    	for (int i = 0; i < n; ++i)
    		if (p[i] == pt)
    			return 3;	// 3: 点在多边形的顶点上
    	for (int i = 0; i < n; ++i)
    	{
    		Line v = Line(p[i], p[(i + 1) % n]);
    		if (Point_on_seg(pt, v))
    			return 2;	// 2: 点在多边形的边上
    	}
    	int num = 0;
    	for (int i = 0; i < n; ++i)
    	{
    		int j = (i + 1) % n;
    		int c = sgn(Cross(pt - p[j], p[i] - p[j]));
    		int u = sgn(p[i].y - pt.y);
    		int v = sgn(p[j].y - pt.y);
    		if (c > 0 && u < 0 && v >= 0)
    			++num;
    		if (c < 0 && u >= 0 && v < 0)
    			--num;
    	}
    	return num != 0;	// 1: 点在内部；0: 点在外部
    }

  2. 多边形的面积

    1
    2
    3
    4
    5
    6
    7
    double Polygon_area(Point* p, int n)
    {
    	double area = 0;
    	for (int i = 0; i < n; ++i)
    		area += Cross(p[i], p[(i + 1) % n]);
    	return area / 2;
    }

  3. 求多边形的重心

    1
    2
    3
    4
    5
    6
    7
    8
    9
    Point Polygon_center(Point* p, int n)
    {
    	Point ans(0, 0);
    	if (Polygon_area(p, n) == 0)
    		return ans;
    	for (int i = 0; i < n; ++i)
    		ans = ans + (p[i] + p[(i + 1) % n]) * Cross(p[i], p[(i + 1) % n]);
    	return ans / Polygon_area(p, n) / 6;
    }

凸包

凸多边形：所有内角大小都在 $[0, \ \pi]$ 范围内的 简单多边形。

凸包：在平面上能包含所有给定点的最小凸多边形叫做凸包。

定义：对于给定集合 X，所有包含 X 的凸集的交集 S 被称为 X 的 凸包。

IMG_9302

凸包用最小的周长围住了给定的所有点。如果一个凹多边形围住了所有的点，它的周长一定不是最小，如下图。根据三角不等式，凸多边形在周长上一定是最优的。

求法

常用的求法有 Graham 扫描法和 Andrew 算法，这里主要介绍 Andrew 算法。

Andrew

该算法的时间复杂度为 $O(n \log n)$，其中 n 为待求凸包点集的大小，同时复杂度的瓶颈也在于对所有点坐标的双关键字排序。

首先把所有点以横坐标为第一关键字，纵坐标为第二关键字排序。

显然排序后最小的元素和最大的元素一定在凸包上。而且因为是凸多边形，我们如果从一个点出发逆时针走，轨迹总是“左拐”的，一旦出现右拐，就说明这一段不在凸包上。因此我们可以用一个单调栈来维护上下凸壳。

因为从左向右看，上下凸壳所旋转的方向不同，为了让单调栈起作用，我们首先 升序枚举 求出下凸壳，然后 降序 求出上凸壳。

求凸壳时，一旦发现即将进栈的点（P）和栈顶的两个点（$S_{1}, \ S_{2}$，其中 $S_{1}$ 为栈顶）行进的方向向右旋转，即叉积小于 0：$\overrightarrow{S_{2}S_{1}} \ \times \ \overrightarrow{S_{1}P} \ < \ 0$，则弹出栈顶，回到上一步，继续检测，直到 $\overrightarrow{S_{2}S_{1}} \ \times \ \overrightarrow{S_{1}P} \ \geq \ 0$ 或者栈内仅剩一个元素为止。

通常情况下不需要保留位于凸包边上的点，因此上面一段中 $\overrightarrow{S_{2}S_{1}} \ \times \ \overrightarrow{S_{1}P} \ < \ 0$ 这个条件中的 $<$ 可以视情况改为 $\leq$，同时后面一个条件应改为 $>$。

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
// stk[] 是整型，存的是下标
// p[] 存储向量或点
tp = 0;                       // 初始化栈
std::sort(p + 1, p + 1 + n);  // 对点进行排序
stk[++tp] = 1;
//栈内添加第一个元素，且不更新 used，使得 1 在最后封闭凸包时也对单调栈更新
for (int i = 2; i <= n; ++i)
{
	while (tp >= 2 && (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)	// * 操作符被重载为叉积
		used[stk[tp--]] = 0;
	used[i] = 1;  // used 表示在凸壳上
	stk[++tp] = i;
}
int tmp = tp;  // tmp 表示下凸壳大小
for (int i = n - 1; i > 0; --i)
if (!used[i])
{
	//      ↓求上凸壳时不影响下凸壳
	while (tp > tmp && (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)
		used[stk[tp--]] = 0;
	used[i] = 1;
	stk[++tp] = i;
}
for (int i = 1; i <= tp; ++i)  // 复制到新数组中去
h[i] = p[stk[i]];
int ans = tp - 1;

根据上面的代码，最后凸包上有 ans 个元素（额外存储了 1 号点，因此 h 数组中有 ans + 1 个元素），并且按逆时针方向排序。周长就是 $\sum_{i=1}^{ans}|\overrightarrow{h_{i}h_{i+1}}|$

About this Post

This post is written by OwlllOvO, licensed under CC BY-NC 4.0.

#C++#计算几何