如何确定 2D 点是否在多边形内？
https://stackoverflow.com/questions/217578
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian题
我正在尝试创建一个 快速地 多边形内的 2D 点算法，用于命中测试（例如 Polygon.contains(p:Point)）。对于有效技术的建议将不胜感激。
解决方案
对于图形，我宁愿不喜欢整数。许多系统使用整数进行UI绘制（毕竟像素是整数），但macOS例如使用float来表示一切。 macOS只知道点，一个点可以转换为一个像素，但根据显示器的分辨率，它可能会转换为其他像素。在视网膜屏幕上，半个点（0.5 / 0.5）是像素。尽管如此，我从未注意到macOS UI比其他UI明显慢。毕竟3D API（OpenGL或Direct3D）也适用于浮点数和现代图形库，经常利用GPU加速。
现在你说速度是你主要考虑的问题，好吧，让我们加快速度。在运行任何复杂算法之前，首先进行简单测试。在多边形周围创建一个轴对齐的边界框。这非常简单，快速，并且已经可以安全地进行大量计算。这是如何运作的？迭代多边形的所有点，找到X和Y的最小值/最大值。
E.g。你有分数(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, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
// Definitely not within the polygon!
}
这是第一个针对任何一点运行的测试。正如您所看到的，此测试速度非常快，但也非常粗糙。要处理边界矩形内的点，我们需要更复杂的算法。有几种方法可以计算出来。如果多边形可以有孔或者总是是实心的，那么哪种方法也起作用。以下是实体（一个凸面，一个凹面）的示例：
这是一个有洞的人：
绿色的中间有一个洞！
最简单的算法，可以处理上述所有三种情况并且仍然非常快，命名为光线投射。该算法的想法非常简单：从多边形外部的任何位置绘制虚拟光线到您的点，并计算它撞击多边形一侧的频率。如果命中数是偶数，则它在多边形之外，如果它是奇数，则它在内部。
绕组数算法将是一种替代方案，它对于非常接近多边形线的点更准确，但它也慢得多。由于有限的浮点精度和圆角问题，射线投射可能因太靠近多边形的点而失败，但实际上这几乎不是问题，就好像一个点靠近一侧，它通常在视觉上甚至不可能观众识别它是否已经在里面或者仍然在外面。
你还有上面的边框，还记得吗？只需在边界框外选取一个点，并将其用作光线的起点。例如。点(Xmin - e/p.y)肯定在多边形之外。
但是什么是e？好吧，v1x1/v1y1（实际上是epsilon）为边界框提供了一些填充。正如我所说，如果我们太靠近多边形线，光线跟踪就会失败。由于边界框可能等于多边形（如果多边形是一个轴对齐的矩形，边界框等于多边形本身！），我们需要一些填充来使这个安全，就是这样。你应该选择多大的v1x2/v1y2？不太大。它取决于您用于绘图的坐标系比例。如果您的像素步长为1.0，那么只需选择1.0（但0.1也可以使用）
现在我们有了射线的开始和结束坐标，问题从<！>转移; 是多边形 <！>中的点;到<！>“光线与多边形边相交的频率 <！>”;因此，我们不能像以前一样使用多边形点，现在我们需要实际的边。一个方总是由两点来定义。
side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:
你需要对所有方面测试光线。将光线视为矢量，将每一边视为矢量。射线必须完全击中每一侧或从不击中每一侧。它不能两次击中同一侧。 2D空间中的两条线总是恰好相交一次，除非它们是平行的，在这种情况下它们永远不会相交。然而，由于向量的长度有限，两个向量可能不是平行的，但仍然不会相交，因为它们太短而不能相互见面。
// 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
}
到目前为止，但是你如何测试两个向量是否相交？这里有一些C代码（未经过测试），应该可以解决这个问题：
#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;
}
输入值是向量1（v2x1/v2y1和v2x2/v2y2）和向量2（YES和NO）的两个端点。所以你有2个向量，4个点，8个坐标。 == 0.0f和< epsilon很清楚。 <=>增加交叉点，<=>什么都不做。
COLLINEAR怎么样？这意味着两个矢量位于相同的无限线上，取决于位置和长度，它们根本不相交或者它们在无数个点中相交。我不完全确定如何处理这种情况，我不认为它是交叉点。嗯，这种情况在实践中相当罕见，因为浮点舍入错误;更好的代码可能不会测试<=>，而是测试像<=>这样的东西，其中epsilon是一个相当小的数字。
如果你需要测试更多的点，你可以通过将多边形边的线性方程标准形式保存在内存中来加快整个过程，所以你不必每次都重新计算这些。这将在每次测试中为您节省两次浮点乘法和三次浮点减法，以换取在内存中存储每个多边形侧的三个浮点值。这是典型的内存与计算时间的权衡。
最后但并非最不重要：如果您可以使用3D硬件来解决问题，那么有一个有趣的选择。让GPU为您完成所有工作。创建一个屏幕外的绘画表面。用黑色完全填充。现在让OpenGL或Direct3D绘制多边形（或者甚至是所有多边形，如果你只想测试点是否在任何一个内，但你不关心哪一个）并用不同的方法填充多边形颜色，例如白色。要检查点是否在多边形内，请从绘图表面获取此点的颜色。这只是一个O（1）内存提取。
当然，只有绘图表面不必很大时，此方法才可用。如果它不能适合GPU内存，这种方法比在CPU上运行要慢。如果它必须是巨大的并且您的GPU支持现代着色器，您仍然可以通过将上面显示的光线投射实现为GPU着色器来使用GPU，这绝对是可能的。对于更大数量的多边形或大量的测试点，这将得到回报，考虑一些GPU将能够并行测试64到256个点。但请注意，将数据从CPU传输到GPU并返回总是很昂贵，因此只需针对几个简单的多边形测试几个点，其中点或多边形是动态的并且会经常变化，GPU方法很少会付出代价关闭。
其他提示
我认为下面的代码是最好的解决方案（摘自 这里):
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;
}
论据
- 反转:多边形中的顶点数。是否在最后重复第一个顶点已在上面提到的文章中讨论过。
- 顶点，垂直:包含多边形顶点的 x 和 y 坐标的数组。
- testx，暴躁的:测试点的 X 和 y 坐标。
它既短又高效，并且适用于凸多边形和凹多边形。正如之前所建议的，您应该首先检查边界矩形并单独处理多边形孔。
这背后的想法非常简单。作者是这样描述的：
我从测试点水平运行一条半无限射线（增加 x，固定 y），并计算它穿过的边数。在每次交叉处，光线都会在内部和外部之间切换。这称为乔丹曲线定理。
每当水平射线穿过任何边缘时，变量 c 就会从 0 切换到 1，以及从 1 切换到 0。所以基本上它是跟踪交叉的边数是偶数还是奇数。0 表示偶数，1 表示奇数。
这是 nirg给出的答案的C＃版本，来自这位RPI教授。请注意，使用该RPI源代码需要归因。
顶部添加了边界框检查。但是，正如James Brown指出的那样，主代码几乎与边界框检查本身一样快，因此边界框检查实际上可以减慢整体操作，如果您检查的大多数点都在边界框内。因此，您可以将边界框退出，或者如果它们不经常更改形状，则可以预先计算多边形的边界框。
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( q.X, minX );
maxX = Math.Max( q.X, maxX );
minY = Math.Min( q.Y, minY );
maxY = Math.Max( q.Y, maxY );
}
if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > 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 > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
{
inside = !inside;
}
}
return inside;
}
以下是M. Katz基于Nirg方法的答案的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(q.x, minX);
maxX = Math.max(q.x, maxX);
minY = Math.min(q.y, minY);
maxY = Math.max(q.y, maxY);
}
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
}
var i = 0, j = polygon.length - 1;
for (i, j; i < polygon.length; j = i++) {
if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
isInside = !isInside;
}
}
return isInside;
}
计算点p与每个多边形顶点之间的方向角度和。如果总取向角度是360度，则该点在内部。如果总数为0，则该点在外面。
我更喜欢这种方法，因为它更稳健，更少依赖于数值精度。
计算交叉点数均匀度的方法是有限的，因为在计算交叉点数量时你可以“击中”一个顶点。
EDIT：By The Way，这种方法适用于凹凸多边形。
编辑：我最近在这个主题上发现了一整页维基百科文章。bobobobo引用的 Eric Haines文章非常出色。特别有趣的是比较算法性能的表格;角度求和方法与其他方法相比非常糟糕。同样有趣的是优化，例如使用查找网格将多边形进一步细分为<！>“;在<！>中;和<！>“; out <！>”;即使在具有<！> gt的多边形上，扇区也可以使测试速度极快1000方面。
无论如何，这是早期的，但我的投票是转到<！>“交叉点<！>”;方法，这几乎是Mecki所描述的。然而，我发现它最为自豪地由David Bourke描述和编纂。我喜欢不需要真正的三角函数，它适用于凸面和凹面，并且随着边数的增加，它的表现也相当不错。
顺便说一句，这是Eric Haines的文章中关于随机多边形测试的性能表之一。
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
这个问题非常有趣。我有另一个可行的想法，不同于这篇文章的其他答案。
这个想法是使用角度之和来决定目标是在内部还是外部。如果目标位于区域内，则目标和每两个边界点形成的角度之和将为360.如果目标位于外部，则总和将不是360.角度具有方向。如果角度向后，角度是负角度。这就像计算蜿蜒数字一样。
请参阅此图片以基本了解该想法：
我的算法假设顺时针是正方向。这是一个潜在的输入：
[[-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
快速版的由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 = (i.y > y) != (j.y > y)
let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
if a && b { c = !c }
j = i
}
return c
}
}
当我在 Michael Stonebraker 担任研究员时，我做了一些工作 - 你知道，这位教授想出了 Ingres ， PostgreSQL 等。
我们意识到最快的方法是首先做一个边界框，因为它的速度非常快。如果它在边界框之外，它就在外面。否则，你会做更努力的工作......
如果你想要一个很好的算法，请查看开源项目PostgreSQL源代码以进行地理工作......
我想指出，我们从来没有深入了解右手与左手（也可以表达为<！>“;内部<！>”vs <！>;在<！>之外的问题..
更新
BKB的链接提供了大量合理的算法。我正在研究地球科学问题，因此需要一个在纬度/经度上工作的解决方案，它有一个特殊的手性问题 - 是较小区域内还是较大区域内的区域？答案是<！>“方向<！>”;顶点问题很重要 - 它是左手或右手，这样你可以指出任何区域为<！>“;在<！>内;任何给定的多边形因此，我的工作使用了该页面上列举的解决方案三。
另外，我的作品在行<！>上使用了单独的函数<！>;测试
...自从有人问：我们发现当顶点数量超过一定数量时，边界框测试最好 - 如果需要，在进行更长时间的测试之前做一个非常快速的测试...边界框由简单地取最大x，最小x，最大y和最小y并将它们组合在一起制作一个盒子的四个点......
对于接下来的人的另一个提示：我们做了所有更复杂的和<！>“调光<！>”;在网格空间中计算所有在平面上的正点，然后重新投射回<！>“真实的<！>”;经度/纬度，从而避免了当经过一条经线180和处理极区时可能出现的缠绕错误。工作得很好！
非常喜欢Nirg发布并由bobobobo编辑的解决方案。我只是让它变得友好，并且对我的使用更加清晰：
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;
}
David Segond的答案几乎是标准的一般答案，而Richard T是最常见的优化，尽管其他一些也是如此。其他强有力的优化基于不太通用的解决方案。例如，如果要检查具有大量点的相同多边形，对多边形进行三角测量可以大大加快速度，因为有许多非常快速的TIN搜索算法。另一个是如果多边形和点在低分辨率的有限平面上，比如屏幕显示，您可以将多边形绘制到给定颜色的内存映射显示缓冲区中，并检查给定像素的颜色以查看它是否位于在多边形中。
与许多优化一样，这些优化基于特定情况而非一般情况，并根据摊销时间而非单次使用产生收益。
在这个领域工作，我发现Joeseph O'Rourkes的计算几何在C'ISBN 0-521-44034-3中是一个很好的帮助。
简单的解决方案是将多边形划分为三角形并按照解释的方式测试三角形这里
如果你的多边形是 CONVEX ，那么可能会有更好的方法。将多边形视为无限线的集合。每条线将空间分成两部分。对于每一点，很容易说它是在线的一侧还是另一侧。如果一个点位于所有线的同一侧，则它位于多边形内。
我意识到这是旧的，但这是一个在Cocoa中实现的光线投射算法，以防任何人感兴趣。不确定这是最有效的做事方式，但它可以帮助别人。
- (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;
}
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");
}
}
nirg答案的C＃版本在这里：我将只分享代码。这可能会节省一些时间。
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;
}
没有什么比问题的归纳定义更美妙了。为了完整起见，你在prolog中有一个版本，它也可以澄清 ray casting 背后的背景：
基于 http：//中简单算法的模拟www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
一些助手谓词：
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]).
给出2点A和B（线（A，B））的线的等式是：
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
线的旋转方向很重要 设置为时钟方式的边界和反时钟方式的孔。 我们要检查点（X，Y），即测试点是否在左边 我们行的半平面（这是一个品味的问题，它也可能是 右边，但也必须改变边界线的方向 这种情况），这是从点到右（或左）投射光线 并确认与线的交叉点。我们选择投射 光线在水平方向（再次是味道的问题， 它也可以在垂直方向上完成，具有类似的限制），所以我们有：
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
现在我们需要知道该点是否位于左侧（或右侧） 只有线段，而不是整个平面，所以我们需要 将搜索仅限制为此段，但这很容易 在该段内只有一个点可以更高 垂直轴上的Y比。因为这是一个更强的限制 需要成为第一个检查，所以我们首先只采取那些线 满足此要求，然后检查其可能性。由约旦 曲线定理投射到多边形的任何光线必须在a处相交 偶数行。所以我们完成了，我们将光线投射到 正确，然后每次它与一条线相交，切换它的状态。 但是在我们的实现中，我们要检查它的长度 满足给定限制的解决方案包并决定 坚持它。对于多边形中的每一行，必须这样做。
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).
Java版本：
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端口：
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;
}
}
VBA VERSION：
注意：请记住，如果您的多边形是地图中的一个区域，纬度/经度是Y / X值而不是X / Y（纬度= Y，经度= X），因为根据我的理解，历史意义来自于当经度不是衡量标准的时候。
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(q.X, minX)
maxX = vbMax(q.X, maxX)
minY = vbMin(q.Y, minY)
maxY = vbMax(q.Y, maxY)
Next i
If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > 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 > p.Y) Then
If (polygon(j).Y < p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - 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 < p.Y) Then
If (polygon(j).Y > p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - 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
p.X = <ENTER X VALUE HERE>
p.Y = <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
我已经做了一个Python实现 尼尔格的 由 小码哥发布于 代码:
输入
- 边界点： 构成多边形的节点。
边界框位置： 要过滤的候选点。（在我的实现中，从边界框创建。
（输入是元组列表，格式如下：
[(xcord, ycord), ...])
退货
- 多边形内部的所有点。
def polygon_ray_casting(self, bounding_points, bounding_box_positions):
# Arrays containing the x- and y-coordinates of the polygon's vertices.
vertx = [point[0] for point in bounding_points]
verty = [point[1] for point in bounding_points]
# Number of vertices in the polygon
nvert = len(bounding_points)
# Points that are inside
points_inside = []
# For every candidate position within the bounding box
for idx, pos in enumerate(bounding_box_positions):
testx, testy = (pos[0], pos[1])
c = 0
for i in range(0, nvert):
j = i - 1 if i != 0 else nvert - 1
if( ((verty[i] > testy ) != (verty[j] > testy)) and
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
c += 1
# If odd, that means that we are inside the polygon
if c % 2 == 1:
points_inside.append(pos)
return points_inside
同样，这个想法取自 这里
在C中没有使用光线投射的多边形测试中的一个点。它可以用于重叠区域（自交叉），请参阅use_holes参数。
/* 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];
}
注意：这是一个不太理想的方法，因为它包含大量对atan2f的调用，但是读取此线程的开发人员可能会感兴趣（在我的测试中，它比使用行交集方法慢〜23倍） ）。
对于检测多边形上的命中，我们需要测试两件事：
- 如果Point位于多边形区域内。 （可以通过射线投射算法完成）
- 如果Point位于多边形边界上（可以通过用于折线（线）上的点检测的相同算法完成）。
在确定点是否在复杂多边形内部。本文提供了一种解决问题的简便方法，因此上述案例不需要特殊处理。
您可以通过检查通过将所需点连接到多边形顶点而形成的区域是否与多边形本身的区域匹配来执行此操作。
或者您可以检查从您的点到每对两个连续多边形顶点到您的检查点的内角之和是否总和为360，但我觉得第一个选项更快，因为它不涉及分裂或三角函数逆的计算。
我不知道如果你的多边形里面有一个洞会发生什么，但在我看来，主要想法可以适应这种情况
您也可以在数学社区发布问题。我打赌他们有百万种方法可以做到这一点
如果您正在寻找一个java脚本库，那么Polygon类的javascript google maps v3扩展可以检测一个点是否位于其中。
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
使用 qt 时（Qt 4.3+） ，可以使用QPolygon的功能 containsPoint
答案取决于你是否有简单或复杂的多边形。简单多边形不得有任何线段交叉点。所以他们可以有洞，但线不能相互交叉。复杂区域可以具有线交叉 - 因此它们可以具有重叠区域，或者仅通过单个点彼此接触的区域。
对于简单多边形，最好的算法是Ray cast（交叉数）算法。对于复杂多边形，此算法不会检测重叠区域内的点。因此，对于复杂多边形，您必须使用绕组数算法。
这是一篇优秀的文章，其中包含两种算法的C实现。我尝试了它们并且效果很好。
nirg解决方案的Scala版本（假设边界矩形预检是单独完成的）：
def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {
val length = polygon.length
@tailrec
def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
if (i == length)
tracker
else {
val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
}
}
oddIntersections(0, length - 1, tracker = false)
}
这是golang版本的@nirg答案（灵感来自@@ m-katz的C＃代码）
func isPointInPolygon(polygon []point, testp point) bool {
minX := polygon[0].X
maxX := polygon[0].X
minY := polygon[0].Y
maxY := polygon[0].Y
for _, p := range polygon {
minX = min(p.X, minX)
maxX = max(p.X, maxX)
minY = min(p.Y, minY)
maxY = max(p.Y, maxY)
}
if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
return false
}
inside := false
j := len(polygon) - 1
for i := 0; i < len(polygon); i++ {
if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
inside = !inside
}
j = i
}
return inside
}
令人惊讶的是没有人提前提出这个问题，但对于需要数据库的实用主义者来说：MongoDB对Geo查询提供了很好的支持，包括这个。
您正在寻找的是：
db.neighborhoods.findOne（{geometry：{$ geoIntersects：{$ geometry：{ type：<！> quot; Point <！> quot;，coordinates：[<！> quot; longitude <！> quot;，<！> quot; latitude <！> quot; ]}}} }）
Neighborhoods是以标准GeoJson格式存储一个或多个多边形的集合。如果查询返回null，则它不会相交，否则就是。
这里有很好的记录： https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
分类在330不规则多边形网格中的超过6,000个点的性能不到一分钟，根本没有优化，包括用各自的多边形更新文档的时间。
