计算几何入门 | 李博的博客
0%

计算几何入门

发表于 更新于 分类于 阅读次数： Valine：

计算几何基础

单位向量 $(Unit\ vector)$

​ 对于任意向量 $\vec a$ ，不论方向如何，若其大小为单位长度，则称其为 $\vec a$ 方向上的单位向量 $(Unit\ vector)$ 。单位向量通常被记为 $\vec u$ 。
特殊地，三维笛卡尔坐标系上的三个基向量 $\vec i=(1,0,0),\vec j=(0,1,0),\vec k=(0,0,1)$ 都是单位向量。
点或向量的结构体如下：

1
2
3
4
5
6
7
8
9
// sgn返回x经过eps处理的符号，负数返回-1，正数返回1，x的绝对值如果足够小，就返回0。
const double eps = 1e-8;
int sgn(double x) { return x < -eps ? -1 : x > eps ? 1 : 0; }
struct Point{
    double x,y;
    Point(){}
    Point(double x,double y):x(x),y(y){}
};
typedef Point Vector;

向量相关的运算

向量与向量的加法

​ 向量的加法满足平行四边形法则和三角形法则。具体地，两个向量 $\vec a$ 和 $\vec b$ 相加，得到的是另一个向量。这个向量可以表示为 $\vec a$ 和 $\vec b$ 的起点重合后，以它们为邻边构成的平行四边形的一条对角线（以共同的起点为起点的那一条，见下图左），或者表示为将的 $\vec a$ 终点和 $\vec b$ 的起点重合后，从 $\vec a$ 的起点指向 $\vec b$ 的终点的向量：

1
2
3
Vector operator + (Vector A,Vector B){
    return Vector(A.x+B.x,A.y+B.y);
}

向量与向量的减法

​ 两个向量 $\vec a$ 和 $\vec b$ 的相减，则可以看成是向量 $\vec a$ 加上一个与 $\vec b$ 大小相等，方向相反的向量。又或者，$\vec a$ 和 $\vec b$ 的相减得到的向量可以表示为 $\vec a$ 和 $\vec b$ 的起点重合后，从 $\vec b$ 的终点指向 $\vec a$ 的终点的向量：

1
2
3
Vector operator - (Point A,Point B){
    return Vector(A.x-B.x,A.y-B.y);
}

向量与数的乘法

1
2
3
Vector operator * (Vector A,double p){
    return Vector(A.x*p,A.y*p);
}

向量与数的除法

1
2
3
Vector operator / (Vector A,double p){
    return Vector(A.x/p,A.y/p);
}

$Dot\ product$

  • ${\vec a}\cdot{\vec b}=|\vec a||\vec b|\cos{\theta}$
    • 这里$\ |\vec a|\ $表示$\vec a$的模（长度），$\theta$表示两个向量之间的角度。
  • 运算结果：标量
  • 用处：
    • 判断正交[^1]（点积为零）
    • 计算两向量的夹角 $cos\beta=\dfrac{ {\vec a}\cdot{\vec b} }{|\vec a||\vec b|}$
1
2
3
4
5
6
7
8
9
10
//点积
double Dot(Vector A,Vector B){
    return A.x*B.x+A.y*B.y;
}
double Length(Vector A){
    return sqrt(Dot(A,A));
}
double Angle(Vector A,Vector B){
    return acos(Dot(A,B)/Length(A)/Length(B));//用连除代替乘法，防止乘法溢出
}

$Cross\ product$

  • 运算结果：向量

  • 用处：计算面积或判断点的位置关系等。

  • 使用右手定则确定叉积的方向

  • 叉积 ${\vec a} \times {\vec b}$（垂直方向、紫色）随着向量 $\vec a$（蓝色）和 $\vec b$（红色）的夹角变化。 叉积垂直于两个向量，模长在两者平行时为零、在两者垂直时达到最大值 $|\vec a||\vec b|$。

1
2
3
4
5
6
7
8
//叉积
double Cross(Vector A,Vector B){
    return (A.x*B.y-A.y*B.x);
}
//三角形面积的二倍的叉乘公式
double Area2(Point A,Point B,Point C){
    return Cross(B-A,C-A);
}

向量的旋转

将向量 $\vec a=(x,y)$ 绕原点逆时针旋转$\alpha$度。

假设模长$d=\sqrt{x^2+y^2}$
那么$\vec a=(x,y)=(d\cos{\beta},d\sin{\beta})$
旋转$\alpha$度就是：
$$
\begin{eqnarray}
\vec{\grave a} &=& (d\cos{(\beta+\alpha)},d\sin{(\beta+\alpha)}) \\
&=& (d(\cos{\beta}\cdot\cos{\alpha}-\sin{\beta}\cdot\sin{\alpha}),
d(\sin{\beta}\cdot\cos{\alpha}+\cos{\beta}\cdot\sin{\alpha})) \\
&&将x=d\cos{\beta}, y=d\sin{\beta}代进去\\
&=& (x\cos{\alpha}-y\sin{\alpha},x\sin{\alpha}+y\cos{\alpha})
\end{eqnarray}
$$

1
2
3
4
//将向量p绕原点逆时针旋转a度
Point rotate(Point p, double a) {
    return Point(p.x*cos(a) - p.y*sin(a), p.x*sin(a) + p.y*cos(a));
}

例题

例题1

给定一个简单多边形，求其面积
输入：多边形（定点按照逆时针顺序排列）
输出：面积S

化繁为简，多边形一般分解为多个三角形。
首先考虑给定一个三角形，如何计算面积？

  1. 在解析几何里，$\triangle ABC$的面积可以通过如下方式求的：
    1. 获得点坐标
    2. 获得边长$a,b,c$
    3. 海伦公式求得面积$S=\sqrt{s(s-a)(s-b)(s-c)}$，其中$s=\dfrac{a+b+c}{2}$
  2. 在计算几何里，我们知道$\triangle ABC$的面积就是$\vec {AB}$和$\vec {AC}$叉积的绝对值的一半。

显然，利用叉积来计算面积的精度误差远小于第一种方法。
下面函数就可以计算多边形的有向面积了。适合凸多边形[^2]和凹多边形。

1
2
3
4
5
6
7
8
//多边形的有向面积
double PolygonArea(Point* p,int n){
    double area=0;
    for(int i=1;i<n-1;i++){
        area+=Cross(p[i]-p[0],p[i+1]-p[0]);
    }
    return area/2;
}

例题2

例题： https://codeforces.com/contest/498/problem/A
题意：问你你家到你学校的平面上有很多直线（给出一般式的系数），问你至少穿过多少条直线才能到达学校？
代码题解： https://blog.csdn.net/intmainhhh/article/details/96175992

判断点$A(x_1,y_1)$和$B(x_2,y_2)$是否在直线$f(x,y)=ax+by+c=0$的两侧。

  • 直接将点$A$和$B$带入$f(x,y)$，若$f(x1,y1)\times f(x2,y2) < 0$，则两点位于直线两侧。
  • 显然，$f(x,y)=0$时说明$(x,y)$在直线上。

二维几何基础算法部件

点在线段之间

  1. 首先需要满足点在线段所在直线
  2. 然后点的 $x$ 在线段端点之间
1
2
3
bool PointOnSegment(Point p, Point a, Point b) {
    return !sgn(Cross(p-a, b-a)) && sgn(Dot(p-a, p-b)) <= 0;
}

求两直线交点

前提是两直线不平行（叉积为0）
用参数方程表示直线：两直线分别为$P+t\vec v$ 和$Q+t\vec w$
设向量$u=P-Q$，交点在第一条直线的参数为$t_1$，第二条直线上的参数为$t_2$，两直线连立可以解的：
$t_1=\dfrac{Cross(w,u)}{Cross(v,w)}$，$t_2=\dfrac{Cross(v,u)}{Cross(v,w)}$
计算$t_1$过程如下，$t_2$同理
$$
\begin{eqnarray}
P+t_1\vec v &=& Q+t_2\vec w \\
P-Q &=& \vec u = t_2\vec w-t_1\vec v \\
\vec u\times\vec w &=& t_2\vec w\times \vec w-t_1\vec v\times\vec w \\
\vec u\times\vec w &=& -t_1\vec v\times\vec w \\
t_1 &=& \dfrac{Cross(w,u)}{Cross(v,w)}
\end{eqnarray}
$$

1
2
3
4
5
6
7
//求两直线交点
//调用前请确保P+tv和Q+tw有唯一交点，当且仅当Cross(v,w)非0
Point GetLineIntersection(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
//点到直线的距离
double DistanceToLine(Point P,Point A,Point B){
    Vector v1=B-A,v2=P-A;
    return fabs(Cross(v1,v2))/Length(v1);
}

点到线段的距离

显然和点到直线的距离不一样，这个需要判断点到线段所在直线的垂足是否在线段上。
自己画个图就明白了。

1
2
3
4
5
6
7
8
//点到线段的距离
double DistanceToSegment(Point P,Point A,Point B){
    if(A==B) return Length(P-A);
    Vector v1=B-A,v2=P-A,v3=P-B;
    if(sgn(Dot(v1,v2))<0) return Length(v2);
    else if(sgn(Dot(v1,v3))>0) return Length(v3);
    else return fabs(Cross(v1,v2))/Length(v1);
}

点在直线上的投影

设$Q=A+t_0\vec v$，$\vec v$是$\vec{AB}$
因为$PQ$ $\bot$ $AB$，两个向量点积为0，因此：$Dot(\vec v,P-(A+t_0\vec v))=0$
根据分配律(中学书上有证明，自己画个图做投影就能看出来)有：$Dot(\vec v,P-A)-t_0\times Dot(\vec v,\vec v)=0$
因此：$t_0=\dfrac{Dot(\vec v,P-A)}{Dot(\vec v,\vec v)}$

1
2
3
4
5
//点在直线上的投影
Point GetLineProjection(Point P,Point A,Point B){
    Vector v=B-A;
    return A+v*(Dot(v,P-A)/Dot(v,v));
}

线段相交判定

首先规定：“规范相交”指线段相交，并且交点不是线段端点。“非规范相交”可以在端点。

1
2
3
4
5
6
7
8
9
10
//判断两线段严格相交
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
    double c1=Cross(a2-a1,b1-a1), c2=Cross(a2-a1,b2-a1),
           c3=Cross(b2-b1,a1-b1), c4=Cross(b2-b1,a2-b1);
    return sgn(c1)*sgn(c2)<0 && sgn(c3)*sgn(c4)<0;
}
//判断点是否在线段上
bool PointOnSegment(Point P,Point a1,Point a2){
    return sgn(Cross(a1-P,a2-P))==0 && sgn(Dot(a1-P,a2-P))<=0;
}

点在多边形内判定

主要有两种方法：1，射线法。2，转角法。
假设多边形顶点按照逆时针顺序排列。

  1. 射线法：判断点在多边形内：从该点做一条水平向右的射线，统计射线与多边形相交的情况，若相交次数为偶数，则说明该点在形外，否则在形内。为了便于交点在定点或射线与某些边重合时的判断，可以将每条边看成左开右闭的线段，即若交点为左端点就不计算。

    一般用转角法比较简单，感兴趣的自己根据上面的步骤实现代码。

  2. 转角法：看多边形相对于这个点转了多少度。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    int PointInPolygon(Point p,Point*poly,int n)
    {
        int wn=0;
        for(int i=0;i<n;++i){
            if(PointOnSegment(p,poly[i],poly[(i+1)%n])) return -1;//在边界上
            int k=sgn( Cross(poly[(i+1)%n]-poly[i], p-poly[i] ) );//判断点在边的哪侧,1:左,-1:右
            int d1=sgn( poly[i].y-p.y );
            int d2=sgn( poly[(i+1)%n].y-p.y );
            if(k>0 && d1<=0 && d2>0) wn++;//左开右闭，没有k==0的舍去重边
            if(k<0 && d2<=0 && d1>0) wn--;
        }
        if(wn!=0) return 1;//内部
        return 0;//外部
    }

与圆有关的相关计算

自学

凸包

凸包就是把所有点包起来并且面积最小的凸多边形[^2]。

常用的求解凸包的算法有$Graham$算法和$Andrew$算法。

这里仅介绍基于水平排序[^3]的$Andrew$算法。先按照水平序排序，删除重复点后得到序列$p_1,p_2,···$，然后把$p_1,p_2$放到凸包中，从$p_3$开始，当新点在凸包”前进“方向的左边时继续，否则一次删除最近加入凸包的点，直到新点在左边。

然后看个完整的过程

1
2
3
bool operator < (const Point& a,const Point& b){
    return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
//当精度要求高时，需要用sgn三态函数进行比较
int ConvexHull(Point* p,int n,Point* ch){
    sort(p,p+n);//先比较x坐标，再比较y坐标
    n=unique(p,p+n)-p;//去重，这个算法如果有重复的点就会出错
    int m=0;
    for(int i=0;i<n;++i){
        while(m>1 && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-1])<=0 ) --m;//如果不希望在凸包边上有点，可以将<=改为<
        ch[m++]=p[i];
    }
    int k=m;
    for(int i=n-2;i>=0;i--){
        while(m>k && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-1])<=0) --m;
        ch[m++]=p[i];
    }
    if(n>1) m--;
    return m;
}

半平面交

自学

平面区域

自学

小技巧和小陷阱

  1. double hypot(double x,double y)用来求$\sqrt{x^2+y^2}$

  2. atan2(y,x)函数用来求向量(x,y)的极角（单位：弧度）。三角函数都很慢，不可大量使用。

  3. 有些题目明确说明了坐标是整型，那么我们一定要把点的结构体中的坐标类型改成int，否则精度上容易出错。

  4. 因为被计算机表示浮点数的方式所限制，CPU在进行浮点数计算时会出现误差。如执行0.1 + 0.2 == 0.3结果往往为false，在四则运算中，加减法对精度的影响较小，而乘法对精度的影响更大，除法的对精度的影响最大。所以，在设计算法时，为了提高最终结果的精度，要尽量减少计算的数量，尤其是乘法和除法的数量。
    浮点数与浮点数之间不能直接比较，要引入一个eps常量。eps是epsilon（$\epsilon$）的简写，在数学中往往代表任意小的量。在对浮点数进行大小比较时，如果他们的差的绝对值小于这个量，那么我们就认为他们是相等的，从而避免了浮点数精度误差对浮点数比较的影响。eps在大部分题目时取1e-8就够了，但要根据题目实际的内容进行调整。

    1
    2
    3
    // sgn返回x经过eps处理的符号，负数返回-1，正数返回1，x的绝对值如果足够小，就返回0。
    const double eps = 1e-8;
    int sgn(double x) { return x < -eps ? -1 : x > eps ? 1 : 0; }

  5. 用scanf输入浮点数时，double的占位符是%lf，但是浮点数double在printf系列函数中的标准占位符是%f而不是%lf，使用时最好使用前者，因为虽然后者在大部分的计算机和编译器中能得到正确结果，但在有些情况下会出错（比如在POJ上）。

  6. 当提供给C语言中的标准库函数double sqrt (double x)的x为负值时，sqrt会返回nan，输出时会显示成nan或-1.#IND00（根据系统的不同）。在进行计算几何编程时，经常有对接近零的数进行开方的情况，如果输入的数是一个极小的负数，那么sqrt会返回nan这个错误的结果，导致输出错误。解决的方法就是将sqrt包装一下，在每次开方前进行判断。double mysqrt(double x){return max(0.0, sqrt(x))};

  7. 大部分的标程的输出是不会输出负零的，有时这样的结果是错误的，所以在没有Special Judge的题目要求四舍五入时，不要忘记对负零进行特殊判断。

涉及的定义

[^3]: 水平排序：$Andrew$算法水平排序是按照 $x$ 从小到大排序(如果 $x$ 相同，按照 $y$ 从小到大排序)，在 $Graham$ 中先排$y$
[^2]: 凸多边形：把凸多边形的任意一条边向两边无限延长为一条直线时，其他各边都在此直线的同旁，那么这个多边形叫做凸多边形。
[^1]: 正交：是线性代数的概念，是垂直这一直观概念的推广。作为一个形容词，只有在一个确定的内积空间中才有意义。若内积空间中两向量的内积为0，则称它们是正交的。如果能够定义向量间的夹角，则正交可以直观的理解为垂直。

如果对您有帮助，请我喝杯奶茶？