ACM算法总结平面几何

最新推荐文章于 2022-07-13 10:22:19 发布

瓴龍

最新推荐文章于 2022-07-13 10:22:19 发布

阅读量580

点赞数 2

分类专栏： ACM算法总结文章标签：几何

本文链接：https://blog.csdn.net/dragonylee/article/details/103103761

版权

ACM算法总结专栏收录该内容

23 篇文章 6 订阅

订阅专栏

基础知识

像一些点、向量、点乘、叉乘、圆等等知识，应该都是最基本的平面几何知识，这里就不再赘述，主要考虑如何实现这些平面几何的运算。

首先是定义一些常量和宏定义：

#define RAD(x) (pi*x/180)
#define DEG(x) (x*180/pi)
const double eps=1e-8;
const double pi=acos(-1);
int dcmp(double x) {return fabs(x)<eps?0:(x<0?-1:1);}

其中dcmp是一个三态函数，用于判断一个浮点数与0的关系，排除了精度误差。

然后是点和向量的定义（注意虽然在几何中点和向量是本质不同的东西，但是它们的形式都是一样的）：

typedef struct point
{
    double x,y;
    point(double x=0,double y=0):x(x),y(y) {}
    friend ostream & operator << (ostream &,point a) {cout<<fixed<<setprecision(3)<<a.x<<' '<<a.y; return cout;}

    point operator + (point b) {return point(x+b.x,y+b.y);}
    point operator - (point b) {return point(x-b.x,y-b.y);}
    point operator * (double p) {return point(x*p,y*p);}
    point operator / (double p) {return point(x/p,y/p);}
    double operator * (point b) {return x*b.x+y*b.y;}       //点乘、内积
    double operator ^ (point b) {return x*b.y-y*b.x;}       //叉乘、外积
    bool operator < (const point &b) const {return x<b.x || (x==b.x && y<b.y);}
}vec;

其中叉乘的意义需要特别注意一下，如果向量 $\times b$ 是一个正数，那么由右手螺旋定则可以知道， $b$ 在 $a$ 的逆时针方向；里面的友元函数是为了方便格式化输出，调试用。

接下来是圆的定义：

struct circle
{
    point c; double r;
    circle(point c,double r):c(c),r(r) {}
    friend ostream & operator << (ostream &,circle a) {cout<<fixed<<setprecision(3)<<a.c<<'\t'<<a.r;return cout;}
};

接下来是直线的定义：

struct line
{
    point p; vec v; double ang;
    line(point p=point(0,0),vec v=vec(0,1)):p(p),v(v) {ang=atan2(v.y,v.x);}
    bool operator < (const line &l) const {return ang<l.ang;}
};

注意atan2的取值范围是0到 $\pi$ ，对于直线来说刚好（射线就不行）。

接下来是多边形和**直线集（或者称为多面集）**的定义：

typedef vector<point> polygon;
typedef vector<line> lineset;

注意这里的多边形顶点一定要按照某一时针方向储存！（不够智能？不然鬼知道怎么给这些点排序）

基础运算

向量长度：向量的二范数，也就是这个向量自己点乘自己然后开方。

double length(vec a) {return sqrt(a*a);}

向量夹角：两个向量的夹角，范围为 $[0,\pi]$ ，可以用余弦定理求得。

double angle(vec a,vec b) {return acos(a*b/length(a)/length(b));}

三角形面积：这个分为两种表示形式，一种是两个向量的形式，一种是三个顶点的形式。

double tri_area(vec a,vec b) {return fabs(a^b)/2;}
double tri_area(point a,point b,point c) {return fabs((b-a)^(c-a))/2;}

向量旋转：将一个向量逆时针旋转某一度数，我们知道旋转可以通过左乘一个矩阵来实现。

vec vec_rotate(vec a,double x) {return vec(a.x*cos(x)-a.y*sin(x),a.x*sin(x)+a.y*cos(x));}

两直线交点：两条直线的交点，要保证必须有交点，即两直线不平行。可以通过参数方程求解。

point line_intersect(line a,line b)
{
    double t=((b.p-a.p)^b.v)/(a.v^b.v);
    return a.p+(a.v*t);
}

点到直线距离：一个点到一条直线的距离。用三角形面积公式求解。

double dis_to_line(point a,line b)
{
    point c=b.p+b.v-a;
    return fabs((c^b.v)/length(b.v));
}

点在直线的投影：点在直线上的投影点。参数方程结合内积求解。

point line_projection(point a,line b)
{
    double t=((a*b.v)-(b.p*b.v))/(b.v*b.v);
    return b.p+b.v*t;
}

点关于直线的对称：一个点关于一条直线的对称点。算出投影后延长即可。

point line_symmetry(point a,line b)
{
    point c=line_projection(a,b);
    return a+((c-a)*2);
}

多边形周长：一个多边形的周长，计算每边边长求和即可。

double polygon_perimeter(polygon a)
{
    double ans=0;
    REP(i,0,a.size()-2) ans+=length(a[i]-a[i+1]);
    return ans+length(a[a.size()-1]-a[0]);
}

多边形面积：一个多边形的面积，拆分成三角形求解。（下面代码对任意多边形（凸、凹）适用）

double polygon_area(polygon a)
{
    double ans=0;
    REP(i,1,a.size()-2) ans+=((a[i]-a[0])^(a[i+1]-a[0]));
    return fabs(ans/2);
}

三角形外接圆：三角形的外接圆，圆心坐标可以用两点之间距离公式推导出来。

circle circumscribed_circle(point p1,point p2,point p3)
{
    double bx=p2.x-p1.x,by=p2.y-p1.y;
    double cx=p3.x-p1.x,cy=p3.y-p1.y;
    double d=(bx*cy-by*cx)*2;
    double ax=(cy*(bx*bx+by*by)-by*(cx*cx+cy*cy))/d+p1.x;
    double ay=(bx*(cx*cx+cy*cy)-cx*(bx*bx+by*by))/d+p1.y;
    point p(ax,ay);
    return circle(p,length(p1-p));
}

三角形内切圆：三角形的内切圆。

circle inscribed_circle(point p1,point p2,point p3)
{
    double a=length(p2-p3);
    double b=length(p3-p1);
    double c=length(p1-p2);
    point p=(p1*a+p2*b+p3*c)/(a+b+c);
    return circle(p,dis_to_line(p,line(p1,p2-p1)));
}

凸包

在实向量空间中，对于给定的集合 $S$ ，所有包含 $S$ 的凸集的交集定义为这个集合的凸包。详细内容可以参看《凸优化》中凸集相关部分的描述。这里只给出几条性质（~~好像没什么用~~）：

凸包内任意两点的凸组合仍然在凸包内；
一个集合是凸的当且仅当它的凸包是它本身；

对于二维平面点集，求解凸包的常用方法是Graham扫描法（复杂度 $n\log n$ ）。具体做法是首先按照横坐标对点集排序，然后用栈维护一个下凸序列，左右扫两次就可以了。代码如下：

polygon convex_hull(polygon a)
{
    sort(a.begin(),a.end());
    int m=0,n=a.size();
    point *q=new point[n+5];
    REP(i,0,n-1)
    {
        while(m>1 && ((q[m-1]-q[m-2])^(a[i]-q[m-2]))<=0) m--;
        q[m++]=a[i];
    }
    int k=m;
    REP_(i,n-2,0)
    {
        while(m>k && ((q[m-1]-q[m-2])^(a[i]-q[m-2]))<=0) m--;
        q[m++]=a[i];
    }
    if(n>1) m--;
    polygon ans;
    REP(i,0,m-1) ans.push_back(q[i]);
    return ans;
}

旋转卡壳

这是一种思维，本质上是挖掘了凸多边形点边距的规律（先增后减，极值单方向变化）。

旋转卡壳有16种读音，但我认为最正确的一种应该是—— $旋(xu\acute{a}n)转(zhu\check{a}n)卡(qi\check{a})壳(k\acute{e})$ 。

旋转卡壳可以用来解决类似最远点对、最远点边距、最小面积矩形覆盖等等（先取凸包），下图是旋转卡壳求最远点边距的示意图：

当单方向枚举边时，对应最优解的点的变化也是单方向的，并且枚举完所有边之后，最优解的点也刚好旋转了一圈。可以看出这样的算法复杂度是O(n)的。

其它的算法类似：最远点对旋转时可以用下一个点的变化方向来判断该谁旋转。比如说当前处理到p[i]和p[j]，其中p是凸包，那么就根据(p[i+1]-p[i])和(p[j+1]-p[j])的夹角来判断就行了（把图画出来就很好理解）。

求点集的最远点对代码如下：

double max_point_dis(polygon a)
{
    polygon p=convex_hull(a);
    int n=p.size(),i=0,j=0,si,sj;
    if(n==2) return length(p[0]-p[1]);
    REP(k,0,n-1)
    {
        if(!(p[i]<p[k])) i=k;
        if(p[j]<p[k]) j=k;
    }
    double ans=0; si=i,sj=j;
    while(i!=sj || j!=si)
    {
        ans=max(ans,length(p[i]-p[j]));
        (((p[(i+1)%n]-p[i])^(p[(j+1)%n]-p[j]))<0)?(i=(i+1)%n):(j=(j+1)%n);
    }
    return ans;
}

对于最小面积或周长矩形覆盖问题其实是一样的，容易想到（其实就是感觉）最优解一定会有一条边和凸包的某条边重合，于是只要枚举这些边，其实整个矩形就已经定下来了。这里随着边的变化有三个点跟着同方向移动，枚举边求相应最小值即可。

半平面交

给定多个半平面（即直线的某一侧无限大区域），求所有半平面的交。基本思路是将所有直线按照极角排序，然后维护一个双端队列（因为同一个极角对应了两个直线方向，所以双端都要维护），加入新的半平面时通过比较这条直线和上一个交点的关系去维护队列。

我们规定直线的左侧是合法的半平面，最后输出所有半平面交出来的凸多边形（一定是凸的，因为半平面本身就是一个凸集，而有限个凸集的交一定是凸集），代码如下：

bool on_left(point p,line L) {return (L.v^(p-L.p))>0;}
polygon half_plane_intersect(lineset L)
{
    sort(L.begin(),L.end());
    int l=0,r=0,n=L.size();
    point *p=new point[n];
    line *q=new line[n];
    q[0]=L[0];
    REP(i,1,n-1)
    {
        while(l<r && !on_left(p[r-1],L[i])) r--;
        while(l<r && !on_left(p[l],L[i])) l++;
        q[++r]=L[i];
        if(fabs(q[r].v^q[r-1].v)<eps)
        {
            r--;
            if(on_left(L[i].p,q[r])) q[r]=L[i];
        }
        if(l<r) p[r-1]=line_intersect(q[r-1],q[r]);
    }
    while(l<r && !on_left(p[r-1],q[l])) r--;
    if(r-l<=1) return (polygon)0;
    p[r]=line_intersect(q[r],q[l]);
    polygon ans;
    REP(i,l,r) ans.push_back(p[i]);
    return ans;
}

如果题目不保证最后交出来的一定是一个有限的凸多边形，那么我们可以添加一个坐标很大的矩形框，最后再把它删除就可以了。

最小圆覆盖

找到一个最小半径的圆，覆盖平面上所有给定的点。

做法是使用随机增量法，大致思路是在添加一个新的点时更新当前圆的圆心和半径，如果新增的点在当前圆内，则什么也不做；如果新增的点在当前圆外，则以新增的点为圆心，双重枚举之前的所有点，找到一个最小的圆。由于最小圆覆盖的圆一定经过某两个点或者某三个点，所以算法正确性得到保证，另外可以证明在随机化的条件下，时间复杂度的期望是O(n) 。

代码：

circle min_circle(polygon p)
{
    random_shuffle(p.begin(),p.end());
    int n=p.size();
    point c=p[0];
    double r=0;
    REP(i,0,n-1)
    {
        if(dcmp(length(c-p[i])-r)<=0) continue;
        c=p[i],r=0;
        REP(j,0,i-1)
        {
            if(dcmp(length(c-p[j])-r)<=0) continue;
            c=(p[i]+p[j])/2; r=length(c-p[i]);
            REP(k,0,j-1)
            {
                if(dcmp(length(c-p[k])-r)<=0) continue;
                circle x=circumscribed_circle(p[i],p[j],p[k]);
                c=x.c,r=x.r;
            }
        }
    }
    return circle(c,r);
}

扫描线

首先介绍离散化的思想，一般来说有些数值太大无法直接使用数组下标表示，我们就从小到大给每个数分配一个id，然后用这个id表示这个数。每次可以使用二分查找在logn的时间内查询某个数的id。代码实现如下：

sort(a+1,a+1+n);
len=unique(a+1,a+n+1)-a-1;

查询id：

int ID(int x) {return lower_bound(a+1,a+len+1,x)-a;}

扫描线，顾名思义就是一条线在图形上扫描，但是肯定是按照一定的顺序和关键点去扫描的，比较经典的应用是计算平面内多个矩形的覆盖面积。大致思路如下：

将所有横坐标离散化，纵坐标当做扫描线的关键点（也即每次递增扫下一个纵坐标），用线段树维护相邻两个横坐标之间的线段出现的次数，每经过一个矩形的下边沿，相应所有线段出现次数加一，反之经过一个上边沿，出现次数减一；对于相邻两个纵坐标统计有效边沿长度和，累计答案。

这里的线段树的使用有一些技巧，因为我们需要的是某一个区间中所有出现次数大于0的结点的长度和，所以lazy标签可以不用pushdown，然后如果当前区间的lazy大于0，其长度和可以直接计算；否则，其长度和为其左右儿子的长度和之和。

下面是一道经典例题代码：（HDU1542 Atlantis）

#include <bits/stdc++.h>
#define mem(a,b) memset(a,b,sizeof(a))
#define REP(i,a,b) for(int i=(a);i<=(int)(b);i++)
#define REP_(i,a,b) for(int i=(a);i>=(b);i--)
using namespace std;
typedef long long LL;

const int maxn=1005;
int lazy[maxn],n,T,len;
double id[maxn],s[maxn];
struct edge
{
    double l,r,y;
    int flag;
    bool operator < (const edge &a) const {return y<a.y;}
}e[maxn];
int ID(double x) {return lower_bound(id+1,id+len+1,x)-id;}

#define chl k<<1
#define chr k<<1|1
#define mid ((l+r)>>1)

void push_up(int k,int l,int r)
{
    if(lazy[k]>0) s[k]=id[r+1]-id[l];
    else s[k]=s[chl]+s[chr];
}

void update(int k,int l,int r,int ll,int rr,int flag)
{
    if(r<ll || l>rr) return;
    if(l>=ll && r<=rr)
    {
        lazy[k]+=flag;
        push_up(k,l,r);
        return;
    }
    update(chl,l,mid,ll,rr,flag);
    update(chr,mid+1,r,ll,rr,flag);
    push_up(k,l,r);
}

int main()
{
    //freopen("input.txt","r",stdin);
    while(~scanf("%d",&n))
    {
        if(!n) break;
        REP(i,1,n)
        {
            double x,y,xx,yy;
            scanf("%lf%lf%lf%lf",&x,&y,&xx,&yy);
            id[2*i-1]=x; id[2*i]=xx;
            e[2*i-1]=(edge){x,xx,yy,-1};
            e[2*i]=(edge){x,xx,y,1};
        }
        n*=2;
        sort(e+1,e+n+1);
        sort(id+1,id+n+1);
        mem(s,0); mem(lazy,0);
        len=unique(id+1,id+n+1)-id-1;

        update(1,1,n-1,ID(e[1].l),ID(e[1].r)-1,e[1].flag);
        double ans=0;
        REP(i,2,n)
        {
            ans+=(e[i].y-e[i-1].y)*s[1];
            update(1,1,n-1,ID(e[i].l),ID(e[i].r)-1,e[i].flag);
        }
        printf("Test case #%d\nTotal explored area: %.2f\n\n",++T,ans);
    }

    return 0;
}