计算几何_三角剖分 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;
}

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

Code:

#include<unordered_set>
#include<unordered_map>
#include<functional>
#include<algorithm>
#include<string.h>
#include<iostream>
#include<iterator>
#include<cstring>
#include<numeric> 
#include<cstdio>
#include<vector>
#include<bitset>
#include<queue>
#include<stack>
#include<cmath>
#include<set>
#include<map>
#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;
}