GIS点在多边形内算法.doc

1.判断点在多边形内外的简单算法 - 改进弧长法今天学图形学的时候发现了一个求多边形内外的超简单算法，当时觉得非常惊喜，后来晚上上完选修的时候在走廊遇到bug，bug也是很惊喜地感慨道：居然有甘简单既办法都捻唔到！遂将其写下，供大家分享，希望不会太火星。 这个算法是源自计算机图形学基础教程（孙家广，清华大学出版社），在该书的48-49页，名字可称为“改进的弧长法”。该算法只需O(1)的附加空间，时间复杂度为O(n)，但系数很小；最大的优点是具有很高的精度，只需做乘法和减法，若针对整数坐标则完全没有精度问题。而且实现起来也非常简单，比转角法和射线法都要好写且不易出错。 首先从该收中摘抄一段弧长法的介绍：“弧长法要求多边形是有向多边形，一般规定沿多边形的正向，边的左侧为多边形的内侧域。以被测点为圆心作单位圆，将全部有向边向单位圆作径向投影，并计算其中单位圆上弧长的代数和。若代数和为0，则点在多边形外部；若代数和为2则点在多边形内部；若代数和为，则点在多边形上。” 按书上的这个介绍，其实弧长法就是转角法。但它的改进方法比较厉害：将坐标原点平移到被测点P，这个新坐标系将平面划分为4个象限，对每个多边形顶点P ，只考虑其所在的象限，然后按邻接顺序访问多边形的各个顶点P，分析P和Pi+1，有下列三种情况：（1）Pi+1在P的下一象限。此时弧长和加/2；（2）Pi+1在P的上一象限。此时弧长和减/2；（3）Pi+1在Pi的相对象限。首先计算f=yi+1*x-xi+1*y（叉积），若f=0，则点在多边形上；若f0，弧长和加。 最后对算出的代数和和上述的情况一样判断即可。 实现的时候还有两点要注意，第一个是若P的某个坐标为0时，一律当正号处理；第二点是若被测点和多边形的顶点重合时要特殊处理。 以上就是书上讲解的内容，其实还存在一个问题。那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错。如边(3,0)-(-3,0)，按以上的处理，象限分别是第一和第二，这样会使代数和加/2，有可能导致最后结果是被测点在多边形外。而实际上被测点是在多边形上（该边穿过该点）。 对于这点，我的处理办法是：每次算P和Pi+1时，就计算叉积和点积，判断该点是否在该边上，是则判断结束，否则继续上述过程。这样牺牲了时间，但保证了正确性。 具体实现的时候，由于只需知道当前点和上一点的象限位置，所以附加空间只需O(1)。实现的时候可以把上述的“/2”改成1，“”改成2，这样便可以完全使用整数进行计算。不必考虑顶点的顺序，逆时针和顺时针都可以处理，只是最后的代数和符号不同而已。整个算法编写起来非常容易。 针对以上算法，我写了一个代码，拿ZOJ 1081 Points Within进行测试，顺利Accepted。这证明该算法的正确性还是可以保障的。以下附上我的代码：/ ZOJ 1081 , 改进弧长法判点在形内形外#include #include const int MAX = 101 ;struct point int x , y ; pMAX ;int main() int n , m , i , sum , t1 , t2 , f , prob = 0 ; point t ; while ( scanf( %d , &n ) , n ) if( prob + ) printf ( n ); printf ( Problem %d:n , prob ) ; scanf ( %d , &m ) ; for ( i = 0 ; i n ; i + ) scanf ( %d%d , &p.x , &p.y ) ; pn = p0 ; while ( m - ) scanf ( %d%d , &t.x , &t.y ); for ( i = 0 ; i =0 ? ( p0.y=0?0:3 ) : ( p0.y=0?1:2 ); / 计算象限 for ( sum = 0 , i = 1 ; i = n ; i + ) if ( !p.x & !p.y ) break ;/ 被测点为多边形顶点 f = p.y * pi-1.x - p.x * pi-1.y ; / 计算叉积 if ( !f & pi-1.x*p.x = 0 & pi-1.y*pi.y =0 ? ( p.y=0?0:3 ) : ( p.y=0?1:2 ) ; / 计算象限 if ( t2 = ( t1 + 1 ) % 4 ) sum += 1 ; / 情况1 else if ( t2 = ( t1 + 3 ) % 4 ) sum -= 1 ; / 情况2 else if ( t2 = ( t1 + 2 ) % 4 ) / 情况3 if ( f 0 ) sum += 2 ; else sum -= 2; t1 = t2 ; if ( i=n | sum ) printf( Withinn ) ; else printf(Outsiden ) ; for ( i = 0 ; i 0) rc-min_x=rc-max_x=vl0.x;rc-min_y=rc-max_y=vl0.y; else rc-min_x=rc-min_y=rc-max_x=rc-max_y=0;/*=0?noverticesatall*/ for(i=1;i if(vli.xmin_x)rc-min_x=vli.x; if(vli.ymin_y)rc-min_y=vli.y; if(vli.xrc-max_x)rc-max_x=vli.x; if(vli.yrc-max_y)rc-max_y=vli.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)相交; 以下是引用片段： Copycode/*p,qisonthesameoflinel*/ staticintis_same(constvertex_t*l_start,constvertex_t*l_end,/*linel*/ constvertex_t*p, constvertex_t*q) doubledx=l_end-x-l_start-x; doubledy=l_end-y-l_start-y; doubledx1=p-x-l_start-x; doubledy1=p-y-l_start-y; doubledx2=q-x-l_end-x; doubledy2=q-y-l_end-y; return(dx*dy1-dy*dx1)*(dx*dy2-dy*dx2)0?1:0); /*2linesegments(s1,s2)areintersect?*/ staticintis_intersect(constvertex_t*s1_start,constvertex_t*s1_end, constvertex_t*s2_start,constvertex_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)内的程序： 以下是引用片段： Copycodeintpt_in_poly(constvertex_t*vl,intnp,/*polygonvlwithnpvertices*/ constvertex_t*v) inti,j,k1,k2,c; rect_trc; vertex_tw; if(npxxrc.max_x|v-yyrc.max_y) return0; /*Setahorizontalbeaml(*v,w)fromvtotheultraright*/ w.x=rc.max_x+DBL_EPSILON; w.y=v-y; c=0;/*Intersectionpointscounter*/ for(i=0;i j=(i+1)%np; if(is_intersect(vl+i,vl+j,v,&w) c+; elseif(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+; if(k2InPoly: q = ); PrintPoint(q); putchar(n); /* Shift so that q is the origin. Note this destroys the polygon. This is done for pedogical clarity. */ for( i = 0; i n; i+ ) for( d = 0; d DIM; d+ ) Pid = Pid - qd; /* For each edge e=(i-1,i), see if crosses ray. */ for( i = 0; i 0 ) & ( Pi1Y 0 ) & ( Pi Y 0 ) != ( Pi1Y 0 ) ) /* e straddles ray, so compute intersection with ray. */ x = (PiX * (double)Pi1Y - Pi1X * (double)PiY) / (double)(Pi1Y - PiY); /* printf(straddles: x = %gt, x); */ /* crosses ray if strictly positive intersection. */ if (x 0) Rcross+; /* printf(Right cross=%dt, Rcross); */ /* if e straddles the x-axis when reversed. */ /* if( ( ( Pi Y = 0 ) ) | ( ( Pi1Y = 0 ) ) ) */ if ( ( PiY 0 ) != ( Pi1Y 0 ) ) /* e straddles ray, so compute intersection with ray. */ x = (PiX * Pi1Y - Pi1X * PiY) / (double)(Pi1Y - PiY); /* printf(straddles: x = %gt, x); */ /* crosses ray if strictly positive intersection. */ if (x 0) Lcross+; /* printf(Left cross=%dn, Lcross); */ /* q on the edge if left and right cross are not the same parity. */ if( ( Rcross % 2 ) != (Lcross % 2 ) ) return e; /* q inside iff an odd number of crossings. */ if( (Rcross % 2) = 1 ) return i; else return o; void PrintPoint( tPointi p ) int i; putchar(); for ( i = 0; i DIM; i+ ) printf(%d, pi); if ( i != DIM-1 ) putchar(,); putchar(); /* Reads in the coordinates of the vertices of a polygon from stdin, puts them into P, and returns n, the number of vertices. Formatting conventions: etc. */ int ReadPoly( tPolygoni P ) int i, n; do printf( Input the number of vertices:n); scanf( %d, &n ); if ( n = PMAX ) break; printf(Error in read_poly: too many points; max is %dn, PMAX); while ( 1 ); printf( Polygon:n ); printf( i x yn); for ( i = 0; i n; i+ ) scanf( %d %d, &Pi0, &Pi1 ); printf(%3d%4d%4dn, i, Pi0, Pi1); printf(n = %3d vertices readn,n); putchar(n); return n; void PrintPoly( int n, tPolygoni P ) int i; printf(Polygon:n); printf( i x yn); for( i = 0; i n; i+ ) printf(%3d%4d%4dn, i, Pi0, Pi1); /* Copy polygon a to b (overwriting b). */ void Copy( tPolygoni a, tPolygoni b, int n ) int i, j; for ( i=0; i n; i+) for ( j = 0; j DIM; j+ ) bij = aij; bool EqPoint( tPointi a, tPointi b ) int i; for ( i = 0; i DIM; i+ ) if ( ai != bi) return FALSE; return TRUE; 1.已知点point(x,y)和多边形Polygon（x1,y1;x2,y2;.xn,yn;）；2.以point为起点，以无穷远为终点作平行于X轴的直线line(x,y; -,y)；3.循环取得(for(i=0;in;i+)多边形的每一条边side(xi,yi;xi+1,yi+1)，且判断是否平行于X轴，如果平行continue，否则，i+；4.同时判断point(x,y)是否在side上，如果是，则返回1(点在多边形上)，否则继续下面的判断；5.判断线side与line是否有交点，如果有则count+，否则，i+。6.判断交点的总数，如果为奇数则返回0（点在多边形内），偶数则返回2（点在多边形外）。代码：/*射线法判断点q与多边形polygon的位置关系，要求polygon为简单多边形，顶点逆时针排列如果点在多边形内：返回0如果点在多边形边上：返回1如果点在多边形外：返回2*/const double INFINITY = 1e10;const double ESP = 1e-5;const int MAX_N = 1000;struct Point double x, y;struct LineSegment Point pt1, pt2;typedef vector 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) );/判断线段是否包含点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 ) &( ( point.y - line.pt1.y ) * ( point.y - line.pt2.y ) = 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) = 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);/判断点在多边形内bool InPolygon(const Polygon& polygon, Point point)int n = polygon.size();int count = 0;LineSegment line;line.pt1 = point;line.pt2.y = point.y;line.pt2.x = - INFINITY;for( int i = 0; i n; i+ ) /得到多边形的一条边LineSegment side;side.pt1 = polygoni;side.pt2 = polygon(i + 1) % n;if( IsOnline(point, side) ) return1 ;/如果side平行x轴则不作考虑if( fabs(side.pt1.y - side.pt2.y) side.pt2.y ) count+; else if( IsOnline(side.pt2, line) ) if( side.pt2.y side.pt1.y ) count+; else if( Intersect(line, side) ) count+;if ( count % 2 = 1 ) return 0;else return 2;4. 用射线相交来判断点是否在多边形内的算法点在边上，射线通过顶点，射线通过边的情况都考虑了：Public Lx(20) As DoublePublic Ly(20) As DoublePublic Ln As Integer定义多边形的坐标，边数Public Px As DoublePublic Py As Double定义点的坐标该过程判断点是否在多边形内，点在边上判断为在多边形内，方法为射线相交法，返回1为在多边形内，0为在多边形外Public Function pinshp()Dim lc(20) As IntegerDim ls(20) As Integerpinshp = 0inline = 0For z = 0 To Lnzn = z + 1If zn Ln Then zn = 0点在边上判断If Ly(zn) = Ly(z) And Py = Ly(z) ThenIf Px Lx(z) And Px Lx(zn) And Px Ly(z) Then inline = 1End IfIf inline = 1 Then Exit ForIf Ly(zn) Ly(z) Then If Lx(z) + (Lx(zn) - Lx(z) * (Py - Ly(z) / (Ly(zn) - Ly(z) = Px Then inline = 1: Exit For点在边上判断lc(z) = Gcnn(Lx(z), Ly(z), Lx(zn), Ly(zn) 判断是否相交，线段端点在延长线上为不相交ls(z) = Gsid(Lx(z), Ly(z), Lx(zn), Ly(zn) 线段端点是否在延长线上，如在，判断线段在哪侧pinshp = pinshp + lc(z)Next zIf inline = 1 Then pinshp = 1: Exit Function 如改If inline = 1 Then pinshp = 0则点在边上判断为在多边形外For z = 0 To Lnzn = z - 1If zn 0 And ls(zn) 0 ThenIf ls(z) ls(zn) Then pinshp = pinshp + 1End IfNext zpinshp = pinshp And 1End FunctionPublic Function Gcnn(zx1, zy1, zx2, zy2)Gcnn = 0If zx1 = Px And zy1 = Py Then Gcnn = 0: Exit FunctionIf zx2 = Px And zy2 = Py Then Gcnn = 0: Exit FunctionIf zy1 Py And zy2 Py Then Gcnn = 0: Exit FunctionIf zy1 Py And zy2 Px Then Gcnn = 1End FunctionPublic Function Gsid(zx1, zy1, zx2, zy2)Gsid = 0If zx1 = Px And zy1 = Py ThenIf zy2 = Px And zy2 = Py ThenIf zy1 Py Then Gsid = 1 Else Gsid = 2End IfEnd Function5.判断点在多边形内，射线算法 1. 已知点point(x,