[模板]计算几何

Posted on Edited on In

[TOC]

二维

起手式(待补充):

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
typedef long long ll;
typedef long double ld;
const double PI = 3.14159265358979; 
const int maxn = 1e5+10;
struct Point{
	ld x,y;
	Point(ld x = 0,ld y = 0):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 * (ld w){return Point(x*w,y*w);}
	Point operator / (ld w){return Point(x/w,y/w);}
	Point Rotate_90_clock(){return Point(y,-x);}//顺时针旋转90°
	ld len2(){return x*x + y*y;}//长度的平方
	ld len(){return(sqrt(x*x + y*y));}
    bool operator < (const Point & A) const{
		return (x == A.x?y < A.y:x < A.x);
	}
}point[maxn];
struct Circle{
	Point O;
	ld R;
	Circle(Point O = Point(0,0),ld R = 0):O(O),R(R){}
	ld R2(){return R*R;}
};
typedef Point Vector;
ld Dot(Vector A,Vector B){
	return A.x*B.x + A.y*B.y;
}//点积
ld Cross(Vector A,Vector B){
	return A.x*B.y - A.y*B.x;
}//二维叉积
Point Line_inter(Point P0,Point V0,Point P1,Point V1){
	ld t = Cross(P1 - P0,V1)/Cross(V0,V1);
	return P0 + V0 * t;
}//返回两向量交点.P和V分别为向量起点终点*
Circle p2circle(Point A,Point B){
	Circle tans = Circle((A+B)/2,0);
	tans.R = (tans.O - A).len();
	return tans;
}//两点的外接圆
Circle p2circle(Point A,Point B,Point C){
	Circle tans = Circle(Line_inter((A+B)/2,(B - A).Rotate_90_clock(),(A+C)/2,(C - A).Rotate_90_clock()),0);
	tans.R = (tans.O - A).len();
	return tans;
}//三点共圆
typedef Point Vector;
inline double to_rad(double deg){
	return deg*(PI/180);
}
Vector Rotate(Vector A,double rad){//逆时针 
	return Vector(A.x*cos(rad) - A.y*sin(rad),A.x*sin(rad) + A.y*cos(rad));
}

点到直线及线段距离

1
2
3
4
5
6
7
8
9
10
11
ld LinePointDist(Point p,Point a,Point b){//点到直线
	Vector v1 = b-a,v2 = p-a;
	return fabs(Cross(v1,v2) / len(v1));
}
ld SegPointDist(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) < eps)return len(v2);
	if(Dot(v1,v3) > eps)return len(v3);
	return LinePointDist(p,a,b);
}

按角度旋转点

设点\(A(x,y)\)绕点\(R(x_{r},y_{r})\)逆时针旋转\(\theta\)弧度后得到点\(A’(x',y')\), 则有: \[ x' = (x - x_r)\cdot\cos(\theta) - (y - y_r)\cdot\sin(\theta) + x_r\\ y' = (x - x_r)\cdot\sin(\theta) + (y - y_r)\cdot\cos(\theta) + y_r \]

1
2
3
4
5
void rotate(ld & x,ld & y,ld theta,ld rx,ld ry){
    ld nx = (x - rx)*cos(theta) - (y - ry)*sin(theta) + rx;
    ld ny = (x - rx)*sin(theta) + (y - ry)*cos(theta) + ry;
    x = nx,y = ny;
}

凸包

输入顶点, 返回凸包. raw为原顶点数组,hull为凸包顶点数组

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
vector<Point> ConvexHull(vector<Point>raw){
    if(raw.size() == 1)return raw;
	sort(raw.begin(),raw.end());
	int n = raw.size(),idx = 0;
	vector<Point>hull(2*n+10);
	for(int i = 0;i < n;i++){
		while(idx > 1 and Cross(hull[idx-1] - hull[idx-2],raw[i] - hull[idx-1]) <= 0)idx--;
		hull[idx++] = raw[i];
	}
	int k = idx;
	for(int i = n-2;i >= 0;i--){
		while(idx > k and Cross(hull[idx-1] - hull[idx-2],raw[i] - hull[idx-1]) <= 0)idx--;
		hull[idx++] = raw[i];
	}
	hull.resize(idx-1);
	return hull;
}

闵可夫斯基和

点集\(A,B\)的闵可夫斯基和为\(C = \{a+b|a\in A,b \in B\}\)

img

复杂度\(O(n)\)

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
vector<Point> minkowski(vector<Point> & C1,vector<Point> & C2){
	int n = (int)C1.size(),m = (int)C2.size();
	vector<Point>s1(n),s2(m),A(n+m+10);
	for(int i = 0;i < n-1;i++)s1[i] = C1[i+1] - C1[i];
	s1[n-1] = C1[0] - C1[n-1];
	for(int i = 0;i < m-1;i++)s2[i] = C2[i+1] - C2[i];
	s2[m-1] = C2[0] - C2[m-1];
	int tot = 0;
	A[tot] = C1[0] + C2[0];
	int p1 = 0,p2 = 0;
	while(p1 < n and p2 < m){
		tot++;
		A[tot] = A[tot-1] + (Cross(s1[p1],s2[p2]) >= eps?s1[p1++]:s2[p2++]);
	}
	while(p1 < n){
		++tot;
		A[tot] = A[tot-1] + s1[p1++];
	}
	while(p2 <= m){
		++tot;
		A[tot] = A[tot-1] + s2[p2++];
	}
	A.resize(tot+1);
	return A;
}

多边形面积

_Point为多边形点集,n为顶点数

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

最小圆覆盖

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double ld;
const int maxn = 1e5+10;
struct Point{
	ld x,y;
	Point(ld x = 0,ld y = 0):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 * (ld w){return Point(x*w,y*w);}
	Point operator / (ld w){return Point(x/w,y/w);}
	Point Rotate_90_clock(){return Point(y,-x);}//顺时针旋转90°
	ld len2(){return x*x + y*y;}//长度的平方
	ld len(){return(sqrt(x*x + y*y));}
}point[maxn];
struct Circle{
	Point O;
	ld R;
	Circle(Point O = Point(0,0),ld R = 0):O(O),R(R){}
	ld R2(){return R*R;}
};
typedef Point Vector;
ld Dot(Vector A,Vector B){
	return A.x*B.x + A.y*B.y;
}//点积
ld Cross(Vector A,Vector B){
	return A.x*B.y - A.y*B.x;
}//二维叉积
Point Line_inter(Point P0,Point V0,Point P1,Point V1){
	ld t = Cross(P1 - P0,V1)/Cross(V0,V1);
	return P0 + V0 * t;
}//返回两向量交点.P和V分别为向量起点终点*
Circle p2circle(Point A,Point B){
	Circle tans = Circle((A+B)/2,0);
	tans.R = (tans.O - A).len();
	return tans;
}//两点的外接圆
Circle p2circle(Point A,Point B,Point C){
	Circle tans = Circle(Line_inter((A+B)/2,(B - A).Rotate_90_clock(),(A+C)/2,(C - A).Rotate_90_clock()),0);
	tans.R = (tans.O - A).len();
	return tans;
}//三点共圆
Circle Min_Cir_Cover(Point * point,int n){
	random_shuffle(point+1,point+n+1);
	Circle res;
	for(int i = 1;i <= n;i++){
		if((point[i] - res.O).len2() > res.R2()){
			res = Circle(point[i],0);
			for(int j = 1;j < i;j++){
				if((point[j] - res.O).len2() > res.R2()){
					res = p2circle(point[i],point[j]);
					for(int k = 1;k < j;k++){
						if((point[k] - res.O).len2() > res.R2()){
							res = p2circle(point[i],point[j],point[k]);
						}
					}
				}
			}
		}
	}
	return res;
}//最小圆覆盖,输入点集数组和点的数量,返回一个圆

判断点在凸包内部及边界

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
bool Left(Point P,Point A,Point B) {
	return Cross(B - A,P - A) > 0;
}
bool Right(Point P,Point A,Point B) {
	return Cross(P - A,B - A) > 0;
}
bool Point_in_Convex(Point P,vector<Point> & C){
	int n = C.size();
	if (Left(P,C[0],C[n - 1]) or Right(P,C[0],C[1]))return 0;
	int l = 1, r = n - 1;
	while (l < r - 1) {
		int m = (l + r) / 2;
		if (Left(P, C[0], C[m])) {
			l = m;
		}
		else {
			r = m;
		}
	}
	if (!Right(P, C[l], C[r]))return 1;
	return 0;
}

三维

起手式:

鸽子biss

叉乘写为矩阵形式

\[ \vec a \times \vec b = A^* \vec b = \left( \begin{array}{cccc} 0 & -z_1 & y_1\\ z_1 & 0 & -x_1\\ -y_1 & x_1 & 0 \end{array} \right) \left( \begin{array}{cccc} x_2\\ y_2\\ z_2 \end{array} \right) \]

\(A^*\)称为dual martix

最小球覆盖

来自上交红书《ACM国际大学生程序设计竞赛-算法与实现》

下标从0开始

复杂度: \(O(n)\)

输入:

​ npoint 全局变量, 点的个数

​ pt 全局变量, 点的坐标

调用 smallest_ball()

输出:

​ res 全局变量, 球心坐标

​ radius 全局变量, 球的半径

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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
//
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
const int maxn = 2e5+10;
double eps = 1e-9;
struct Tpoint{
	double x,y,z;
	Tpoint(double x = 0,double y = 0,double z = 0):x(x),y(y),z(z){}
};
int npoint,nouter;
Tpoint pt[maxn],outer[4],res;
double radius,tmp;
inline double dist(Tpoint p1,Tpoint p2){
	double dx = p1.x - p2.x,dy = p1.y - p2.y,dz = p1.z - p2.z;
	return (dx*dx + dy*dy + dz*dz);
}
inline double dot(Tpoint p1,Tpoint p2){
	return p1.x*p2.x + p1.y*p2.y + p1.z*p2.z;
}
void ball(){
	Tpoint q[3];double m[3][3], sol[3],L[3],det;
	res.x = res.y = res.z = radius = 0;
	switch(nouter){
		case 1: res = outer[0];break;
		case 2:
			res.x = (outer[0].x+outer[1].x)/2;
			res.y = (outer[0].y+outer[1].y)/2;
			res.z = (outer[0].z+outer[1].z)/2;
			radius = dist(res,outer[0]);
			break;
		case 3:
			for(int i = 0;i < 2;i++){
				q[i].x = outer[i+1].x - outer[0].x;
				q[i].y = outer[i+1].y - outer[0].y;
				q[i].z = outer[i+1].z - outer[0].z;
			}
			for(int i = 0;i < 2;i++)
				for(int j = 0;j < 2;j++)m[i][j] = dot(q[i],q[j])*2;
			for(int i = 0;i < 2;i++)sol[i] = dot(q[i],q[i]);
			if(fabs(det = m[0][0]*m[1][1] - m[0][1]*m[1][0])<eps)return;
			L[0] = (sol[0]*m[1][1] - sol[1]*m[0][1])/det;
			L[1] = (sol[1]*m[0][0] - sol[0]*m[1][0])/det;
			res.x = outer[0].x + q[0].x*L[0] + q[1].x*L[1];
			res.y = outer[0].y + q[0].y*L[0] + q[1].y*L[1];
			res.z = outer[0].z + q[0].z*L[0] + q[1].z*L[1];
			radius = dist(res,outer[0]);
			break;
		case 4:
			for(int i = 0;i < 3;i++){
				q[i].x = outer[i+1].x - outer[0].x;
				q[i].y = outer[i+1].y - outer[0].y;
				q[i].z = outer[i+1].z - outer[0].z;
				sol[i] = dot(q[i],q[i]);
			}
			for(int i = 0;i < 3;i++)
				for(int j = 0;j < 3;j++)m[i][j] = dot(q[i],q[j])*2;
			det = m[0][0] * m[1][1] * m[2][2]
				+ m[0][1] * m[1][2] * m[2][0]
				+ m[0][2] * m[2][1] * m[1][0]
				- m[0][2] * m[1][1] * m[2][0]
				- m[0][1] * m[1][0] * m[2][2]
				- m[0][0] * m[1][2] * m[2][1];
			if(fabs(det) < eps)return;
			for(int j = 0;j < 3;j++){
				for(int i = 0;i < 3;i++)m[i][j] = sol[i];
				L[j] = (m[0][0] * m[1][1] * m[2][2]
					+ m[0][1] * m[1][2] * m[2][0]
					+ m[0][2] * m[2][1] * m[1][0]
					- m[0][2] * m[1][1] * m[2][0]
					- m[0][1] * m[1][0] * m[2][2]
					- m[0][0] * m[1][2] * m[2][1]
					)/det;
				for(int i = 0;i < 3;i++)m[i][j] = dot(q[i],q[j])*2;
			}
			res = outer[0];
			for(int i = 0;i < 3;i++){
				res.x += q[i].x * L[i];
				res.y += q[i].y * L[i];
				res.z += q[i].z * L[i];
			}
			radius = dist(res,outer[0]);
	}
}
void minball(int n){
	ball();
	if(nouter < 4){
		for(int i = 0;i < n;i++){
			if(dist(res,pt[i]) - radius > eps){
				outer[nouter] = pt[i];
				++nouter;
				minball(i);
				--nouter;
				if(i > 0){
					Tpoint Tt = pt[i];
					memmove(&pt[1],&pt[0],sizeof(Tpoint)*i);
					pt[0] = Tt;
				}
			}
		}
	}
}
double smallest_ball(){
	radius = -1;
	for(int i = 0;i < npoint;i++){
		if(dist(res,pt[i]) - radius > eps){
			nouter = 1;
			outer[0] = pt[i];
			minball(i);
		}
	}
	return sqrt(radius);
}
signed main(){
	cin >> npoint;
	for(int i = 0;i < npoint;i++){
		int x,y,z;
		cin >> x >> y >> z;
		pt[i] = Tpoint(x,y,z);
	}
	printf("%.10lf",smallest_ball());
	return 0;
}