GIS点在多边形内算法.docx

《GIS点在多边形内算法.docx》由会员分享，可在线阅读，更多相关《GIS点在多边形内算法.docx（31页珍藏版）》请在优知文库上搜索。

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

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

3、:(3) Pi+1在Pi的相对象限。首先计算f=yi+l*xri+l*y(叉积)，若f=0,则点在多边形上：若f0,弧长和加11。最终对算出的代数和和上述的状况一样推断即可。实现的时候还有两点要留意,第个是若P的某个坐标为0时,一律当正号处理：其次点是若被测点和多边形的顶点重合时要特别处理。以上就是书上讲解的内容，其实还存在一个问题。那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错。如边(3,0)-(-3,0),按以上的处理，象限分别是第一和其次，这样会使代数和加112,有可能导致最终结果是被测点在多边形外。而事实上被测点是在多边形上(该边穿过该点)。对于这点，我的处理方法

4、是：每次算P和Pi+1时，就计算叉积和点积，推断该点是否在该边上，是则推断结束，否则接着上述过程。这样牺牲了时间，但保证了正确性。详细实现的时候，由于只需知道当前点和上一点的象限位置，所以附加空间只需0(1)。实现的时候可以把上述的“JT/2”改成1,“丸”改成2,这样便可以完全运用整数进行计完。不必考虑顶点的依次，逆时针和顺时针都可以处理，只是最终的代数和符号不同而已。整个免法编写起来特别简洁。针对以上算法,我写了一个代码，拿ZOJlO81PointSWithin进行测试，顺当ACCePted。这证明该算法的正确性还是可以保障的。以下附上我的代码：/ZOJ1081,改进弧长法判点在形内形外i

5、ncludeincludeconstintMAX=101;structpointintx,y;pMX;intmain()(intn,m,i,sum,tl,t2,f,prob=0;pointt;while(scanf(ld”,&n),n)(if(prob+)printf(n);printf(Problemd:n”,prob);scanf(%d,&m);for(i=0;in;i+)scanf(%d,&p.x,&p.y):pn=p0;while(m-)(scanf(c%d”,&t.x,&t.y);f(i=0;i=0?(p0.y=0?0:3):(pO.y=O?l:2);/计算象限for(sum=0,i

6、=1;i=n;i+)(if(!口点y)b向/被测点为多边形顶点f=p.y*pi-l.X-p.X*pi-l.y;/计算又积if(!f&pi-l.x*p.X=0&pi-l.y*pi.y=0?(p.y=0?0:3):(p.y=0?l:2);/计算象限if(t2=(tl+1)%4)sum+=1;/状况1elseif(t2=(tl3)4)am-=1；/状况2elseif(t2=(tl+2)%4)/状况3i11C0sumc*sum-=2;tl=t2;if(i=n；)sun)Prinlf(,Vithin11)；elsePrintf(Outsiden);f(i=O;i0)(rc-minx=rc-maxx=vl

7、0.x;rc-miny=rc-max_y=vlO.y;elserc-minx=rc-min_y=rc-maxx=rc-maxy=0;*=0?noverticesatall*/for(i=l;i(if(vli.xif(vli.yif(vli.yrc-min_x)rc-miny)rc-max_x)rc-maxy)rc-min_xrc-minyrc-max_xrc-maxy=vli.x;=vli.y;=vli.x;=vli.y;当点满意落在多边形外包矩形内的条件，要进一步推断点(V)是否在多边形(vl：np)内。本程序采纳射线法，由待测试点(V)水平引出一条射线B(v,w),计算B及VI边线的交点数

8、目，记为c,依据奇内偶外原则(C为奇数说明V在Vl内，否则V不在Vl内)推断点是否在多边形内。详细原理就不多说。为计算线段间是否存在交点，引入下面的函数:(1) is_same推断2(p、q)个点是(1)否(O)在直线l(l_start,1.end)的同侧;(2) is_intersect用来推断2条线段(不是直线)si、s2是否(0)相交;以下是引用片段：Copycode*p,qisonthesameofline1*/staticintis_same(constvertex_t*l_start,constvertext*lend,*line1*/constvertext*q)(doubled

9、x=lend-x-lstart-x;doubledy=l_end-y-l-start-y;doubledxl=p-x-lstart-x:doubledyl=-y-l-start-y:constP.vertex_t*doubledx2=q-x-l_end-x;doubledy2=q-y-lend-y;return(dx*dy1-dy*dx1)*(dx*dy2dy*dx2)0?1:0);*21inesegments(slts2)areintersect?*/staticintisintersect(constvertext*sistart,constvertex_t*sl_end,constver

10、tex_t*s2_start,constvertex_t*s2_end)(return(issame(sistart,siend,s2start,s2end)=0&is_same(s2_start,s2_end,S1.Start,S1.end)=0)?1:0下面的函数pt_in_poly就是推断点(v)是(1)否(0)在多边形（vl：np）内的程序:以下是引用片段：Copycodeintptinpoly(constvertext*vl,intnp,*polygonvlwithnpvertices*/constvertext*v)(inti,j,kl,k2,c;rect_trc;vertextw;if(npxxrc.max_xv-yyrc.maxy)returnO;*Setahorizontalbeam1(*v,w)fromvtotheultraright*/w.x=rc.max_x+DB1._EPSI1.ON;w.y=v-y;c=0;/*Intersectionpointscounter*/for(i=0;ij=(i+l)%np:if(isintersect(vl+i,vl+j,v,&w)(c+;elseif(vli.y=w.y)(kl=(np+i-l)%np;while(kl!=i&vlkl.y=w.y)kl=(np+kl-l)%np;