GIS点在多边形内算法.docx

1.推断点在多边形内外的简洁算法一改进弧长法今日学图形学的时候发觉了个求多边形内外的超简洁算法，当时觉得特别惊喜，后来晚上上完选修的时候在走廊遇到bug,bug也是很惊喜地感慨道：尽然有甘简洁既方法都捻唔到！遂将其写下，供大家共享，希望不会太火星。这个算法是源自计算机图形学基础教程(孙家广，清华高校出版社)，在该书的48-49页，名字可称为“改进的弧长法”。该算法只需O(I)的附加空间，时间困难度为0(n),但系数很小；最大的优点是具有很高的精度，只需做乘法和减法，若针对整数坐标则完全没有精度问题。而且实现起来也特别简洁，比转角法和射线法都要好写且不易出错。首先从该收中摘抄段弧长法的介绍：“弧长法要求多边形是有向多边形，一般规定沿多边形的正向，边的左侧为多边形的内侧域。以被测点为圆心作单位圆，将全部有向边向单位圆作径向投影，并计算其中单位圆上弧长的代数和。若代数和为0,则点在多边形外部；若代数和为2兀则点在多边形内部；若代数和为11,则点在多边形上。”按书上的这个介绍,其实弧长法就是转角法。但它的改进方法比较厉害：将坐标原点平移到被测点P，这个新坐标系将平面划分为4个象限，对每个多边形顶点P,只考虑其所在的象限，然后按邻接依次访问多边形的各个顶点P，分析P和Pi+1,有下列三种状况：(1)Pi+1在P的下一象限。此时弧长和加112;(2) Pi+1在P的上一象限。此时弧长和减112:(3) Pi+1在Pi的相对象限。首先计算f=yi+l*xri+l*y(叉积)，若f=0,则点在多边形上：若f<0,弧长和减丸：若f>0,弧长和加11。最终对算出的代数和和上述的状况一样推断即可。实现的时候还有两点要留意,第个是若P的某个坐标为0时,一律当正号处理：其次点是若被测点和多边形的顶点重合时要特别处理。以上就是书上讲解的内容，其实还存在一个问题。那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错。如边(3,0)-(-3,0),按以上的处理，象限分别是第一和其次，这样会使代数和加112,有可能导致最终结果是被测点在多边形外。而事实上被测点是在多边形上(该边穿过该点)。对于这点，我的处理方法是：每次算P和Pi+1时，就计算叉积和点积，推断该点是否在该边上，是则推断结束，否则接着上述过程。这样牺牲了时间，但保证了正确性。详细实现的时候，由于只需知道当前点和上一点的象限位置，所以附加空间只需0(1)。实现的时候可以把上述的“JT/2”改成1,“丸”改成2,这样便可以完全运用整数进行计完。不必考虑顶点的依次，逆时针和顺时针都可以处理，只是最终的代数和符号不同而已。整个免法编写起来特别简洁。针对以上算法,我写了一个代码，拿ZOJlO81PointSWithin进行测试，顺当ACCePted。这证明该算法的正确性还是可以保障的。以下附上我的代码：/ZOJ1081,改进弧长法判点在形内形外include<stdio.h>include<math.h>constintMAX=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("Problem'd:n”,prob);scanf("%d",&m);for(i=0;i<n;i+)scanf("%d,&p.x,&p.y):pn=p0;while(m-)(scanf(''c%d”,&t.x,&t.y);f(i=0;i<=n;i÷*-).-=tx,Py-=t.y;/坐标平移tl=p0.x>=0?(p0.y>=0?0:3):(pO.y>=O?l:2);/计算象限for(sum=0,i=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)break;/点在边上t2=p.x>=0?(p.y>=0?0:3):(p.y>=0?l:2);/计算象限if(t2=(tl+1)%4)sum+=1;/状况1elseif(t2=(tl÷3)4)am-=1；/状况2elseif(t2=(tl+2)%4)/状况3i11C0sumc*sum-=2;tl=t2;if(i<=n；)sun)Prinlf(,Vithin11")；elsePrintf("Outsiden");f(i=O;i<=n;i-H-)p.x=tx,p.y+=t.y;/复原坐标)return0;2*C本文是采纳射线法推断点是否在多边形内的C语言程序。多年前，诘我自己实现了这样个算法。但是随着时间的推移，我确定重写这忘个代码。参考周培德的计完几何一书，结合我的实践和阅历，q我信任，在这个算法的实现上，这是你迄今为止遇到的最优的代码。到这是个C语言的小算法的实现程序，原来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察6下来竟然个符合须要的也没有。我对自己高校读书时写的代码没右有信念，所以，确定重新写一个，并把它放到这里，以飨读者。也图增加下B1.OG的点击量。开首先定义点结构如下：以下是引用片段：豺CoPycode*Vertexstructure*/typedefstruct（doublex,y;vertex_t;本算法里所指的多边形，是指由一系列点序列组成的封闭简洁多边形。它的首尾点可以是或不是同个点（不强制要求首尾点是同一个点）。这样的多边形可以是随意形态的，包括多条边在一条肯定直线上。因此，定义多边形结构如下：以下是引用片段：Copycode*Vertexliststructure-polygon*/typedefstructintnumvertices;*Numberofverticesinlist*/vertex_t*vertex;*Vertexarraypointer*/vertexlistt;为加快判别速度，首先计算多边形的外包矩形(recjt),推断点是否落在外包矩形内，只有满意落在外包矩形内的条件的点，才进入下一步的计和。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent,代码如下：以下是引用片段：Copycode*boundingrectangletype*/typedefstruct(doubleminx,miny,maxx,maxy;rect_t;*getsextentofvertices*/voidVerticesgetextent(constvertext*vl,intnp,*invertices*/rectt*rc*outextent*/)inti;if(np>0)(rc->minx=rc->maxx=vl0.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.x<if(vli.y<if(vli.x>if(vli.y>rc->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边线的交点数目，记为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)(doubledx=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,constvertex_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(np<3)returnO;vertices_get_extent(vl,np,&rc);if(v->x<rc.min_xIIv->x>rc.max_xv->y<rc.minyjv->y>rc.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;