我怎样才能确定一个二维点是否在一个多边形？
我正尝试在多边形algorithm中创build一个快速的 2D点，用于命中testing（例如Polygon.contains(p:Point) ）。 有效的技术build议，将不胜感激。
对于graphics，我宁愿不喜欢整数。 许多系统使用整数来绘制UI（毕竟像素是整数），但是例如macOS使用float来表示一切。 macOS只知道点和一个点可以转换为一个像素，但取决于显示器的分辨率，它可能会转化为别的东西。 在视网膜屏幕上半点（0.5 / 0.5）是像素。 不过，我从来没有注意到，macOS用户界面比其他用户界面慢得多。 毕竟，3D API（OpenGL或Direct3D）也可以与浮点数和现代graphics库一起工作，通常可以利用GPU加速。
现在你说速度是你最关心的问题，好吧，让我们去追求速度。 在运行任何复杂的algorithm之前，先做一个简单的testing。 在多边形周围创build一个轴alignment的边界框 。 这很容易，快速，已经可以保证你很多的计算。 这是如何工作的？ 遍历多边形的所有点，findX和Y的最小/最大值。
例如你有分(9/1), (4/3), (2/7), (8/2), (3/6) 。 这意味着Xmin是2，Xmax是9，Ymin是1，Ymax是7.具有两个边（2/1）和（9/7）的矩形外部的点不能在多边形内。
// p is your point, px is the x coord, py is the y coord if (px < Xmin || px > Xmax || py < Ymin || py > Ymax) { // Definitely not within the polygon! }
这是任何一点运行的第一个testing。 正如你所看到的，这个testing速度非常快，但也很粗糙。 要处理边界矩形内的点，我们需要更复杂的algorithm。 有几种方法可以计算出来。 哪种方法的作用还取决于如果多边形可以有孔或将永远是固体的事实。 这里是一些固体的例子（一个凸，一个凹）：
这里有一个洞：
绿色的中间有一个洞！
最简单的algorithm，可以处理以上所有三种情况，并且仍然非常快，被命名为光线投射 。 algorithm的思想非常简单：从多边形以外的任意位置绘制虚拟光线到点，并计算它碰撞多边形边的频率。 如果命中的数量是偶数，那么它在多边形之外，如果它是奇数，它就在里面。
蜿蜒数algorithm是一种替代方法，对于多边形线非常丢失的点更准确，但也更慢。 由于浮点精度和舍入问题的限制，光线投射可能会失败，因为浮点精度和舍入问题有限，但是实际上这并不是一个问题，就好像一个点靠近一边，通常在视觉上甚至不可能观众辨认它是否已经在里面或者仍然在外面。
你还有上面的边框，记得吗？ 只需在边界框外select一个点，并将其用作射线的起点。 例如，点(Xmin - e/py)确实在多边形之外。
但是什么是e ？ 那么， e （实际上epsilon）给边界框一些填充 。 正如我所说，如果我们开始太靠近多边形线，射线跟踪失败。 由于边界框可能等于多边形（如果多边形是一个轴alignment的矩形，边界框等于多边形本身！），我们需要一些填充来使这个安全，就这样。 你应该select多大的e ？ 不太大。 这取决于您用于绘制的坐标系统比例。 如果你的像素步长是1.0，那就select1.0（但是0.1也可以）
现在我们有了起点和终点坐标的射线，问题从“ 是多边形内的点 ”转移到“ 射线与多边形相交的频率 ”。 因此，我们不能像以前一样使用多边形点，现在我们需要实际的边。 一边总是由两点定义。
side 1: (X1/Y1)-(X2/Y2) side 2: (X2/Y2)-(X3/Y3) side 3: (X3/Y3)-(X4/Y4) :
你需要testing四面八方的光线。 考虑射线是vector，每一边都是vector。 射线必须完全击中每一侧或永远不会。 它不能两次击中同一方。 二维空间中的两条线总是只相交一次，除非它们是平行的，在这种情况下它们从不相交。 然而，由于vector的长度有限，两个vector可能并不平行，但仍然不会相交，因为它们太短而不能彼此相遇。
// Test the ray against all sides int intersections = 0; for (side = 0; side < numberOfSides; side++) { // Test if current side intersects with ray. // If yes, intersections++; } if ((intersections & 1) == 1) { // Inside of polygon } else { // Outside of polygon }
到目前为止，但是如何testing两个vector相交？ 这里有一些C代码（未testing），应该做的伎俩：
#define NO 0 #define YES 1 #define COLLINEAR 2 int areIntersecting( float v1x1, float v1y1, float v1x2, float v1y2, float v2x1, float v2y1, float v2x2, float v2y2 ) { float d1, d2; float a1, a2, b1, b2, c1, c2; // Convert vector 1 to a line (line 1) of infinite length. // We want the line in linear equation standard form: A*x + B*y + C = 0 // See: http://en.wikipedia.org/wiki/Linear_equation a1 = v1y2 - v1y1; b1 = v1x1 - v1x2; c1 = (v1x2 * v1y1) - (v1x1 * v1y2); // Every point (x,y), that solves the equation above, is on the line, // every point that does not solve it, is not. The equation will have a // positive result if it is on one side of the line and a negative one // if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector // 2 into the equation above. d1 = (a1 * v2x1) + (b1 * v2y1) + c1; d2 = (a1 * v2x2) + (b1 * v2y2) + c1; // If d1 and d2 both have the same sign, they are both on the same side // of our line 1 and in that case no intersection is possible. Careful, // 0 is a special case, that's why we don't test ">=" and "<=", // but "<" and ">". if (d1 > 0 && d2 > 0) return NO; if (d1 < 0 && d2 < 0) return NO; // The fact that vector 2 intersected the infinite line 1 above doesn't // mean it also intersects the vector 1. Vector 1 is only a subset of that // infinite line 1, so it may have intersected that line before the vector // started or after it ended. To know for sure, we have to repeat the // the same test the other way round. We start by calculating the // infinite line 2 in linear equation standard form. a2 = v2y2 - v2y1; b2 = v2x1 - v2x2; c2 = (v2x2 * v2y1) - (v2x1 * v2y2); // Calculate d1 and d2 again, this time using points of vector 1. d1 = (a2 * v1x1) + (b2 * v1y1) + c2; d2 = (a2 * v1x2) + (b2 * v1y2) + c2; // Again, if both have the same sign (and neither one is 0), // no intersection is possible. if (d1 > 0 && d2 > 0) return NO; if (d1 < 0 && d2 < 0) return NO; // If we get here, only two possibilities are left. Either the two // vectors intersect in exactly one point or they are collinear, which // means they intersect in any number of points from zero to infinite. if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR; // If they are not collinear, they must intersect in exactly one point. return YES; }
input值是向量1（ v1x1/v1y1和v1x2/v1y2 ）和向量2（ v2x1/v2y1和v2x2/v2y2 ）的两个端点 。 所以你有2个向量，4点，8个坐标。 YES和NO清楚。 YES增加十字路口， NO什么也不做。
COLLINEAR呢？ 这意味着两个vector位于相同的无限长线上，取决于位置和长度，它们根本不相交，或者相交于无穷多个点。 我不是很确定如何处理这种情况，我不会把它算作交叉。 那么由于浮点舍入错误，这种情况在实践中是非常less见的; 更好的代码可能不会testing== 0.0f ，而是像< epsilon ，其中epsilon是一个相当小的数字。
如果您需要testing更多点数，通过将多边形的线性方程标准forms保留在内存中，您肯定可以加速整个过程，所以您不必每次都重新计算这些点。 这将为您节省两次浮点乘法和每次testing的三次浮点减法，以换取每个多边形在存储器中存储三个浮点值。 这是一个典型的记忆与计算时间的折衷。
最后但并非最不重要的：如果您可以使用3D硬件来解决问题，那么有一个有趣的select。 让GPU为你做所有的工作。 创build一个屏幕外的绘画表面。 用黑色填充它。 现在让OpenGL或者Direct3D绘制你的多边形（如果你只是想testing点是否在其中任何一个，但是你不关心哪一个），然后用不同的多边形填充多边形（甚至所有的多边形）颜色，如白色。 要检查点是否在多边形内，从绘图表面获取此点的颜色。 这只是一个O（1）的内存提取。
当然，这种方法只有在您的绘图表面不必很大时才可用。 如果它不能适应GPU内存，这种方法比在CPU上执行速度慢。 如果它必须是巨大的，并且你的GPU支持现代着色器，你仍然可以通过实现上面所示的光线投射作为GPU着色器来使用GPU，这绝对是可能的。 对于更多数量的多边形或大量的点来testing，这将会得到回报，考虑一些GPU将能够并行testing64到256点。 不过要注意的是，将数据从CPU传输到GPU并返回是非常昂贵的，所以为了仅仅对几个简单的多边形进行几点testing，其中点或多边形是dynamic的并且将频繁变化，GPU方法将很less支付closures。
我认为下面这段代码是最好的解决scheme（从这里取得 ）：
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy) { int i, j, c = 0; for (i = 0, j = nvert-1; i < nvert; j = i++) { if ( ((verty[i]>testy) != (verty[j]>testy)) && (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ) c = !c; } return c; }
参数
- nvert ：多边形中的顶点数。 在上面提到的文章中讨论了是否重复第一个顶点。
- vertx，verty ：包含多边形顶点的x和y坐标的数组。
- testx，testy ：testing点的X坐标和Y坐标。
它既简短又高效，适用于凸多边形和凹多边形。 如前所述，您应该首先检查边界矩形，然后分别处理多边形的孔。
这背后的想法很简单。 作者描述如下：
我从testing点水平运行半无限射线（增加x，固定y），并计算它横过的边数。 在每个十字路口，射线在内部和外部之间切换。 这被称为Jordan曲线定理。
每当水平光线穿越任何边缘时，variablesc从0切换到1，从1切换到0。 所以基本上它是跟踪交叉的边数是偶数还是奇数。 0表示偶数，1表示奇数。
这是来自RPI教授的 nirg给出的答案的C＃版本。 请注意，使用该RPI源代码需要归因。
边界框检查已添加在顶部。 但是，正如詹姆斯·布朗所指出的那样，主代码几乎和边界框自检一样快，所以边界框检查实际上可以减缓整体操作，在大多数检查点位于边界框内的情况下。 所以你可以离开边界框检查出来，或者替代方法是预先计算你的多边形的边界框，如果他们不经常改变形状的话。
public bool IsPointInPolygon( Point p, Point[] polygon ) { double minX = polygon[ 0 ].X; double maxX = polygon[ 0 ].X; double minY = polygon[ 0 ].Y; double maxY = polygon[ 0 ].Y; for ( int i = 1 ; i < polygon.Length ; i++ ) { Point q = polygon[ i ]; minX = Math.Min( qX, minX ); maxX = Math.Max( qX, maxX ); minY = Math.Min( qY, minY ); maxY = Math.Max( qY, maxY ); } if ( pX < minX || pX > maxX || pY < minY || pY > maxY ) { return false; } // http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html bool inside = false; for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ ) { if ( ( polygon[ i ].Y > pY ) != ( polygon[ j ].Y > pY ) && pX < ( polygon[ j ].X - polygon[ i ].X ) * ( pY - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X ) { inside = !inside; } } return inside; }
这是一个基于Nirg方法的M. Katz答案的JavaScript变体：
function pointIsInPoly(p, polygon) { var isInside = false; var minX = polygon[0].x, maxX = polygon[0].x; var minY = polygon[0].y, maxY = polygon[0].y; for (var n = 1; n < polygon.length; n++) { var q = polygon[n]; minX = Math.min(qx, minX); maxX = Math.max(qx, maxX); minY = Math.min(qy, minY); maxY = Math.max(qy, maxY); } if (px < minX || px > maxX || py < minY || py > maxY) { return false; } var i = 0, j = polygon.length - 1; for (i, j; i < polygon.length; j = i++) { if ( (polygon[i].y > py) != (polygon[j].y > py) && px < (polygon[j].x - polygon[i].x) * (py - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) { isInside = !isInside; } } return isInside; }
计算点p与每个多边形顶点之间的angular度取向总和。 如果总的方位angular是360度，那么点就在里面。 如果总数为0，则该点在外部。
我更喜欢这种方法，因为它更稳健，对数值精度的依赖性更小。
计算交叉点数量的均匀性的方法是有限的，因为在计算交叉点数量的过程中可以“碰撞”顶点。
编辑：顺便说一下，这种方法适用于凹和多边形的多边形。
编辑：我最近发现了关于这个主题的整个维基百科文章 。
Bob Bobobo引用的Eric Haines文章非常出色。 特别有趣的是比较algorithm性能的表格; angular度求和方法相对于其他方法来说是非常糟糕的。 另外有趣的是，像使用查找网格进一步将多边形细分为“入”和“出”扇区的优化可以使testing甚至在具有> 1000侧的多边形上快得多。
无论如何，现在已经过时了，但是我的投票是“交叉”的方法，这几乎是Mecki描述的我所想的。 然而，我发现它最被David Bourke形容和编纂 。 我喜欢没有真正需要的三angular函数，它适用于凸和凹，并且在边数增加的情况下performance相当好。
顺便说一句，下面是Eric Haines关于兴趣的文章中的一个performance表格，在随机多边形上进行testing。
number of edges per polygon 3 4 10 100 1000 MacMartin 2.9 3.2 5.9 50.6 485 Crossings 3.1 3.4 6.8 60.0 624 Triangle Fan+edge sort 1.1 1.8 6.5 77.6 787 Triangle Fan 1.2 2.1 7.3 85.4 865 Barycentric 2.1 3.8 13.8 160.7 1665 Angle Summation 56.2 70.4 153.6 1403.8 14693 Grid (100x100) 1.5 1.5 1.6 2.1 9.8 Grid (20x20) 1.7 1.7 1.9 5.7 42.2 Bins (100) 1.8 1.9 2.7 15.1 117 Bins (20) 2.1 2.2 3.7 26.3 278
这个问题太有意思了。 我有另一个可行的想法不同于这篇文章的其他答案。
这个想法是使用angular度的总和来决定目标是内部还是外部。 如果目标在一个区域内，则目标与每两个边界点的angular度forms总和将是360.如果目标在外部，则总和将不是360.angular度具有方向。 如果angular度向后，angular度是负的。 这就像计算绕组数量一样 。
请参考这张图片，对这个想法有一个基本的了解： 
我的algorithm假定顺时针是正方向。 这是一个潜在的投入：
[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]
以下是实现该想法的Python代码：
def isInside(self, border, target): degree = 0 for i in range(len(border) - 1): a = border[i] b = border[i + 1] # calculate distance of vector A = getDistance(a[0], a[1], b[0], b[1]); B = getDistance(target[0], target[1], a[0], a[1]) C = getDistance(target[0], target[1], b[0], b[1]) # calculate direction of vector ta_x = a[0] - target[0] ta_y = a[1] - target[1] tb_x = b[0] - target[0] tb_y = b[1] - target[1] cross = tb_y * ta_x - tb_x * ta_y clockwise = cross < 0 # calculate sum of angles if(clockwise): degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C))) else: degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C))) if(abs(round(degree) - 360) <= 3): return True return False
当我是Michael Stonebraker下的研究员的时候，我做了一些工作 – 你知道，有Ingres ， PostgreSQL等等的教授。
我们意识到，最快的方法是首先做一个边界框，因为它是超快的。 如果它在边界框之外，则在外面。 否则，你做更艰苦的工作…
如果您想要一个很好的algorithm，请参阅开源项目PostgreSQL的源代码以获取地理工作…
我想指出的是，我们从来没有对左撇子和左撇子的洞察力（也可以expression为“内部”与“外部”的问题。
UPDATE
BKB的链接提供了许多合理的algorithm。 我在研究地球科学的问题，因此需要一个经纬度的解决scheme，它有一个特殊的问题 – 是在较小的区域还是在较大的区域？ 答案是，顶点的“方向”很重要 – 它可以是左手或右手，这样就可以将任何一个区域指定为“在”任何给定的多边形“内部”。 因此，我的工作使用了该页面上列举的三个解决scheme。
另外，我的工作使用了“在线”testing的单独function。
…因为有人问：我们发现，当vertices的数量超过了某个数字时，bounding boxtesting是最好的 – 在做必要的更长的testing之前做一个非常快速的testing…边界框是通过简单的最大的x，最小的x，最大的y和最小的y，把它们放在一起做出一个盒子的四个点。
另一个提示是：我们在网格空间中完成了所有更复杂的“光调”计算，然后重新投影到“真实”的经度/纬度上，从而避免了可能的错误当一条交叉的经线180和处理极地时，环绕着它。 工作很好！
真的很喜欢Nirg发布的由bobobobo编辑的解决scheme。 我只是让它更友好，更清晰易读
function insidePoly(poly, pointx, pointy) { var i, j; var inside = false; for (i = 0, j = poly.length - 1; i < poly.length; j = i++) { if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside; } return inside; }
由nirg答案的迅捷版本：
extension CGPoint { func isInsidePolygon(vertices: [CGPoint]) -> Bool { guard !vertices.isEmpty else { return false } var j = vertices.last!, c = false for i in vertices { let a = (iy > y) != (jy > y) let b = (x < (jx - ix) * (y - iy) / (jy - iy) + ix) if a && b { c = !c } j = i } return c } }
David Segond的答案几乎是标准的一般答案，理查德T是最常见的优化，尽pipe其他方面也是如此。 其他强大的优化是基于较不通用的解决scheme。 例如，如果要检查具有多个点的相同多边形，则对多边形进行三angular化可以大大加快速度，因为有许多非常快的TINsearchalgorithm。 另一种情况是，如果多边形和点在低分辨率的有限平面上（比如屏幕显示），则可以将多边形绘制到存储器映射的显示缓冲区中，并以给定的颜色进行绘制，然后检查给定像素的颜色，看它是否存在在多边形。
像许多优化一样，这些基于特定的而不是一般的情况，并且基于分摊时间而不是一次性使用来获得有益效果。
在这个领域工作，我发现Joeseph O'Rourkes'计算几何在C'ISBN 0-521-44034-3是一个很大的帮助。
微不足道的解决办法是将多边形分割成三angular形，然后按照这里所述的方法testing三angular形
如果你的多边形是CONVEX，那么可能会有更好的方法。 将多边形视为无限行的集合。 每条线将空间分成两部分。 对于每一个点来说，如果在线的一侧或另一侧容易说出来。 如果一个点在所有线的同一侧，那么它在多边形内。
我意识到这是旧的，但这里是在Cocoa中实施的光线投射algorithm，以防任何人感兴趣。 不知道这是做事的最有效的方式，但它可能会帮助别人。
- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point { NSBezierPath *currentPath = [path bezierPathByFlatteningPath]; BOOL result; float aggregateX = 0; //I use these to calculate the centroid of the shape float aggregateY = 0; NSPoint firstPoint[1]; [currentPath elementAtIndex:0 associatedPoints:firstPoint]; float olderX = firstPoint[0].x; float olderY = firstPoint[0].y; NSPoint interPoint; int noOfIntersections = 0; for (int n = 0; n < [currentPath elementCount]; n++) { NSPoint points[1]; [currentPath elementAtIndex:n associatedPoints:points]; aggregateX += points[0].x; aggregateY += points[0].y; } for (int n = 0; n < [currentPath elementCount]; n++) { NSPoint points[1]; [currentPath elementAtIndex:n associatedPoints:points]; //line equations in Ax + By = C form float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y; float _B_FOO = point.x - (aggregateX/[currentPath elementCount]); float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y); float _A_BAR = olderY - points[0].y; float _B_BAR = points[0].x - olderX; float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY); float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO); if (det != 0) { //intersection points with the edges float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det; float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det; interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint); if (olderX <= points[0].x) { //doesn't matter in which direction the ray goes, so I send it right-ward. if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) { noOfIntersections++; } } else { if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) { noOfIntersections++; } } } olderX = points[0].x; olderY = points[0].y; } if (noOfIntersections % 2 == 0) { result = FALSE; } else { result = TRUE; } return result; }
对于testing点的示例方法，用nirg的答案的Obj-C版本。 Nirg的回答对我来说效果不错。
- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test { NSUInteger nvert = [vertices count]; NSInteger i, j, c = 0; CGPoint verti, vertj; for (i = 0, j = nvert-1; i < nvert; j = i++) { verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue]; vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue]; if (( (verti.y > test.y) != (vertj.y > test.y) ) && ( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) ) c = !c; } return (c ? YES : NO); } - (void)testPoint { NSArray *polygonVertices = [NSArray arrayWithObjects: [NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)], [NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)], [NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)], [NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)], [NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)], [NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)], nil ]; CGPoint tappedPoint = CGPointMake(23.0, 70.0); if ([self isPointInPolygon:polygonVertices point:tappedPoint]) { NSLog(@"YES"); } else { NSLog(@"NO"); } }
我也认为360总结只适用于凸多边形，但事实并非如此。
This site has a nice diagram showing exactly this, and a good explanation on hit testing: Gamasutra – Crashing into the New Year: Collision Detection
C# version of nirg's answer is here: I'll just share the code. It may save someone some time.
public static bool IsPointInPolygon(IList<Point> polygon, Point testPoint) { bool result = false; int j = polygon.Count() - 1; for (int i = 0; i < polygon.Count(); i++) { if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) { if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) { result = !result; } } j = i; } return result; }
Java Version:
public class Geocode { private float latitude; private float longitude; public Geocode() { } public Geocode(float latitude, float longitude) { this.latitude = latitude; this.longitude = longitude; } public float getLatitude() { return latitude; } public void setLatitude(float latitude) { this.latitude = latitude; } public float getLongitude() { return longitude; } public void setLongitude(float longitude) { this.longitude = longitude; } } public class GeoPolygon { private ArrayList<Geocode> points; public GeoPolygon() { this.points = new ArrayList<Geocode>(); } public GeoPolygon(ArrayList<Geocode> points) { this.points = points; } public GeoPolygon add(Geocode geo) { points.add(geo); return this; } public boolean inside(Geocode geo) { int i, j; boolean c = false; for (i = 0, j = points.size() - 1; i < points.size(); j = i++) { if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) && (geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude())) c = !c; } return c; } }
.Net port:
static void Main(string[] args) { Console.Write("Hola"); List<double> vertx = new List<double>(); List<double> verty = new List<double>(); int i, j, c = 0; vertx.Add(1); vertx.Add(2); vertx.Add(1); vertx.Add(4); vertx.Add(4); vertx.Add(1); verty.Add(1); verty.Add(2); verty.Add(4); verty.Add(4); verty.Add(1); verty.Add(1); int nvert = 6; //Vértices del poligono double testx = 2; double testy = 5; for (i = 0, j = nvert - 1; i < nvert; j = i++) { if (((verty[i] > testy) != (verty[j] > testy)) && (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i])) c = 1; } }
There is nothing more beutiful than an inductive definition of a problem. For the sake of completeness here you have a version in prolog which might also clarify the thoughs behind ray casting :
Based on the simulation of simplicity algorithm in http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
Some helper predicates:
exor(A,B):- \+A,B;A,\+B. in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)). inside(false). inside(_,[_|[]]). inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) + X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]). get_line(_,_,[]). get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).
The equation of a line given 2 points A and B (Line(A,B)) is:
(YB-YA) Y - YA = ------- * (X - XA) (XB-YB)
It is important that the direction of rotation for the line is setted to clock-wise for boundaries and anti-clock-wise for holes. We are going to check whether the point (X,Y), ie the tested point is at the left half-plane of our line (it is a matter of taste, it could also be the right side, but also the direction of boundaries lines has to be changed in that case), this is to project the ray from the point to the right (or left) and acknowledge the intersection with the line. We have chosen to project the ray in the horizontal direction (again it is a matter of taste, it could also be done in vertical with similar restrictions), so we have:
(XB-XA) X < ------- * (Y - YA) + XA (YB-YA)
Now we need to know if the point is at the left (or right) side of the line segment only, not the entire plane, so we need to restrict the search only to this segment, but this is easy since to be inside the segment only one point in the line can be higher than Y in the vertical axis. As this is a stronger restriction it needs to be the first to check, so we take first only those lines meeting this requirement and then check its possition. By the Jordan Curve theorem any ray projected to a polygon must intersect at an even number of lines. So we are done, we will throw the ray to the right and then everytime it intersects a line, toggle its state. However in our implementation we are goint to check the lenght of the bag of solutions meeting the given restrictions and decide the innership upon it. for each line in the polygon this have to be done.
is_left_half_plane(_,[],[],_). is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] = [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA)); is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test). in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon), in_range(Y,YA,YB). all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line), in_y_range_at_poly(Coordinate,Line,Polygon), Lines). traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count). % This is the entry point predicate inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).
VBA VERSION:
Note: Remember that if your polygon is an area within a map that Latitude/Longitude are Y/X values as opposed to X/Y (Latitude = Y, Longitude = X) due to from what I understand are historical implications from way back when Longitude was not a measurement.
CLASS MODULE: CPoint
Private pXValue As Double Private pYValue As Double '''''X Value Property''''' Public Property Get X() As Double X = pXValue End Property Public Property Let X(Value As Double) pXValue = Value End Property '''''Y Value Property''''' Public Property Get Y() As Double Y = pYValue End Property Public Property Let Y(Value As Double) pYValue = Value End Property
MODULE:
Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean Dim i As Integer Dim j As Integer Dim q As Object Dim minX As Double Dim maxX As Double Dim minY As Double Dim maxY As Double minX = polygon(0).X maxX = polygon(0).X minY = polygon(0).Y maxY = polygon(0).Y For i = 1 To UBound(polygon) Set q = polygon(i) minX = vbMin(qX, minX) maxX = vbMax(qX, maxX) minY = vbMin(qY, minY) maxY = vbMax(qY, maxY) Next i If pX < minX Or pX > maxX Or pY < minY Or pY > maxY Then isPointInPolygon = False Exit Function End If ' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html isPointInPolygon = False i = 0 j = UBound(polygon) Do While i < UBound(polygon) + 1 If (polygon(i).Y > pY) Then If (polygon(j).Y < pY) Then If pX < (polygon(j).X - polygon(i).X) * (pY - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then isPointInPolygon = True Exit Function End If End If ElseIf (polygon(i).Y < pY) Then If (polygon(j).Y > pY) Then If pX < (polygon(j).X - polygon(i).X) * (pY - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then isPointInPolygon = True Exit Function End If End If End If j = i i = i + 1 Loop End Function Function vbMax(n1, n2) As Double vbMax = IIf(n1 > n2, n1, n2) End Function Function vbMin(n1, n2) As Double vbMin = IIf(n1 > n2, n2, n1) End Function Sub TestPointInPolygon() Dim i As Integer Dim InPolygon As Boolean ' MARKER Object Dim p As CPoint Set p = New CPoint pX = <ENTER X VALUE HERE> pY = <ENTER Y VALUE HERE> ' POLYGON OBJECT Dim polygon() As CPoint ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1 For i = 0 To <ENTER VALUE HERE> 'Same value as above Set polygon(i) = New CPoint polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through Next i InPolygon = isPointInPolygon(p, polygon) MsgBox InPolygon End Sub
Heres a point in polygon test in C that isn't using ray-casting. And it can work for overlapping areas (self intersections), see the use_holes argument.
/* math lib (defined below) */ static float dot_v2v2(const float a[2], const float b[2]); static float angle_signed_v2v2(const float v1[2], const float v2[2]); static void copy_v2_v2(float r[2], const float a[2]); /* intersection function */ bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr, const bool use_holes) { /* we do the angle rule, define that all added angles should be about zero or (2 * PI) */ float angletot = 0.0; float fp1[2], fp2[2]; unsigned int i; const float *p1, *p2; p1 = verts[nr - 1]; /* first vector */ fp1[0] = p1[0] - pt[0]; fp1[1] = p1[1] - pt[1]; for (i = 0; i < nr; i++) { p2 = verts[i]; /* second vector */ fp2[0] = p2[0] - pt[0]; fp2[1] = p2[1] - pt[1]; /* dot and angle and cross */ angletot += angle_signed_v2v2(fp1, fp2); /* circulate */ copy_v2_v2(fp1, fp2); p1 = p2; } angletot = fabsf(angletot); if (use_holes) { const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f); angletot -= nested * (float)(M_PI * 2.0); return (angletot > 4.0f) != ((int)nested % 2); } else { return (angletot > 4.0f); } } /* math lib */ static float dot_v2v2(const float a[2], const float b[2]) { return a[0] * b[0] + a[1] * b[1]; } static float angle_signed_v2v2(const float v1[2], const float v2[2]) { const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]); return atan2f(perp_dot, dot_v2v2(v1, v2)); } static void copy_v2_v2(float r[2], const float a[2]) { r[0] = a[0]; r[1] = a[1]; }
Note: this is one of the less optimal methods since it includes a lot of calls to atan2f , but it may be of interest to developers reading this thread (in my tests its ~23x slower then using the line intersection method).
For Detecting hit on Polygon we need to test two things:
- If Point is inside polygon area. (can be accomplished by Ray-Casting Algorithm)
- If Point is on the polygon border(can be accomplished by same algorithm which is used for point detection on polyline(line)).
To deal with the following special cases in Ray casting algorithm :
- The ray overlaps one of the polygon's side.
- The point is inside of the polygon and the ray passes through a vertex of the polygon.
- The point is outside of the polygon and the ray just touches one of the polygon's angle.
Check Determining Whether A Point Is Inside A Complex Polygon . The article provides an easy way to resolve them so there will be no special treatment required for the above cases.
You can do this by checking if the area formed by connecting the desired point to the vertices of your polygon matches the area of the polygon itself.
Or you could check if the sum of the inner angles from your point to each pair of two consecutive polygon vertices to your check point sums to 360, but I have the feeling that the first option is quicker because it doesn't involve divisions nor calculations of inverse of trigonometric functions.
I don't know what happens if your polygon has a hole inside it but it seems to me that the main idea can be adapted to this situation
You can as well post the question in a math community. I bet they have one million ways of doing that
If you are looking for a java-script library there's a javascript google maps v3 extension for the Polygon class to detect whether or not a point resides within it.
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3); var isWithinPolygon = polygon.containsLatLng(40, -90);
Google Extention Github
When using qt (Qt 4.3+), one can use QPolygon's function containsPoint
The answer depends on if you have the simple or complex polygons. Simple polygons must not have any line segment intersections. So they can have the holes but lines can't cross each other. Complex regions can have the line intersections – so they can have the overlapping regions, or regions that touch each other just by a single point.
For simple polygons the best algorithm is Ray casting (Crossing number) algorithm. For complex polygons, this algorithm doesn't detect points that are inside the overlapping regions. So for complex polygons you have to use Winding number algorithm.
Here is an excellent article with C implementation of both algorithms. I tried them and they work well.
This only works for convex shapes, but Minkowski Portal Refinement, and GJK are also great options for testing if a point is in a polygon. You use minkowski subtraction to subtract the point from the polygon, then run those algorithms to see if the polygon contains the origin.
Also, interestingly, you can describe your shapes a bit more implicitly using support functions which take a direction vector as input and spit out the farthest point along that vector. This allows you to describe any convex shape.. curved, made out of polygons, or mixed. You can also do operations to combine the results of simple support functions to make more complex shapes.
More info: http://xenocollide.snethen.com/mpr2d.html
Also, game programming gems 7 talks about how to do this in 3d (:
