【笔记篇】最良心的计算几何学习笔记（二）
依然放上本文的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\n",s.x,s.y);
}
puts("END OF OUTPUT");
}
新课OvO
这次学学多边形...
多边形就是多个点首尾顺次相接围成的图形.
多边形分为凹多边形和凸多边形.
大家见到的比较美观的多边形都是凸多边形, 比如正方形一类的.
它们的特点是比较可爱 任取一条边, 整个图形都在这条边所在直线的一侧.
而凹多边形则是不可爱满足这一点的多边形.
多边形的形式非常的多, 我们通常用一个储存了所有顶点坐标的数组来表示多边形.
关于多边形有如下常见的问题:
判断点是否在多边形内部
这是一个非常经典的问题, 导致出现了很多种解法.
射线法(交叉法)
以该点为起点, 做一条射线(一般方向取水平向右), 计算与多边形交点的奇偶性.
可以看到 多边形内部的点为起点做的射线与多边形的交点个数是奇数, 而以外部的点为起点做的射线与多边形的交点个数是偶数.
想法虽然很美好, 但是很显然, 这样做如果交点恰好是某个顶点就会出现问题.
这样我们就需要新的定义.
一种标准的做法是, 在左侧或者下边的点属于内部，在右侧或者顶边的点在外部.(就是计左不计右或计下不计上辣~)
这样有公共的边界线段的两个多边形就不能拥有公共点(起码我们不视为这样)
在这里我们的交叉法有一个专用的边缘交叉规则:
- 一条向上的边包括它的起始点, 不包括它的终止点.(记下不记上)
- 一条向下的边不包括它的起始点, 包括它的终止点.(记上不记下)
- 不计所有水平的边. (既不上也不下的特殊情况不予考虑)
- 交点必须严格出现在P的右侧. (右边界属于外部, 左边界属于内部)
举个栗子, 交点是顶点的情况大约有下面几种:
根据边缘交叉规则, 绿色的点才是满足条件的交点.
我们可以很容易的写出算法流程:
- 设点的坐标为\((x_0,y_0)\),
- 遍历多边形的每一条边e
- 看线段e是否跨立直线\(y=y_0\).
- 根据向上还是向下分别判断是否符合#1 #2 #3
- 求出e与直线\(y=y_0\)交点的横坐标\(x\), 跟\(x_0\)比较, 如果\(x>x_0\)(满足#4), 则\(ans++\).
- 判断\(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)\).
因为是求交点所以日常联立.. 根据必修二所学知识, 我们直线方程可以设两点式.
然后将①代入②, 随便(???)化一波式子得到
或者什么其他式子也可以..只要能算就行..然后就可以了..
这个算法比较简单好用, 而且速度也不错. 时间复杂度是\(O(n)\)的.. 也比较常用.
绕数法
这个射线法已经是比较好用多了, 但是有一个缺陷..
多边形还有一种分类是简单和非简单...
非简单多边形长得就比凹多边形更为难受...
它们竟然会自我重叠...
就像这样...(图中的数字表示顶点的顺序)
这个多边形就很"不简单".. 因为黄色的区域自我重叠了..
我知道这样非常的不好看, 但是就是有这样的图形我也很绝望啊←_←
然后再用射线法的时候就会出现一些问题..
你看这样一统计.. 交点个数是偶数 所以在外面?! 这样就出问题了...
那我们该怎么办?? 办法总比问题多, 我们可以再想有没有别的方法嘛...
有一种神奇的操作, 我们可以计算多边形绕某个点绕过的角度...
虽然画的有点乱,但是还是可以跟着紫色箭头看一下的..
先从这个点向每个顶点做射线, 然后从1绕一圈绕回1. 看转过的角度是不是0.
这个点在外部, 当且仅当绕过的角度是0.(其实就是没绕嘛, 这么就好理解了)
很明显图中的角度是\(2\pi\)而不是0, 所以在内部.
在外部的点就是0了.像这样:
可以看到, 黄色区域的点也可以被判断到内部(不过绕了两圈)
那我们如何统计呢?
因为都是向量的夹角, 我们又想到了点积. 我们最后要求的就是
但是这里我们用到了一个昂贵的函数arccos, 这样运行会变慢, 而且精度会受到影响.
或许最初这也就是绕数法不如射线法常用的原因?(据说可能会慢20倍?)
所以我们要考虑优化.
我也不知道谁想出了一种天才的优化方式..
这种方式借鉴了射线法, 依然做一条射线, 经过上行边时, 计数器+1; 经过下行边时, 计数器-1, 最后看计数器是否为0, 为0则点在多边形外, 否则点在多边形内.
举个栗子:
图上数得很清楚了 自己看看就能懂..
但是还有一些细节问题..
比如边界问题.. 见上面的边缘交叉规则即可..
在这里我们没有必要求出交点的坐标, 我们只需要判断交点会不会在定点的右边,
所以我们汇总一下一共有这几种情况:
对于上行边来说, 边的起始点指向定点的向量要在边的逆时针方向; 对于下行边来说, 则要在顺时针方向, 这样我们一点积就出来了..
还是给出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\)之外的话就不需要考虑了(因为是凸多边形).
然后我们在\(1\sim8\)之间二分, 每次判断顺逆时针关系.
就可以在\(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找好了...
最后发现自己其实一天只学了判断点在多边形内..
心情复杂...
