《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
11、;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;k2=(i+l)%np;whi1e(k2!=
12、i&vlk2.y=w.y)k2=(k2+l)%np;if(kl!=k2&is_same(v,&w,vl+kl,vl+k2)=0)c+;if(k2InPoly:q=);PrintPoint(q)：putchar(,);*Shiftsothatqistheorigin.Notethisdestroysthepolygon.ThisisdoneforpedogicalcIarity.*/for(i=O;in;i+)(for(d=O;dDIM;d+)Pid=Pid-qd;)*Foreachedgee=(i-l,i),seeifcrossesray.*/for(i=O;iO)&(PilYO)&(PiYO
13、)!=(PilYO)*estraddlesray,socomputeintersectionwithray.*/x=(PiX*(double)PilY-PilX*(double)PiY)/(double)(PilY-PiY)：*printf(,/straddles:x=%gt”,x);*/*crossesrayifstrictlypositiveintersection.*/if(x0)Rcross+;)*printf(wRightCrOSS=%dt”,Rcross)；*/*ifestraddlesthex-axiswhenreversed.*/*if(PiY=O)H(PilY=O)*/if(
14、PiYO)!=(piiO)*estraddlesray,socomputeintersectionwithray.*/x=(PX*PilY-PilX*PiY)/(double)(PilY-PiY);*Printf(straddles:x=%gt”,);*/*crossesrayifstrictlypositiveintersection.*/if(x0)1.cross+;)*printf(w1.eftCroSS=%dn”,1.cross);*/)*qontheedgeifleftandrightcrossarenotthesameparity.*/if(Rcross%2)(1.cross%2)
15、returnanoddnumbero*qinsideiffcrossings.*/if(Rcross%2)returni;elsereturno;)voidPrintPoint(tPointip)putchar(,()；for(i=0;iDIM;i+PrintfP)；if(i!=DIM-I)putcharC,)putchar(,);*Readsinthecoordinatesoftheverticesofapolygonfromstdin,putsthemintoP,andreturnsn,thenumberofvertices.Formattingconventions:etc.*/intR
16、eadPoly(tPolygoniP)(inti,n;doprintf(Inputthenumberofvertices:nz,);scanf(%d,&n);if(n=PMX)break;Printf(Errorinreadpoly:toomanypoints;maxis%dn”,PMAX);)while1)；printf(Polygon:n);rintf(iXyn);for(i=0;in;i+)(scanf(%d%d,&Pi0,&Pi1)；printf(*%3d%4d%4dni,Pi0,Pi1);)printf(wn=%3dverticesreadn”,n);utchar(,n,);retu
17、rnn;)voidPrintPoly(intn,tPoIygoniP)(inti;Printf(Polygon:n);printf(*ixyn*);for(i=O;in;i+)printf(w%3d%4d%4dnz，,i,Pi0,P*Copypolygonatob(overwritingb).*/voidCopy(tPolygonia,tPoIygonib,intn)(inti,j;for(i=0;in;i+)for(j=0;jDIM;j+)bij=aij;)boolEqPoint(tPointia,tPointib)inti;for(i=O;iDIM;i+aibi)returnFA1.SE;
18、returnTRUE;)1. 已知点POint(X,y)和多边形POIygOn(xl,yl;x2,y2;.xn,yn;)；2. 以Point为起点，以无穷远为终点作平行于X轴的直线Iine(X,y;-8,y)；1.1. 循环取得(for(i=0;in;i+)多边形的每一条边Side(Xi,yi;xi+l,yi+l),且推断是否平行于X轴,假如平行continue,否则，i+；4. 同时推断Point(x,y)是否在Side上，假如是,则返回1(点在多边形上)，否则接着下面的推断；5. 推断线Side及Iine是否有交点，假如有则count+,否则，i+o6. 推断交点的总数，假如为奇数则返回0
19、(点在多边形内)，偶数则返回2(点在多边形外)。代码：/*射线法推断点q及多边形polygon的位置关系，要求p。Iygon为简洁多边形，顶点逆时针排列假如点在多边形内：返回O假如点在多边形边上：返回1假如点在多边形外：返回2constdoubleINFINITY=IelO;constdoub1eESP=le-5:constintMAX_N=1000;structPointdoublex,y;struct1.ineSegmentPointptl,pt2;typedefvectorPolygon;/计算叉乘POPlPOP2idoubleMultiply(Pointpl,Pointp2lPoint
20、p)return(pl.x-p.x)*(p2.y-p.y)-(p2.x-p.x)*(pl.y-p.y);/推断线段是否包含点pointboolIsOnline(Pointpoint,1.ineSegmentline)return(fabs(Multiply(1ine.ptl,1ine.pt2,point)ESP)&(point,x-line.ptl.x)*(point,x-line.pt2.x)=0)(point,y-line,ptI.y)*(point,y-line.pt2.y)=nin(1.2.ptl.x,1.2.pt2.x)&(max(1.2.ptl.x,1.2.pt2.x)=min(1
21、.l.pt1.x,1.I.pt2.x)&(max(1.I.ptl.y,1.I.pt2.y)=min(1.2.pt1.y,1.2.pt2.y)&(max(1.2.ptl.y,1.2.pt2.y)=min(1.l.pt1.y,1.I.pt2.y)&(Multiply(1.2.ptl,1.I.pt2,1.I.ptl)*Multiply(1.I.pt2,1.2.pt2,1.l.ptl)=0)&(Multiply(1.1.ptl,1.2.pt2,1.2.ptl)*Multiply(1.2.pt2,1.1.Pt2,1.2.ptl)=0)/推断点在多边形内boolInPolygon(constPOlygOn
22、&polygon,Pointpoint)intn=polygon,size();intcount=O;1.ineSegmentline;1ine.ptl=point;line.pt2.y=point,y;line.pt2.x=-INFINITY：for(inti=0：in;i+)/得到多边形的一条边1.ineSegmentside;side,ptl=polygoni;side.pt2=polygon(i+1)%n;if(IsOnline(point,side)returnl;/假如Side平行X轴则不作考虑if(fabs(side.ptl.y-side.pt2.y)side.pt2.y)cou
23、nt+;elseif(IsOnline(side.pt2,1ine)if(side.pt2.yside.ptl.y)count+;elseif(Intersectdine,side)count+;if(count%2=I)return0;elsereturn2;4.用射线相交来推断点是否在多边形内的算法点在边上，射线通过顶点，射线通过边的状况都考虑了:Public1.x(20)sDoub1ePublic1.y(20)AsDoublePublic1.nsInteger定义多边形的坐标，边数PublicPxAsDoublePublicPysDouble定义点的坐标该过程推断点是否在多边形内，点在边
24、上推断为在多边形内，方法为射线相交法，返回1为在多边形内，O为在多边形外PublicFunctionpinshp()DimIc(20)AsIntegerDimIs(20)AsIntegerpinshp=0inline=0Forz=0To1.nzn=z+1Ifzn1.nThenzn=0点在边上推断If1.y(zn)=1.y(z)AndPy=1.y(z)ThenIfPx1.x(z)AndPx1.x(zn)AndPx1.y(z)Theninline=1EndIfIfinline=1ThenExitForIf1.y(zn)O1.y(z)ThenIf1.x(z)+(1.x(zn)-1.x(z)*(Py-
25、1.y(z)/(1.y(zn)-1.y(z)=PxTheninline=1:ExitFor点在边上推断lc(z)=Gcnn(1.x(z),1.y(z),1.x(zn),1.y(zn)推断是否相交,线段端点在延长线上为不相交ls(z)=Gsid(1.x(z),1.y(z),1.x(zn)t1.y(zn)线段端点是否在延长线上，如在，推断线段在哪侧pinshp=pinshp+lc(z)NextzIfinline=1Thenpinshp=1:ExitFunction如改Ifinline=1Thenpinshp=0则点在边上推断为在多边形夕卜Forz=0To1.nzn=z-1Ifzn0AndIs(zn
26、)0ThenIfls(z)Ols(zn)Thenpinshp=pinshp+1EndIfNextzpinshp=pinshpAnd1EndFunctionPublicFunctionGcnn(zxl,zyl,zx2,zy2)Gcnn=0Ifzxl=PxAndzyl=PyThenGcnn=0:ExitFunctionIfzx2=PxAndzy2=PyThenGcnn=0:ExitFunctionIfzylPyAndzy2PyThenGcnn=0:ExitFunctionIfzylPyAndzy2PxThenGcnn=1EndFunctionzyl,zx2,zy2)Pub1icFunctionGs
27、id(zxl,Gsid=0Ifzxl=PxAndzyl=PyThenIfzy2=PxAndzy2=PyThenIfzylPyThenGsid=1ElseGsid=2EndIfEndFunction5.推断点在多边形内，射线算法1. 已知点POint(x,y)和多边形POlygon(xl,yl;x2,y2：.xn,yn;)；2. 以point为起点，以无穷远为终点作平行于X轴的直线Iine(X,y;-8,y)：1.1. .循环取得(for(i=0;in;i+)多边形的每一条边Side(Xi,yi;xi+l,yi+l),且推断是否平行于X轴，假如平行continue,否则，i+；4. 同时推断Po
28、int(x,y)是否在Side上，假如是，则返回1(点在多边形上)，否则接着下面的推断；5. 推断线Side及Iine是否有交点，假如有则count+,否则，i+。6.推断交点的总数，假如为奇数则返回O（点在多边形内)，偶数则返回2(点在多边形外)。代码：/*射线法推断点q及多边形polygon的位置关系，要求polygon为简洁多边形，顶点逆时针排列假如点在多边形内：返回0假如点在多边形边上：返回1假如点在多边形外：返回2constdoubleINFINITY=IelO;constdoub1eESP=le-5：constintMAX_N=1000;structPointdoublex,y;s
29、truct1.ineSegmentPointptl,pt2;typedefvectorPolygon;/计算叉乘IPOPlIIPOP21doubleMultiply(Pointpl,Pointp2,Pointp)return(pl.x-p.x)*(p2.y-p.y)-(p2.x-p.x)*(pl.y-py)；/推断线段是否包含点pointboolIsOnline(Pointpoint,1.ineSegmentline)return(fabs(Multiply(line.ptl,line,pt2,point)ESP)&(point,x-line.ptl.x)*(point,x-1ine.pt2.
30、x)=0)(point,y-line.ptl.y)*(point,y-1ine.pt2.y)=min(1.2.ptl.x,1.2.pt2.x)&(max(1.2.ptl.x,1.2.pt2.x)=min(1.l.pt1.x,1.I.pt2.x)&(max(1.I.ptl.y,1.I.pt2.y)=min(1.2.ptl.y,1.2.pt2.y)&(max(1.2.ptl.y,1.2.pt2.y)=min(1.l.ptl.y,1.I.pt2.y)&(Multiply(1.2.ptl,1.I.pt2,1.I.ptl)*Multily(1.l.pt2,1.2.pt2,1.l.ptl)=0)&(Mul
31、tiply(1.I.ptl,1.2.pt2,1.2.pt!)*Multiply(l2.pt2,1.I.pt2,1.2.ptl)=0)/推断点在多边形内boo)InPoIygon(constPoIygon&polygon,Pointpoint)intn=polygon,size();intcount=O;1.ineSegmentline;line,ptl=point;line.pt2.y=point,y;line.pt2.x=-INFINITY;for(inti=O;in;i+)/得到多边形的一条边1.ineSegmentside;side,ptl=polygoni;side.pt2=polygon(i+1)%n;if(IsOnline(point,side)return1;/假如side平行x轴则不作考虑if(fabs(side.ptl.y-side.pt2.y)side.pt2.y)count+;elseif(IsOnline(side.pt2,1ine)if(side.pt2.yside.ptl.y)count+;elseif(Intersectdine,side)count+;if(count%2=1)return0：)elsereturn2;
