判斷點在多邊形內(nèi)的多種寫法

1、判斷點在多邊形內(nèi)的多種寫法(射線算法)  (2010-10-09 17:04:24)轉(zhuǎn)載標簽： 計算幾何 射線法 雜談分類： 經(jīng)驗總結(jié)*                射線算法一               *1.  &

2、#160;      已知點point(x,y)和多邊形Polygon（x1,y1;x2,y2;.xn,yn;）；2.         以point為起點，以無窮遠為終點作平行于X軸的直線line(x,y; -,y)；3.         循環(huán)取得(for(i=0;i<n;i+)多邊形的每一條邊side(xi,yi;xi+1,yi+1)，且判斷是否平行

3、于X軸，如果平行continue，否則，i+；4.         同時判斷point(x,y)是否在side上，如果是，則返回1(點在多邊形上)，否則繼續(xù)下面的判斷；5.         判斷線side與line是否有交點，如果有則count+，否則，i+。6.         判斷交點的總數(shù)，如果為奇數(shù)則返回0（點在多邊形內(nèi)

4、），偶數(shù)則返回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) );    / 判斷點在多邊形內(nèi)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;           *本文是采用射線法判斷點是否在多邊形內(nèi)的C語言程序。多年前，我自己實現(xiàn)了這樣一個算法。但是隨著時間的推移，我決定重寫這個代碼。參考周培德的計算幾何一書，結(jié)合我的實踐和經(jīng)驗，我相信，在這個算法的實現(xiàn)上，這是你迄今為止遇到的最優(yōu)的代碼。這是個C語言的小算法的實現(xiàn)程序，本來不想放到這里?？墒牵斘易约阂獙崿F(xiàn)這樣一個算法的時候，想在網(wǎng)上找個現(xiàn)成的，考察下來竟然一個符合需要的也沒有。我對自己大學(xué)讀書時寫的代碼沒有信心，所以，決定重新寫一個，并把它放到這里，以饗讀者

20、。也增加一下BLOG的點擊量。首先定義點結(jié)構(gòu)如下： 以下是引用片段：typedef structdouble x, y; vertex_t;本算法里所指的多邊形，是指由一系列點序列組成的封閉簡單多邊形。它的首尾點可以是或不是同一個點（不強制要求首尾點是同一個點）。這樣的多邊形可以是任意形狀的，包括多條邊在一條絕對直線上。因此，定義多邊形結(jié)構(gòu)如下： 以下是引用片段：typedef structint num_vertices;vertex_t *vertex; vertexlist_t;為加快判別速度，首先計算多邊形的外包矩形（rect_t），判斷點是否落在外包矩形內(nèi)，只有滿足落在外包矩形內(nèi)的條

21、件的點，才進入下一步的計算。為此，引入外包矩形結(jié)構(gòu)rect_t和求點集合的外包矩形內(nèi)的方法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;當點滿足落在多邊形外包矩形內(nèi)的條件，要進一步判斷點（v）是否在多邊形（

23、vl：np）內(nèi)。本程序采用射線法，由待測試點（v）水平引出一條射線B（v，w），計算B與vl邊線的交點數(shù)目，記為c，根據(jù)奇內(nèi)偶外原則（c為奇數(shù)說明v在vl內(nèi)，否則v不在vl內(nèi)）判斷點是否在多邊形內(nèi)。具體原理就不多說。為計算線段間是否存在交點，引入下面的函數(shù)：（1）is_same判斷2（p、q）個點是（1）否（0）在直線l（l_start，l_end）的同側(cè)；（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;下面的函數(shù)pt_in_poly就是判斷點（v）是（1）否（0）在多邊形（vl：np）內(nèi)的程序： 以下是

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;判斷點在多邊形內(nèi)的多種寫法(射線算法)  (2010-10-09 17:04:24)轉(zhuǎn)載標簽： 計算幾何 射線法 雜談分類： 經(jīng)驗總結(jié)*                射線算法一           &

29、#160;   *1.         已知點point(x,y)和多邊形Polygon（x1,y1;x2,y2;.xn,yn;）；2.         以point為起點，以無窮遠為終點作平行于X軸的直線line(x,y; -,y)；3.         循環(huán)取得(for(i=0;i<n;i

30、+)多邊形的每一條邊side(xi,yi;xi+1,yi+1)，且判斷是否平行于X軸，如果平行continue，否則，i+；4.         同時判斷point(x,y)是否在side上，如果是，則返回1(點在多邊形上)，否則繼續(xù)下面的判斷；5.         判斷線side與line是否有交點，如果有則count+，否則，i+。6.      

31、60;  判斷交點的總數(shù)，如果為奇數(shù)則返回0（點在多邊形內(nèi)），偶數(shù)則返回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、;判斷點在多邊形內(nèi)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;  射線算法二               *本文是采用射線法判斷點是否在多邊形內(nèi)的C語言程序。多年前，我自己實現(xiàn)了這樣一個算法。但是隨著時間的推移，我決定重寫這個代碼。參考周培德的計算幾何一書，結(jié)合我的實踐和經(jīng)驗，我相信，在這個算法的實現(xiàn)上，這是你迄今為止遇到的最優(yōu)的代碼。這是個C語言的小算法的實現(xiàn)程序，本來不想放到這里?？墒牵斘易约阂獙崿F(xiàn)這樣一個算法的時候，想在網(wǎng)上找個現(xiàn)成的，考察下來竟然一個符合需要的也沒有。我對

47、自己大學(xué)讀書時寫的代碼沒有信心，所以，決定重新寫一個，并把它放到這里，以饗讀者。也增加一下BLOG的點擊量。首先定義點結(jié)構(gòu)如下： 以下是引用片段：typedef structdouble x, y; vertex_t;本算法里所指的多邊形，是指由一系列點序列組成的封閉簡單多邊形。它的首尾點可以是或不是同一個點（不強制要求首尾點是同一個點）。這樣的多邊形可以是任意形狀的，包括多條邊在一條絕對直線上。因此，定義多邊形結(jié)構(gòu)如下： 以下是引用片段：typedef structint num_vertices;vertex_t *vertex; vertexlist_t;為加快判別速度，首先計算多邊形的

48、外包矩形（rect_t），判斷點是否落在外包矩形內(nèi)，只有滿足落在外包矩形內(nèi)的條件的點，才進入下一步的計算。為此，引入外包矩形結(jié)構(gòu)rect_t和求點集合的外包矩形內(nèi)的方法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;當點滿足落在多邊形外包矩形內(nèi)的條件，要進一步判斷點（v）是否在多邊形（vl：np）內(nèi)。本程序采用射線法，由待測試點（v）水平引出一條射線B（v，w），計算B與vl邊線的交點數(shù)目，記為c，根據(jù)奇內(nèi)偶外原則（c為奇數(shù)說明v在vl內(nèi)，否則v不在vl內(nèi)）判斷點是否在多邊形內(nèi)。具體原理就不多說。為計算線段間是否存在交點，引入下面的函數(shù)：（1）is_same判斷2（p、q）個點是（1）否（0）在直線l（l_start，l_end）的同側(cè)；（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;下面的函數(shù)pt_in_p

53、oly就是判斷點（v）是（1）否（0）在多邊形（vl：np）內(nèi)的程序： 以下是引用片段：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軟件開發(fā)中，經(jīng)常要用到一些幾何的算法，比如三角網(wǎng)構(gòu)建，多邊形的剖分，點，線，面之間的關(guān)系。而點與多邊形關(guān)系的判斷是一項非常重要的基礎(chǔ)工作。在點與多邊形關(guān)系的判斷中，經(jīng)常用到的方法是射線法和夾角和方法，其中射線法能夠針對帶島的多邊形進行判斷，而夾角和方法就顯得無能為力。射線法的基本思想是：從待判斷的點向某一個方向引射線，計算和多邊形交點的個數(shù)，如果個數(shù)是偶數(shù)或者0則點在多邊形外，如果是奇數(shù)，則在多邊形內(nèi)。這個只是最基本的判別情況，還有

56、一些復(fù)雜的情況需要特殊處理：（射線經(jīng)過頂點）：當射線經(jīng)過頂點時，判斷就會出現(xiàn)異常情況，現(xiàn)在規(guī)定，線段的兩個端點，相對于另一個端點在上面的頂點稱為上端點，下面是下端點，如果經(jīng)過下端點，則認為邊和射線不相交。（點在邊上）：這種情況也不能用交點個數(shù)的奇偶性來判斷了，要快速地判斷這個點是否在邊上。射線法改進：傳統(tǒng)的射線法一開始就直接計算點和多邊形的交點個數(shù)，這樣的話，會花費大量的時間來作拓撲關(guān)系的判斷。改進的算法是首先利用多邊形的最小外接矩形迅速排出掉不在MBR內(nèi)的點，然后利用交點個數(shù)的奇偶性判斷：下面的函數(shù)是射線和邊關(guān)系以及交點個數(shù)判斷：cpp view plaincopy1. /射線和線

57、段的關(guān)系 :相交返回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.     /交點在射線右側(cè)，則不相交  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.     /穿過下端點也不計數(shù)  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.     /如果點不在多邊形的最小外接矩形中，則一定不在多邊形內(nèi)  4.     MyEnvelope env; &#