判断点在多边形内的多种写法

1、判断点在多边形内的多种写法(射线算法)  (2010-10-09 17:04:24)转载标签： 计算几何 射线法 杂谈分类： 经验总结*                射线算法一               *1.  &

2、#160;      已知点point(x,y)和多边形Polygon（x1,y1;x2,y2;.xn,yn;）；2.         以point为起点，以无穷远为终点作平行于X轴的直线line(x,y; -,y)；3.         循环取得(for(i=0;i<n;i+)多边形的每一条边side(xi,yi;xi+1,yi+1)，且判断是否平行

3、于X轴，如果平行continue，否则，i+；4.         同时判断point(x,y)是否在side上，如果是，则返回1(点在多边形上)，否则继续下面的判断；5.         判断线side与line是否有交点，如果有则count+，否则，i+。6.         判断交点的总数，如果为奇数则返回0（点在多边形内

4、），偶数则返回2（点在多边形外）。 代码： const double INFINITY = 1e10;const double ESP = 1e-5;const int MAX_N = 1000; struct Point          double x, y;     struct LineSegment         

5、Point pt1, pt2;    typedef vector<Point> Polygon; / 计算叉乘 |P0P1| × |P0P2|double Multiply(Point p1, Point p2, Point p0)       return ( (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y) );  

6、  / 判断线段是否包含点pointbool IsOnline(Point point, LineSegment line)         return( ( fabs(Multiply(line.pt1, line.pt2, point) < ESP ) &&     ( ( point.x - line.pt1.x ) * ( point.x - line.pt2.x ) <= 0 ) &am

7、p;&     ( ( point.y - line.pt1.y ) * ( point.y - line.pt2.y ) <= 0 ) );    / 判断线段相交bool Intersect(LineSegment L1, LineSegment L2)         return( (max(L1.pt1.x, L1.pt2.x) >= min(L2.pt1.x, L2.pt2.

8、x) &&     (max(L2.pt1.x, L2.pt2.x) >= min(L1.pt1.x, L1.pt2.x) &&     (max(L1.pt1.y, L1.pt2.y) >= min(L2.pt1.y, L2.pt2.y) &&     (max(L2.pt1.y, L2.pt2.y) >= min(L1.pt1.y, L1.pt2.y) &&&#

9、160;    (Multiply(L2.pt1, L1.pt2, L1.pt1) * Multiply(L1.pt2, L2.pt2, L1.pt1) >= 0) &&     (Multiply(L1.pt1, L2.pt2, L2.pt1) * Multiply(L2.pt2, L1.pt2, L2.pt1) >= 0) );    / 判断点在多边形内bool InPolygon(const Polygon&a

10、mp; polygon, Point point)         int n = polygon.size();     int count = 0;     LineSegment line;     line.pt1 = point;     line.pt2.y = point.y; 

11、60;   line.pt2.x = - INFINITY;      for( int i = 0; i < n; i+ )        / 得到多边形的一条边           LineSegment side;       

12、   side.pt1 = polygoni;          side.pt2 = polygon(i + 1) % n;           if( IsOnline(point, side) )             &

13、#160;            return1 ;              / 如果side平行x轴则不作考虑          if( fabs(side.pt1.y - side.pt2.y) < ESP )

14、0;                        continue;                       if( IsO

15、nline(side.pt1, line) )                         if( side.pt1.y > side.pt2.y ) count+;              

16、60;      else if( IsOnline(side.pt2, line) )                       if( side.pt2.y > side.pt1.y ) count+;        

17、             else if( Intersect(line, side) )                       count+;       

18、60;             if ( count % 2 = 1 ) return 0;else return 2;         *                射线算法二   &

19、#160;           *本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的计算几何一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者

20、。也增加一下BLOG的点击量。首先定义点结构如下： 以下是引用片段：typedef structdouble x, y; vertex_t;本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点（不强制要求首尾点是同一个点）。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下： 以下是引用片段：typedef structint num_vertices;vertex_t *vertex; vertexlist_t;为加快判别速度，首先计算多边形的外包矩形（rect_t），判断点是否落在外包矩形内，只有满足落在外包矩形内的条

21、件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下： 以下是引用片段：typedef structdouble min_x, min_y, max_x, max_y; rect_t;void vertices_get_extent (const vertex_t* vl, int np,rect_t* rc )int i;if (np > 0)rc->min_x = rc->max_x = vl0.x; rc->min_y = rc->max_y = vl0.y;elserc-

22、>min_x = rc->min_y = rc->max_x = rc->max_y = 0;for(i=1; iif(vli.x < rc->min_x) rc->min_x = vli.x;if(vli.y < rc->min_y) rc->min_y = vli.y;if(vli.x > rc->max_x) rc->max_x = vli.x;if(vli.y > rc->max_y) rc->max_y = vli.y;当点满足落在多边形外包矩形内的条件，要进一步判断点（v）是否在多边形（

23、vl：np）内。本程序采用射线法，由待测试点（v）水平引出一条射线B（v，w），计算B与vl边线的交点数目，记为c，根据奇内偶外原则（c为奇数说明v在vl内，否则v不在vl内）判断点是否在多边形内。具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：（1）is_same判断2（p、q）个点是（1）否（0）在直线l（l_start，l_end）的同侧；（2）is_intersect用来判断2条线段（不是直线）s1、s2是（1）否（0）相交； 以下是引用片段：static int is_same(const vertex_t* l_start, const vertex_t* l_end,

24、const vertex_t* p,const vertex_t* q)double dx = l_end->x - l_start->x;double dy = l_end->y - l_start->y;double dx1= p->x - l_start->x;double dy1= p->y - l_start->y;double dx2= q->x - l_end->x;double dy2= q->y - l_end->y;return (dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) >

25、0? 1 : 0);static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,const vertex_t* s2_start, const vertex_t* s2_end)return (is_same(s1_start, s1_end, s2_start, s2_end)=0 &&is_same(s2_start, s2_end, s1_start, s1_end)=0)? 1: 0;下面的函数pt_in_poly就是判断点（v）是（1）否（0）在多边形（vl：np）内的程序： 以下是

26、引用片段：int pt_in_poly ( const vertex_t* vl, int np,const vertex_t* v)int i, j, k1, k2, c;rect_t rc;vertex_t w;if (np < 3)return 0;vertices_get_extent(vl, np, &rc);if (v->x < rc.min_x | v->x > rc.max_x | v->y < rc.min_y | v->y > rc.max_y)return 0;w.x = rc.max_x + DBL_EPSI

27、LON;w.y = v->y;c = 0;for(i=0; ij = (i+1) % np;if(is_intersect(vl+i, vl+j, v, &w)C+;else if(vli.y=w.y)k1 = (np+i-1)%np;while(k1!=i && vlk1.y=w.y)k1 = (np+k1-1)%np;k2 = (i+1)%np;while(k2!=i && vlk2.y=w.y)k2 = (k2+1)%np;if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)=0)c+

28、;if(k2 <= i)break;i = k2;return c%2;判断点在多边形内的多种写法(射线算法)  (2010-10-09 17:04:24)转载标签： 计算几何 射线法 杂谈分类： 经验总结*                射线算法一           &

29、#160;   *1.         已知点point(x,y)和多边形Polygon（x1,y1;x2,y2;.xn,yn;）；2.         以point为起点，以无穷远为终点作平行于X轴的直线line(x,y; -,y)；3.         循环取得(for(i=0;i<n;i

30、+)多边形的每一条边side(xi,yi;xi+1,yi+1)，且判断是否平行于X轴，如果平行continue，否则，i+；4.         同时判断point(x,y)是否在side上，如果是，则返回1(点在多边形上)，否则继续下面的判断；5.         判断线side与line是否有交点，如果有则count+，否则，i+。6.      

31、60;  判断交点的总数，如果为奇数则返回0（点在多边形内），偶数则返回2（点在多边形外）。 代码： const double INFINITY = 1e10;const double ESP = 1e-5;const int MAX_N = 1000; struct Point          double x, y;     struct LineSegment  

32、60;      Point pt1, pt2;    typedef vector<Point> Polygon; / 计算叉乘 |P0P1| × |P0P2|double Multiply(Point p1, Point p2, Point p0)       return ( (p1.x - p0.x) * (p2.y - p0.y) - (p2.x

33、 - p0.x) * (p1.y - p0.y) );    / 判断线段是否包含点pointbool IsOnline(Point point, LineSegment line)         return( ( fabs(Multiply(line.pt1, line.pt2, point) < ESP ) &&     ( ( point.x - line.pt1.x ) *

34、 ( point.x - line.pt2.x ) <= 0 ) &&     ( ( point.y - line.pt1.y ) * ( point.y - line.pt2.y ) <= 0 ) );    / 判断线段相交bool Intersect(LineSegment L1, LineSegment L2)         return( (max(L1.pt1.x

35、, L1.pt2.x) >= min(L2.pt1.x, L2.pt2.x) &&     (max(L2.pt1.x, L2.pt2.x) >= min(L1.pt1.x, L1.pt2.x) &&     (max(L1.pt1.y, L1.pt2.y) >= min(L2.pt1.y, L2.pt2.y) &&     (max(L2.pt1.y, L2.pt2.y) >

36、;= min(L1.pt1.y, L1.pt2.y) &&     (Multiply(L2.pt1, L1.pt2, L1.pt1) * Multiply(L1.pt2, L2.pt2, L1.pt1) >= 0) &&     (Multiply(L1.pt1, L2.pt2, L2.pt1) * Multiply(L2.pt2, L1.pt2, L2.pt1) >= 0) );    / 

37、;判断点在多边形内bool InPolygon(const Polygon& polygon, Point point)         int n = polygon.size();     int count = 0;     LineSegment line;     line.pt1 = point;   

38、60; line.pt2.y = point.y;     line.pt2.x = - INFINITY;      for( int i = 0; i < n; i+ )        / 得到多边形的一条边           LineSegment side;

39、60;         side.pt1 = polygoni;          side.pt2 = polygon(i + 1) % n;           if( IsOnline(point, side) )      

40、0;                   return1 ;              / 如果side平行x轴则不作考虑          if( fabs(

41、side.pt1.y - side.pt2.y) < ESP )                         continue;                 

42、0;     if( IsOnline(side.pt1, line) )                         if( side.pt1.y > side.pt2.y ) count+;        

43、             else if( IsOnline(side.pt2, line) )                       if( side.pt2.y > side.pt1.y ) count+; 

44、60;                   else if( Intersect(line, side) )                       count+; 

45、                    if ( count % 2 = 1 ) return 0;else return 2;         *             

46、60;  射线算法二               *本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的计算几何一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对

47、自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。首先定义点结构如下： 以下是引用片段：typedef structdouble x, y; vertex_t;本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点（不强制要求首尾点是同一个点）。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下： 以下是引用片段：typedef structint num_vertices;vertex_t *vertex; vertexlist_t;为加快判别速度，首先计算多边形的

48、外包矩形（rect_t），判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下： 以下是引用片段：typedef structdouble min_x, min_y, max_x, max_y; rect_t;void vertices_get_extent (const vertex_t* vl, int np,rect_t* rc )int i;if (np > 0)rc->min_x = rc->max_x = vl0.x; rc-&

49、gt;min_y = rc->max_y = vl0.y;elserc->min_x = rc->min_y = rc->max_x = rc->max_y = 0;for(i=1; iif(vli.x < rc->min_x) rc->min_x = vli.x;if(vli.y < rc->min_y) rc->min_y = vli.y;if(vli.x > rc->max_x) rc->max_x = vli.x;if(vli.y > rc->max_y) rc->max_y = vl

50、i.y;当点满足落在多边形外包矩形内的条件，要进一步判断点（v）是否在多边形（vl：np）内。本程序采用射线法，由待测试点（v）水平引出一条射线B（v，w），计算B与vl边线的交点数目，记为c，根据奇内偶外原则（c为奇数说明v在vl内，否则v不在vl内）判断点是否在多边形内。具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：（1）is_same判断2（p、q）个点是（1）否（0）在直线l（l_start，l_end）的同侧；（2）is_intersect用来判断2条线段（不是直线）s1、s2是（1）否（0）相交； 以下是引用片段：static int is_same(const ve

51、rtex_t* l_start, const vertex_t* l_end,const vertex_t* p,const vertex_t* q)double dx = l_end->x - l_start->x;double dy = l_end->y - l_start->y;double dx1= p->x - l_start->x;double dy1= p->y - l_start->y;double dx2= q->x - l_end->x;double dy2= q->y - l_end->y;retur

52、n (dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,const vertex_t* s2_start, const vertex_t* s2_end)return (is_same(s1_start, s1_end, s2_start, s2_end)=0 &&is_same(s2_start, s2_end, s1_start, s1_end)=0)? 1: 0;下面的函数pt_in_p

53、oly就是判断点（v）是（1）否（0）在多边形（vl：np）内的程序： 以下是引用片段：int pt_in_poly ( const vertex_t* vl, int np,const vertex_t* v)int i, j, k1, k2, c;rect_t rc;vertex_t w;if (np < 3)return 0;vertices_get_extent(vl, np, &rc);if (v->x < rc.min_x | v->x > rc.max_x | v->y < rc.min_y | v->y > rc.m

54、ax_y)return 0;w.x = rc.max_x + DBL_EPSILON;w.y = v->y;c = 0;for(i=0; ij = (i+1) % np;if(is_intersect(vl+i, vl+j, v, &w)C+;else if(vli.y=w.y)k1 = (np+i-1)%np;while(k1!=i && vlk1.y=w.y)k1 = (np+k1-1)%np;k2 = (i+1)%np;while(k2!=i && vlk2.y=w.y)k2 = (k2+1)%np;if(k1 != k2 &&

55、; is_same(v, &w, vl+k1, vl+k2)=0)c+;if(k2 <= i)break;i = k2;return c%2;在GIS软件开发中，经常要用到一些几何的算法，比如三角网构建，多边形的剖分，点，线，面之间的关系。而点与多边形关系的判断是一项非常重要的基础工作。在点与多边形关系的判断中，经常用到的方法是射线法和夹角和方法，其中射线法能够针对带岛的多边形进行判断，而夹角和方法就显得无能为力。射线法的基本思想是：从待判断的点向某一个方向引射线，计算和多边形交点的个数，如果个数是偶数或者0则点在多边形外，如果是奇数，则在多边形内。这个只是最基本的判别情况，还有

56、一些复杂的情况需要特殊处理：（射线经过顶点）：当射线经过顶点时，判断就会出现异常情况，现在规定，线段的两个端点，相对于另一个端点在上面的顶点称为上端点，下面是下端点，如果经过下端点，则认为边和射线不相交。（点在边上）：这种情况也不能用交点个数的奇偶性来判断了，要快速地判断这个点是否在边上。射线法改进：传统的射线法一开始就直接计算点和多边形的交点个数，这样的话，会花费大量的时间来作拓扑关系的判断。改进的算法是首先利用多边形的最小外接矩形迅速排出掉不在MBR内的点，然后利用交点个数的奇偶性判断：下面的函数是射线和边关系以及交点个数判断：cpp view plaincopy1. /射线和线

57、段的关系 :相交返回1，不相交返回0，射线起点在线段上返回-1  2. int IsIntersectAnt(double x,double y,double X1,double Y1,double X2,double Y2)  3.   4.     /计算线段的最小和最大坐标值  5.     double minX,maxX,minY,max

58、Y;  6.     minX = X1;  7.     maxX = X2;  8.     if (minX > maxX)  9.       10.         minX

59、 = X2;  11.         maxX = X1;  12.       13.     minY = Y1;  14.     maxY = Y2;  15.     

60、if (minY > maxY)  16.       17.         minY = Y2;  18.         maxY = Y1;  19.       20. &#

61、160; 21.     /射线与边无交点的快速判断  22.     if (y<minY | y>maxY | x<minX)  23.       24.         return 0;  25.   

62、    26.   27.     /如果是水平线段，在线段上返回-1，否则返回0  28.     if (fabs(maxY - minY) < eps)  29.       30.         return

63、60;(x >= minX && x <= maxX)? (-1):0;  31.       32.   33.     /计算射线与边所在直线的交点的横坐标  34.     double x0 = X1 + (double)(y -&#

64、160;Y1)*(X2 - X1)/(Y2 - Y1);  35.       36.     /交点在射线右侧，则不相交  37.     if (x0 > x)  38.       39.      

65、0;  return 0;  40.       41.     /交点和射线起点相同  42.     if (fabs(x-x0)< eps)  43.       44.         retur

66、n -1;  45.       46.     /穿过下端点也不计数  47.     if (fabs(y-minY) < eps)  48.       49.         return 0;&

67、#160; 50.       51.     return 1;  52.   53.   上面的是计算射线和一条边的情况，对于一个多边形所以只要逐个判断射线和它的边就可以了cpp view plaincopy1. int MyPolygon:PointInPolygon(const MyPoint& poPoint)  2.   3.     /如果点不在多边形的最小外接矩形中，则一定不在多边形内  4.     MyEnvelope env; &#