远航休息栈
首页 > 笔记 > 计算几何学习笔记

计算几何学习笔记

原理比较少，主要涉及实现。

综述

计算几何computational geometry，研究几何模型和数据处理的学科，探讨几何形体的计算机表示。分析和综合，研究如何灵活、有效的建立几何形体的数学模型以及在计算机中更好地存储和管理这些模型数据。

基础定义

1.精度控制

因为double（long double）存在精度误差，所以不能直接用(a.x==b.x&&a.y==b.y)这种方式比较两个点（向量）是否相等

const double eps=1e-7;
int dcmp(double x)
{
	if (fabs(x)<eps) return 0;
	else
		if (x>0)
			return 1;
		else
			return -1;
}

比较实数的时候先做差再用dcmp比较

2.计算\(π\)

const double pi=acos(-1.0);

3.定义向量

按照我的理解，我觉得这里的向量可以看做一个方向的标记，它可以表示从原点向x,y移动，坐标系中所有的点都这样移动，它就可以表示直线的方向了。

struct Vector
{
	double x,y;
	Vector (double a=0,double b=0)
	{
		x=a,y=b;
	}
};

4.定义点

因为它和向量一样，就可以偷个懒

typedef Vector Point;

5.定义直线

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

struct Line
{
	Point p;
	Vector v;
	Line(Point P=Point(0,0),Vector V=Vector(0,0))
	{
		p=P,v=V;
	}
};

向量

1.四则运算

  • 加法：横纵坐标分别相加
  • 减法：横纵坐标分别相减
  • 数乘：横纵坐标分别乘以乘数
  • 除法：横纵坐标分别除以除数
  • 相等：横纵左边分别相等
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 b) {return Vector(a.x*b  ,a.y*b  );}
Vector operator / (Vector a,double b) {return Vector(a.x/b  ,a.y/b  );}
bool   operator ==(Vector a,Vector b) {return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}

2.旋转

我们让向量V逆时针旋转\(\theta\)度，则

images

推导公式

\[sin(a+b)=sina*cosb+sinb*cosa\]

\[cos(a+b)=cosa*cosb-sina*sinb\]

推导略

结果\((x*cos\theta-y*sin\theta,x*sin\theta+y*cos\theta)\)

Vector Rotate(Vector a,double rad)
{
	return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}

3.点积

\(a \cdot b\)的集合意义为a在b上的投影长度乘以b的长度

\( a \cdot b=|a||b|cos<a,b> \)

其中\(<a,b>\)表示ab之间的夹角

坐标表示

\[a=(x_1,y_1),b=(x_2,y_2),a \cdot b=(x_1 * x_2,y_1 * y_2)\]

证明

images

我们现在有两个向量,V1,V2,现在要求\(V_1 \cdot V_2\)

先从V2向V1做一条垂线,设它的长度为h,则

\(V_1 \cdot V_2=h \times |V1|\)

\(\because\)    \(V_1=\frac{y_1}{\cos \alpha}\,V_2=\frac{y_2}{\cos \beta} \)

\(\therefore\)    \(V_1 \cdot V_2=\frac{y_2}{\cos \beta } \times \cos (\beta – \alpha ) \times \frac{y_1}{\cos \alpha }\)

\(V_1 \cdot V_2=\frac{y_1y_2}{\cos \beta \cos \alpha } \times \cos (\beta – \alpha )\)

\(\because\)    \(\cos (\beta – \alpha )= \cos \alpha \cos \beta + \sin \alpha \sin \beta \)

\(\therefore\)    \(V_1 \cdot V_2=y_1y_2 + \sin \alpha \sin \beta \frac{y_1y_2}{\cos \beta \cos \alpha }\)

\(\because\)    \(\frac{\sin \alpha }{\cos \alpha }=\tan \alpha \)

\(\therefore\)    \(V_1 \cdot V_2=y_1y_2 + \tan \alpha \tan \beta y_1y_2\)

\(\therefore\)    \(V_1 \cdot V_2=y_1y_2 + x_1x_2\)

\(证毕\)

double Dot(Vector a,Vector b)
{
    return a.x*b.x+a.y*b.y;
}

点积的性质

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

点积的应用

求向量的模长

double Len(Vector a){return sqrt(Dot(a,a));}

求向量之间的夹角

double Angle(Vector a,Vector b){return acos(Dot(a,b)/Len(a)/Len(b));}

4.叉积

两个二维向量的叉积为一个标量

\(a \times b\)的几何意义就是ab所形成的平行四边形的有向面积。

坐标表示

\[a=(x_1,y_1),b=(x_2,y_2),a \times b=x_1 * y_2 – x_2 * y_1\]

double Cross(Vector a,Vector b)
{
    return a.x*b.y-a.y*b.x;
}

叉积的性质

顺着第一个向量V1看，假如V2在左边，叉积>0，否则<0，且\(Cross(V_1,V_2)=-Cross(V_2,V_1)\)，若两个同向，则叉积为0

images

poj 2318
poj 2398

5.单位向量与单位法向量

\(\frac{a}{|a|}\)为a的单位向量，方向相同，长度为1

与a垂直的是a的单位法向量

Vector Normal(Vector a)
{
    double L=Len(a);
    return Vector(-a.y/L,a.x/L);
}

6.判断向量关系

结合点积和叉积的结果，可以更好的判断两个向量的关系。

images

第一个数是点积，第二个数是叉积

7.其他

两向量垂直，点积为0

两向量共线（平行），叉积为0

点与直线

直线可以用直线上一点\(P_0\)和一个方向向量V表示。对于直线上所有的点P，\(P=P_0+tV\)，其中t叫作参数。

如果已知了直线上两点A,B，则V就为B-A，所以参数方程就为\(A+(B-A)t\)

参数方程可以供射线、线段使用，射线的t>0，线段的t在0~1之间。这样，直线的内容就可以应用到射线和线段了

1.判断直线关系

平行：两条直线各任选两个点组成两个向量平行（叉积为0）

相交：只需要判断是否平行即可（叉积为0 ）

重合：在平行的基础上，在两条直线上各选一个点组成一个向量在去与前两个判平行（叉积为0）

poj 1269

2.直线交点

设直线分别为\(P+tV\)和\(Q+tW\)且设向量\(u=P-Q\)，设交点在直线1的参数为\(t_1\)，交点在直线2的参数为\(t_2\)

这样可以列一个方程，解得

\(t_1=\frac{cross(w,u)}{cross(v,w)},t_2=\frac{cross(v,u)}{cross(v,w)}\)

代码（调用前确保有唯一交点，且cross(v,w)非零）

Point line_cross(Line a,Line b)
{
	Vector u=a.p-b.p;
	double t=cross(b.v,u)/cross(a.v,b.v);
	return a.p+t*a.v;
}

4.点到直线的距离

这个可以直接叉积求出，距离为\(h=\frac{V_1 \times V_2}{|V_1|}\)，即\(平行四边形的高=\frac{平行四边形的面积}{平行四边形的底}\)

代码：

double Dis_To_Line(Point p,Line l)
{
	return fabs(Cross(l.v-l.p,p-l.p)/Len(l.v-l.p));
}

5.点到线段的距离

分两种情况：一种是点到线段的垂线在线段上，这种情况同点到直线的距离，否则是点到点的距离

images

 

我们设投影点为Q，如果在AB上，就是P到直线AB的距离，不在AB上，假如在射线BA上，就为PA距离，不在就是PB距离。判断Q的位置可以用点积来进行。

double Dis_To_Seg(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(v);
	else return fabs(Cross(v,w)/Len(v));
}

6.点到直线的投影

还是上图，我们设\(直线AB为A+tV，V=\overrightarrow{AB}\)，设Q的参数为\(t_0\)，则\(Q=A+t_0V\)。根据\(PQ \perp AB\)，所以两个向量的点积应该为0。

所以\(V \cdot (P-(A+t_0V))=0\)，所以\(V \cdot (P-A) – t_0 * V \cdot V\)，就可以解了。

Point Line_Pro(Point p,Point a,Point b)
{
	Vector v=b-a;
	return a+v*(Dot(v,p-a)/Dot(v,v));
}

7.判断直线和线段是否相交

利用叉积的性质

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

bool Line_Seg_Cross(Line l,Point a,Point b)
{
	Vector v=l.p-l.v,w=a-l.p,u=b-l.p;
	return dcmp(Cross(v,w))!=dcmp(Cross(v,u));
}

poj 1269

8.判断线段和线段是否相交

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

bool Seg_Cross(Point a1,Point a2,Point b1,Point b2)
{
	double c1=Cross(a2-a1,b1-b1),c2=Cross(a2-a1,b2-a1),
	       c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);
	return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}

poj 2653
poj 2826
poj 1410

三角形

1.三角形面积

double Area_T(point a, point b,point c)
{
	return Cross(b-a,c-a)*0.5;
}

2.求三角形的外心

首先判断是否三点共线，共线的话即为线段的中点

否则任求两边的垂直平分线的交点即可

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)))
	{
        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 Line_cross(Line(p,v),Line(q,w));
}

3.求三角形的重心

三角形的重心为三条中线的交点

经过推导可得G=(A+B+C)/3

Point TG(Point a,Point b,Point c)
{
    return (a+b+c)/3;
}

多边形

1.多边形的面积

假如多边形的凸的，我们直接找一点进行三角剖分就可以了。

但是因为叉积是有向的，所以凹多边形可以直接跑。

double Area_P(Point ploy[],int num)
{
	double area=0;
	for (int i=2;i<num;i++)
		area+=Cross(ploy[i]-ploy[1],ploy[i+1]-ploy[1]);
	return area/2;
}

poj 1654

2.点是否在多边形内

从目标点出发引一条射线，看这条射线和多边形所有边的交点数目。如果有奇数个交点，则说明在内部，如果有偶数个交点，则说明在外部。

一般选的射线与x轴平行。

具体的话，判断一下该点在边的那边，在左边的话，就是大于第一个点的y，小于第二个点的y；右边是小于第一个点的y，大于第二个点的y

注意判断重合的情况。

bool Point_In_Ploy(Point p,Point poly[],int num)
{
	int ans=0,k,d1,d2;
	for (int i=1;i<=num;i++)
	{
		if (!dcmp(Dis_To_Seg(p,poly[i],poly[i%num+1]))) return 0;
		k=dcmp(Cross(poly[i%num+1]-poly[i],p-poly[i]));
		d1=dcmp(poly[i].y-p.y);
		d2=dcmp(poly[i%num+1].y-p.y);
		if (k>0&&d1<=0&&d2>0) ++ans;
		if (k<0&&d1>=0&&d2<0) --ans;
	}
	if (ans) return 1;
	return 0;
}

poj 1584

凸包

1.定义

凸包指的是把给定点包围在内部的、面积最小的凸多边形。

2.计算流程

常用Graham扫描法

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

通过判断叉积是否==0可以控制是否有三点共线的情况

最终栈中的点就是凸包中的点

需要判断，输入不能有重复的点

还有就是坐标不能只按照x或y坐标排序，否则处理不了所有点都在一条直线上的情况

时间复杂度O(nlogn)，排序慢

COGS 896
pku 1113
pku 1228
pku 1873
pku 3348

3.凸包应用：旋转卡壳

我们接着讲一讲旋转卡壳算法。（那个字念 qia 。。。搞了好久才发现读都读错了。）

旋转卡壳可以用于求凸包的直径、宽度，两个不相交凸包间的最大距离和最小距离等。虽然算法的思想不难理解，但是实现起来真的很容易让人“卡壳”。

images

其实简单来说就是用一对平行线“卡”住凸包进行旋转。

被一对卡壳正好卡住的对应点对称为对踵点，对锺点的具体定义不好说，不过从图上还是比较好理解的。

可以证明对踵点的个数不超过\( \frac{3N}{2}\)个 也就是说对踵点的个数是O(N)的

对踵点的个数也是我们下面解决问题时间复杂度的保证。

卡壳呢，具体来说有两种情况：

1.

images

一种是这样，两个平行线正好卡着两个点；

2.

images

一种是这样，分别卡着一条边和一个点。

而第二种情况在实现中比较容易处理，这里就只研究第二种情况。

在第二种情况中 我们可以看到 一个对踵点和对应边之间的距离比其他点要大

也就是一个对踵点和对应边所形成的三角形是最大的 下面我们会据此得到对踵点的简化求法。

images

下面给出一个官方的写法：

Compute the polygon’s extreme points in the y direction. Call them ymin and ymax. Construct two horizontal lines of support through ymin and ymax. Since this is already an anti-podal pair, compute the distance, and keep as maximum. Rotate the lines until one is flush with an edge of the polygon. A new anti-podal pair is determined. Compute the new distance, compare to old maximum, and update if necessary. Repeat steps 3 and 4 until the anti-podal pair considered is (ymin,ymax) again. Output the pair(s) determining the maximum as the diameter pair(s).

images

要是真的按这个实现起来就麻烦到吐了。。

根据上面的第二种情况，我们可以得到下面的方法：

如果qa，qb是凸包上最远两点，必然可以分别过qa，qb画出一对平行线。通过旋转这对平行线，我们可以让它和凸包上的一条边重合,如图中蓝色直线，可以注意到，qa是凸包上离p和qb所在直线最远的点。于是我们的思路就是枚举凸包上的所有边，对每一条边找出凸包上离该边最远的顶点，计算这个顶点到该边两个端点的距离，并记录最大的值。

直观上这是一个O(n2)的算法，和直接枚举任意两个顶点一样了。

然而我们可以发现 凸包上的点依次与对应边产生的距离成单峰函数（如下图：）

images

这个性质就很神了。

根据这个凸包的特性，我们注意到逆时针枚举边的时候，最远点的变化也是逆时针的，这样就可以不用从头计算最远点，而可以紧接着上一次的最远点继续计算。于是我们得到了O(n)的算法。这就是所谓的“旋转”吧！

利用旋转卡壳，我们可以在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(Dis_To_Line(stack[now],stack[i],stack[i%top+1])-Dis_To_Line(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;
}

其实是比较好理解的。

poj 2187
poj 2079
bzoj 1069
bzoj 1185

半平面交

1.什么是半平面？

顾名思义，半平面就是指平面的一半，我们知道，一条直线可以将平面分为两个部分，那么这两个部分就叫做两个半平面。

2.半平面怎么表示呢？

二维坐标系下，直线可以表示为ax + by + c = 0，那么两个半平面则可以表示为ax + by + c >= 0 和ax + by + c < 0，这就是半平面的表示方法。

3.半平面的交是啥啊？

其实就是一个方程组，让你画出满足若干个式子的坐标系上的区域（类似于线性规划的可行域），方程组就是由类似于上面的这些不等式组成的。

4.半平面交可以干什么？

半平面交虽然说是半平面的问题，但它其实就是关于直线的问题。一个一个的半平面其实就是一个一个有方向的直线而已。

半平面交的一个重要应用就是求多边形的核 。 多边形的核又是啥？  它是平面简单多边形的核是该多边形内部的一个点集，该点集中任意一点与多边形边界上一点的连线都处于这个多边形内部。就是一个在一个房子里面放一个摄像 头，能将所有的地方监视到的放摄像头的地点的集合即为多边形的核。经常会遇到让你判定一个多边形是否有核的问题。

5.如何求半平面交呢

我们把初始答案设成整个平面，然后逐渐加入每一个半平面，维护现在的半平面交。在实现中，我们一般都用四个人工半平面来框一个很大的矩形，这样加入半平面就相当于切割这个矩形。

切割非常简单，用叉积判断点的位置就可以了

Point poly[N],npoly[N];int cnt,ncnt;
void Crate_Square()
{
	cnt=0;
	poly[++cnt]=Point(INF,INF);
	poly[++cnt]=Point(INF,-INF);
	poly[++cnt]=Point(-INF,-INF);
	poly[++cnt]=Point(-INF,INF);
}
void Cut_Ploy(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 (Line_Seg_Cross(a,b,c,d))
			npoly[++ncnt]=Line_cross(Line(a,b-a),Line(c,d-c));
	}
	cnt=ncnt;
	for (int i=1;i<=ncnt;i++)
		poly[i]=npoly[i];
}

pku 3130
pku 1474
pku 1279
pku 3525
pku 3384

后记

个人认为，计算几何是一种非常优美的展现几何的一种方式，它把繁重的解析几何用优美简洁的向量、点来表示出来。今后还需多加练习。