计算几何

发表于 | 分类于

计算几何

解决几何问题的有力工具

OI中的几何知识考查通常是需要大规模运算的问题。如果对于每道题都去现成的推公式，必将会浪费大量的时间。对此，可以解决一切几何问题的解析几何便十分常用。而将解析几何模式化的工具就是向量。

向量的基本运算

向量的表示：$\vec a $ 向量a

向量的加减运算:满足平行四边形法则

$\vec {AB}+\vec {BC}=\vec{AC}$

$(a,b)+(c,d)=(a+c,b+c)$

$\vec{AB}-\vec{AC}=\vec{CB}$

$(a,b)-(c,d)=(a-b,c-d)$

向量的点乘、数量积

$\vec x · \vec y=|\vec x||\vec y|cos<\vec x,\vec y>$

$(a,b)·(c,d)=ac+bd$

几何意义为x向量在y向量上投影的与y向量模长的乘积

Attention:这是一个数量

向量的叉积

$\vec x ×\vec y=|\vec x||\vec y|sin<\vec x,\vec y>$

$(a,b)×(c,d)=ad-bc$

Attention:这是一个向量

同时，两个点相减是向量，相加无意义

两向量垂直，则点积为$0$，共线则叉积为$\vec 0$

code

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
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.y-a.y*b.x;	
   }
   double dot(P a,P b)
   {
       return a.x*b.x+a.y*b.y;
   }
   bool ver(P a,P b)//判断是否垂直 
   {
       return dot(a,b)==0; 
   }
   bool col(P a,P b)//判断是否共线
   {
       return cross(a,b)==0;
   }
   double dis(P a,P b)
   {
       return sqrt(dot(b-a,b-a));
   }

注意精度

凸包

定义：凸包是可以包裹平面上某个点集的最小凸多边形。

性质：凸包的顶点一定在点集内。

算法

Granham算法的变种：Andrew（更快更好记）

思路：先对点按x，y排序，然后进行从左到右扫描先求出下凸包，再求出上凸包。

Code：cogs896圈奶牛

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

题目有待补充

半平面交

定义：求几个半平面相交的部分。

算法：排序增量算法

表示：规定一条向量的左边是这条向量对应的半平面

According to：经典论文：

求n个半平面的交有三种做法： 第一种就是用每个平面去切割已有的凸多边形，复杂度O(n^2)。

第二种就是传说中的分治算法。将n个半平面分成两个部分，分别求完交之后再将两个相交的区域求交集。由于交出来的都是凸多边形，利用凸多边形的交可以在O(n)时间内完成的性质，将复杂度降为O(nlogn)。

第三种就是ZZY大牛的那篇论文提到的他自创的排序增量算法。

但是他的那种做法还是有些复杂，在网上找到evalls写的一个很优美的版本：

将所有半平面按极角排序，对于极角相同的，选择性的保留一个。$O(nlogn)$

step2. 使用一个双端队列(deque)，加入最开始2个半平面。

step3. 每次考虑一个新的半平面：

a.while deque顶端的两个半平面的交点在当前半平面外:删除deque顶端的半平面

b.while deque底部的两个半平面的交点在当前半平面外:删除deque底部的半平面

c.将新半平面加入deque顶端step4.删除两端多余的半平面。

具体方法是：

a.while deque顶端的两个半平面的交点在底部半平面外:删除deque顶端的半平面

b.while deque底部的两个半平面的交点在顶端半平面外:删除deque底部的半平面重复a,b直到不能删除为止。

step5:计算出deque顶端和底部的交点即可。 这个算法描述的非常清晰。当初写的时候有两个地方想的不太明白：step 1如何选择性的保留一个。step3如何判断交点在半平面外。

其实这两个问题都可以用叉积来解决。首先根据给定的两点顺序规定好极角序。

假定两点$o1,o2$的输入方向是顺时针，那么另一点P是否在其平面内只要判断o1P这个向量是否在$o1,o2$这个向量的右手边即可。

对于相同角度的两个半平面$(a_1a_2，b_1b_2)$，可以看$a_1b_1$这个向量是否在$a_1 a_2$这个向量的右手边，每次都要选择更靠近右手边的那个半平面。

直线的向量表示

$line_{AB}=Point_{A}+k×\vec {AB}$

这样表示的好处是，如果有两条直线有交点的话，那么直线$P+k_1\vec v,Q+k_2\vec w$
的交点在两条直线上的参数就是:
$$\vec u=P-Q$$
$$k_1=\frac {cross(w,u)} {cross(v,w)} $$
$$k_2=\frac {cross(u,v)} {cross(v,w)} $$

推导略。

例题：
Cqoi2006凸多边形
求几个凸多边形的交集：裸的半平面交（代码中半平面交部分见下题）

ZJOI2008瞭望塔

题面较长，简单描述一下：大概是要求一些点，这些点之间的连线下面都是山，然后要在一个点上选一个最低高度使得这个点可以看到所有点。

题解：对于相邻的两个点构造向量，求半平面交。
半平面交求得的区域就是可以建瞭望塔的区域，然后就是求得到的区域的边界与下方折线的边界的最近距离。
这个最短距离只可能出现在上下边界的顶点上，所以计算出来取最小值即可。

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
#include<iostream>  
#include<cstdio>  
#include<cmath>  
#include<cstdlib>  
#include<cstring>  
#include<algorithm>  
#define F(i,j,n) for(int i=j;i<=n;i++)  
#define D(i,j,n) for(int i=j;i>=n;i--)  
#define ll long long  
#define maxn 1005  
#define eps 1e-8  
#define inf 1000000000  
using namespace std;  
int n,cnt,tot;  
double ans=1e60;  
struct P{double x,y;}p[maxn],a[maxn];  
struct L{P a,b;double angle;}l[maxn],q[maxn];  
inline int read()  
{  
    int x=0,f=1;char ch=getchar();  
    while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}  
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}  
    return x*f;  
}  
inline P operator -(P a,P b){return (P){a.x-b.x,a.y-b.y};}  
inline double operator *(P a,P b){return a.x*b.y-a.y*b.x;}  
inline bool operator <(L a,L b)  
{  
    if (fabs(a.angle-b.angle)<eps) return (a.b-a.a)*(b.b-a.a)>0;  
    else return a.angle<b.angle;  
}  
inline P inter(L l1,L l2)  
{  
    double k1=(l2.b-l1.a)*(l1.b-l1.a),k2=(l1.b-l1.a)*(l2.a-l1.a),t=k1/(k1+k2);  
    return (P){l2.b.x+(l2.a.x-l2.b.x)*t,l2.b.y+(l2.a.y-l2.b.y)*t};  
}  
inline bool judge(L l1,L l2,L t)  
{  
    P p=inter(l1,l2);  
    return (t.b-t.a)*(p-t.a)<0;  
}  
inline void pre()  
{  
    p[0].x=p[1].x;p[0].y=inf;  
    p[n+1].x=p[n].x;p[n+1].y=inf;  
    F(i,0,n) l[++cnt]=(L){p[i],p[i+1]};  
    F(i,1,cnt) l[i].angle=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x);  
    sort(l+1,l+cnt+1);  
}  
inline void hpi()  
{  
    int top=0,bot=1;  
    tot=1;  
    F(i,2,cnt)  
    {  
        if (fabs(l[i].angle-l[i-1].angle)>=eps) tot++;  
        l[tot]=l[i];  
    }  
    cnt=tot;  
    q[++top]=l[1];q[++top]=l[2];  
    F(i,3,cnt)  
    {  
        while (bot<top&&judge(q[top-1],q[top],l[i])) top--;  
        while (bot<top&&judge(q[bot+1],q[bot],l[i])) bot++;  
        q[++top]=l[i];  
    }  
    while (bot<top&&judge(q[top-1],q[top],q[bot])) top--;  
    while (bot<top&&judge(q[bot+1],q[bot],q[top])) bot++;  
    tot=0;  
    F(i,bot,top-1) a[++tot]=inter(q[i],q[i+1]);  
}  
inline void getans()  
{  
    F(k,1,tot) F(i,1,n-1)  
    {  
        P t=(P){a[k].x,-1};  
        if (a[k].x>=p[i].x&&a[k].x<=p[i+1].x)  
            ans=min(ans,a[k].y-inter((L){p[i],p[i+1]},(L){t,a[k]}).y);  
    }  
    F(k,1,n) F(i,1,tot-1)  
    {  
        P t=(P){p[k].x,-1};  
        if (p[k].x>=a[i].x&&p[k].x<=a[i+1].x)  
            ans=min(ans,inter((L){a[i],a[i+1]},(L){t,p[k]}).y-p[k].y);  
    }  
}  
int main()  
{  
    n=read();  
    F(i,1,n) p[i].x=read();  
    F(i,1,n) p[i].y=read();  
    pre();  
    hpi();  
    getans();  
    printf("%.3lf\n",ans);  
    return 0;  
}

旋转卡壳

思路就是固定一条边，然后逆时针旋转着计算答案。因为边逆时针枚举，所以贡献答案的点也是逆时针旋转的，只要对于需要计算的凸多边形来操作就行了。

POJ2185

给出n个点，求平面内的最远点对。

题解：
首先答案肯定在凸包上。先求出凸包之后，只需要求凸包的直径，这时采用旋转卡壳来求。

code

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
#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std; 
const int N=1e5+1e3+7;
const double eps=1e-12;
struct P{
    double x,y;
}p[N],ch[N];
int n,m;
double ans;
P operator -(P a,P b)
{
    P ret;
    ret.x=a.x-b.x;
    ret.y=a.y-b.y;
    return ret;
}
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 dot(b-a,b-a);
}
bool cmp(P a,P b)
{
    return (a.x<b.x)||(a.x==b.x&&a.y<=b.y);
}
void pre()
{
    sort(p+1,p+n+1,cmp);
    for(int i=1;i<=n;i++)
    {
        while(m>1&&cross(ch[m]-ch[m-1],p[i]-ch[m-1])<=0)
            m--;
        ch[++m]=p[i];
    }
    int k=m;
    for(int i=n-1;i>=1;i--)
    {
        while(m>k&&cross(ch[m]-ch[m-1],p[i]-ch[m-1])<=0)
            m--;
        ch[++m]=p[i];
    }
    if(n>1)
        m--;
}
void GR()
{
    ch[m+1]=ch[1];
    int now=2;
    for(int i=1;i<=m;i++)
    {
        while(cross((ch[i+1]-ch[i]),(ch[now]-ch[i]))<cross((ch[i+1]-ch[i]),(ch[now+1]-ch[i])))
        {
            now++;
            if(now==m+1)
                now=1;
        }
        ans=max(ans,dis(ch[now],ch[i]));
    }
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        scanf("%lf%lf",&p[i].x,&p[i].y);
    pre();
    GR();
    printf("%d",(int)(ans));
}

例题见后续博客