@dxbdly 2021-08-29T05:38:38.000000Z 字数 9076 阅读 235

G2020 19寒假集训Day6

信息学——19寒假集训

G2020 19寒假集训Day6

今天是G2020寒假集训Day6，今天的主题是~~毒瘤的~~ 计算几何

数学不是很好，有写得不对或写得不是很清楚的请见谅。

$Update\,2020.10.25$ ：有些东西的确写得不准确啊……凑合着看吧，咕了。

一. 向量

$1.$ 定义：

向量的定义:

向量表示的是一个有大小和方向的量，
在平面坐标系下它与点一样也用两个数来表示。
这两个数的实际含义是将这个向量的起点移至原点后终点的坐标。

我们用 $\vec{v}$ 来表示一个向量。其实我们可以将向量看成点的另一种形式。

$2.$ 向量的基本运算

关于向量的四则运算，还是比较简单的

直接放代码啦~~~

//The code is from dxbdly
struct Vec{
    double x,y;
    Vec(double X=0,double Y=0):x(X),y(Y){}
    Vec operator+(const Vec &k)const {return Vec(x+k.x,y+k.y);}
    Vec operator-(const Vec &k)const {return Vec(x-k.x,y-k.y);}
    Vec operator*(const int &k)const {return Vec(x*k,y*k);}
    Vec operator/(const int &k)const {return Vec(x/k,y/k);}
}

应该能看懂吧，下面看几个比较特殊的。

$3.$ 向量的点积

设向量 $\vec{u}(x_1,y_1)$ ，向量 $\vec{v}(x_2,y_2)$

则两向量的点积：

$\vec{u}*\vec{v}=x_1*x_2+y_1*y_2$

点积的几何意义：

若两向量 $\vec{u}$ 与 $\vec{v}$ 的夹角为 $\Theta$ ，他们的点积为：

$\vec{u}*\vec{v}=|\vec{u}||\vec{v}|\cos\Theta$

使用点积判断两向量夹角情况：

且 两 向 量 垂 直 点 积 为 正 点 积 为 点 积 为 负

$\begin{cases} \Theta<90^o &\text{点积为正} \\ \text{$\Theta=90^o$且两向量垂直} &\text{点积为0} \\ \Theta>90^o & \text{点积为负} \end{cases}$

（PS：点积满足交换律）

代码：

//The code is from dxbdly
inline double dj(Vec a,Vec b)//点积
{
    return a.x*b.x+a.y*b.y;
}

$4.$ 向量的叉积

叉积与点积比较相似。

设向量 $\vec{u}(x_1,y_1)$ ，向量 $\vec{v}(x_2,y_2)$

两向量的叉积：

$\vec{u}*\vec{v}=x_1*y_2-x_2*y_1$

叉积的几何意义：

由 与 构 成 的 平 行 四 边 形 面 积

$\vec{u}*\vec{v}=\text{由$\vec{u}$与$\vec{v}$构成的平行四边形面积}$

所以我们容易求出两向量构成的三角形面积。

$S_\Delta=\vec{u}*\vec{v}$

使用叉积判断两向量位置关系：

在 左 边 与 方 向 相 同 在 右 边

$\begin{cases} \text{$\vec{u}$在$\vec{v}$左边} & \vec{v}*\vec{u}>0\\ \text{$\vec{u}$与$\vec{v}$方向相同} & \vec{v}*\vec{u}=0\\ \text{$\vec{u}$在$\vec{v}$右边} & \vec{v}*\vec{u}<0 \end{cases}$

补充：
判断线与向量的位置关系只需在直线上任取一向量
看两向量之间的位置关系即可。

（PS：注意叉积不满足交换律！！！）

代码：

//The code is from dxbdly
inline double cj(Vec a,Vec b)
{
    return a.x*b.y-a.y*b.x;
}

$5.$ 向量的极角

对于一个向量 $\vec{u}(x,y)$ 我们可以通过 $atan2$ 函数( $c++$ 内有)得到他的极角 $(\vec{u}.at)$ 。

$\vec{u}.at=atan2(y,x)$

$6.$ 向量的旋转

设向量 $\vec{u}(x,y)$ 旋转了 $\Theta^o$

旋转公式：

$\vec{u'}=Vec(x\cos\Theta-y\cos\Theta,x \sin \Theta +y\sin\Theta)$

矩阵形式:

$\begin{bmatrix} x' \\ y' \\ \end{bmatrix}= \begin{bmatrix} \cos\Theta&-\sin\Theta \\ \sin\Theta&\cos\Theta \\ \end{bmatrix} \begin{bmatrix} x \\ y \\ \end{bmatrix}$

二.凸包

1.凸包定义与性质：

定义：

凸包：平面凸包是指覆盖平面上n个点的最小的凸多边形。
上凸壳：凸包的上半部分
下凸壳：凸包的下半部分。

性质：

凸壳的凸性：对于一个凸壳，所有线段的斜率满足单调性。

2.Graham算法

题面

算法本质：

Graham扫描算法维护一个凸壳
通过不断在凸壳中加入新的点和去除影响凸性的点 
最后形成凸包

算法流程：

1. 排序
2. 扫描

我们先将所有向量按极角序排序。

然后维护一个存点的栈要求栈内斜率单调(满足凸性)。

我们每次将一个点加入栈中，判断是否会与前面冲突。

如果有冲突将前面的点删除即可，最后留下的点便是凸包。

判断的方法么:

//The code is from dxbdly
inline bool check(Vec a1,Vec a2,Vec b1,Vec b2)
{
    return (a2-a1)*(b2-b1)<0.0;
}

总代码：(不是我写的……)

//The code is not from dxbdly
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
int n,top;
double ans;
struct data{
   double x,y;
}p[10001],s[10001];
double dis(data a,data b)
{return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
double mul(data p1,data p2,data p0)
{return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);}
inline bool cmp(data a,data b)
{
    if(mul(a,b,p[0])==0)return dis(a,p[0])<dis(b,p[0]);
    return mul(a,b,p[0])>0;
}
void graham()
{
    top=2;data t;
    int k=0;
    for(int i=1;i<n;i++)
        if((p[k].y>p[i].y)||(p[k].y==p[i].y&&p[k].x>p[i].x))k=i;
    t=p[0];p[0]=p[k];p[k]=t;
    sort(p+1,p+n,cmp);
    s[0]=p[0],s[1]=p[1],s[2]=p[2];
    for(int i=3;i<n;i++)
    {
        while(top&&mul(p[i],s[top],s[top-1])>=0)
           top--;
        s[++top]=p[i];
    }
    s[++top]=p[0];
    for(int i=0;i<top;i++)
        ans+=dis(s[i],s[i+1]);
}
int main()
{
    //freopen("fc.in","r",stdin);
    //freopen("fc.out","w",stdout);
    scanf("%d",&n);
    for(int i=0;i<n;i++)
        scanf("%lf%lf",&p[i].x,&p[i].y);
    graham();
    printf("%.2lf",ans);
    return 0;
}

3.Andrew算法

算法思路：

1.按x坐标为第一优先级，按y坐标为第二优先级排序。
2.分别维护上下凸壳，将两个凸壳合并形成凸包。

直接放代码：

//The code is from dxbdly
#include<bits/stdc++.h>
#define int long long
#define INF 1e12
using namespace std;
inline int read()
{
    int x=0;
    bool f=0;
    char c=getchar();
    while (!isdigit(c))
          f|=(c=='-'),c=getchar();
    while (isdigit(c))
          x=(x<<3)+(x<<1)+(c^48),c=getchar();
    return f?-x:x;
}
int n,top;
struct Vec{
    double x,y;
}point[100005];
Vec st[100005];
double ans; 
inline bool operator <(Vec a,Vec b)
{
    return a.x!=b.x?a.x<b.x:a.y<b.y; 
}
inline Vec operator -(Vec a,Vec b)
{
    return (Vec){b.x-a.x,b.y-a.y};
}
inline double operator *(Vec a,Vec b) 
{
    return a.x*b.y-a.y*b.x;
}
inline bool check(Vec a1,Vec a2,Vec b1,Vec b2)
{
    return (a1-a2)*(b1-b2)<0.0;
}
inline double get_dis(Vec a,Vec b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline void get_up()
{
    top=0,st[++top]=point[1],st[++top]=point[2];
    for (register int i=3;i<=n;++i)
        {
        while (top>1&&check(st[top-1],st[top],st[top],point[i]))
              top--;
        st[++top]=point[i];
        }
    for (register int i=1;i<top;++i)
        ans+=get_dis(st[i],st[i+1]);
}
inline void get_down()
{
    top=0,st[++top]=point[n],st[++top]=point[n-1];
    for (register int i=n-2;i>=1;--i)
        {
        while (top>1&&check(st[top-1],st[top],st[top],point[i]))
              top--;
        st[++top]=point[i];
        }
    for (register int i=1;i<top;++i)
        ans+=get_dis(st[i],st[i+1]);    
}
signed main(){
    n=read();
    for (register int i=1;i<=n;++i)
        scanf("%lf%lf",&point[i].x,&point[i].y);
    sort(point+1,point+n+1);
    get_up(),get_down();
    printf("%0.2lf",ans);
    return 0;
}

三.半平面交

题面

1.定义：

半平面：由一条向量将平面分为两半，其中一半为半平面。
半平面交：所有半平面的最大相交部分组成的多边形。

2.求法：

将线转换为两个向量表示，考虑与凸包相似的求法，将线段按向量极角序排序，维护一个栈，使其满足凸性。

就可以了吗？？？

答案是肯定不行的，比如下面这种情况:

你会发现，加入新的线(蓝线)后与之前的首尾线段都有冲突！

所以我们应该用单调队列维护

将半平面交的组成线段确定之后，我们需要求出这些线段的交点。

求交点的话直接用叉积搞一下即可，不多赘述（~~其实我不会~~)

代码：

//The code is from dxbdly
inline Vec get_inter(Line a,Line b)
{
    double k1=(b.b-a.a)*(a.b-a.a),k2=(a.b-a.a)*(b.a-a.a);
    double cnt=k1/(k1+k2);
    return (Vec){b.b.x+(b.a.x-b.b.x)*cnt,b.b.y+(b.a.y-b.b.y)*cnt};
}

我们还要会一件事：求任意 $n$ 边形面积

对于求任意n边形面积问题
我们采用分割法，将其分割成n-1个三角形。
对于每个三角形，我们利用叉积的几何意义求出其面积。
将n-1个三角形面积相加即可。

算法流程：

1.按极角序排序。
2.类Graham算法求出半平面交的边集合。
3.求出所有交点坐标。
3.用叉积求出面积。

其中还有一些细节，在这里列举一些

(其实计算几何就是细节多)

1.对于平行的情况我们只取最左侧的边。
2.注意边界的处理。

具体的话可以根据代码看。

//The code is from dxbdly
#include<bits/stdc++.h>
using namespace std;
inline int read()
{
    int x=0;
    char c=getchar();
    bool f=0;
    while (!isdigit(c))
          f=f|(c=='-'),c=getchar();
    while (isdigit(c))
          x=(x<<3)+(x<<1)+(c^48),c=getchar();
    return f?-x:x;
}
int n,m,tot,len;
struct Vec{
    double x,y;
}point[1505],ans[1505];
struct Line{
    Vec a,b;
    double at;
}line[1505],q[1505];
inline Vec operator -(Vec a,Vec b) 
{
    return (Vec){a.x-b.x,a.y-b.y};
}
inline double operator *(Vec a,Vec b)
{
    return (a.x*b.y)-(a.y*b.x);
}
inline bool operator <(Line a,Line b)
{
    return a.at!=b.at?a.at<b.at:(a.b-a.a)*(b.b-a.a)>0;
}
inline Vec get_inter(Line a,Line b)
{
    double k1=(b.b-a.a)*(a.b-a.a),k2=(a.b-a.a)*(b.a-a.a);
    double cnt=k1/(k1+k2);
    return (Vec){b.b.x+(b.a.x-b.b.x)*cnt,b.b.y+(b.a.y-b.b.y)*cnt};
}
inline bool check(Line a,Line b,Line c)
{
    Vec inter=get_inter(a,b);
    return (c.b-c.a)*(inter-c.a)<0;
}
inline void work()
{
    int head=1,tail=0,cnt=0;
    for (register int i=1;i<=tot;++i)
        {
        if (line[i].at!=line[i-1].at)
           cnt++;
        line[cnt]=line[i];
        }
    tot=cnt;
    q[++tail]=line[1],q[++tail]=line[2];
    for (register int i=3;i<=tot;++i)
        {
        while (head<tail&&check(q[tail-1],q[tail],line[i]))
              tail--;
        while (head<tail&&check(q[head+1],q[head],line[i]))
              head++;
        q[++tail]=line[i];
        }
    while (head<tail&&check(q[tail-1],q[tail],q[head]))
          tail--;
    while (head<tail&&check(q[head+1],q[head],q[tail]))
          head++;
    q[tail+1]=q[head];
    for (register int i=head;i<=tail;++i)   
        ans[++len]=get_inter(q[i],q[i+1]);
}
inline double get_ans()
{
    if (len<3)
       return 0;
    double sum=0;
    ans[++len]=ans[1];
    for (register int i=1;i<=len;++i)
        sum+=ans[i]*ans[i+1];
    return fabs(sum)/2;
}
int main(){
    n=read();
    for (register int i=1;i<=n;++i)
        {
        m=read();   
        for (register int j=1;j<=m;++j)
            scanf("%lf%lf",&point[j].x,&point[j].y);
        point[m+1]=point[1];
        for (register int j=1;j<=m;++j)
            ++tot,line[tot].a=point[j],line[tot].b=point[j+1];
        }
    for (register int i=1;i<=tot;++i)
        line[i].at=atan2(line[i].b.y-line[i].a.y,line[i].b.x-line[i].a.x);
    sort(line+1,line+tot+1);
    work();
    printf("%0.3f",get_ans());
    return 0;
}

4.旋转卡壳

~~这是一个有16种读音的算法~~

旋(xuan2)转(zhuan4)卡(qia3)壳(qiao4)

题面

是求凸包（或多边形）最长点对距离的方法。

算法思路不难，写起来真的卡壳……

算法原理：

过凸包上两点用两条平行线将凸包“卡”住
每次旋转两条平行线，使得始终有一条线贴着凸包的一条边。
设这样的两个点为对踵点，则答案一定在这些对踵点中。

至于实现：

构造出凸包。
枚举一条凸包上的边(模拟旋转)
找到其对应的对踵点。

那么如何寻找对应的对踵点尤为重要。

$1^o$ 暴力

对于每一条边枚举所有点，然后找到最优点

时间复杂度 $O(n^2)$

$2^o$ 旋转卡壳

观察发现。

我们可以发现凸包上的边与每个点的三角形面积成单峰函数

所以我们在逆时针枚举边的时候，对应的对踵点也会逆时针变化，所以我们每次从上一点开始即可。

时间复杂度 $O(n)$

总结一下算法流程：

1.构造出凸包。
2.枚举凸包上的边
3.从上一次的答案开始枚举对踵点
4.判断是否更新。

代码：（我没调出来……，所以不是我的）

//The code is not from dxbdly
#include<bits/stdc++.h>
#define eps 1e-8
#define inf 1000000000
#define pa pair<int,int>
#define ll long long 
using namespace std;
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;
}
int n,top;
double ans;
double sqr(double x)
{
    return x*x;
}
struct P{
    double x,y;
    P(){}
    P(double _x,double _y):x(_x),y(_y){}
    friend P operator +(P a,P b){
        return P(a.x+b.x,a.y+b.y);
    }
    friend P operator -(P a,P b){
        return P(a.x-b.x,a.y-b.y);
    }
    friend double operator*(P a,P b){
        return a.x*b.y-a.y*b.x;
    }
    friend double operator/(P a,P b){
        return a.x*b.x+a.y*b.y;
    }
    friend bool operator==(P a,P b){
        return fabs(a.x-b.x)<eps&&fabs(a.y-b.y)<eps;
    }
    friend bool operator!=(P a,P b){
        return !(a==b);
    }
    friend bool operator<(P a,P b){
        if(fabs(a.y-b.y)<eps)return a.x<b.x;
        return a.y<b.y;
    }
    friend double dis2(P a){
        return sqr(a.x)+sqr(a.y);
    }
    friend void print(P a){
        printf("%.2lf %.2lf\n",a.x,a.y);
    }
}p[50005],q[50005];
bool cmp(P a,P b)
{
    if(fabs((b-p[1])*(a-p[1]))<eps)return dis2(a-p[1])<dis2(b-p[1]);
    return (a-p[1])*(b-p[1])>0;
}
void graham()
{
    for(int i=1;i<=n;i++)
        if(p[i]<p[1])swap(p[i],p[1]);
    sort(p+2,p+n+1,cmp);
    q[++top]=p[1];q[++top]=p[2];
    for(int i=3;i<=n;i++)
    {
        while((q[top]-q[top-1])*(p[i]-q[top-1])<eps&&top>1)top--;
        q[++top]=p[i];
    }
}
void RC()
{
    q[top+1]=q[1];
    int now=2;
    for(int i=1;i<=top;i++)
    {
        while((q[i+1]-q[i])*(q[now]-q[i])<(q[i+1]-q[i])*(q[now+1]-q[i]))
        {
            now++;
            if(now==top+1)now=1;
        }
        ans=max(ans,dis2(q[now]-q[i]));
    }
}
int main()
{
    n=read();
    for(int i=1;i<=n;i++)
        p[i].x=read(),p[i].y=read();
    graham();
    RC();
    printf("%d",(int)ans);
    return 0;
}

这里面有些代码没调出来，所以不是我写的，到时候调出来了再更新吧，先凑合一下。

G2020 19寒假集训Day6

一. 向量

1.定义：

2. 向量的基本运算

3. 向量的点积

4. 向量的叉积

5.向量的极角

6.向量的旋转

二.凸包

1.凸包定义与性质：

2.Graham算法

3.Andrew算法

三.半平面交

1.定义：

2.求法：

4.旋转卡壳

1^o 暴力

2^o 旋转卡壳

内容目录

选择主题

$1.$ 定义：

$2.$ 向量的基本运算

$3.$ 向量的点积

$4.$ 向量的叉积

$5.$ 向量的极角

$6.$ 向量的旋转

$1^o$ 暴力

$2^o$ 旋转卡壳