GIS点在多边形内算法#知识参考

1、1 .在多边形内外判断点的简单算法弧长法的改进今天学习图学的时候，发现了求多边形内外的超单纯算法。 那时非常吃惊，晚上选择结束的时候在走廊上遇到了臭虫，臭虫也变得惊人地感慨深刻，然后写下来，让大家共享，希望不要太火星。该算法来源于计算机图形学基础教程 (孙家广，清华大学出版社)，本书第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时，减去弧长和，如果f0，则加上弧长和。关于最后计算出的代数，可以按照与上述情况相同的方式进行判断。在实现时，第一，p的某坐标为0的情况下，必须注意一律为正的第二点，在被测定点与多边形的顶点重叠的情况下，需要特殊的处理。以上是书中说明的内容，其实还有一个问题。 如果多边形的边位于坐标轴上，且两个顶点分别位于原点的两侧，则会发生错误。 如边(3，0 )-(-3，0 )，在以上处理中，象限分别是第一和第二，将代数和/2相加，最终测量点有可能在多边形之外。 实际上被测量点在多

4、边形上(这一边通过这一点)。对于这一点，我的应对方法是每次计算p和Pi 1时计算叉积和点积，判断该点是否在那里，否则判断为结束，否则继续上述过程。 这样牺牲了时间，但保证了正确性。具体地说，因为只知道当前点和前点象限的位置，所以附加空间只需要O(1)。 实现时，通过将上述的“/2”变更为1，将“”变更为2，就可以完全用整数计算。 不需要考虑顶点的顺序，能逆时针和顺时针处理的只是最后的代数和符号不同。 写整个算法很容易。对于以上的算法，我测试了ZOJ 1081 Points Within，写了一个顺利接受Accepted的代码。 这证明了该算法的正确性得到了保障。下面附上我的代码。/ZOJ 10

5、81，改善弧长法的判定点在形的内外#include#includeconst int MAX=101；结构点 int x，y； TRint main ()举止int n、m、I、和、t1、t2、f、prob=0；点t；while (扫描( % d ，n )，n )举止if (打印)打印(n )；打印( problem 360n ，prob )扫描( % d ，m )；for (i=0； PS； I )扫描( % d % d ，p.x，p.y )；pn=p0；while (m - )举止扫描( % d % d ，t.x，t.y )；for (i=0； i=n； i ) p.x -=t.x、p.y

6、 -=t.y； /坐标位移t1=p0.x=0？ (p0.y=0？ (p 0.y=0？ 1:2； 计算/象限for (和=0，i=1； i=n； PS )举止PK (！ p.x！ p.y )中断； /被测量点是多边形的顶点f=p.y * pi-1.x - p.x * pi-1.y； /计算叉积PK (！ FP I-1 .x * p.x=0p I-1 .y * p I .y=0)中断； /点在边上t2=p.x=0？ (p.y=0？ (p.y=0)？ 1:2； 计算/象限PS (T2=(t11 ) %4)和=1； /情况1else if (t2=(t1 3 ) % 4 ) sum -=1； /情况2

7、else if (t2=(t1 2 ) % 4 ) /情况3举止PC (f0)和=2； else sum -=2；以下t1=t2；以下if (I=n| sum )打印( withinn )； else printf(Outsiden )；for (i=0； i=n； i ) p.x =t.x、p.y =t.y； /恢复坐标以下以下返回0；以下2.c语言中的实现点多边形内算法本稿是用放射线法来判断点是否在多边形内的c语言程序。 多年前，我自己实现了这种算法。 但是随着时间的流逝，我决定重写这个代码。 借鉴周培德的计算几何本，结合我的实践和经验，相信这个算法的实现，这是你见过的最好的代码。这是c语

8、言的小算法的实现程序，本来不想放在这里。 但是，我自己要实现这样的算法时，想在互联网上找现成的，考察也不符合需求。 我对自己在大学学习时写的代码没有自信，所以决定再写一遍，把它放在这里，向读者提出请求。 也请增加博客的点击数。首先，将点结构定义如下参照剪辑如下所示复制代码/* vertex结构* /typedefstruct举止双倍，y；vertex_t；本算法中所谓的多边形，是指由一系列的点列构成的闭合的单纯的多边形。 首尾点可以是同一点，也可以不是同一点(首尾点不一定是同一点)。 这种多边形可以是任意形状，包括多个边缘在绝对直线上。 因此，多边形结构的定义如下参照剪辑如下所示复印代码/*

9、vertexliststructure-polygon * /typedefstruct举止intnum_vertices； /*numberoffverticesinlist*/vertex_t*vertex； /*Vertexarraypointer*/vertexlist_t；为了提高判别速度，首先计算多边形的外接矩形(rect_t )，判断点是否收敛于外接矩形内，只对满足收敛于外接矩形内的条件的点进行下一个计算。 为此，引入用于确定外接矩形结构rect_t和点集合的外接矩形中的方法vertices_get_extent。 代码如下所示。参照剪辑如下所示复印代码/* boundingrc

10、tangletype * /typedefstruct举止双最小值_ x，最小值_ y，最大值_ x，最大值_ y；rect_t；/*getsextentofvertices*/voidvvertices _ get _ extent (const vertex _ t * VL，intnp，/*invertices*/rect_t*rc/*outextent*/)举止inti；将PS (PK ) )rc-min_x=rc-max_x=vl0.x； rc-min_y=rc-max_y=vl0.y；elseRC-min _ x=RC-min _ y=RC-max _ x=RC-max _ y=0

11、； /*=0？ noverticesatall*/以下for(i=1； 举止PS (PS _ x ) PS-PS _ x=PS.x；PS (PS _ y ) PS _ y=PS；y；PS (PS ).xrc-max _ x=VL I .x；PS (PS ).yrc-max _ y ) RC-max _ y=VL I .y；以下以下当点满足收敛在多边形外接矩形内的条件时，还判断点(v )是否在多边形(vl:np )内。 在本程序中，使用放射线法，从被试验点(v )水平引出一条放射线B(v，w )，计算b和vl的边界线的交点数，标记为c，根据奇内偶外原则(c是奇数，v不在vl内，否则v不在vl内)

12、，判定点是否在多边形内具体原理不言而喻。 要计算线段之间是否有交点，请引入以下函数(1)is_same确定第二点是否与直线l(l_start，l_end )相同(2)is_intersect用于判断两条线段(不是直线) s1、s2是(1)还是(0)。参照剪辑如下所示复印代码/* p，qisonthesameoflinel*/staticintis _ same (const vertex _ t * l _ start，constvertex_t*l_end，/*linel*/constvertex_t*pconstvertex_t*q )举止doubleledx=l_end-x-l_star

13、t-x；doubleledy=l_end-y-l_start-y；dx1=p-x-l _ star t-x；doubleledy1=p-y-l_start-y；双倍dx2=q-x-l _ end-x；doubleledy2=q-y-l_end-y；回复(dx * dy1-dy * dx1 ) * (dx * dy2-dy * dx2 ) 0？ 1:0；以下/* 2线段(s 1，s2)areintersect？ */staticintis _ intersect (const vertex _ t * S1 _ start，constvertex_t*s1_end )constvertex_t*

14、s2_start，constvertex_t*s2_end )举止返回(is _ same (s1_ start，s1_end，s2_start，s2_end)=0AA (S2 _ start，s2_end，s1_start，s1_end)=0？ 1:0；以下以下函数pt_in_poly是判断点(v )是否在多边形(vl:np )内的程序参照剪辑如下所示copycodeintpt _ in _ poly (const vertex _ t * VL，intnp，/*polygonvlwithnpvertices*/constvertex_t*v )举止inti，j，k1，k2，c；rect_t

15、rc；vertex_tw；PS (PS3)返回0；vertices_get_extent(vl、np、rc )if(v-xxrc.max_x|v-yyrc.max_y )返回0；/*Setahorizontalbeaml(*v，w)fromvtotheultraright*/w.x=rc.max_x DBL_EPSILON；w.y=v-y；c=0； /*intersectpointscounter*/for(i=0； 举止j=(i 1)%np；PS (PS，PS，PS，PS，w ) )举止c；以下elseif(vli.y=w.y )举止k1=(np i-1)%np；PS (PK 1！=ivlk1.y=w.yk1=(np k1-1)%np；k2=(i