数学：计算几何二维平面基础

Posted on Edited on

Preface

学校计算几何专题刚刚好讲完，放一下计算几何用到的板子；

Content

Problem

前摇：常见常数，误差比较

由于计算几何常常带有浮点类型的运算经常等号判定会有偏差，所以这里采用误差比较的方式判定相等；

1
2
3
4
5
6
7
8
9
10
11
12
13
typedef double db;
constexpr double eps=1e-6;
const double pi=acos(-1);
int sign(double k){   
    if (k>eps) return 1; 
    else if (k<-eps) return -1; 
    return 0;
}
int dbcmp(double x,double y){
    if(y-x>eps) return 1;
    else if(y-x<-eps) return -1;
    return 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
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
template<typename T>struct point { // 点类
	T x, y;
	bool operator == (const point &a) const {
		return (abs(x - a.x) <= eps && abs(y - a.y) <= eps);
	}
	bool operator < (const point &a) const {
		if (abs(x - a.x) <= eps) return y < a.y;
		return x < a.x;
	}
	T length2() const {
		return (*this) * (*this);
	}
	T length() const {
		return sqrt(length2());
	}
	point unit() {
		double len = length();
		return {x / len, y / len};
	}
	point operator -() const {
		return { -x, -y};
	}
	point operator + (const point &a) const {
		return {x + a.x, y + a.y};
	}
	point operator - (const point &a) const {
		return {x - a.x, y - a.y};
	}
	point operator * (const T k) const {
		return {k * x, k * y};
	}
	point operator / (const T k) const {
		return {x / k, y / k};
	}
	T operator * (const point &a) const {
		return x * a.x + y * a.y;
	}
	T operator ^ (const point &a) const {
		return x * a.y - y * a.x;
	}
	point rotate(const double rad) const {
		double c = cos(rad), s = sin(rad);
		return {x*c - y * s, x*s + y * c};
	}
	point rotate90() const {
		return { -y, x};
	}
	double get_angle (const point &a) const {
		return atan2(y, x);
	}
	friend istream & operator >> (istream&, point &a) {
		scanf("%lf %lf", &a.x, &a.y);
		return cin;
	}
	friend ostream & operator << (ostream&, point &a) {
		printf(" ( %.6lf , %.6lf ) ", a.x, a.y);
		return cout;
	}
	// to-left测试
	int toleft (const point &a) const {
		const auto t = (*this)^a;
		return (t > eps) - (t < -eps);
	}
};

直线模板

1
2
3
4
5
6
7
8
9
template<typename T> struct line { 
	point<T>  p, v;
	bool operator == (const line &a) const {
		return abs(v ^ a.v) <= eps && abs(v ^ (p - a.p)) <= eps;
	}
	int toleft(const point<T> &a) const {
		return v.toleft(a - p);
	}
};

线段模板

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
template<typename T> struct segment {
	point<T> a, b;
	// -1 点在线段端点 | 0 点不在线段上 | 1 点严格在线段上
	int is_on(const point<T> &p) const {
		if (p == a || p == b) return -1;
		return (p - a).toleft(p - b) == 0 && (p - a) * (p - b) < -eps;
	}
    
	// -1 直线经过线段端点 | 0 线段和直线不相交 | 1 线段和直线严格相交
	int is_inter(const line<T> &l) const {
		if (l.toleft(a) == 0 || l.toleft(b) == 0) return -1;
		return l.toleft(a) != l.toleft(b);
	}
    
	// -1 在某一线段端点处相交 | 0 两线段不相交 | 1 两线段严格相交
	int is_inter(const segment<T> &s) const {
		if (is_on(s.a) || is_on(s.b) || s.is_on(a) || s.is_on(b)) 				return -1;
		const line<T> l {a, b - a}, ls {s.a, s.b - s.a};
		return l.toleft(s.a) * l.toleft(s.b) == -1 && 
            	ls.toleft(a) * ls.toleft(b) == -1;
	}

	// 点到线段距离
	long double dis(const point<T> &p) const {
		if ((p - a) * (b - a) < -eps || (p - b) * (a - b) < -eps) 				return min(p.dis(a), p.dis(b));
		const line<T> l {a, b - a};
		return l.dis(p);
	}

	// 两线段间距离
	long double dis(const segment<T> &s) const {
		if (is_inter(s)) return 0;
		return min({dis(s.a), dis(s.b), s.dis(a), s.dis(b)});
	}
};

普通多边形模板

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
template<typename T> struct polygon {
	std::vector<point<T>> p;//逆时针存储
	inline size_t nxt(const size_t i) const {
		return i == p.size() - 1 ? 0 : i + 1;
	}
	inline size_t pre(const size_t i) const {
		return i == 0 ? p.size() - 1 : i - 1;
	}
	pair<bool, int> winding(const point<T> &a) const {
		int cnt = 0;
		for (size_t i = 0; i < p.size(); i++) {
			point<T> u = p[i], v = p[nxt(i)];
			if (abs((a - u) ^ (a - v)) <= eps && (a - u) * (a - v) <= eps) return {true, 0}; 
             // 点a在线段uv上
			if (abs(u.y - v.y) <= eps) continue;
			line<T> uv = {u, v - u};
			if (u.y < v.y - eps && uv.toleft(a) <= 0) continue; // 方向向上且不在左侧
			if (u.y > v.y + eps && uv.toleft(a) >= 0) continue; // 方向向上且不在右侧
			if (u.y < a.y - eps && v.y >= a.y - eps) cnt++; // 向上穿过
			if (u.y >= a.y - eps && v.y < a.y - eps) cnt--; // 向下穿过
		}
		return {false, cnt};
	}

	// 计算周长
	double perimeter() const {
		double sum = 0;
		for (size_t i = 0; i < p.size(); i++) sum += dist(p[i], p[nxt(i)]);
		return sum;
	}
	// 计算面积
	double area() const {
		double sum = 0;
		for (size_t i = 0; i < p.size(); i++) sum += p[i] ^ p[nxt(i)];
		return abs(sum) / 2;
	}
};

圆模板

1
2
3
4
5
6
7
8
9
template<typename T> struct circle { 
	point<T>  o;
	T r;
	int is_in(const point<T> &q) const {
		point<T> qo=q-o;
		return (qo.length()-r)<=-eps;
	}
};

点与点之间（向量与向量）之间的关系

极角排序

利用库函数极角排序

先将所有点的极角算出来，然后按极角的值进行排序,不是将 atan2(y,x) 函数写在排序的比较函数里；

注意，atan2(double y,double x) 函数返回的是原点至点 $$( x , y )$$的方位角，即与与x轴正半轴的夹角，范围是$$(-\pi,\pi \ ]$$；

在极角排序中如果返回的极角值相同则按与原点距离升序；

相对于外积排序这个速度更快，但是可能有精度问题，如果一直$$WA$$就考虑把这个换掉；

注意：使用atan2（）应该先传y再传x；

1
2
3
4
5
6
7
8
9
10
11
bool cmp(point a, point b) {
    if(dcmp(atan2(a.y, a.x) - atan2(b.y, b.x)) == 0) 
        return a.x < b.x;
    return atan2(a.y, a.x) < atan2(b.y, b.x);
}//极点为原点时直接套用即可
point<double> o;
bool cmp(point a, point b) {
    if(dcmp(atan2(a.y-o.y, a.x-o.x) - atan2(b.y-o.y, b.x-o.x)) == 0) 
        return a.x < b.x;
    return atan2(a.y-o.y, a.x-o.x) < atan2(b.y-o.y, b.x-o.x);
} //改变api可以自定义极点

排序结果如图所示: $$A-B-C-D-E-F-G-H-I$$ ；

利用外积极角排序

一般利用外积可能会出现转一圈小于自己的情况，因为右手定则只能判断某两个的相对情况，所以我们利用象限作为排序的第一关键字就可以解决这个问题了；

相对于前者精度更优但是速度一般，如果题目有容错就不要写这个了；

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
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
void Convex(vector<point<db>>&Convex_A,vector<point<db>>A){
	sort(A.begin(),A.end(),[](point<db>&l,point<db>&r){
		return l.x==r.x?l.y<r.y:l.x<r.x;
	});
	for(int i=1;i<A.size();i++) A[i]=A[i]-A[0];
	sort(A.begin()+1,A.end(),[](const point<db>&l,const point<db>&r){
		return sign(l^r)==0?l.length()<r.length():sign(l^r)>0;
	});
	for(int i=1;i<A.size();i++) A[i]=A[i]+A[0];
	Convex_A.resize(A.size());
	Convex_A[0]=A[0];
	Convex_A[1]=A[1];
	int cnt=2;
	for(int i=2;i<A.size();i++){
		while(cnt>1&&(sign((Convex_A[cnt-1]-Convex_A[cnt-2])^(A[i]-Convex_A[cnt-1]))<=0)) cnt--;
		Convex_A[cnt++]=A[i];
	}
	Convex_A.resize(cnt);
	return ;
}

int inquadrant(const point<db> &v){
	if(v.x>=0&&v.y>0) return 1;
	else if(v.x<0&&v.y>=0) return 2;
	else if(v.x<=0&&v.y<0) return 3;
	else if(v.x>0&&v.y<=0) return 4;
	else return 0;
}

const int maxn=1e5+5;
point<double> p[maxn];
point<double> cvx[maxn];
int tot,n;
bool cmp(point<db> l,point<db> r){
    return (sign(l^r)==0)
	?(l.length()<r.length())
	:sign(l^r)>0;
}
signed main ()
{
	int n;
	cin>>n;
	for(int i=0;i<n;i++) cin>>p[i];
	sort(p,p+n);
	for(int i=1;i<n;i++) p[i]=p[i]-p[0];
	sort(p+1,p+n,cmp);
	for(int i=1;i<n;i++) p[i]=p[i]+p[0];
	cvx[0]=p[0];
	cvx[1]=p[1];
	tot=2;
	for(int i=2;i<n;i++){
		while(tot>1&&(sign((cvx[tot-1]-cvx[tot-2])^(p[i]-cvx[tot-1]))<=0)) tot--;
		cvx[tot++]=p[i];
	}
	cvx[tot]=cvx[0];
    db ccc=0;
    for(int i=0;i<tot;i++){
        ccc+=(cvx[i+1]-cvx[i]).length();
    }
    //printf("%.2lf",ccc);
	int m;
	cin>>m;
	while(m--){
		point<db> v;
		cin>>v;
		if(tot==1) {
			if(v.x!=cvx[0].x||v.y!=cvx[0].y) cout<<"NO"<<endl;
			else cout<<"YES"<<endl;
		}
		else if(tot==2){
			segment<db> l={cvx[0],cvx[1]};
			if(l.is_on(v)==0) cout<<"NO"<<endl;
			else cout<<"YES"<<endl;
		}
		else{
			point<db> oa=cvx[1]-cvx[0];
			point<db> ov=v-cvx[0];
			point<db> ob=cvx[tot-1]-cvx[0];
			if(sign(oa^ov)<0) cout<<"NO"<<endl;
			else if(sign(oa^ov)==0){
				segment<db> l={cvx[0],cvx[1]};
				if(l.is_on(v)==0) cout<<"NO"<<endl;
				else cout<<"YES"<<endl;
			}
			else if(sign(ob^ov)>0) cout<<"NO"<<endl;
			else if(sign(ob^ov)==0){
				segment<db> l={cvx[0],cvx[tot-1]};
				if(l.is_on(v)==0) cout<<"NO"<<endl;
				else cout<<"YES"<<endl;
			}
			else{
				int l=1,r=tot-1,mid=0;
				while(l<r){
					mid=(l+r)/2;
					point om=cvx[mid]-cvx[0];
					if(sign(om^ov)==0) l=r=mid;
					else if(sign(om^ov)<0) r=mid;
					else l=mid+1;
				}
				point lr=cvx[l]-cvx[l-1];
				point lv=cvx[l]-v;
				if(sign(lr^lv)<=0) cout<<"YES"<<endl;
				else cout<<"NO"<<endl;
			}
		}
	}
	return 0;
}

点定位（$$PIP$$问题）

点定位问题大致分为三种方法：外积二分法，光线投影法，$$Winding \ Number$$（回转数法）；

其中，外积二分基于凸多边形的三角剖分，常与点集的凸包联系起来，通过确定极点并二分极角序确定点的位置；光线投影法属于扫描线算法的一种；回转数法适用于任意多边形，但是其经典形式存在着一定的精度问题；

平面最近点对

判断点与直线（射线，线段）的位置关系

$$ToLeftTest$$ : 判断点是否在直线的左侧

在上面的点模板和直线模板中也有用到，所以也可以引用上面的结构体内联的函数；

$$ToLeftTest$$ 就是利用向量的外积判断正负，外积为正说明在右手转过的角小于 $$\pi$$ ,点在直线左侧，大于$$\pi$$ 则说明点在直线右侧，如果等于零，在直线两点和该点构成的三角形面积为零，因此在直线上；

1
2
3
int to_left(const line<double> &l, const point<double> &p) {
	return l.v.toleft(p - l.p);
}

点到直线距离

在直线上取两点可以转化为三角形的高，利用外积求出面积，利用距离公式求出底；

1
2
3
double dist(const line<double> &l, const point<double> &p) {
	return abs(l.v ^ (p - l.p)) / l.v.length();
}

点到直线的投影点

利用内积关于投影的性质，将向量运算转化为模运算；

1
2
3
point<double> projection(const line<double> &l, const point<double> &p) {
	return l.p + l.v * ((l.v * (p - l.p)) / (l.v * l.v));
}

两直线交点

1
2
3
point<double> intersection(const line<double> &a, const line<double> &b) {
	return a.p + a.v * ((b.v ^ (a.p - b.p)) / (a.v ^ b.v));
}

判断直线与直线（射线、线段）的关系

半平面以及半平面交

判断圆与点、线、圆的关系

两圆的位置关系：相等，相离，内含，外切，内切，相交；

倒不如说考虑问题应该按照这个顺序想，因为前面的特殊情况往往蕴含着不一样的答案；

1
2
3
4
5
6
7
8
9
10
11
12
//两圆关系：5相离，4外切，3相交，2内切，1内含 
int relationcircle(circle k1,circle k2){
    db d=k1.o.dis(k2.o);
    db l=fabs(k1.r-k2.r);
    if(sign(d)==0&&sign(l)==0) return 0;
    else if(sign(d-k1.r-k2.r)>0)return 5;
    else if(sign(d-k1.r-k2.r)==0)return 4;
    else if(sign(d-k1.r-k2.r)<0 && sign(d-l)>0)return 3;
    else if(sign(d-l)==0)return 2;
    else if(sign(d-l)<0)return 1; 
    else exit(0);// to check another error
}

圆与圆的面积交问题

很明显在最朴素相交情况下，两圆面积交可以表示为两个弓形的面积之和，这样我们就转化为容斥原理求解即可；

注意我们这里使用$$Heron$$公式求解三角形的面积，尽量绕开对反三角的求解；

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int relationcircle(circle k1,circle k2){
    db d=k1.o.dis(k2.o);
    db l=fabs(k1.r-k2.r);
    if(sign(d)==0&&sign(l)==0) return 0;
    else if(sign(d-k1.r-k2.r)>0)return 5;
    else if(sign(d-k1.r-k2.r)==0)return 4;
    else if(sign(d-k1.r-k2.r)<0 && sign(d-l)>0)return 3;
    else if(sign(d-l)==0)return 2;
    else if(sign(d-l)<0)return 1; 
    else exit(0);
}
db circlearea(circle k1,circle k2){
    int rel=relationcircle(k1,k2);
    if(rel>=4)return 0.0;
    if(rel<=2)return min(k1.area(),k2.area());
    db d=k1.o.dis(k2.o);
    db hf=(k1.r+k2.r+d)/2.0;
    db ss=2*sqrt(hf*(hf-k1.r)*(hf-k2.r)*(hf-d));//Heron
    db a1=acos((k1.r*k1.r+d*d-k2.r*k2.r)/(2.0*k1.r*d));
    a1=a1*k1.r*k1.r;
    db a2=acos((k2.r*k2.r+d*d-k1.r*k1.r)/(2.0*k2.r*d));
    a2=a2*k2.r*k2.r;
    return a1+a2-ss;
}

圆与多边形面积交问题

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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
//#pragma GCC optimize(2)
#include<bits/stdc++.h>
#include<iostream>
using namespace std;
typedef double db;
constexpr double eps=1e-5;
const double pi=acos(-1);
int sign(double k){   
    if (k>eps) return 1; 
    else if (k<-eps) return -1; 
    return 0;
}
int dbcmp(double x,double y){
    if(y-x>eps) return 1;
    else if(y-x<-eps) return -1;
    return 0;
}
template<typename T>struct point { // 点类
	T x, y;
	bool operator == (const point &a) const {
		return (abs(x - a.x) <= eps && abs(y - a.y) <= eps);
	}
	bool operator < (const point &a) const {
		if (abs(x - a.x) <= eps) return y < a.y;
		return x < a.x;
	}
	T length2() const {
		return (*this) * (*this);
	}
	T length() const {
		return sqrt(length2());
	}
	point unit() {
		double len = length();
		return {x / len, y / len};
	}
	point operator -() const {
		return { -x, -y};
	}
	point operator + (const point &a) const {
		return {x + a.x, y + a.y};
	}
	point operator - (const point &a) const {
		return {x - a.x, y - a.y};
	}
	point operator * (const T k) const {
		return {k * x, k * y};
	}
	point operator / (const T k) const {
		return {x / k, y / k};
	}
	T operator * (const point &a) const {
		return x * a.x + y * a.y;
	}
	T operator ^ (const point &a) const {
		return x * a.y - y * a.x;
	}
	point rotate(const double rad) const {
		double c = cos(rad), s = sin(rad);
		return {x*c - y * s, x*s + y * c};
	}
	point rotate90() const {
		return { -y, x};
	}
	double get_angle (const point &a) const {
		return atan2(y, x);
	}
	friend istream & operator >> (istream&, point &a) {
		scanf("%lf %lf", &a.x, &a.y);
		return cin;
	}
	friend ostream & operator << (ostream&, point &a) {
		printf(" ( %.6lf , %.6lf ) ", a.x, a.y);
		return cout;
	}
	// to-left测试
	int toleft (const point &a) const {
		const auto t = (*this)^a;
		return (t > eps) - (t < -eps);
	}
};


template<typename T> struct line { 
	point<T>  p, v;
	bool operator == (const line &a) const {
		return abs(v ^ a.v) <= eps && abs(v ^ (p - a.p)) <= eps;
	}
	int toleft(const point<T> &a) const {
		return v.toleft(a - p);
	}
	int is_on(const point<T> &q) const {
		point<T> pq=p-q;
		return abs((pq^v))<=eps;
	}
};
template<typename T> struct segment {
	point<T> a, b;
	// -1 点在线段端点 | 0 点不在线段上 | 1 点严格在线段上
	int is_on(const point<T> &p) const {
		if (p == a || p == b) return -1;
		return (p - a).toleft(p - b) == 0 && (p - a) * (p - b) < -eps;
	}
    
	// -1 直线经过线段端点 | 0 线段和直线不相交 | 1 线段和直线严格相交
	int is_inter(const line<T> &l) const {
		if (l.toleft(a) == 0 || l.toleft(b) == 0) return -1;
		return l.toleft(a) != l.toleft(b);
	}
    
	// -1 在某一线段端点处相交 | 0 两线段不相交 | 1 两线段严格相交
	int is_inter(const segment<T> &s) const {
		if (is_on(s.a) || is_on(s.b) || s.is_on(a) || s.is_on(b)) 				
		return -1;
		const line<T> l {a, b - a}, ls {s.a, s.b - s.a};
		return l.toleft(s.a) * l.toleft(s.b) == -1 && 
            	ls.toleft(a) * ls.toleft(b) == -1;
	}

	// 点到线段距离
	double dis(const point<T> &p) const {
		if ((p - a) * (b - a) < -eps || (p - b) * (a - b) < -eps) 				
		return min(p.dis(a), p.dis(b));
		const line<T> l {a, b - a};
		return l.dis(p);
	}

	// 两线段间距离
	double dis(const segment<T> &s) const {
		if (is_inter(s)) return 0;
		return min({dis(s.a), dis(s.b), s.dis(a), s.dis(b)});
	}
};
template<typename T> struct circle { 
	point<T>  o;
	T r;
	int is_in(const point<T> &q) const {
		point<T> qo=q-o;
		return (qo.length()-r)<-eps;
	}
};

int getSegCircleIntersection(line<db> L,circle<db> C,
vector<point<db>>& sol) {
	double a=L.v.x;
	double b=L.p.x-C.o.x;
	double c=L.v.y;
	double d=L.p.y-C.o.y;
	double e=a*a+c*c;
	double f=2*(a*b + c*d);
	double g=b*b+d*d-C.r*C.r;
	double delta=f*f-4*e*g;
	double t1,t2;
	int ans = 0;
	if(sign(delta)<0) return 0;
	if(sign(delta)==0) {
		t1=t2=-f/(2*e);
		if(sign(t1)>=0&&sign(t1-1)<=0){
			ans++;
			sol.push_back(L.p+L.v*t1);
		}
		return ans;
	}
	t1=(-f-sqrt(delta))/(2*e);
	t2=(-f+sqrt(delta))/(2*e);
	if(t1>t2) swap(t1,t2);
	if(sign(t1)>=0&&sign(t1-1)<=0){
		ans++;
		sol.push_back(L.p+L.v*t1);
	}
	if(sign(t2)>=0&&sign(t2-1)<= 0) {
		ans++;
		sol.push_back(L.p+L.v*t2);
	}
	return ans;
}
double TriangleArea(point<db>A,point<db> B,point<db> C) {
	return (double)(fabs((B-A)^(C-A)) )/ 2;
}
double Angle(point<db> A, point<db> B) {
	if(sign(A^B)==0) return 0;
	return (double)acos(A*B / A.length() / B.length());
}
double IntersectionArea(circle<db> C,point<db> A,point<db> B) {
	if((A==B)|| (B ==C.o)||(C.o==A)) return 0;
	line<db> L={A,B-A};
	int cnt=0;
	bool inA, inB;
	if(inA=C.is_in(A)) cnt++;
	if(inB=C.is_in(B)) cnt++;
	if(cnt == 2) return TriangleArea(C.o, A, B);
	if(cnt == 1) {
		vector<point<db>> q;
		getSegCircleIntersection(L, C, q);
		if(inB) swap(A, B);
		double theta = Angle(q[0]-C.o, B-C.o);
		return C.r*C.r*theta/2 + TriangleArea(C.o, A, q[0]);
	}

	vector<point<db>> q;
	int sz = getSegCircleIntersection(L, C, q);
	if(sz <= 1) {
		double theta = Angle(A-C.o, B-C.o);
		return C.r*C.r*theta/2;
	}

	double theta = Angle(C.o-A, C.o-q[0]) + Angle(C.o-B,C.o-q[1]);
	return C.r*C.r*theta/2 + TriangleArea(C.o,q[0], q[1]);
}
void solve(){
    int n;
    db k;
    cin>>n>>k;
    vector<point<db>> poly(n);
	for(int i=0;i<n;i++) cin>>poly[i];
	point<db> T,S;
	cin>>T>>S;
	// point<db> X=S*(1/(k+1))+T*(k/(k+1));
	// point<db> Y=S*(1/(1-k))+T*(k/(k-1));
	point<db> M=S*(1/(1-k*k))-T*((k*k)/(1-k*k));
	db rr=(S-T).length()*k/(1-k*k);
	circle<db> C={M,rr};

	double ans=0;
	for(int i=0;i<n;i++){
		int sgn;
		if( sign(((poly[i]-C.o)^(poly[(i+1)%n]-C.o)))>0) sgn = 1;
		else sgn=-1;
		ans+=sgn*IntersectionArea(C, poly[i], poly[(i+1)%n]);
	}
	printf("%e\n",fabs(ans));
	return ;
}
int main ()
{
    //ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    #ifdef local
    freopen("in.txt","r",stdin);
    freopen("out1.txt","w",stdout);
    #endif // local

    int t;
    cin>>t;
    while(t--) solve();
    return 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
36
int get_tangents(circle A, circle B, pnode *a, pnode *b){
    int cnt = 0;        //存切点用
    if(dcmp(A.r - B.r) < 0)
    {
        swap(A, B);
        swap(a, b);
    }
    double d = sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));     //圆心距
    double rdiff = A.r - B.r;      //两圆半径差
    double rsum  = A.r + B.r;       //两圆半径和
    if(dcmp(d - rdiff) < 0)
        return 0;        //1.内含
    double base = atan2(B.y - A.y, B.x - A.x);      //向量AB的极角
    if(dcmp(d) == 0) return -1;        //2.重合
    if(dcmp(d - rdiff) == 0)
    {       //3.内切
        a[cnt] = b[cnt] = A.point(base);
        cnt++;
        return 1;
    }
    double ang = acos((A.r - B.r) / d);
    a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;      //4.相交（外切、外离的外公切线也在此求出）
    a[cnt] = A.point(base - ang); b[cnt] = B.point(base - ang); cnt++;      //两条外公切线的切点
    if(dcmp(d - rsum) == 0)
    {        //5.外切
        a[cnt] = b[cnt] = A.point(base);
        cnt++;
    }
    else if(dcmp(d - rsum) > 0)
    {      //6.外离
        double ang = acos((A.r + B.r) / d);
        a[cnt] = A.point(base + ang); b[cnt] = B.point(PI + base + ang); cnt++；
        a[cnt] = A.point(base - ang); b[cnt] = B.point(PI + base - ang); cnt++;
    }
    return cnt;
}

圆与直线的反演变换

最小圆覆盖

$$Ἀπολλώνιος$$圆

凸集问题

动态凸包

Minkowski Sum

来测下板子啊

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void Minkowski(vector<point<db>>&s,vector<point<db>> a,vector<point<db>> b){
    int n=a.size();
    int m=b.size();
    vector<point<db>> s1,s2;
    for(int i=0;i<n;i++) s1.push_back(a[(i+1)%n]-a[i]);
    for(int i=0;i<m;i++) s2.push_back(b[(i+1)%m]-b[i]);
    int i=0,j=0,cnt=0;
    s.push_back(a[0]+b[0]);
    while(i<n&&j<m){
        if(sign(s1[i]^s2[j])>=0) s.push_back(s[cnt++]+s1[i++]);
        else s.push_back(s[cnt++]+s2[j++]);
    }
    while(i<n) s.push_back(s[cnt++]+s1[i++]);
    while(j<m) s.push_back(s[cnt++]+s2[j++]);
}

基本几何量的计算

三角形的有向面积（逆时针为正）

我们知道，两个向量的外积表示某个平行四边形的面积，那么我们只需要计算三角形两条边的向量外积的一半就是这个三角形的有向面积，实际面积即为有向距离的绝对值；

1
2
3
4
5
6
double triarea(Point p, Point q, Point s){
     return 
       0.5*(p.x * q.y - p.y * q.x
       + q.x * s.y - q.y * s.x
       + s.x * p.y - s.y * p.x);
 }

三角形$$Heron$$公式（三个顶点不知情况下）

Pick定理与整点问题

计算几何与数据结构

KD-Tree

Example

Remark

1