「计算几何」计算几何从入门到入土

这里我们来介绍一下计算几何的相关知识qwq

计算几何是【真】初中内容qwq

我们先来看看计算几何最重要的问题:精度问题

精度控制

计算几何经常牵扯到浮点数的运算,所以就会产生精度误差,因此我们需要设置一个 epseps(偏差值),一般取 10710^{-7}101010^{-10} 之间

比如我们比较两个 double 类型的变量是否相等,我们就可以这样写

1
2
3
inline bool xxx (double first, double second) {
return abs(first - second) < eps;
}

点是计算几何中常用的东西惹qwq

我们一般通过表示坐标的方法进行描述

1
2
3
4
5
struct Node {
double posx, posy;

Node (double a = 0, double b = 0) : posx (a) , posy (b) {}
}

两点距离

两点距离是小学内容吧qwq

A(x0,y0)A(x_0, y_0)B(x1,y1)B(x_1,y_1) 的距离是 (x0x1)2+(y0y1)2\sqrt{(x_0-x_1)^2+(y_0-y_1)^2}

代码实现:

1
2
3
inline double dis (const Node &a, const Node &b) {
return sqrt ( (a.posx - b.posx) * (a.posx - b.posx) + (a.posy - b.posy) * (a.posy - b.posy) );
}

拓展

对于 nn 维的两个点 A(a1,a2,,an)A(a_1, a_2,\cdots,a_n)B(b1,b2,,bn)B(b_1, b_2,\cdots,b_n) 的距离为

dis=i=1n(aibi)2\text{dis}=\sqrt{\sum^{n}_{i=1}(a_i-b_i)^2}
1
/* 代码略 */

中点

对于点 A(x0,y0)A(x_0, y_0) , B(x1,y1)B(x_1,y_1) 的中点坐标为 P(x0+x12,y0+y12)P\left({x_0+x_1\over 2},{y_0+y_1\over 2}\right)

1
2
3
4
5
6
Node mid (const Node &a, const Node &b) {
Node ans;
ans.posx = (a.posx + b.posx) / 2;
ans.posy = (a.posy + b.posy) / 2;
return ans;
}

向量

有大小有方向的量qwq

向量(vector)是一个有大小和方向的量,在几何中,它被表示为带箭头的线段。向量可以用起点终点的坐标来表示

img

如图就是几个向量

向量 ww 可以表示为 (12)\begin{pmatrix}-1\\-2\end{pmatrix}

通常用一对数 (a,b)(a,b) 来表示 (ab)\begin{pmatrix}a\\b\end{pmatrix}

1
2
3
4
struct vector {
double x,y;
vector (double _x = 0, double _y = 0) : x(_x) , y(_y) {}
}

不要和 stl 中的容器弄混惹qwq

向量的模

向量的模,即向量的长度

对于 a=(xy)a=\begin{pmatrix}x\\y\end{pmatrix} 它的长度为 a=x2+y2|a|=\sqrt{x^2+y^2}

1
2
3
double length (vector a) {
return sqrt (a.x * a.x + a.y * a.y);
}

向量的运算

我们要注意的是

  1. pointpoint=vector\text{point} - \text{point}=\text{vector}
  2. point±vector=point\text{point}\pm\text{vector}=\text{point}
1
2
3
4
vector operator + (vector a, vector b) { return vector (a.x + b.x, a.y + b.y); }
vector operator - (vector a, vector b) { return vector (a.x - b.x, a.y - b.y); }
vector operator * (vector a, double p) { return vector (a.x * p, a.y * p); }
vector operator / (vector a. double p) { return vector (a.x / p, a.y / p); }

向量的点积

aba \cdot b 的几何意义为 aabb 上的投影长度乘以 bb 的模长

ab=abcosθa\cdot b=|a||b|\cos \theta , 其中 θ\thetaaa , bb 之间的夹角

两个向量 a,ba,b 的夹角 θ\theta 定义为从 aa 旋转到 bb 所经过的角度,逆时针为正,顺时针为负。

a=(x0y0),b=(x1y1)ab=x0×x1+y0×y1a=\begin{pmatrix}x_0\\y_0\end{pmatrix},b=\begin{pmatrix}x_1\\y_1\end{pmatrix}\Rightarrow a\cdot b=x_0\times x_1+y_0\times y_1

点积的应用

  1. aba \bot b ,有 ab=0a\cdot b=0,反之亦然
  2. 求两个向量的夹角,若 ab>0a\cdot b >0,夹角为锐角,ab<0a\cdot b<0 为钝角

二维叉积

两个向量的叉积是一个标量,a×ba\times b 的几何意义为他们所形成的平行四边形的 有向面积

a=(x0y0),b=(x1y1)a×b=x0y1x1y0a=\begin{pmatrix}x_0\\y_0\end{pmatrix}, b=\begin{pmatrix}x_1\\y_1\end{pmatrix} \Rightarrow a\times b=x_0y_1-x_1y_0
1
2
3
double cross (vector a, vector b) {
return a.x * b.y - a.y * b.x;
}

叉积的作用

直观理解,假如 bbaa 的左边,则有向面积为正,假如在右边则为负。假如 b,ab,a 共线,则叉积为 00

所以叉积可以判断平行qwq

向量的旋转

a=(xy)a=\begin{pmatrix}x\\y\end{pmatrix} 可以看成是 x×(10)+y×(01)x\times \begin{pmatrix}1\\0\end{pmatrix}+y\times \begin{pmatrix}0\\1\end{pmatrix}

将其旋转 θ\theta °,

一些常用的方法

点到直线的距离

利用叉积求面积,然后除以平行四边形的底边长,得到平行四边形的高即点到直线的距离

img

【真】小学内容qwq

1
2
3
4
double point_vector_dist (point p, point a, point b) {
vector v = p - a; vector u = b - a;
return fabs (cross (v, u)) / length (u);
}

点到线段的距离

img

比点到直线的距离稍微复杂

如果平行四边形的高在区域之外的话就不合理,这时候需要计算点到距离较近的端点的距离

1
2
3
4
5
6
7
double point_dists (point p, point a, point b) {
if (a == b) return len(p - a);
vector v1 = b - a, v2 = p - a, v3 = p - b;
if (dot(v1,v2) < 0) return len(v2);
else if (dot(v1,v3) > 0) return len(v3);
return fabs(cross(v1,v2)) / len(v1);
}

判断两个线段是否相交

跨立实验

判断一条线段的两端是否在另一条线段的两侧(两个端点与另一线段的叉积乘积为负)。需要正反判断两侧

dcmp(double a) 是我们自己定义的一个函数,函数原型:

1
2
3
4
int dcmp (double x) {
if (fabs(x) < eps) return 0;
else return x < 0 ? -1 : 1;
}

Code:

1
2
3
4
5
bool 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 dcmp(c1) * dcmp(c2) < 0 && dcmp(d1) * dcmp(d2) < 0;
}

判断点是否在多边形内部

经典问题

我们来看一个简单的栗子:

img

这个画图软件已经标记出来了多边形的内部,我们来看看在程序中如何辨别:

我们从这个点随便引一条射线,如果他与多边形有奇数个交点,就说明它在多边形内部

如果他与多边形有偶数个交点,就说明他在多边形外部

注意:这个多边形可能不是凸多边形

img

我们来看一个复杂的栗子:

img

1
2
3
4
5
6
7
8
9
10
11
12
13
int point_in (point p, point *a, int n) {  
int wn = 0, k, d1, d2;
for (register int i = 1; i <= n; i++) {
if (dcmp (dists(p, a[i], a[(i + 1 - 1) % n + 1])) == 0) return -1;//判断点是否在多边形的边界上
k = dcmp(cross(a[ (i + 1 - 1) % n + 1 ] - a[i], p - a[i]));
d1 = dcmp(a[i].y - p.y);
d2 = dcmp(a[ (i + 1 - 1) % n + 1 ].y - p.y);
if (k > 0 && d1 <= 0 && d2 > 0) wn++;
if (k < 0 && d2 <= 0 && d1 > 0) wn--;
}
if (wn) return 1;
else return 0;
}

求任意多边形面积

任取一个点进行三角剖分,用叉积求三角形的面积。因为叉积是有向面积,所以任意多边形都使用

img

虽然我们加了很多没用的面积,但是我们会减掉(容斥),这就求出来面积了

注意最后取绝对值

这里注意一个十分zz的错误,我们最后的面积是要除以 22 的(惨痛的经历.wa

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

求多边形的重心

同样方法将多边形三角剖分

算出每个三角形的重心,套用质点组的重心公式即可

质点组重心公式 三个点 A,B,CA,B,C

x=ma×xa+mb×xb+mc×xcma+mb+mcx={m_a\times x_a+m_b\times x_b+m_c\times x_c\over m_a+m_b+m_c}y=ma×ya+mb×yb+mc×ycma+mb+mcy={m_a\times y_a+m_b\times y_b+m_c\times y_c\over m_a+m_b+m_c}

mm 表示权,三角形的有向面积

这个还没有用到qwq(应该没有良心出题人会出这种题目吧qwq

凸包

凸包(Convex Hull\text{Convex Hull})是一个计算几何(图形学)中的概念。

在一个实数向量空间 VV 中,对于给定集合 XX ,所有包含 XX 的凸集的交集 SS 被称为 XX 的凸包。

我们使用通俗易懂的语言描述一下

给定二维平面上的点集,凸包就是将最外层的点连接起来构成的凸多边型,它能包含点集中所有的点

如图所示

img

比如我们看这个图,这个图的凸包大概是这样的:

img

凸包最常见的应用是求平面上距离最远的两个点,一般有两种算法来计算包含n个点的点集的凸包。这种两种算法都是按逆时针方向的顺序输出凸包

Graham 算法

时间复杂度为 O(nlog2n)O(n\log_2 n)

伪代码:

Graham(Q)\text {Graham}(Q)

  1. let p0 be the point in Q with the minimun y-coordinate\text {let } p_0 \text { be the point in } Q\text { with the minimun }y \text {-coordinate}
  2. let<p1,p2,,pn> be the remaining points in Qsort by polar angle  in counterclockwise order around p0\begin{matrix}\text {let}<p_1, p_2,\cdots,p_n>\text { be the remaining points in } Q\\\text {sort by polar angle }\\\text { in counterclockwise order around }p_0\end{matrix}
  3. push(p0,Stack)\text {push}(p_0,\text {Stack})
  4. push(p1,Stack)\text {push}(p_1, \text {Stack})
  5. push(p2,Stack)\text {push}(p_2,\text {Stack})
  6. for i in 3n then\text {for }i\text { in } 3 \to n \text { then}
  7.                   pop(Stack)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{pop}(\text{Stack})
  8.          push (pi,Stack)\ \ \ \ \ \ \ \ \ \text{push }(p_i,\text{Stack})
  9. return Stack\text{return Stack}

这显然是不显然的(???

我们来画几张较为显然的图来讲解一下qwq

img1

现在平面上有若干点组成的点集 GG

我们江他按照 AA 极角排序:

为什么按照A呢qwq

因为凸包中坐标极值的点一定在凸包内x,yx,y 的 最大/最小 值),比如图中的 A,C,FA,C,F 一定在凸包中

这里我们就选择最左边的点惹qwq

2

3

我们江 AA 和按极角排序最小的点 DD 连线

建立一个栈,栈内只有两个元素 A,DA,D

4

我们连第一小和第二小的点 DFDFFFADAD 的左侧,于是使点 FF 入栈

我们按照这个过程进行下去,我们发现连接 FGFG 之后 GGFDFD 的右侧,我们就将 FF 出栈,GG 入栈

5

按照以上步骤:

6

7

这样就完成求这个点集的凸包了

完成时栈的内容:[A,D,C,B][A,D,C,B]

这个正好就是凸包的点,我们将他倒序输出就是顺时针凸包节点

代码实现

from:Menci-Blog
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
// 求凸包用的点
int n;
Pt a[MAXN + 1];

// 凸包极角排序的比较函数
inline bool compare(const Pt &a, const Pt &b)
{
// 两个向量
Vec va = a - ::a[1], vb = b - ::a[1];
double t = cross(va, vb);
if (!dcmp(t)) return t > 0; // OA -> OB 是逆时针,则 A 极角序在先
else return va.norm() < vb.norm(); // norm 较小的长度较小
}

struct Poly
{
std::vector<Pt> pts;shixian

// 求凸包(Convex),结果储存在自身 pts 中
void convex()
{
// 找出最左下角的点
int id = 1;
for (int i = 1; i <= n; i++)
{
if (a[i].x < a[id].x || (a[i].x == a[id].x && a[i].y < a[id].y)) id = i;
}
if (id != 1) std::swap(a[1], a[id]);

// 排序
std::sort(a + 2, a + n + 1, &compare);

// 极角序扫描
pts.push_back(a[1]);
for (int i = 2; i <= n; i++)
{
// 比较,如果最后一个点需要被删掉则弹出(pop_back)
while (pts.size() >= 2 && cross(pts.back() - pts[pts.size() - 2], a[i] - pts.back()) <= 0) pts.pop_back();
pts.push_back(a[i]);
}
}
};

作业

【真】板子题

最大土地面积

题目

我们可以 O(n4)O(n^4) 暴力枚举

等等....我们好像学过另外一种求四边形面积的方法qwq

我们发现可以枚举一条边,再以该边为对角线求出左右两边的三角形最大值

我们又发现,好像最大的四边形的顶点都在凸包上!

为什么会变成这样呢……两件快乐的优化重合在一起。而这两份优化,又给我带来更多的分数。得到的,本该是像梦境一般AC的时间……但是,为什么,我会被他们打死呢?

首先,凸包是肯定要求的。

受到之前的启发,我们仍然考虑枚举一条对角线,设现在枚举到的顶点为 i,ji, j

设另两个点为 a,ba, b 。我们先来看 bb

首先,我们不难得出如果确定 i,ji,j ,那 bbijij 的距离一定是单峰的。那如果 jj 开始逆时针转动,峰显然也会逆时针转动(可以想象是整个凸包转过来了)

img

于是我们让 bb 逆时针跑就行了

aa 也同理

我们发现 bb 对每个点只扫了一遍,对于每个 ii 复杂度 O(n)O(n) ,同理 j,bj,b 都是,所以总复杂度 O(n2)O(n^2)

旋转卡壳

其实无论是求凸包的最大直径也好,还是凸包上距离最远的两个点也好,都可以使用万能方法---枚举暴力来求解,但是枚举暴力的时间复杂度为 O(n2)O(n^2),而接下来要介绍的卡壳法能够时间复杂度降到 O(n)O(n)

顾名思义,就是通过一对平行线来卡住凸包上的两个相对的点,然后通过旋转这对平行线能够得到凸包上距离最远的一对点

先来看一个直观图:

如果这样直接去实现不太好实现

我们可以转化一下,就是这一对平行线上在旋转的过程中我们可以让其中一条与凸包的一条边重合,此时另一条线的顶点是距这条边最远的点

我们可以依次枚举每条边求得距离这条最远的顶点,然后计算这个顶点到边两个端点的距离并取较大的一个,枚举完成时得到的最大距离即为直径长度

但是这样算的话因为枚举了所有边,每次又要枚举所有顶点,因此时间复杂度和暴力枚举一样。但是可以用下图观察到,当枚举的边逆时针旋转时,最远点也是跟着逆时针变化,这样我们可以不用每次枚举所有的顶点,直接从上次的最远点开始继续计算即可,此时复杂度为 O(n)O(n)。自己可以在纸上画下图模拟一下,就明白了

img

计算点到直线的距离时没必要直接计算,我们知道最远的顶点和这条组成的三角形面积一定是最大的,所以可以直接计算三角形面积,而三角形面积可以用叉积来计算

作业

板子

最大土地面积 - 旋转卡壳


End.

坚持原创技术分享,您的支持将鼓励我继续创作!