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

计算几何入门

计算几何基础

单位向量 $(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,但是浮点数doubleprintf系列函数中的标准占位符是%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,则称它们是正交的。如果能够定义向量间的夹角,则正交可以直观的理解为垂直。

如果对您有帮助,请我喝杯奶茶?