计算几何(Computational Geometry),是一系列使用计算机解决几何问题的算法。与解析几何相比,计算几何更适合计算机运算,精度较高,运算速度较快,并且易于编写。
本文只包含二维计算几何在 OI 中的部分应用。
浮点误差
程序设计中,考虑到浮点数 double
有精度误差,在比较时,通常允许一定的误差,即对于两个数 、,如果 ,则认为 。一般根据题目要求,(代码中命名为 EPS
)取一个较小值,如 。
const double EPS = 1e-8;
// 带误差比较,返回 x 是否等于 y
inline bool dcmp(double x, double y = 0)
{
return fabs(x - y) <= EPS;
}
向量
向量(vector)是一个有大小和方向的量,在几何中,它被表示为带箭头的线段。向量可以用起点和终点的坐标来表示 —— 从点 到点 的向量表示为 。
向量的书写,两个大写字母上加一个箭头(表示方向),或者加粗的小写字母 。
向量没有位置,即向量可以在平面内任意平移而保持其本身的性质不变 —— 所以我们可以将向量平移至其起点与原点 重合,用终点坐标表示这个向量,下文称此方法为坐标表示法,除特殊说明外,下文表示向量均使用此方法。
两个向量的夹角 定义为从 旋转到 所经过的角度,逆时针为正,顺时针为负。
向量范数(norm),这里简单定义向量 的范数为 ,即长度的平方 。这样定义的好处是,可以方便地比较两个向量的长度,而不需要开平方。
与向量相对的概念是标量(scalar),标量只有大小没有方向。
向量的基本运算
两个向量相加,遵循平行四边形定则(或三角形定则),。
两个向量相减与相加类似,除此之外,点可以加减一个向量,表示将点按照向量所指向的方向进行平移。
向量的数乘,将一个向量乘或除以 ,表示将向量所对应的有向线段长度变为原来的 倍或 分之一。
可以发现,向量和点的运算是相似的,在程序中,可以使用同一种类型来表示向量和点:
/*
* 向量(Vector)或点
*
* 使用原点到一个点 (x, y) 的有向线段表示向量
* 从点 A 到点 B 的向量表示为 A - B
*/
typedef struct Vec
{
double x, y;
Vec(double x = 0, double y = 0) : x(x), y(y) {}
// 相加
Vec operator+(const Vec &v) const
{
return Vec(x + v.x, y + v.y);
}
// 相减
Vec operator-(const Vec &v) const
{
return Vec(x - v.x, y - v.y);
}
// 数乘(伸长、缩短)
Vec operator*(double d) const
{
return Vec(x * d, y * d);
}
Vec operator/(const double d) const
{
return Vec(x / d, y / d);
}
// 范数,用来比较长度,等于长度的平方
double norm() const
{
return x * x + y * y;
}
} Pt;
为了方便阅读,我们定义类型点 Pt
为向量 Vec
的别名,在代码中区别「点」和「向量」。同时,应记住点与向量的运算中的类型:
- 向量 向量 向量
- 点 向量 点
- 向量 数 向量
向量乘法
向量的乘法分为两种 —— 点乘(Dot)和叉乘(Cross)。
点乘
向量的点乘表示为 ,设 ,则
点乘的结果是一个标量,它表示 在 所在直线上的投影长度与 长度的乘积,即
从上式中可以看出,因为 ,所以 ,即点乘满足交换律。
叉乘
向量的叉乘表示为 ,设 ,则
叉乘的结果是一个标量,它表示以 与 为邻边的平行四边形的(有向)面积,即
这个平行四边形可以以 为底,底边上的高即为 。
结果的符号与 相同,即当 在 逆时针方向时,结果为正。
如果需要求以两向量为两边的三角形面积,可以求出平行四边形面积后除以二。
代码实现
// 点乘
double dot(const Vec &a, const Vec &b)
{
return a.x * b.x + a.y * b.y;
}
// 叉乘
double cross(const Vec &a, const Vec &b)
{
return a.x * b.y - a.y * b.x;
}
直线与线段
用直线上的两个点表示直线,用线段的两端点表示线段。
直线上的两个点组成的向量,被称为这条直线的一个方向向量。
点在直线上的判定
如果点 在直线上,则以 、 为邻边的平行四边形不存在(面积为零),即 。
bool include(const Pt &p) const
{
return dcmp(cross(a - p, b - p));
}
判断两直线关系
- 如果一条直线上的两个点都在另一条实现上,则两条直线重合;
- 如果不重合的两条直线的方向向量叉乘得零(无法组成平行四边形),则两条直线平行。
// 两直线关系(交点个数)
// 0 表示平行(无交点)
// 1 表示相交(一个交点)
// -1 表示重合(无数个交点)
static int relation(const Line &a, const Line &b)
{
if (a.include(b.a) && a.include(b.b)) return -1;
else if (dcmp(cross(a.b - a.a, b.b - b.a))) return 0;
else return 1;
}
求两直线交点
在判断相交后,可以求两直线 、 的交点。设其交点为 ,分别过点 、 作 的垂线,连接 、、、。作向量 ,只需要求出 ,即可得到 。
因为
所以
所以
即
用叉乘计算面积,求出向量 即可,注意计算面积时的方向。
// 求两直线交点(需要保证两直线有交点)
static Pt intersect(const Line &a, const Line &b)
{
double s1 = cross(b.a - a.a, b.b - a.a), s2 = cross(b.b - a.b, b.a - a.b);
return a.a + (a.b - a.a) * s1 / (s1 + s2);
}
点在线段上的判定
如果点 在线段 所在的直线上,则以 、 为邻边的平行四边形不存在(面积为零),即 。
此时还需判断 是否在 与 之间。
- 如果 ,则表示 ();
- 如果 ,则表示 ();
- 如果 ,则表示 或 ,即 与 或 重合。
后两种情况表明点 在 上。
// 线段包含点(点在线段上)
bool include(const Pt &p) const
{
// PA × PB = 0:PA 与 PB 共线,即点在线段所在的直线上
// PA · PB < 0:PA 与 PB 方向不同(A 和 B 分别在 P 的两边),如果 PA · PB = 0 则 P = A 或 P = B
return dcmp(cross(a - p, b - p)) && dot(a - p, b - p) <= 0;
}
线段相交 / 求交点
求对应所在直线的交点,判断是否在线段上即可。
多边形
多边形(polygon)通常用顺时针或逆时针排列的顶点来表示。
struct Poly
{
std::vector<Pt> pts;
};
点在多边形内的判定
判断点是否在多边形内部,常使用射线法,即从该点引一条平行于 轴的射线,如果它与多边形有奇数个交点,则该点在多边形内,否则在多边形外。
考虑相邻两个顶点 与 ,判断由 引出的水平向右的射线有没有与线段 相交。作向量 ,分别过 、 作 轴的垂线,垂足分别为 、,设
如果 在 的顺时针方向上,即 ,则需要 在 上方, 在 下方,即
注意到 可以为 ,即射线穿过 被认为与当前边相交,因为可能会有这种情况
此时不应认为 在多边形内,如果因为 被认为被射线穿过了一次,则会误判。解决方案是对于每一条边,只将射线穿过上方点(也可以取下方点)的情况下认为射线与边相交,这样 会被统计两次(或零次),不会影响判断。
如果 在 的逆时针方向上,即 ,则需要 在 下方, 在 上方,即 。
代码实现
struct Poly
{
std::vector<Pt> pts;
bool include(const Pt &p) const
{
int cnt = 0;
// 判断与每条边有没有交点
for (size_t i = 0; i < pts.size(); i++)
{
// 枚举相邻的每两个点
const Pt &a = pts[i], &b = pts[(i + 1) % pts.size()];
// 如果点 P 在边 AB 上
if (Seg(a, b).include(p)) return true;
// 详见图
double d1 = a.y - p.y, d2 = b.y - p.y, tmp = cross(a - p, b - p);
if ((tmp >= 0 && d1 >= 0 && d2 < 0) || (tmp <= 0 && d1 < 0 && d2 >= 0)) cnt++;
}
// 奇数的交点
return cnt % 2 == 1;
}
};
多边形面积
按顺序枚举每两个相邻的顶点,将每两个相邻的顶点与原点围成的三角形面积累加起来。因为所有不包含在多边形内的区域都被一正一负计算了偶数次,所以最终结果只包含多边形内的部分。
如图,浅蓝色部分被计算了两次,点 左下深色部分,在多边形内的部分被计算了三次,不在多边形内的部分被计算了四次,其它在多边形内部的部分被计算了一次。
如果按照逆时针顺序枚举,则结果为正,否则结果为负。
// 多边形面积(有向面积)
double area() const
{
double res = 0;
for (size_t i = 0; i < pts.size(); i++)
{
// 枚举每两个点
const Pt &a = pts[i], &b = pts[(i + 1) % pts.size()];
res += cross(a, b);
}
return res / 2;
}
凸包
在平面内,对于给定的一些点,包含这些点的最小凸多边形被称为凸包(Convex Hull)。
求凸包的一种方法是 Graham's Scan,算法的流程为:
- 选出横坐标最小(横坐标相同时纵座标最小)的点,作为极点,该点一定在凸包上;
- 对其它点按照与极点的极角排序,极角相同的按照与极点的距离排序;
- 用栈维护当前在凸包上的点,按顺序考虑每一个点是否可以替换栈顶的点。
这里极角是指,过极点 作水平向右的向量 ,则点 的极角为 。比较两个点 、 的极角,如果 经过逆时针旋转到 ,即 ,则 极角序较小。
枚举每一个点,考虑新加入一个点的影响:
当前已被加入到凸包中的点为 (极点)、、、,考虑新加入点 的影响 —— 相对于 向顺时针方向旋转了(即 ),此时这些边一定无法组成一个凸多边形。所以要将 点从凸包中删除,加入 点。这种情况我们称 的连线不合法。
加入一个点时,应当考虑与栈顶点的连线是否合法,如果不合法,则将栈顶点删除,继续考虑栈顶点,直到栈中只剩一个点或连线合法。
// 求凸包用的点
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;
// 求凸包(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]);
}
}
};
完整代码
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <vector>
const int MAXN = 100000;
const double EPS = 1e-8;
// 带误差比较
inline bool dcmp(double x, double y = 0)
{
return fabs(x - y) <= EPS;
}
/*
* 向量(Vector)或点
*
* 使用原点到一个点 (x, y) 的有向线段表示向量
* 从点 A 到点 B 的向量表示为 A - B
*/
typedef struct Vec
{
double x, y;
Vec(double x = 0, double y = 0) : x(x), y(y) {}
// 相加
Vec operator+(const Vec &v) const
{
return Vec(x + v.x, y + v.y);
}
// 相减
Vec operator-(const Vec &v) const
{
return Vec(x - v.x, y - v.y);
}
// 数乘(伸长、缩短)
Vec operator*(double d) const
{
return Vec(x * d, y * d);
}
Vec operator/(const double d) const
{
return Vec(x / d, y / d);
}
// 范数,用来比较长度,等于长度的平方
double norm() const
{
return x * x + y * y;
}
} Pt;
// 点乘
double dot(const Vec &a, const Vec &b)
{
return a.x * b.x + a.y * b.y;
}
// 叉乘
double cross(const Vec &a, const Vec &b)
{
return a.x * b.y - a.y * b.x;
}
// 线段(Segment),用两个点表示
struct Seg
{
Pt a, b;
Seg(const Pt &a, const Pt &b) : a(a), b(b) {}
// 线段包含点(点在线段上)
bool include(const Pt &p)
{
// PA × PB = 0:PA 与 PB 共线,即点在线段所在的直线上
// PA · PB = 0:PA 与 PB 方向不同(A 和 B 分别在 P 的两边),如果 PA · PB = 0 则 P = A 或 P = B
return dcmp(cross(a - p, b - p)) && dot(a - p, b - p) <= 0;
}
};
// 直线,用两个点表示
struct Line
{
Pt a, b;
Line() {} // 提供一个不需要参数的构造函数
Line(const Pt &a, const Pt &b) : a(a), b(b) {}
bool include(const Pt &p) const
{
return dcmp(cross(a - p, b - p));
}
// 两直线关系(交点个数)
// 0 表示平行(无交点)
// 1 表示相交(一个交点)
// -1 表示重合(无数个交点)
static int relation(const Line &a, const Line &b)
{
if (a.include(b.a) && a.include(b.b)) return -1;
else if (dcmp(cross(a.b - a.a, b.b - b.a))) return 0;
else return 1;
}
// 求两直线交点(需要保证两直线有交点)
static Pt intersect(const Line &a, const Line &b)
{
double s1 = cross(b.a - a.a, b.b - a.a), s2 = cross(b.b - a.b, b.a - a.b);
return a.a + (a.b - a.a) * s1 / (s1 + s2);
}
};
// 求凸包用的点
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;
bool include(const Pt &p) const
{
int cnt = 0;
// 判断与每条边有没有交点
for (size_t i = 0; i < pts.size(); i++)
{
// 枚举相邻的每两个点
const Pt &a = pts[i], &b = pts[(i + 1) % pts.size()];
// 如果点 P 在边 AB 上
if (Seg(a, b).include(p)) return true;
// 详见图
double d1 = a.y - p.y, d2 = b.y - p.y, tmp = cross(a - p, b - p);
if ((tmp >= 0 && d1 >= 0 && d2 < 0) || (tmp <= 0 && d1 < 0 && d2 >= 0)) cnt++;
}
// 奇数的交点
return cnt % 2 == 1;
}
// 多边形面积(有向面积)
double area() const
{
double res = 0;
for (size_t i = 0; i < pts.size(); i++)
{
// 枚举每两个点
const Pt &a = pts[i], &b = pts[(i + 1) % pts.size()];
res += cross(a, b);
}
return res / 2;
}
// 求凸包(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(a[i] - pts[pts.size() - 2], pts.back() - pts[pts.size() - 2]) >= 0) pts.pop_back();
pts.push_back(a[i]);
}
}
};