计算几何_三角剖分 POJ3675 望远镜

三角剖分

求一个圆和一个多边形的面积交.

前置知识：

用向量求圆和直线的交点.

POJ3675 望远镜

做法：

将多边形按照逆时针方向，圆心向所有顶点连边，剖分成多个三角形，然后分类讨论求三角形和圆的面积交即可。

1.A和B都在圆内

$S_{\cap }=S_{ABC}$

2.A和B都在圆外

$\theta =arccos(\frac{\vec{a}\cdot \vec{b}}{|\vec{a}|\cdot |\vec{b}|})$

$S_{\cap }=\frac{\theta r^2}{2}$

3.A在圆内B在圆外

$S=扇形(C{p}_{b}B)+S_{△AC{p}_b}$

4.A在圆外B在圆内

此情况与上一情况相同

5.A和B都在圆外的特殊情况

$S=S_{扇形AC{p}_a}+S_{扇形BC{p}_b}+S_{△p_{a}C{p}_b}$

所以我们先求出直线AB与圆的两个交点和圆心到线段AB的最短距离，然后进行讨论。

1.如果a和b距离圆心的距离都小于r，则说明是第一种情况

2.如果a和b的距离都大于r且圆心到ab线段的距离大于r说明是第二种情况

3.然后依次讨论da和mind的长度来确定是否是第三种或者第四种情况

4否则为最后一种情况。

double circle_triangle(Poi a,Poi b,Cir c={{0,0},r}){ //求圆c和三角形的面积交
	auto da=dist(c.p,a),db=dist(c.p,b);
	if(sign(c.r-da)>=0&&sign(c.r-db)>=0) return a*b/2;
	if(!sign(a*b)) return  0.0;
	 //直线ab和圆的交点
	auto mind=poi_to_segment(a,b,c.p);
	// debug(mind);
	vector<Poi> intersection=get_circle_line_intersection(a,b,c);
	Poi pa,pb;
	// if(intersection.size()<2) return 0.0; 
	if(intersection.size()==2)
		pa=intersection[0],pb=intersection[1];
	if(sign(c.r-mind)<=0)	return sector_area(a,b,c);
	if(sign(c.r-da) >= 0) {//a在圆内，b在圆外
		return (a-c.p)*(pb-c.p)/2.0+sector_area(pb,b,c);
	}
	if(sign(c.r-db) >= 0)
		return sector_area(a,pa,c)+(pa-c.p)*(b-c.p)/2.0;
	return sector_area(a,pa,c)+sector_area(pb,b,c)+(pa-c.p)*(pb-c.p)/2.0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

在这里顺便放一下我的计算几何模板:

Code:

#include
#include
#include
#include
#include
#include
#include
#include
#include 
#include
#include
#include
#include
#include
#include
#include
#include
#define x first
#define y second
using namespace std;
//================================DEBUG
#define FASTIO cin.tie(nullptr) -> sync_with_stdio(false)
#define debug(a) cout << #a": " << a << endl;
#define rep(i, ll, rr) for(int i = ll; i <= rr; ++ i)
#define per(i, rr, ll) for(int i = rr; i >= ll; -- i)
//================================IO
typedef long long LL; typedef unsigned long long ULL; typedef long double LD;
inline LL read(){LL s=0,w=1;char ch=getchar();for(;!isdigit(ch); ch = getchar())if(ch == '-') w = -1; for (; isdigit(ch);ch = getchar())s=(s<<1)+(s<<3)+(ch^48);return s*w;}
inline void print(LL x,int op=10){if(!x){putchar('0');if(op)putchar(op);return;}char F[40];LL tmp=x>0?x:-x;if(x<0)putchar('-');int cnt=0;while(tmp>0){F[cnt++]=tmp%10+'0';tmp/=10;}while(cnt>0)putchar(F[--cnt]);if(op)putchar(op);}
inline void print128(__int128_t x){if(x < 0) {putchar('-');x = -x;}if(x/10) print128(x/10);putchar(x%10+'0');}
template <typename T>void read(T &x){x=0;int f=1;char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+(ch^48);ch=getchar();}x*=f;return;}
template <typename T>void write(T x){if(x<0){putchar('-');x=-x;}if(x>9)write(x/10);putchar(x%10+'0');return;}      
//================================HABIT
LL fpower(LL a,LL b,LL mod) {LL ans = 1; while(b){ if(b & 1) ans = ans * (a % mod) % mod; a = a % mod * (a % mod) % mod; b >>= 1;} return ans; } 
LL Mod(LL a,LL mod){return (a%mod+mod)%mod;}
LL gcd(LL a,LL b) {return b?gcd(b,a%b):a;}
int mov[8][2]={1,0,0,1,-1,0,0,-1,1,1,-1,-1,1,-1,-1,1};
//================================DEFINE
// #define int LL
#define double long double
typedef pair<int,int> pii;
typedef pair<double,double> pdd;
const double eps=1e-9,PI=acos(-1);
//================================GEOMETRY
int sign(double x) {if(fabs(x)<eps) return 0;return x>0?1:-1;}
struct Poi{
	double x,y;
	Poi operator-(Poi b){return {x-b.x,y-b.y};}
	Poi operator+(Poi b){return {x+b.x,y+b.y};}
	Poi operator*(double k){return {x*k,y*k};}
	Poi operator/(double k){return {x/k,y/k};}
	Poi norm(){double d=sqrt(x*x+y*y);return {x/d,y/d};}
	double operator*(Poi b){return x*b.y-y*b.x;}
	double operator&(Poi b){return x*b.x+y*b.y;}
	bool operator==(Poi b){return sign(x-b.x)==0&&sign(y-b.y)==0;}
	bool operator<(Poi b){return sign(x-b.x)<0||(sign(x-b.x)==0&&sign(y-b.y)<0);}
};
double cross(Poi a,Poi b){return a.x*b.y-a.y*b.x;}
double area(Poi a,Poi b,Poi c){return cross({b.x-a.x,b.y-a.y},{c.x-a.x,c.y-a.y});}
double dist(Poi a,Poi b){double dx=a.x-b.x;double dy=a.y-b.y;return sqrt(dx*dx+dy*dy);}
struct Cir{Poi p;double r;};
struct Line{Poi st,ed;};
double get_angle(const Line &a){ //获得直线的角度
	return atan2(a.ed.y-a.st.y,a.ed.x-a.st.x);
}
//直线按照角度的排序函数
bool cmp(Line &a,Line &b){  double A=get_angle(a),B=get_angle(b); if(sign(A-B)==0) return sign(area(a.st,a.ed,b.ed))<0;return A<B;}
//求直线p+kv和直线q+kw的交点
Poi get_line_intersection(Poi p,Poi v,Poi q,Poi w){ auto u=p-q; double t=cross(w,u)/cross(v,w); return {p.x+v.x*t,p.y+v.y*t};}
//两条线的交点
Poi get_line_intersection(Line a,Line b){ return get_line_intersection(a.st,a.ed-a.st,b.st,b.ed-b.st);}
//bc的交点是否再a的右侧
bool on_right(Line a,Line b,Line c){ auto jiao=get_line_intersection(b,c);return sign(area(a.st,a.ed,jiao))<=0;}
//将一个点顺时针旋转d度
Poi rotate(Poi a,double b){return {a.x*cos(b)+a.y*sin(b), -a.x*sin(b)+a.y*cos(b)};}
//获取中垂线
Line get_perpendicular_bisector(Poi a,Poi b){return {(a+b)/2,rotate(b-a,PI/2.0)};}
//三点确定圆
Cir get_cir(Poi a,Poi b,Poi c){auto u=get_perpendicular_bisector(a,b),v=get_perpendicular_bisector(a,c);auto p=get_line_intersection(u.st,u.ed,v.st,v.ed);return {p, dist(p,a)};}
/*random_suffle(p+1,p+1+n);  点随机化*/
//================================
const int N=200010,M=N*2,mod=1e9+7;
vector<Poi> p;
Poi m,c;int n;
double r;

double len(Poi a){return sqrt(a&a);}

bool on_segment(Poi p,Poi a,Poi b){ //判断c是否在线段ab上
	return !sign((p-a)*(p-b)) && sign((p-a)&(p-b))<=0;
}

vector<Poi> get_circle_line_intersection(Poi a,Poi b,Cir c={{0,0},r}){ //线段ab和圆c的交点
	vector<Poi> ans;
	auto e=get_line_intersection(a, b-a, c.p, rotate(b-a,PI/2)); //弦与中垂线的交点
	auto d=dist(c.p, e); //弦心距
	if(!on_segment(e,a,b)) d=min(dist(c.p,a), dist(c.p, b)); 
	if(sign(c.r-d)<=0) return ans;
	auto len=sqrt(c.r*c.r-dist(c.p, e)*dist(c.p, e));
	Poi pa=e+(a-b).norm()*len,pb=e+((b-a).norm()*len);
	ans.push_back(pa);
	ans.push_back(pb);
	return ans;
}

double poi_to_segment(Poi a,Poi b,Poi c={0,0}){ //点到线段的距离
	auto e=get_line_intersection(a, b-a, c, rotate(b-a,PI/2)); //弦与中垂线的交点
	auto d=dist(c, e); //弦心距
	if(!on_segment(e,a,b)) d=min(dist(c,a), dist(c, b)); 
	return d;
}

double sector_area(Poi a,Poi b,Cir c){ //c为圆，acb扇形面积
	auto angle=acos((a&b)/len(a)/len(b));
	if(sign(a*b)<0) angle=-angle;
	return c.r*c.r*angle/2.0;
}

double circle_triangle(Poi a,Poi b,Cir c={{0,0},r}){ //求圆c和三角形的面积交
	auto da=dist(c.p,a),db=dist(c.p,b);
	if(sign(c.r-da)>=0&&sign(c.r-db)>=0) return a*b/2;
	if(!sign(a*b)) return  0.0;
	 //直线ab和圆的交点
	auto mind=poi_to_segment(a,b,c.p);
	// debug(mind);
	vector<Poi> intersection=get_circle_line_intersection(a,b,c);
	Poi pa,pb;
	// if(intersection.size()<2) return 0.0; 
	if(intersection.size()==2)
		pa=intersection[0],pb=intersection[1];
	if(sign(c.r-mind)<=0)	return sector_area(a,b,c);
	if(sign(c.r-da) >= 0) {//a在圆内，b在圆外
		return (a-c.p)*(pb-c.p)/2.0+sector_area(pb,b,c);
	}
	if(sign(c.r-db) >= 0)
		return sector_area(a,pa,c)+(pa-c.p)*(b-c.p)/2.0;
	return sector_area(a,pa,c)+sector_area(pb,b,c)+(pa-c.p)*(pb-c.p)/2.0;
}

double area(){
	double res=0;
	for(int i=0;i<n;++i){
		res += circle_triangle(p[i], p[(i+1)%n]);
	}
	return fabs(res);
}

void solve(){
	p.clear();
	p.resize(n);
	for(auto &u:p){
		scanf("%Lf%Lf",&u.x,&u.y);
	}
	printf("%.2Lf\n",area());
}
//===============================
signed main(){
	// int _=1;
	while(scanf("%Lf%d",&r,&n)!=-1)	
		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
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