详谈判断点在多边形内的七种方法(最全面) hdu1756 hrbust1429 为例

来源:互联网 发布:淘宝业务流程分析 编辑:程序博客网 时间:2024/06/05 23:05

这几天在学计算几何,学到点定位的判断点在多边形内,书上提到了三种方法,但是有些方法的代码不全。于是网上找了找,又发现更多判断的方法,一时兴起决定学习一下,看看到底有多少种,结果一个大坑。。。
网上好多介绍的不详细(特别是转角法,最后还是google出来的),而且有些方法叫不同的名字,有点难搞啊,花了我一天多的时间。。TAT

话不多说,下面分享一下。有些方法我会介绍清楚但不会画图详解,希望大家自己画图体会。

涉及到的题目
hdu1756
就是裸题,多边形点顺时针给出(有的算法会强调给出顺序)

射线法

时间复杂度:O(n) 适用范围:任意多边形
个人认为是非常不错的算法(不需考虑精度误差和多边形点给出的顺序),可以作为第一选择。

算法思想:
以被测点Q为端点,向任意方向作射线(一般水平向右作射线),统计该射线与多边形的交点数。如果为奇数,Q在多边形内;如果为偶数,Q在多边形外。计数的时候会有一些特殊情况,如图

图片已经把特殊情况和算法实现说的很清楚了,下面我直接贴代码,具体可看代码注释。

const double eps = 1e-6;const double PI = acos(-1);//三态函数,判断两个double在eps精度下的大小关系int dcmp(double x){    if(fabs(x)<eps) return 0;    else        return x<0?-1:1;}//判断点Q是否在P1和P2的线段上bool OnSegment(Point P1,Point P2,Point Q){    //前一个判断点Q在P1P2直线上 后一个判断在P1P2范围上    return dcmp((P1-Q)^(P2-Q))==0&&dcmp((P1-Q)*(P2-Q))<=0;}//判断点P在多边形内-射线法bool InPolygon(Point P){    bool flag = false; //相当于计数    Point P1,P2; //多边形一条边的两个顶点    for(int i=1,j=n;i<=n;j=i++)    {        //polygon[]是给出多边形的顶点        P1 = polygon[i];        P2 = polygon[j];        if(OnSegment(P1,P2,P)) return true; //点在多边形一条边上        //前一个判断min(P1.y,P2.y)<P.y<=max(P1.y,P2.y)        //这个判断代码我觉得写的很精妙 我网上看的 应该是大神模版        //后一个判断被测点 在 射线与边交点 的左边        if( (dcmp(P1.y-P.y)>0 != dcmp(P2.y-P.y)>0) && dcmp(P.x - (P.y-P1.y)*(P1.x-P2.x)/(P1.y-P2.y)-P1.x)<0)            flag = !flag;    }    return flag;}

hdu1756-射线法代码

#include <bits/stdc++.h>using namespace std;typedef long long ll;const double eps = 1e-6;const double PI = acos(-1);int n,m;struct Point{    double x,y;    Point(double x=0,double y=0):x(x),y(y){}    //向量+    Point operator +(const Point &b)const    {        return Point(x+b.x,y+b.y);    }    //向量-    Point operator -(const Point &b)const    {        return Point(x-b.x,y-b.y);    }    //点积    double operator *(const Point &b)const    {        return x*b.x + y*b.y;    }    //叉积    //P^Q>0,P在Q的顺时针方向;<0,P在Q的逆时针方向;=0,P,Q共线,可能同向或反向    double operator ^(const Point &b)const    {        return x*b.y - b.x*y;    }}polygon[105];typedef Point Vector;//三态函数,判断两个double在eps精度下的大小关系int dcmp(double x){    if(fabs(x)<eps) return 0;    else        return x<0?-1:1;}//判断点Q是否在P1和P2的线段上bool OnSegment(Point P1,Point P2,Point Q){    return dcmp((P1-Q)^(P2-Q))==0&&dcmp((P1-Q)*(P2-Q))<=0;}//判断点P在多边形内-射线法bool InPolygon(Point P){    bool flag = false;    Point P1,P2; //多边形一条边的两个顶点    for(int i=1,j=n;i<=n;j=i++)    {        P1 = polygon[i];        P2 = polygon[j];        if(OnSegment(P1,P2,P)) return true; //点在多边形一条边上        if( (dcmp(P1.y-P.y)>0 != dcmp(P2.y-P.y)>0) && dcmp(P.x - (P.y-P1.y)*(P1.x-P2.x)/(P1.y-P2.y)-P1.x)<0)            flag = !flag;    }    return flag;}int main(){    while(~scanf("%d",&n))    {        for(int i=1;i<=n;i++) scanf("%lf %lf",&polygon[i].x,&polygon[i].y);        Point test;        scanf("%d",&m);        while(m--)        {            scanf("%lf %lf",&test.x,&test.y);            if(InPolygon(test)) printf("Yes\n");            else printf("No\n");        }    }    return 0;}

角度和判断法

时间复杂度:O(n) 适用范围:任意多边形
感觉这个方法和之后要介绍的转角法类似,个人感觉转角法就是这个方法的优化变形。个人非常不推荐这个算法,最好就是不用。这个算法对精度的要求很高(会造成很大精度误差),不强调多边形点给出顺序,我用这个算法没过hdu1756。

算法思想:
连接被测点与多边形所有顶点所形成的所有角的角度和在精度范围内等于2π则该点在多边形内,否则在多边形外。
不推荐,所以就不画图了。

const double eps = 1e-6;const double PI = acos(-1);//三态函数,判断两个double在eps精度下的大小关系int dcmp(double x){    if(fabs(x)<eps) return 0;    else        return x<0?-1:1;}//判断点Q是否在P1和P2的线段上bool OnSegment(Point P1,Point P2,Point Q){    //前一个判断点Q在P1P2直线上 后一个判断在P1P2范围上    return dcmp((P1-Q)^(P2-Q))==0&&dcmp((P1-Q)*(P2-Q))<=0;}//判断点P在多边形内-角度和判断法//精度要求高 最好不用bool InPolygon(Point P){    double angle = 0;    Point P1,P2; //多边形一条边的两个顶点    Vector V1,V2; //以被测点为原点 P1 P2与P形成的向量    for(int i=1,j=n;i<=n;j=i++)    {        P1 = polygon[i];        P2 = polygon[j];        if(OnSegment(P1,P2,P)) return true; //点在多边形一条边上        V1 = P1-P;        V2 = P2-P;        double res = atan2(V2.y,V2.x)-atan2(V1.y,V1.x);        res = abs(res);        if(dcmp(res-PI)>0) res = 2*PI-res;        angle += res;    }    return dcmp(2*PI-angle)==0;}

atan2(y,x)计算y/x的反tan值,在[π,+π]范围内。
hdu1756的该方法的代码没过,就不贴了。

转角法

时间复杂度:O(n) 适用范围:任意多边形
个人感觉是O(n)算法的第二推荐,该算法本来对精度要求较高,之后会有一个改进让其不用考虑精度误差,不过该算法要强调多边形点给出的顺序。
一般博客都以多边形正向即逆时针介绍,我这里也主要介绍逆时针,但hdu1756是顺时针给出,我会在括号中介绍一下顺时针(其实本质是一样的),顺时针具体在代码注释中提一下。不会画图,希望大家看了思想自己画图体会一下。

算法思想:
转角法非常简单,按照多边形顶点逆时针顺序,从P点到顶点Vi分别做连线,其中αi为Vi和Vi+1之间的夹角。其中α角度逆时针为正,顺时针为负,这样所有到顶点做连线之间夹角和为(环绕数)0,这点P在多边形外部,否则在内部。(感觉和角度和判断法本质一样,加了个方向)
(顺时针就是角度顺时针为正,逆时针为负)

直接环绕数的推导会需要用到反三角函数,这样即会耗时又会造成较大的精度误差,所以这里有一个优化。
从P点向右做射线R,如果边从射线R下方跨到上方,那么穿越+1,如果从上方跨到下方,则是-1。最终和为wn环绕数。如下图所示:
(写博客的时候才发现其实本质不就是射线吗,不过理解后代码会感觉写的比射线简单)

这种方法不必去计算射线和边的交点,但需要判断点P和边的左右关系,而且对于方向向上和向下的边的判断规则不同。对于方向向上的边,如果穿过射线,那么P是在有向边的左侧;方向向下的边如果穿过射线,那么P在有向边的右边(意思是说判断点P与边的关系,而不是相对坐标系内位置)。

这里有一点要注意,如下图

这里射线经过了BCCD并穿过C点,但是要计算两次,一次+1一次-1,代码中会有体现。

const double eps = 1e-6;const double PI = acos(-1);//三态函数,判断两个double在eps精度下的大小关系int dcmp(double x){    if(fabs(x)<eps) return 0;    else        return x<0?-1:1;}//判断点Q是否在P1和P2的线段上bool OnSegment(Point P1,Point P2,Point Q){    //前一个判断点Q在P1P2直线上 后一个判断在P1P2范围上    return dcmp((P1-Q)^(P2-Q))==0&&dcmp((P1-Q)*(P2-Q))<=0;}//判断点P在多边形内-转角法(多边形点顺时针给出)bool InPolygon(Point P){    int wn = 0;    Point P1,P2; //多边形一条边的两个顶点    Vector V1,V2; //被测点与P1分别于P2形成的向量    for(int i=1,j=n;i<=n;j=i++)    {        P1 = polygon[j];  //顺时针中前一点        P2 = polygon[i];  //顺时针中后一点        if(OnSegment(P1,P2,P)) return true; //点在多边形一条边上        V1 = P2-P1;        V2 = P-P1;        int k = dcmp(V1^V2); //用于判断被测点在有向边的左右        int d1 = dcmp(P1.y-P.y); //用于判断向上还是向下穿过        int d2 = dcmp(P2.y-P.y);        //V1在V2的顺时针方向即测试点在有向边左边 并且有向边向上穿过        if(k>0 && d1<=0&&d2>0) wn--;          //V1在V2的逆时针方向即测试点在有向边右边 并且有向边向下穿过        if(k<0 && d1>0&&d2<=0) wn++;  //上面wn+和wn- 一个允许起点在射线上另一个允许终点在射线上 最后特殊情况就会算两次//逆时针的wn+和wn-自己画图体会一下    }    //不为0即在多边形内 不管是正还是负    return wn!=0;}

hdu1756-转角法代码

#include <bits/stdc++.h>using namespace std;typedef long long ll;const double eps = 1e-6;const double PI = acos(-1);int n,m;struct Point{    double x,y;    Point(double x=0,double y=0):x(x),y(y){}    //向量+    Point operator +(const Point &b)const    {        return Point(x+b.x,y+b.y);    }    //向量-    Point operator -(const Point &b)const    {        return Point(x-b.x,y-b.y);    }    //点积    double operator *(const Point &b)const    {        return x*b.x + y*b.y;    }    //叉积    //P^Q>0,P在Q的顺时针方向;<0,P在Q的逆时针方向;=0,P,Q共线,可能同向或反向    double operator ^(const Point &b)const    {        return x*b.y - b.x*y;    }}polygon[105];typedef Point Vector;//三态函数,判断两个double在eps精度下的大小关系int dcmp(double x){    if(fabs(x)<eps) return 0;    else        return x<0?-1:1;}//判断点Q是否在P1和P2的线段上bool OnSegment(Point P1,Point P2,Point Q){    return dcmp((P1-Q)^(P2-Q))==0&&dcmp((P1-Q)*(P2-Q))<=0;}//判断点P在多边形内-转角法(多边形点顺时针给出)bool InPolygon(Point P){    int wn = 0;    Point P1,P2; //多边形一条边的两个顶点    Vector V1,V2; //被测点与P1分别于P2形成的向量    for(int i=1,j=n;i<=n;j=i++)    {        P1 = polygon[j];  //顺时针中前一点        P2 = polygon[i];  //顺时针中后一点        if(OnSegment(P1,P2,P)) return true; //点在多边形一条边上        V1 = P2-P1;        V2 = P-P1;        int k = dcmp(V1^V2);        int d1 = dcmp(P1.y-P.y);        int d2 = dcmp(P2.y-P.y);        if(k>0 && d1<=0&&d2>0) wn--;  //测试点在有向边左边 有向边向上穿过        if(k<0 && d1>0&&d2<=0) wn++;  //测试点在有向边右边 有向边向下穿过    }    return wn!=0;}int main(){    while(~scanf("%d",&n))    {        for(int i=1;i<=n;i++) scanf("%lf %lf",&polygon[i].x,&polygon[i].y);        Point test;        scanf("%d",&m);        while(m--)        {            scanf("%lf %lf",&test.x,&test.y);            if(InPolygon(test)) printf("Yes\n");            else printf("No\n");        }    }    return 0;}

改进弧长法

时间复杂度:O(n) 适用范围:任意多边形
该算法感觉是转角法的另一种优化,也解决了传统转角法的精度问题,也要求多边形点给出的顺序。

算法思想:
以被测点O为坐标原点,将平面划分为4个象限,对每个多边形顶点P[i],计算其所在的象限,然后顺序访问多边形的各个顶点P[i],分析P[i]和P[i+1],有下列三种情况:
1. P[i+1]在P[i]的下一象限。此时弧长和加π/2;(代码中+1)
2. P[i+1]在P[i]的上一象限。此时弧长和减π/2;(代码中-1)
3. P[i+1]在Pi的相对象限。利用叉积res=OP[i]xOP[i+1]计算OP[i]与OP[i+1]的关系。
若f=0,OP[i]与OP[i+1]共线,点在多边形边上;若f<0,OP[i]在OP[i+1]逆时针方向,弧长和减π(代码中-2);若f>0,OP[i]在OP[i+1]顺时针方向,弧长和加π(代码中+2)。

有点累,逆时针的代码自己画图推一推吧,

//判断点P在多边形内-改进弧长法(多边形点顺时针给出)//用这个还不如用上一个转角法bool InPolygon(Point P){    int q1,q2,ans=0;    Point P1,P2;    Vector V1,V2;    for(int i=1,j=n;i<=n;j=i++)    {        P1 = polygon[j];        P2 = polygon[i];        V1 = P1-P;        V2 = P2-P;        if(OnSegment(P1,P2,P)) return true;        q1 = V1.x>0 ? (V1.y>0 ? 0:3) : (V1.y>0 ? 1:2);        q2 = V2.x>0 ? (V2.y>0 ? 0:3) : (V2.y>0 ? 1:2);        int g = (q2-q1+4)%4;        if(g==1) ans--; //在上一象限        if(g==3) ans++; //在下一象限        if(g==2) dcmp(V1^V2)>0 ? (ans-=2) : (ans+=2); //在相对象限    }    return ans!=0;}

hdu1756-改进弧长法代码

#include <bits/stdc++.h>using namespace std;typedef long long ll;const double eps = 1e-6;const double PI = acos(-1);int n,m;struct Point{    double x,y;    Point(double x=0,double y=0):x(x),y(y){}    //向量+    Point operator +(const Point &b)const    {        return Point(x+b.x,y+b.y);    }    //向量-    Point operator -(const Point &b)const    {        return Point(x-b.x,y-b.y);    }    //点积    double operator *(const Point &b)const    {        return x*b.x + y*b.y;    }    //叉积    //P^Q>0,P在Q的顺时针方向;<0,P在Q的逆时针方向;=0,P,Q共线,可能同向或反向    double operator ^(const Point &b)const    {        return x*b.y - b.x*y;    }}polygon[105];typedef Point Vector;//三态函数,判断两个double在eps精度下的大小关系int dcmp(double x){    if(fabs(x)<eps) return 0;    else        return x<0?-1:1;}//判断点Q是否在P1和P2的线段上bool OnSegment(Point P1,Point P2,Point Q){    return dcmp((P1-Q)^(P2-Q))==0&&dcmp((P1-Q)*(P2-Q))<=0;}//判断点P在多边形内-改进弧长法bool InPolygon(Point P){    int q1,q2,ans=0;    Point P1,P2;    Vector V1,V2;    for(int i=1,j=n;i<=n;j=i++)    {        P1 = polygon[j];        P2 = polygon[i];        V1 = P1-P;        V2 = P2-P;        if(OnSegment(P1,P2,P)) return true;        q1 = V1.x>0 ? (V1.y>0 ? 0:3) : (V1.y>0 ? 1:2);        q2 = V2.x>0 ? (V2.y>0 ? 0:3) : (V2.y>0 ? 1:2);        int g = (q2-q1+4)%4;        if(g==1) ans--; //在上一象限        if(g==3) ans++; //在下一象限        if(g==2) dcmp(V1^V2)>0 ? (ans-=2) : (ans+=2);    }    return ans!=0;}int main(){    while(~scanf("%d",&n))    {        for(int i=1;i<=n;i++) scanf("%lf %lf",&polygon[i].x,&polygon[i].y);        Point test;        scanf("%d",&m);        while(m--)        {            scanf("%lf %lf",&test.x,&test.y);            if(InPolygon(test)) printf("Yes\n");            else printf("No\n");        }    }    return 0;}

叉积(点线)判断法

时间复杂度:O(n) 适用范围:凸多边形

我觉得该算法只作为了解吧。
对于多边形(正向,即逆时针),如果一个点它的所有有向边的左边,那么这个点一定在多边形内部。利用叉积正好可以判断点与给定边的关系,即点是在边的左边右边还是边上。
这里说一下,(P叉乘Q)P^Q>0说明P在Q的顺时针方向,<0说明P在Q的逆时针方向,=0说明P和Q共线。

没有代码,我觉得只要一下就行了,不推荐。

面积法

时间复杂度:O(n)(但是时间应该会比之前的O(n)的长一点)
适用范围:凸多边形

了解即可,有精度要求,强调多边形点给出的方向(逆时针)。
算法思想:
如果点在多边形内部或者边上,那么点与多边形所有边组成的三角形面积和等于多边形面积。多边形的面积可以用叉积计算即连接坐标原点和各顶点形成向量,所有向量叉积的0.5的和即为多边形面积。不过计算面积是会有一定误差的,需要设置精度的误差范围。

没有代码。

二分

时间复杂度:O(logn) 适用范围:凸多边形
强掉多边形给出的方向。

这种方法以hrbust1429为例。 题目链接
这道题的题意是已知构成凸多边形A的n个点的坐标,和点集B的m个点的坐标,求这B的m个点是否都在凸多边形A内(严格内部,就是点不能在多边形边上)。

算法思想:
1. 选择多边形其中一个点为起点,连接其它点作射线。


2. 判断给定的点是否在所有射线包围的区域之内,即判断给定点是否在最左侧射线的左边,或者在最右侧射线的右边。
3. 如果在射线包围的区域之内,选择构成最两侧的射线的点为left和right,则mid = (left+right)/2,连接给顶点和起点作射线,判断该射线在mid点和起点的哪一边,不断循环,如此用二分法最后求出给定点所在的三角形区域,由此确定了除起点外的一条边。

4. 判断给定点在这条边的左方还是右方,由此判断给定点是否在三角形区域内,也就是是否在多边形内。

//判断点P在多边形内-O(logn) 顺时针给出bool InPolygon(Point P){    if(dcmp((polygon[n]-polygon[1])^(P-polygon[1]))<=0 || dcmp((polygon[2]-polygon[1])^(P-polygon[1]))>=0) //判断在不在区域外        return false;    int l=2,r=n,mid;    //找出所在三角形区域的左边的顶点    while(l<r)    {        mid = (l+r+1)>>1;        if(dcmp((polygon[mid]-polygon[1])^(P-polygon[1]))<=0) l=mid;        else r = mid-1;    }     //判断在多边形外或线上    if(dcmp((polygon[l+1]-polygon[l])^(P-polygon[l]))>=0) return false;    return true;}

hrbust1429

#include <bits/stdc++.h>using namespace std;typedef long long ll;const double eps = 1e-6;const double PI = acos(-1);const int maxn = 1e5+5;int n,m;struct Point{    double x,y;    Point(double x=0,double y=0):x(x),y(y){}    //向量+    Point operator +(const Point &b)const    {        return Point(x+b.x,y+b.y);    }    //向量-    Point operator -(const Point &b)const    {        return Point(x-b.x,y-b.y);    }    //点积    double operator *(const Point &b)const    {        return x*b.x + y*b.y;    }    //叉积    //P^Q>0,P在Q的顺时针方向;<0,P在Q的逆时针方向;=0,P,Q共线,可能同向或反向    double operator ^(const Point &b)const    {        return x*b.y - b.x*y;    }}polygon[maxn];typedef Point Vector;//三态函数,判断两个double在eps精度下的大小关系int dcmp(double x){    if(fabs(x)<eps) return 0;    else        return x<0?-1:1;}//判断点Q是否在P1和P2的线段上bool OnSegment(Point P1,Point P2,Point Q){    return dcmp((P1-Q)^(P2-Q))==0&&dcmp((P1-Q)*(P2-Q))<=0;}//判断点P在多边形内-O(logn) 顺时针给出bool InPolygon(Point P){    if(dcmp((polygon[n]-polygon[1])^(P-polygon[1]))<=0 || dcmp((polygon[2]-polygon[1])^(P-polygon[1]))>=0)        return false;    int l=2,r=n,mid;    while(l<r)    {        mid = (l+r+1)>>1;        if(dcmp((polygon[mid]-polygon[1])^(P-polygon[1]))<=0) l=mid;        else r = mid-1;    }    if(dcmp((polygon[l+1]-polygon[l])^(P-polygon[l]))>=0) return false;    return true;}int main(){    while(~scanf("%d",&n))    {        for(int i=1;i<=n;i++) scanf("%lf %lf",&polygon[i].x,&polygon[i].y);        Point test;        bool flag = true;        scanf("%d",&m);        while(m--)        {            scanf("%lf %lf",&test.x,&test.y);            if(!flag) continue;            if(!InPolygon(test)) flag = false;        }        if(flag) printf("YES\n");        else printf("NO\n");    }    return 0;}

最后总结一下O(n)算法的选择:个人感觉最优先考虑射线,其次是转角,其它的能不用最好不用。

如果对你有所帮助请给个赞,没有也请不要踩。。。QAQ

参考博客:
二分法参考博客

转角法参考博客

面积法 改进弧长法参考博客

原创粉丝点击