GIS点在多边形内算法.docx

1、GIS点在多边形内算法1.判断点在多边形内外的简单算法 - 改进弧长法今天学图形学的时候发现了一个求多边形内外的超简单算法，当时觉得非常惊喜，后来晚上上完选修的时候在走廊遇到bug，bug也是很惊喜地感慨道：居然有甘简单既办法都捻唔到！遂将其写下，供大家分享，希望不会太火星。 这个算法是源自计算机图形学基础教程（孙家广，清华大学出版社），在该书的48-49页，名字可称为“改进的弧长法”。该算法只需O(1)的附加空间，时间复杂度为O(n)，但系数很小；最大的优点是具有很高的精度，只需做乘法和减法，若针对整数坐标则完全没有精度问题。而且实现起来也非常简单，比转角法和射线法都要好写且不易出错。 首先

2、从该收中摘抄一段弧长法的介绍：“弧长法要求多边形是有向多边形，一般规定沿多边形的正向，边的左侧为多边形的内侧域。以被测点为圆心作单位圆，将全部有向边向单位圆作径向投影，并计算其中单位圆上弧长的代数和。若代数和为0，则点在多边形外部；若代数和为2则点在多边形内部；若代数和为，则点在多边形上。” 按书上的这个介绍，其实弧长法就是转角法。但它的改进方法比较厉害：将坐标原点平移到被测点P，这个新坐标系将平面划分为4个象限，对每个多边形顶点P ，只考虑其所在的象限，然后按邻接顺序访问多边形的各个顶点P，分析P和Pi+1，有下列三种情况：（1）Pi+1在P的下一象限。此时弧长和加/2；（2）Pi+1在P的

3、上一象限。此时弧长和减/2；（3）Pi+1在Pi的相对象限。首先计算f=yi+1*x-xi+1*y（叉积），若f=0，则点在多边形上；若f0，弧长和加。 最后对算出的代数和和上述的情况一样判断即可。 实现的时候还有两点要注意，第一个是若P的某个坐标为0时，一律当正号处理；第二点是若被测点和多边形的顶点重合时要特殊处理。 以上就是书上讲解的内容，其实还存在一个问题。那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错。如边(3,0)-(-3,0)，按以上的处理，象限分别是第一和第二，这样会使代数和加/2，有可能导致最后结果是被测点在多边形外。而实际上被测点是在多边形上（该边穿过该

4、点）。 对于这点，我的处理办法是：每次算P和Pi+1时，就计算叉积和点积，判断该点是否在该边上，是则判断结束，否则继续上述过程。这样牺牲了时间，但保证了正确性。 具体实现的时候，由于只需知道当前点和上一点的象限位置，所以附加空间只需O(1)。实现的时候可以把上述的“/2”改成1，“”改成2，这样便可以完全使用整数进行计算。不必考虑顶点的顺序，逆时针和顺时针都可以处理，只是最后的代数和符号不同而已。整个算法编写起来非常容易。 针对以上算法，我写了一个代码，拿ZOJ 1081 Points Within进行测试，顺利Accepted。这证明该算法的正确性还是可以保障的。以下附上我的代码：/ ZOJ

5、 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

6、 ( %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

7、? ( 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)

8、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)内。本程序采用射

9、线法，由待测试点(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,con

10、stvertex_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?*/ stat

11、icintis_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

12、*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,

13、&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 do

14、ne 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 n; i+ ) /* First see if q=(0,0) is a vertex. */ if ( PiX=0 & PiY=0 ) return v; i1 = ( i + n - 1 ) % n; /* printf(e=(%d,%d)t, i1, i); */ /* if e straddles the x-axis. */ /* The commented-out statement is logically equivalent to the one following. */