计算几何入门

计算几何

基础部分

高精度圆周率

const double pi = acos(-1.0);

偏差值

const double eps = 1e-8;

sgn

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

dcmp

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

点和向量

点

struct Point {
	double x, y;
	Point() {}
	Point(double x, double y) : x(x), y(y) {}
};

两点间距离

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

向量

typedef Point Vector;

向量的运算

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; }
	// 等于
};

点积

定义

几何意义
θ：表示向量 b与向量 a的夹角
向量 b在向量 a上的投影长度乘以向量 a的模长

公式

代码

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

应用

1. 判断向量 与向量 的夹角是钝角还是锐角
   若 Dot（A,B）> 0，说明向量 与向量 的夹角是锐角
   若 Dot（A,B）< 0，说明向量 与向量 的夹角是钝角
   若 Dot（A,B）== 0，说明向量 与向量 的夹角是直角
2. 求向量 A的模长

double Len(Vector A) {
	return sqrt(Dot(A, A));
}

求向量 A模长的平方

double Len2(Vector A) {
	return Dot(A, A);
}

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

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

叉积

定义

几何意义

θ：向量 a旋转到向量 b所经过的夹角
|a×b|在数值上等于由向量 a和向量 b构成的平行四边形的面积

代码

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

应用

1. 判断向量 A与向量 B的方向关系
   若 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 ，它们构成的平行四边形的面积
   如下

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

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

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

bool Parallel(Vector A, Vector B) {
	return sgn(Cross(A, B)) == 0;
}

点和线

代码

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);
		}
	}
};

线段的表示

typedef Line Segment;

点和直线的位置关系

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

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的两个端点产生的角是否为钝角

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的距离

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

点在直线上的投影

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;
}

点关于直线的对称点

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的距离

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);
}

两条直线的位置关系

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：相交
}

求两条直线的交点

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);
}

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

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

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

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

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;
}

多边形

代码

Point p[N];

应用

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

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. 多边形的面积

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. 求多边形的重心

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;
}