计算几何

发表于 | 分类于 | 阅读次数:

计算几何一类的题目不仅需要精妙的算法设计，还要注意精度等细节问题，常常是比赛中拉分的题目。

所以毕叔叔说：NOI见到计算几何就弃了吧

##向量
加减 若 $\mathbf{a}=(x_a,y_a),\mathbf{b}=(x_b,y_b)$ 则 $\mathbf{a}+\mathbf{b}=(x_a+x_b,y_a+y_b),\mathbf{a}-\mathbf{b}=(x_a-x_b,y_a-y_b)$
点乘 $\mathbf{a}\cdot\mathbf{b}=|a||b|\cos<\mathbf{a},\mathbf{b}>=(x_ax_b+y_ay_b)$
叉乘 $\mathbf{a}\times\mathbf{b}=(x_ay_b-x_by_a)$
垂直 若两向量的点积为 $0$ ，则两向量垂直
平行 若两向量的叉积为 $0$ ，则两向量平行

1
2
3
4
5
6
7
8
9
10
struct P
{
	double x,y;
	P(double x=0,double y=0):x(x),y(y){};
}
P operator +(P a,P b){return P(a.x+b.x,a.y+b.y);}
P operator -(P a,P b){return P(a.x-b.x,a.y-b.y);}
double dot(P a,P b){return a.x*b.x+a.y*b.y;}
double cross(P a,P b){return a.x*b.y-a.y*b.x;}
double dis(P a,P b){return std::sqrt(cross(b-a,b-a));}

##线段
表示方法：记录两端点，根据需要记录极角

1
2
3
4
5
6
7
struct L
{
	P l,r;
	double ang;
	L(P l=P(),P r=P()):l(l),r(r)
	{ang=std::atan2(r.y-l.y,r.x-l.x);};
};

###求交点

1
2
3
4
5
6
7
P inter(L a,L b)
{
	 double s1=cross(b.l-a.l,a.r-a.l),
	 		s2=cross(a.r-a.l,b.r-a.l);
	 return P((b.l.x*s2+b.r.x*s1)/(s1+s2),
	 		 (b.l.y*s2+b.r.y*s1)/(s1+s2));
}

##圆
表示：记录圆心和半径

1
2
3
4
5
struct C
{
    P cen;double r;
    C(P cen=P(),double r=0):cen(cen),r(r){};
};

两圆间的位置关系

  • 内切
  • 外切
  • 包含
  • 相交
  • 相离

两圆相交的弧度

求两圆的交点计算繁琐，而我们却可以很简单地知道两圆交点对各自的弧度。
若 a 和 b 两圆相交，那么对于

1
2
3
4
5
6
7
8
void inter(C a,C b)
{
    double d=dis(a.cen,b.cen),
           ang0=atan2(a.cen.x-b.cen.x,a.cen.y-b.cen.y),
           x=(a.r*a.r-b.r*b.r+d*d)/(2*d),
           angDet=acos(x/a.r);
    ang0-angDet,ang0+angDet
}

###[BZOJ1043]下落的圆盘

有n个圆盘从天而降，后面落下的可以盖住前面的。求最后形成的封闭区域的周长。看下面这副图, 所有的红色线条的总长度即为所求.
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
/**************************************************************
    Problem: 1043
    User: zhangche0526
    Language: C++
    Result: Accepted
    Time:292 ms
    Memory:1344 kb
****************************************************************/
 
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
const int MAXN=1e3+5;
const double INF=1e60,PI=acos(-1);
 
int n;
 
struct P
{
    double x,y;
    P(double x=0,double y=0):x(x),y(y){};
};
P operator +(P a,P b){return P(a.x+b.x,a.y+b.y);}
P operator -(P a,P b){return P(a.x-b.x,a.y-b.y);}
double cross(P a,P b){return a.x*b.x+a.y*b.y;}
double dis(P a,P b){return std::sqrt(cross(b-a,b-a));}
 
struct C
{
    P cen;double r;
    C(P cen=P(),double r=0):cen(cen),r(r){};
} cir[MAXN];
struct R{double l,r;R(double l=0,double r=0):l(l),r(r){};} ran[MAXN];int rcnt;
bool cmp(const R &a,const R &b){return a.l<b.l;}
inline bool cntn(C a,C b){return a.r>=b.r+dis(a.cen,b.cen);}//a contains b
void inter(C a,C b)
{
    double d=dis(a.cen,b.cen),
           ang0=atan2(a.cen.x-b.cen.x,a.cen.y-b.cen.y),
           x=(a.r*a.r-b.r*b.r+d*d)/(2*d),
           angDet=acos(x/a.r);
    ran[++rcnt]=R(ang0-angDet,ang0+angDet);
}
double calc(int x)
{
    int i;
    for(i=x+1;i<=n;i++)
        if(cntn(cir[i],cir[x])) return 0;
    rcnt=0;
    for(i=x+1;i<=n;i++)
        if(!cntn(cir[x],cir[i])&&cir[x].r+cir[i].r>dis(cir[x].cen,cir[i].cen))
            inter(cir[x],cir[i]);
    double sum=0,now=0;
    for(i=1;i<=rcnt;i++)
    {
        if(ran[i].l<0) ran[i].l+=2*PI;
        if(ran[i].r<0) ran[i].r+=2*PI;
        if(ran[i].l>ran[i].r) ran[++rcnt]=R(0,ran[i].r),ran[i].r=2*PI;
    }
    std::sort(ran+1,ran+rcnt+1,cmp);
    for(i=1;i<=rcnt;i++)
        if(ran[i].l>now)
            sum+=ran[i].l-now,now=ran[i].r;
        else now=std::max(now,ran[i].r);
    sum+=2*PI-now;
    return sum*cir[x].r;
}
 
int main()
{
    int i;
    scanf("%d",&n);
    for(i=1;i<=n;i++)
    {
        double x,y,r;scanf("%lf%lf%lf",&r,&x,&y);
        cir[i]=C(P(x,y),r);
    }
    double ans=0;
    for(i=1;i<=n;i++) ans+=calc(i);
    printf("%.3lf\n",ans);
}

##凸包
给定二维平面上的点集，凸包就是将最外层的点连接起来构成的凸多边型，它能包含点集中所有的点。
直接看裸题：
###[COGS896]圈奶牛

农夫约翰想要建造一个围栏用来围住他的奶牛，可是他资金匮乏。他建造的围栏必须包括他的奶牛喜欢吃草的所有地点。对于给出的这些地点的坐标，计算最短的能够围住这些点的围栏的长度。
地点数 $N\leq10^4$

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
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
 
const int MAXN=1e4+5;
int n;
 
struct P
{
	double x,y;
	P(double x=0,double y=0):x(x),y(y){};
} pt[MAXN];
P operator +(P a,P b){return P(a.x+b.x,a.y+b.y);}
P operator -(P a,P b){return P(a.x-b.x,a.y-b.y);}
double cross(P a,P b){return a.x*b.y-a.y*b.x;}
double dot(P a,P b){return a.x*b.x+a.y*b.y;}
double dis(P a,P b){return sqrt(dot(b-a,b-a));}
 
P ch[MAXN];int ccnt;
bool cmp(const P &a,const P &b)
{return a.x<b.x||(a.x==b.x&&a.y<b.y);}
 
void calCH()
{
	int i;
	std::sort(pt+1,pt+n+1,cmp);
	for(i=1;i<=n;i++)
	{
		while(ccnt>1&&cross(ch[ccnt]-ch[ccnt-1],pt[i]-ch[ccnt-1])<=0)
			ccnt--;
		ch[++ccnt]=pt[i];
	}
	int tmp=ccnt;
	for(i=n-1;i;i--) 
	{
		while(ccnt>tmp&&cross(ch[ccnt]-ch[ccnt-1],pt[i]-ch[ccnt-1])<=0)
			ccnt--;
		ch[++ccnt]=pt[i];
	}
	if(n>1) ccnt--;
}
 
int main()
{
	freopen("fc.in","r",stdin);
	freopen("fc.out","w",stdout);
	int i;
	scanf("%d",&n);
	for(i=1;i<=n;i++)
	{
		double x,y;scanf("%lf%lf",&x,&y);
		pt[i]=P(x,y);
	}
	calCH();
	double ans=0;
	for(i=1;i<=ccnt;i++) ans+=dis(ch[i],ch[i+1]);
	printf("%.2lf\n",ans);
	return 0;
}

###[POJ1113]Wall
花式水题，其实就是求个凸包，只不过每个点是圆，其实我们可以求出凸包的周长加上一个圆的周长即可。

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
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>

const int MAXN=1005;
const double PI=std::acos(-1);

struct P
{
	double x,y;
	P(double x=0,double y=0):x(x),y(y){};
};
P operator +(P a,P b){return P(a.x+b.x,a.y+b.y);}
P operator -(P a,P b){return P(a.x-b.x,a.y-b.y);}
double dot(P a,P b){return a.x*b.x+a.y*b.y;}
double cross(P a,P b){return a.x*b.y-a.y*b.x;}
double dis(P a,P b){return std::sqrt(dot(b-a,b-a));}
bool cmp(const P &a,const P &b)
{return a.x<b.x||(a.x==b.x&&a.y<b.y);}

int N,L;
P pt[MAXN],ch[MAXN];int ccnt;
inline void calCH()
{
	int i;
	std::sort(pt+1,pt+N+1,cmp);
	for(i=1;i<=N;i++)
	{
		while(ccnt>1&&cross(ch[ccnt]-ch[ccnt-1],pt[i]-ch[ccnt-1])<=0)
			ccnt--;
		ch[++ccnt]=pt[i];
	}
	int tmp=ccnt;
	for(i=N-1;i;i--)
	{
		while(ccnt>tmp&&cross(ch[ccnt]-ch[ccnt-1],pt[i]-ch[ccnt-1])<=0)
			ccnt--;
		ch[++ccnt]=pt[i];
	}
	if(N>=2) ccnt--;
}

int main()
{
	int i,x,y;
	scanf("%d%d",&N,&L);
	for(i=1;i<=N;i++) scanf("%d%d",&x,&y),pt[i]=P(x,y);
	calCH();
	double ans=0;
	for(i=1;i<=ccnt;i++) ans+=dis(ch[i],ch[i+1]);
	ans+=2*PI*L;
	if(ans-std::floor(ans)<0.5) ans=std::floor(ans);
	else ans=std::ceil(ans);
	printf("%.0f\n",ans);
	return 0;
}

###[BZOJ1027][JSOI2007]合金

给出 $m$ 种原材料合金，要铸造 $n$ 种需要的合金，求最少需要多少种原材料。

$n,m\leq500$

其实这道题跟凸包算法没关系。

个人认为，这类将一类实际问题转化为抽象问题解决的才是值得做的好题批判一下 ZJ 树的统计那种题。本题其实可以用几何方法做，首先，我们确定：由于第三维可以由前两维确定，所以用平面几何方就可以做。

之后，对于两个原料来说，可以合成的合金就在起线段上；进一步推广：多个原料可以合成的合金就在其凸包总点集中，所以问题就转化为了给出 $n$ 个点，求 $m$ 个点中最少要去多少个围成的多边形能圈住这 $n$ 个点。

然后呢，我们又把几何题转化为了图论题：直接用 Floyd 求最小环。

具体做法是：如果目标点都在线段 $(u,v)$ 的一侧，那我们令 $dis[u][v]=1$ ，否则 $dis[u][v]=\infty$ ，注意初始时 $dis[i][i]=\infty$ ，最后答案就是 $Min{dis[i][i]}$ （由于我们的起点一定要选在环里，而我们不能肯定那个点一定要选，所以要对所有点都跑最小环然后取最小值）

注意要特判所有点都重合或共线的情况。

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
#include<iostream>
#include<cstdio>
#include<cmath>

const int MAXN=505,INF=1e8;
const double eps=1e-10;

int n,m;

struct P
{
	double x,y;
	P(double x=0,double y=0):x(x),y(y){};
}A[MAXN],B[MAXN];

P operator +(P a,P b){return P(a.x+b.x,a.y+b.y);}
P operator -(P a,P b){return P(a.x-b.x,a.y-b.y);}
double cross(P a,P b){return a.x*b.y-a.y*b.x;}

int dis[MAXN][MAXN];

bool col(P a,P b)
{
	int i;
	if(a.x>b.x) std::swap(a,b);
	for(i=1;i<=m;i++)
		if(B[i].x<a.x||B[i].x>b.x) return false;
	if(a.y>b.y) std::swap(a,b);
	for(i=1;i<=m;i++)
		if(B[i].y<a.y||B[i].y>b.y) return false;
	return true;
}

int jud(P a,P b)
{
	int c1=0,c2=0;
	for(int i=1;i<=m;i++)
	{
		double t=cross(b-a,B[i]-a);
		if(t>eps) c1++;
		if(t<-eps) c2++;
		if(c1*c2) return 0;
	}
	if(!c1&&!c2&&col(a,b)){printf("2\n");return -1;}
	if(c1) return 1;
	if(c2) return 2;
	return 3;
}

void floyd()
{
	int ans=INF;
	int i,j,k;
	for(k=1;k<=n;k++)
		for(i=1;i<=n;i++)
			if(dis[i][k]<INF)
				for(j=1;j<=n;j++)
					dis[i][j]=std::min(dis[i][j],dis[i][k]+dis[k][j]);
	for(i=1;i<=n;i++) ans=std::min(ans,dis[i][i]);
	if(ans==INF||ans<=2) printf("-1\n");
	else printf("%d",ans);
}

void solve()
{
	int i,j;
	for(i=1;i<=n;i++)
		for(j=i+1;j<=n;j++)
		{
			int flag=jud(A[i],A[j]);
			if(flag==-1) return;
			if(flag==1) dis[i][j]=1;
			else if(flag==2) dis[j][i]=1;
			else if(flag==3) dis[i][j]=dis[j][i]=1;
		}
	floyd();
}

bool spj()
{
	for(int i=1;i<=m;i++)
		if(fabs(B[i].x-A[1].x)>eps||fabs(B[i].y-A[1].y)>eps)
			return 0;
	printf("1\n");
	return 1;
}

int main()
{
	int i,j;
	scanf("%d%d",&n,&m);
	for(i=1;i<=n;i++) for(j=1;j<=n;j++) dis[i][j]=INF;
	double tmp;
	for(i=1;i<=n;i++) scanf("%lf%lf%lf",&A[i].x,&A[i].y,&tmp);
	for(i=1;i<=m;i++) scanf("%lf%lf%lf",&B[i].x,&B[i].y,&tmp);
	if(spj()) return 0;
	else solve();return 0;
}

###[BZOJ4570][SCOI2016]妖怪

有 $n$ 只妖怪，每只妖怪有 $atk,dnf$ 两种属性，环境对妖怪有影响，在环境 $(a,b)$ 中，一直妖怪可以选择降低自己 $ka$ 点攻击力，提升 $kb$ 点防御力或者，提升自己 $ka$ 点攻击力，降低 $kb$ 点防御力。定义妖怪在某一环境中的战斗力为在此种环境下所能达到的最大 $atk$ 值加上最大 $dfn$ 值。求在最不利的环境下， $n$ 只妖怪中的战斗力最大值。

$n\leq10^6,atk_,dfn\leq10^8$

##半平面交
##旋转卡壳