计算几何 学习笔记

向量、直线与多边形

精度控制

控制精度一般在10−7~10−10之间

const double eps=1e-7;
int dcmp(double x)
{
    if (fabs(x)<eps) return 0;
    else return (x>0)?1:-1;
}
1
2
3
4
5
6
1
2
3
4
5
6

计算π的值

反余弦函数

const double pi=acos(-1.0);
1
1

定义一个向量

struct Vector
{
    double x,y;
    Vector(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
};
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8

定义一个点

二维坐标内的点也可以用向量来表示

typedef Vector Point;
1
1

定义一条直线

可以用直线上的两个点来存储，也可以用一个点和一个方向向量来存储，一般来说怎么方便怎么来。

struct Line
{
    Point p;
    Vector v;
    Line(Point P=Point(0,0),Vector V=Vector(0,0))
    {
        p=P,v=V;
    }
};
1
2
3
4
5
6
7
8
9
1
2
3
4
5
6
7
8
9

向量的四则运算和相等

• 加法：横纵坐标分别相加
• 减法：横纵坐标分别相减
• 数乘：横纵坐标分别乘以乘数
• 除法：横纵坐标分别除以除数
• 相等：横纵左边分别相等

Vector operator + (Vector A,Vector B) {return Vector(A.x+B.x,A.y+B.y);}
Vector operator - (Vector A,Vector B) {return Vector(A.x-B.x,A.y-B.y);}
Vector operator * (Vector A,double a) {return Vector(A.x*a,A.y*a);}
Vector operator / (Vector A,double a) {return Vector(A.x/a,A.y/a);}
bool operator == (Vector A,Vector B) {return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0;}
1
2
3
4
5
1
2
3
4
5

向量的旋转

a=(x,y)可以看成是x*(1,0)+y*(0,1) 
分别旋转两个单位向量 就变成了x*(cosα,sinα)+y*(-sinα,cosα) 
合并后得 (x*cosα-y*sinα,x*sinα+y*cosα)

Vector Rotate(Vector A,double rad)
{
    return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}
1
2
3
4
1
2
3
4

点积

求两个向量的点积

a·b的集合意义为a在b上的投影长度乘以b的模长 
a⋅b=|a||b|cos<a,b>，其中<a,b>为ab间的夹角 
坐标表示 
a=(x1,y1),b=(x2,y2),a⋅b=(x1⋅x2,y1⋅y2)

double Dot(Vector A,Vector B)
{
    return A.x*B.x+A.y*B.y;
}
1
2
3
4
1
2
3
4

点积的性质

• 点积满足交换律
• 两个向量夹角是一个锐角时，点积为正数；为直角时点积为零（因此可用点积判垂直）；为钝角时点积为负数

点积的应用

求向量的模长

double Len(Vector A)
{
    return sqrt(Dot(A,A));
}
1
2
3
4
1
2
3
4

求向量之间的夹角

double Angle(Vector A,Vector B)
{
    return acos(Dot(A,B)/Len(A)/Len(B));
}
1
2
3
4
1
2
3
4

单位向量与单位法向量

a|a|为a的单位向量，与a同向，模为1 
与a垂直的单位向量为a的单位法向量 
求单位法向量

Vector Normal(Vector A)
{
    double L=Len(A);
    return Vector(-A.y/L,A.x/L);
}
1
2
3
4
5
1
2
3
4
5

二维叉积

求两个向量的叉积

两个二维向量的叉积是一个标量，a×b的几何意义为ab所形成的平行四边形的有向面积 
坐标表示 
a=(x1,y1),b=(x2,y2),a×b=x1∗y2−x2∗y1

double Cross(Vector A,Vector B)
{
    return A.x*B.y-A.y*B.x;
}
1
2
3
4
1
2
3
4

叉积的性质

直观理解a×b，如果b在a的左边，叉积为正，右边为负，ab共线即为0 
pku 2318 
pku 2398

判断两条直线是否平行/相交/重合

平行：两条直线各任选两个点组成两个向量平行（叉积为0） 
相交：只需要判断是否平行即可（叉积为0 
重合：在平行的基础上，在两条直线上各选一个点组成一个向量在去与前两个判平行（叉积为0） 
pku 1269

判断直线和线段是否相交

利用叉积的性质 
在直线上任取两个点组成一个向量，再以其中的某一个点为起点到线段的两个端点分别组成两个向量，判断这两个向量和直线上向量的叉积是否同号，同号则不相交，异号则相交 
AB（直线）CD（线段）

double insLS(Point A,Point B,Point C,Point D)
{
    Vector v=B-A,w=C-A,u=D-A;
    return dcmp(Cross(v,w))!=dcmp(Cross(v,u));
}
1
2
3
4
5
1
2
3
4
5

pku 3304

判断两条线段是否相交

和上一个的方法类似，只需要两个线段判断两次即可。同为真才为真。 
pku 2653 
pku 2826 
pku 1410

点到直线的距离

用面积法

double DisTL(Point P,Point A,Point B)
{
    Vector v=B-A,w=P-A;
    return fabs(Cross(v,w)/Len(v));
}
1
2
3
4
5
1
2
3
4
5

点到线段的距离

分两种情况：一种是点到线段的垂线在线段上，这种情况同点到直线的距离，否则是点到点的距离 
判断是哪种情况时用点积的正负判断很方便

double DisTS(Point P,Point A,Point B)
{
    if (A==B) return Len(P-A);
    Vector v=B-A,w=P-A,u=P-B;
    if (dcmp(Dot(v,w))<0) return Len(w);
    else if (dcmp(Dot(v,u))>0) return Len(u);
    else return fabs(Cross(v,w)/Len(v));
}
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8

求直线的交点

根据直线的分点（比值）来求 
同样用面积法

Point GLI(Point P,Vector v,Point Q,Vector w)
{
    Vector u=P-Q;
    double t=Cross(w,u)/Cross(v,w);
    return P+v*t;
}
1
2
3
4
5
6
1
2
3
4
5
6

pku 1269

求三角形的外心

首先判断是否三点共线，共线的话即为线段的中点 
否则任求两边的垂直平分线的交点即可

Point TC(Point A,Point B,Point C)
{
    Point P=(A+B)/2,Q=(A+C)/2;
    Vector v=Rotate(B-A,pi/2),w=Rotate(C-A,pi/2);
    if (dcmp(Cross(v,w))==0)
    {
        if (dcmp(Len(A-B)+Len(B-C)-Len(A-C))==0)
            return (A+C)/2;
        if (dcmp(Len(A-C)+Len(B-C)-Len(A-B))==0)
            return (A+B)/2;
        if (dcmp(Len(A-B)+Len(A-C)-Len(B-C))==0)
            return (B+C)/2;
    }
    return GLI(P,v,Q,w);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

求三角形的重心

三角形的重心为三条中线的交点 
经过推导可得G=(A+B+C)/3

Point TG(Point A,Point B,Point C)
{
    return (A+B+C)/3;
}
1
2
3
4
1
2
3
4

判断点是否在多边形内

射线法 
注意判断点在多边形上的情况 
O(n)

bool PInP(Point P)
{
    int cnt=0,k,d1,d2;
    for (int i=1;i<=n;++i)
    {
        if (dcmp(DisTS(P,poly[i],poly[i%n+1]))==0) return false;
        k=dcmp(Cross(poly[i%n+1]-poly[i],P-poly[i]));
        d1=dcmp(poly[i].y-P.y);
        d2=dcmp(poly[i%n+1].y-P.y);
        if (k>0&&d1<=0&&d2>0) ++cnt;
        if (k<0&&d1>=0&&d2<0) --cnt;
    }
    if (cnt) return true;
    else return false;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

pku 1584

求多边形面积

假如是凸多边形，可以在多边形内任取一点，进行三角剖分 
因为叉积是有向面积，所以做法相同，基准点选在哪个点都可以 
那么我们不妨选第一个点为基准点 
O(n)

double Area()
{
    double area=0;
    for (int i=2;i<n;++i)
        area+=Cross(poly[i]-poly[1],poly[i+1]-poly[1]);
    return area/2;
}
1
2
3
4
5
6
7
1
2
3
4
5
6
7

pku 1654

Pick定理

一个计算点阵中顶点在格点上的多边形面积公式：S=a+b/2-1，其中a表示多边形内部的点数，b表示多边形边界上的点数，s表示多边形的面积。 
pku 2954 
pku 1265

凸包

把给定点包在内部，面积最小的凸多边形 
常用Graham扫描法

• 将点按x坐标第一关键字，y坐标第二关键字排序
• 首先将p1,p2压入栈中
• 判断下一个点和栈中最后一个点组成的向量是否在当前前进方向的左边（求下凸壳），如果是的话就入栈，否则弹栈
• 这样求出了下凸壳，再倒着做一遍就求出来了正凸壳，合起来就是完整的凸包

通过判断叉积是否=0可以控制是否有三点共线的情况 
最终栈中的点就是凸包中的点 
注意输入不能有重复的点！ 
还有就是坐标不能只按照x或y坐标排序，否则处理不了所有点都在一条直线上的情况 
时间复杂度O(nlogn)，瓶颈在排序上

struct Vector
{
    double x,y;
    Vector(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
    bool operator < (const Vector &a) const
    {
        return x<a.x||(x==a.x&&y<a.y);
    }
};
typedef Vector Point;
void Graham()
{
    sort(p+1,p+n+1);
    for (int i=1;i<=n;++i)
    {
        while (top>1&&dcmp(Cross(stack[top]-stack[top-1],p[i]-stack[top-1]))<=0)
            --top;
        stack[++top]=p[i];
    }
    int k=top;
    for (int i=n-1;i>=1;--i)
    {
        while (top>k&&dcmp(Cross(stack[top]-stack[top-1],p[i]-stack[top-1]))<=0)
            --top;
        stack[++top]=p[i];
    }
    if (n>1) --top;
}
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
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

COGS 896 
pku 1113 
pku 1228 
pku 1873 
pku 3348

凸包的应用——旋转卡壳

最远点对问题

给出平面上n个点的坐标，求最远点对 
最远点对一定在凸包上 
某错误结论：凸包上一个点到其它点的距离单峰

被凸包上被一对平行直线卡住的点对叫对踵点 
答案一定出在一对对踵点上 
尝试用通过旋转一对平行直线，枚举所有对踵点

做法一：直接模拟，每次比较哪条直线转到下一条边的角度较小，这种做法精度误差较大 
做法二：按逆时针顺序枚举凸包上一条边，则凸包上到该边所在直线最远的点也单调逆时针旋转，相当于用一条直线卡对面一个顶点

一般来说，做法二比较常用一些

之前让我困惑的是，为什么我们需要求两点之间的最大距离，而单调移动指针的时候保证的是点到直线的最大距离 
一个结论是，假设AB两点之间的距离最大，那么一定存在凸包上一条包含点A的边，并且B是距离这条边最远的点 
这样显得很复杂，其实就是从旋转卡壳最原始的做法考虑，卡住的点一定是对踵点，而等效成一条边和一个点，只有在满足点线距最大的情况下点才有可能成为对踵点，并且在旋转的时候每一个有可能成为对踵点的点都会考虑 
时间复杂度O(n)

double rotating()
{
    if (top==1) return 0;
    if (top==2) return Len(stack[2]-stack[1]);
    int now=1;
    double ans=0;
    for (int i=1;i<=top;++i)
    {
        while (dcmp(DisTL(stack[now],stack[i],stack[i%top+1])-DisTL(stack[now%top+1],stack[i],stack[i%top+1]))<=0)
            now=now%top+1;
        ans=max(ans,Len(stack[now]-stack[i]));
        ans=max(ans,Len(stack[now]-stack[i%top+1]));
    }
    return ans;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

pku 2187 
pku 2079 
bzoj 1069 
bzoj 1185

最小矩形覆盖问题

结论：矩形的某一条边一定和凸包的某一条边重合，并且边界点单调 
bzoj 1185

半平面交

给出若干个半平面，求它们的交 
也就是类似数学必修五的线性规划问题 
半平面交的答案是一个凸壳

首先初始化的时候给它加上一个很大的框 
然后每一次加入一条直线，将直线左边（右边）的点保留，并且计算出来这条直线和当前凸壳的交点加入 
时间复杂度O(n2)

void init()
{
    cnt=0;
    poly[++cnt]=Point(inf,inf);
    poly[++cnt]=Point(inf,-inf);
    poly[++cnt]=Point(-inf,-inf);
    poly[++cnt]=Point(-inf,inf);
}
void halfp(Point A,Point B)
{
    ncnt=0;
    Point C,D;
    for (int i=1;i<=cnt;++i)
    {
        C=poly[i%cnt+1];
        D=poly[(i+1)%cnt+1];
        if (dcmp(Cross(B-A,C-A))<=0)
            npoly[++ncnt]=C;
        if (insLS(A,B,C,D))
            npoly[++ncnt]=GLI(A,B-A,C,D-C);
    }
    cnt=ncnt;
    for (int i=1;i<=cnt;++i)
        poly[i]=npoly[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
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

pku 3130 
pku 1474 
pku 1279 
pku 3525 
pku 3384

最小圆覆盖

给出若干点，求一个最小的圆，覆盖所有的点 
三个点能确定一个圆，答案一定为某一个三角形的外接圆 
有一种非常神奇的复杂度为O(n)的随机增量法解决这个问题

算法流程

• 先对点集random_shuffle处理一下
• 假设已经有了前i-1个点的最小圆覆盖Ci−1，检验点pi是否在Ci−1内，如果在，则Ci=Ci−1，否则，pi一定在圆Ci上，进行下一步
• 包括点pi在内的前j-1个点(j<i)的最小圆覆盖为Cj−1，检验点pj是否在Cj−1内，如果在，则Cj=Cj−1，否则，pj一定在圆Cj上，进行下一步
• 包括点pi和点pj在内的前k-1个点(k<j)的最小圆覆盖为Ck−1，检验点pk是否在Ck−1内，如果在，则Ck=Ck−1，否则，pk一定在圆Ck上。返回

代码

void mcc()
{
    random_shuffle(p+1,p+n+1);
    c=p[1],r=0;
    for (int i=2;i<=n;++i)
        if (dcmp(Len(c-p[i])-r)>0)
        {
            c=p[i],r=0;
            for (int j=1;j<i;++j)
                if (dcmp(Len(c-p[j])-r)>0)
                {
                    c=(p[i]+p[j])/2,r=Len(c-p[i]);
                    for (int k=1;k<j;++k)
                        if (dcmp(Len(c-p[k])-r)>0)
                        {
                            c=TC(p[i],p[j],p[k]);
                            r=Len(c-p[i]);
                        }
                }
        }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

时间复杂度证明

目测复杂度是O(n3)的，实际上是O(n)的

• 首先最后一层循环是O(j)的
• 考虑倒数第二层循环，这个点有3j的概率不在圆内，这个时候要进行一个复杂度为O(j)的循环，所以期望是O(1)，总体复杂度O(i)
• 同理可得，第一层循环的复杂度是O(n)的

hdu 3007 
BZOJ 1336 
BZOJ 1337