【笔记篇】最良心的计算几何学习笔记(二)

依然放上本文的github地址...

作业QwQ

先来说一下上次留下的例题.
poj这道题并没有实数比较模式..
所以被精度势力干翻.
交上去WA掉竟然是因为-0.00和0.00不相等?
根据对拍结果别的地方应该没什么问题了OvO
下面给出并不能AC的"正确"代码:

#include <cstdio>
#include <cmath>
const double eps=1e-8;
inline int dcmp(const double &a){
  	if(fabs(a)<eps) return 0;
	return a<0?-1:1;
}
int errr=0;
inline int gn(int a=0,char c=0,int f=1){
	for(;(c<48||c>57)&&c!='-';c=getchar());
	if(c=='-')f=-1,c=getchar();
	for(;c>47&&c<58;c=getchar())a=a*10+c-48;return a*f;
}
struct vec{
	double x,y;
	vec(double p=0,double q=0):x(p),y(q){}
};
double operator *(const vec &a,const vec &b){
	return a.x*b.y-a.y*b.x;
} //向量运算只需要叉乘
inline vec Cut(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
	errr=0;
	double a1,b1,c1,a2,b2,c2,d;
	a1=p1.y-p0.y; b1=p0.x-p1.x; c1=p0*p1;
	a2=q1.y-q0.y; b2=q0.x-q1.x; c2=q0*q1;
	d=a1*b2-a2*b1;
	if(!dcmp(d)){ //斜率相等
		if(!dcmp(c1/a1-c2/a2)) errr=1; //平行
		else errr=-1;	//重合
		return vec(0,0);
	}
	return vec((b2*c1-b1*c2)/d,(a1*c2-a2*c1)/d);
}
int main(){
	int T=gn(); vec a[4];
	puts("INTERSECTING LINES OUTPUT");
	while(T--){
		for(int i=0;i<4;++i) a[i].x=gn(),a[i].y=gn();
		vec s=Cut(a[0],a[1],a[2],a[3]);
		if(errr<0) puts("NONE"); else if(errr>0) puts("LINE");
		else printf("POINT %.2lf %.2lf
",s.x,s.y);
	}
	puts("END OF OUTPUT");
}

新课OvO

这次学学多边形...
多边形就是多个点首尾顺次相接围成的图形.

多边形分为凹多边形和凸多边形.
大家见到的比较美观的多边形都是凸多边形, 比如正方形一类的.
它们的特点是比较可爱 任取一条边, 整个图形都在这条边所在直线的一侧.
这是一个凸多边形
而凹多边形则是不可爱满足这一点的多边形.
这是一个凹多边形

多边形的形式非常的多, 我们通常用一个储存了所有顶点坐标的数组来表示多边形.
关于多边形有如下常见的问题:

判断点是否在多边形内部

这是一个非常经典的问题, 导致出现了很多种解法.

射线法(交叉法)

以该点为起点, 做一条射线(一般方向取水平向右), 计算与多边形交点的奇偶性.
射线法

可以看到 多边形内部的点为起点做的射线与多边形的交点个数是奇数, 而以外部的点为起点做的射线与多边形的交点个数是偶数.
想法虽然很美好, 但是很显然, 这样做如果交点恰好是某个顶点就会出现问题.
这样我们就需要新的定义.
一种标准的做法是, 在左侧或者下边的点属于内部,在右侧或者顶边的点在外部.(就是计左不计右或计下不计上辣~)
这样有公共的边界线段的两个多边形就不能拥有公共点(起码我们不视为这样)
在这里我们的交叉法有一个专用的边缘交叉规则:

  1. 一条向上的边包括它的起始点, 不包括它的终止点.(记下不记上)
  2. 一条向下的边不包括它的起始点, 包括它的终止点.(记上不记下)
  3. 不计所有水平的边. (既不上也不下的特殊情况不予考虑)
  4. 交点必须严格出现在P的右侧. (右边界属于外部, 左边界属于内部)

举个栗子, 交点是顶点的情况大约有下面几种:
交点是顶点的几种情况

根据边缘交叉规则, 绿色的点才是满足条件的交点.
我们可以很容易的写出算法流程:

  1. 设点的坐标为((x_0,y_0)),
  2. 遍历多边形的每一条边e
  3. 看线段e是否跨立直线(y=y_0).
  4. 根据向上还是向下分别判断是否符合#1 #2 #3
  5. 求出e与直线(y=y_0)交点的横坐标(x), 跟(x_0)比较, 如果​(x>x_0)(满足#4), 则​(ans++).
  6. 判断(ans)的奇偶性, 奇数则在内部, 偶数则在外部.

一道本子……啊呸 板子题.. 稀有的hdu中文题目...
虽然题目描述并不是太严谨, 但是数据水啊, 但是题意很精炼.
刚才分析那么多了直接贴代码吧.(这个题1A了就很爽..

#include <cmath> 
#include <cstdio>
#include <cstring>
const double eps=1e-9;
int dcmp(const double &a){
	if(fabs(a)<eps) return 0;
	return a<0?-1:1;
}
struct point{
	double x,y;
}poly[202],pt;
double operator *(const point &a,const point &b){return a.x*b.y-a.y*b.x;}
//求p0p1两点确定的直线和直线y=y0的交点 
double Cutx(const double &y0,const point &p0,const point &p1){
	return ((y0-p0.y)*p1.x+(p1.y-y0)*p0.x)/(p1.y-p0.y);
}
//交叉法判断点是否在多边形内
//参数p表示点 pts表示多边形的顶点集 pcnt表示顶点数 
bool cn_ptInPoly(const point &p,point *pts,int pcnt){
	pts[pcnt]=pts[0]; int cnt=0;
	for(int i=1;i<=pcnt;++i)
		if(((dcmp(pts[i].y-p.y)>0&&dcmp(pts[i-1].y-p.y)<1)		//#1
			||(dcmp(pts[i].y-p.y)<1&&dcmp(pts[i-1].y-p.y>0)))	//#2(#3)
			&&(dcmp(p.x-Cutx(p.y,pts[i],pts[i-1]))<0))			//#4
			++cnt; //交点计数一次 
	return cnt&1; //交点奇数时在内, 偶数时在外.
}
int main(){
	int n,m;
	while(~scanf("%d",&n)){
		memset(poly,0,sizeof(poly));
		for(int i=0;i<n;++i) scanf("%lf%lf",&poly[i].x,&poly[i].y);
		scanf("%d",&m);
		for(int i=0;i<m;++i) {
			scanf("%lf%lf",&pt.x,&pt.y);
			puts(cn_ptInPoly(pt,poly,n)?"Yes":"No");
		} 
	}
}

这里的话唯一没有说的就是求交点了.
首先经历了一波特判我们知道一定是有且只有一个交点的..
而且纵坐标知道了只需要求横坐标, 就是解析几何化式子的过程了...
我们这里让点(p)的坐标为((x_p,y_p)), 线段两端点分别是((x_0,y_0),(x_1,y_1)).
因为是求交点所以日常联立.. 根据必修二所学知识, 我们直线方程可以设两点式.

[left{egin{matrix} y=y_p\ frac{x_p-x_0}{x_1-x_0}=frac{y-y_0}{y_1-y_0} end{matrix} ight. ]

然后将①代入②, 随便(???)化一波式子得到

[x_p=x=frac{(y_p-y_0)x_1+(y_1-y_p)x_0}{y_1-y_0} ]

或者什么其他式子也可以..只要能算就行..然后就可以了..

这个算法比较简单好用, 而且速度也不错. 时间复杂度是(O(n))的.. 也比较常用.

绕数法

这个射线法已经是比较好用多了, 但是有一个缺陷..
多边形还有一种分类是简单和非简单...
非简单多边形长得就比凹多边形更为难受...
它们竟然会自我重叠...
就像这样...(图中的数字表示顶点的顺序)
一个不"简单"的多边形

这个多边形就很"不简单".. 因为黄色的区域自我重叠了..
我知道这样非常的不好看, 但是就是有这样的图形我也很绝望啊←_←
然后再用射线法的时候就会出现一些问题..
射线法是不对的

你看这样一统计.. 交点个数是偶数 所以在外面?! 这样就出问题了...
那我们该怎么办?? 办法总比问题多, 我们可以再想有没有别的方法嘛...
有一种神奇的操作, 我们可以计算多边形绕某个点绕过的角度...
内-绕数

虽然画的有点乱,但是还是可以跟着紫色箭头看一下的..
先从这个点向每个顶点做射线, 然后从1绕一圈绕回1. 看转过的角度是不是0.
这个点在外部, 当且仅当绕过的角度是0.(其实就是没绕嘛, 这么就好理解了)
很明显图中的角度是(2pi)而不是0, 所以在内部.
在外部的点就是0了.像这样:
外-绕数

可以看到, 黄色区域的点也可以被判断到内部(不过绕了两圈)
黄-绕数

那我们如何统计呢?
因为都是向量的夹角, 我们又想到了点积. 我们最后要求的就是

[sum_{i=0}^{n-1}arccos(frac{vec{PV_i}vec{PV_{i+1}}}{|vec{PV_1}||vec{PV_{i+1}}|}) ]

但是这里我们用到了一个昂贵的函数arccos, 这样运行会变慢, 而且精度会受到影响.
或许最初这也就是绕数法不如射线法常用的原因?(据说可能会慢20倍?)
所以我们要考虑优化.
我也不知道谁想出了一种天才的优化方式..
这种方式借鉴了射线法, 依然做一条射线, 经过上行边时, 计数器+1; 经过下行边时, 计数器-1, 最后看计数器是否为0, 为0则点在多边形外, 否则点在多边形内.

举个栗子:
绕数优化

图上数得很清楚了 自己看看就能懂..
但是还有一些细节问题..
比如边界问题.. 见上面的边缘交叉规则即可..
在这里我们没有必要求出交点的坐标, 我们只需要判断交点会不会在定点的右边,
所以我们汇总一下一共有这几种情况:
#4的几种情况

对于上行边来说, 边的起始点指向定点的向量要在边的逆时针方向; 对于下行边来说, 则要在顺时针方向, 这样我们一点积就出来了..
还是给出hdu1756的代码(又是1A哈哈哈哈

#include <cmath> 
#include <cstdio>
#include <cstring>
const double eps=1e-9;
int dcmp(const double &a){
	if(fabs(a)<eps) return 0;
	return a<0?-1:1;
}
struct point{
	double x,y;
	point(double X=0,double Y=0):x(X),y(Y){}
}poly[202],pt;
point operator -(const point &a,const point &b){return point(a.x-b.x,a.y-b.y);}
double operator *(const point &a,const point &b){return a.x*b.y-a.y*b.x;}
//绕数法确定点是否在多边形内 
bool wn_ptInPoly(const point &p,point *pts,int pcnt){
	pts[pcnt]=pts[0]; int cnt=0;
	for(int i=1;i<=pcnt;++i)
		if((dcmp(pts[i].y-p.y)>0&&dcmp(pts[i-1].y-p.y)<1) 		//#1(#3) 
			&&dcmp((p-pts[i-1])*(pts[i]-pts[i-1]))<0) 			//#4逆时针 
			++cnt; 
		else if((dcmp(pts[i].y-p.y)<1&&dcmp(pts[i-1].y-p.y)>0)	//#2(#3) 
			&&dcmp((p-pts[i-1])*(pts[i]-pts[i-1]))>0)			//#4顺时针
			--cnt;
	return cnt;	//只要不为0 
}
int main(){
	int n,m;
	while(~scanf("%d",&n)){
		memset(poly,0,sizeof(poly));
		for(int i=0;i<n;++i) scanf("%lf%lf",&poly[i].x,&poly[i].y);
		scanf("%d",&m);
		for(int i=0;i<m;++i) {
			scanf("%lf%lf",&pt.x,&pt.y);
			puts(wn_ptInPoly(pt,poly,n)?"Yes":"No");
		} 
	}
}

不过竟然跑了15ms 不知道为什么跟射线法的0ms还有一定差距..
所以简单的多边形还是跑射线法似乎靠谱一点...
不过出于正确性的考虑, 如果不保证多边形简单的话, 还是要用绕数的..

二分法

我们之前说过, 凸多边形非常的美观.
美观的图形就一定要有美观的性质.
于是我们有一种特殊的判断点在多边形内的方法.
凸多边形哟~

我们来看这张图, 第一步, 我们可以做一个"快速排斥实验".
(P.S. 以下均假设点按照逆时针方向给出, 顺时针请自求多福做个对称..)
判断一下(vec {PV_0})(vec{PV_1})(vec{PV_8})的顺逆时针关系.
如果在(angle V_1V_0V_8)之外的话就不需要考虑了(因为是凸多边形).
然后我们在(1sim8)之间二分, 每次判断顺逆时针关系.
就可以在(O(logn))时间内确定是在(angle V_iV_0V_{i+1})内了.
然后我们再(O(1))判断(P)点在(vec {V_iV_{i+1}})左侧还是右侧(用叉积)就好了~
(没找到裸题TAT 那就只能干写了 正确性就不敢说了哟~)

//二分确定点是否在凸多边形内
//P.S. 假设点以逆时针给出 
bool ptInConvexPoly(const point &p,point *pts,int pcnt){
	point po=p-pts[0]; //比较常用,所以存一下 
	if(dcmp(po*(pts[1]-pts[0]))>0
		||dcmp(po*(pts[pcnt-1]-pts[0]))<0)
		return false;  //"快速排斥"一波
	int l=1,r=pcnt-1;
	while(r-l>1){
		int mid=(l+r)>>1;
		if(!dcmp(po*(pts[mid]-pts[0]))<0) //逆时针方向
			l=mid;
		else r=mid;
	}
	return dcmp((p-pts[l])*(pts[l+1]-pts[l]))<0; //在边的左侧 
} 

然后完整的代码依然去github找好了...

最后发现自己其实一天只学了判断点在多边形内..
心情复杂...

原文地址:https://www.cnblogs.com/enzymii/p/8413410.html