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

最新推荐文章于 2023-02-18 02:14:16 发布

Enzymii

最新推荐文章于 2023-02-18 02:14:16 发布

阅读量506

点赞数 1

分类专栏：几何就是个无底洞文章标签：计算几何

本文链接：https://blog.csdn.net/Enzymii/article/details/78965395

版权

几何就是个无底洞专栏收录该内容

7 篇文章 2 订阅

订阅专栏

依然放上本文的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)$ .
因为是求交点所以日常联立.. ~~根据必修二所学知识~~, 我们直线方程可以设两点式.

⎧ ⎩ ⎨ ⎪ ⎪ y = y p x p - x 0 x 1 - x 0 = y - y 0 y 1 - y 0

$\left\{\begin{matrix} y=y_p\\ \frac{x_p-x_0}{x_1-x_0}=\frac{y-y_0}{y_1-y_0} \end{matrix}\right.$
然后将①代入②, 随便(???)化一波式子得到

x p = x = ( y p - y 0 ) x 1 + ( y 1 - y p ) x 0 y 1 - y 0

$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.(其实就是没绕嘛, 这么就好理解了)
很明显图中的角度是 $2\pi$ 而不是0, 所以在内部.
在外部的点就是0了.像这样:
外-绕数

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

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

\sum i = 0 n - 1 a r c c o s (P V i \to P V i + 1 \to | P V 1 \to | | P V i + 1 \to |)

$\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$ 之外的话就不需要考虑了(因为是凸多边形).
然后我们在 $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; //在边的左侧 
}