如何确定一个二维点是否在多边形内?(涉及IT技术)

605

我正在尝试创建一个快速的2D点在多边形内部算法,用于命中测试(例如Polygon.contains(p:Point))。欢迎提供有效技巧建议。


您忘了告诉我们关于左右利手的问题您的看法——这也可以被解释为“内部”与“外部”。 - Richard T
17
是的,我现在意识到这个问题留下了许多未指定的细节,但目前我有点想看看各种回答的多样性。 - Scott Evernden
5
一个拥有90个面的多边形被称为九十边形(或称九角形),一个拥有10,000个面的多边形被称为万边形。 - user263678
1
“最优雅”已经不是目标了,因为我一直在寻找一个“能够工作”的算法。我必须自己想出来:https://dev59.com/GmUq5IYBdhLWcg3wEcVZ#18190354 - davidkonrad
这个库已经实现了它,所以你只需要在 Python 中使用 point.within(polygon),它会返回一个布尔值,表示点是否在多边形内部。 - J Agustin Barrachina
42个回答

881
对于图形,我不喜欢使用整数。许多系统将整数用于UI绘制(毕竟像素是整数),但例如macOS会为所有内容使用浮点数。macOS只知道点,一个点可以转换为一个像素,但是根据显示器分辨率,它可能会转换成其他东西。在视网膜屏幕上,半个点(0.5/0.5)是像素。尽管如此,我从未注意到macOS UI与其他UI相比明显更慢。毕竟,3D API(OpenGL或Direct3D)也使用浮点数,现代图形库往往利用GPU加速。
现在您说速度是您的主要关注点,好的,让我们考虑速度。在运行任何复杂算法之前,首先进行简单测试。在多边形周围创建一个轴对齐的边框框,这非常容易、快速,而且已经可以节省大量计算。如何做到这一点?迭代多边形的所有点,并找到X和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, 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!
}

这是任何点运行的第一个测试。如您所见,此测试速度非常快,但也非常粗糙。为了处理在边界矩形内的点,我们需要更复杂的算法。有几种方法可以计算这个问题。哪种方法有效还取决于多边形是否可能有洞或始终是实心的。以下是实心多边形的示例(一个凸多边形和一个凹多边形): Polygon without hole 这里是一个带有洞的多边形: Polygon with hole 绿色的中间有一个洞!
最简单的算法,可以处理上述所有三种情况,并且仍然非常快,称为射线投射。算法的思想非常简单:从多边形外部的任何位置向您的点绘制虚拟射线,并计算它击中多边形边的次数。如果命中次数是偶数,则在多边形外部,如果是奇数,则在内部。
演示射线如何穿过多边形: Demonstrating how the ray cuts through a polygon 缠绕数算法将是一种替代方案,它对于非常靠近多边形线的点更准确,但速度要慢得多。由于浮点精度有限和舍入问题,射线投射可能会对距离多边形边太近的点失败,但实际上这几乎不是问题,因为如果一个点靠得那么近,观察者通常甚至无法识别它是否已经在内部或仍然在外部。
您还记得上面的边界框吗?只需选择一个在边界框之外的点,并将其用作射线的起点即可。例如,点(Xmin-e/p.y)肯定在多边形外部。
那么,什么是e?好吧,e(实际上是epsilon)给边界框一些填充。正如我所说,如果我们离多边形线太近,则射线跟踪会失败。由于边界框可能等于多边形(如果多边形是轴对齐矩形,则边界框等于多边形本身!),因此我们需要一些填充来使其安全,就是这样。应该选择多大的e?不要太大。这取决于您用于绘图的坐标系比例尺。如果您的像素步长为1.0,则只需选择1.0(但0.1也可以)。
现在,我们有了带有其起始和结束坐标的射线,问题从“点是否在多边形内部”转移到“射线与多边形边相交的次数有多少次”。因此,我们不能像以前那样只使用多边形点,现在我们需要实际的边。一条边始终由两个点定义。
side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:

您需要对所有边进行光线测试。将光线视为向量,每条边也是一个向量。光线必须恰好击中每条边一次或者根本不击中。它不能两次击中同一条边。在二维空间中,两条直线总是恰好相交一次,除非它们是平行的,在这种情况下它们永远不会相交。然而,由于向量具有有限长度,两个向量可能不平行,但仍然永远不会相交,因为它们太短了,永远无法相遇。

// 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的两个端点(v1x1/v1y1和v1x2/v1y2)和向量2的两个端点(v2x1/v2y1和v2x2/v2y2)。所以你有2个向量,4个点,8个坐标。YES和NO很清楚。YES会增加交集,NO不会做任何事情。
那么COLLINEAR呢?它意味着两个向量位于同一无限线上,取决于位置和长度,它们可能根本不相交或在无数个点相交。我不确定如何处理这种情况,无论如何都不会将其视为交集。好吧,这种情况在实践中相当罕见,因为由于浮点舍入误差,更好的代码可能不会测试== 0.0f,而是测试类似于< epsilon的某些东西,其中epsilon是一个相当小的数字。
如果您需要测试更多的点,则可以通过将多边形边的线性方程标准形式存储在内存中来加快整个过程,这样您就不必每次都重新计算。这将为您节省每次测试的两个浮点乘法和三个浮点减法,以换取在内存中存储每个多边形边的三个浮点值。这是典型的内存与计算时间之间的权衡。
最后但并非最不重要的一点:如果您可以使用3D硬件来解决问题,那么有一个有趣的替代方法。只需让GPU为您完成所有工作即可。创建一个屏幕外的绘画表面。将其完全填充为黑色。现在让OpenGL或Direct3D绘制您的多边形(或者甚至是所有多边形,如果您只想测试该点是否在任何多边形中,但您不关心哪个),并用不同的颜色(例如白色)填充多边形。要检查点是否在多边形内,请从绘图表面获取此点的颜色。这只是O(1)的存储器获取。
当然,如果您的绘图表面不能很好地适应GPU内存,则此方法仅可用。如果它无法适应GPU内存,则此方法比在CPU上执行要慢。如果它必须很大且您的GPU支持现代着色器,则仍然可以使用GPU通过将上述射线投射实现为GPU着色器,这绝对是可能的。对于大量的多边形或需要测试的大量点,这将收益,考虑某些GPU将能够并行测试64到256个点。但请注意,从CPU到GPU和反向传输数据始终很昂贵,因此对于仅测试几个点与几个简单多边形的情况,其中点或多边形是动态的并且会经常更改,GPU方法很少会有回报。

34
非常棒的回答。关于使用硬件来实现,我在其他地方也看到过它可能会很慢,因为你需要从显卡获取数据。但我非常赞同减轻CPU负担的原则。有人有任何关于如何在OpenGL中实现这个功能的好参考资料吗? - Gavin
6
+1 因为这很简单!主要问题是如果你的多边形和测试点在网格上对齐(这不是罕见的情况),那么你必须处理“欺骗性”的交点,例如直接穿过多边形的一个点!(产生两个而不是一个的奇偶性)。进入了这个奇怪的领域:https://dev59.com/O3E95IYBdhLWcg3wkeuq。计算机图形学充满了这些特殊情况:理论上简单,实践中复杂。 - Jared Updike
9
@RMorrisey问:您为什么这样认为?我不明白为什么对于凹多边形来说这种方法会失效。当多边形是凹的时候,射线可能会多次离开和重新进入多边形,但最终,如果点在内部则命中计数器将为奇数,如果在外部则为偶数,这也适用于凹多边形。 - Mecki
9
“快速绕数算法”是指http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm上描述的算法,它非常迅速... - S P
21
你使用斜杠来分隔x和y坐标让人感到困惑,它表达的意思是x除以y。更加清晰的方式是使用x,y(即x逗号y)。总体而言,这是一个有用的答案。 - Ash
显示剩余23条评论

625

我认为下面这段代码是最好的解决方案(来自这里):

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: 测试点的 X 和 Y 坐标。

这个算法短小高效,适用于凸多边形和凹多边形。如前所建议,您应该首先检查边界矩形,然后单独处理多边形孔。

其背后的思想非常简单。作者将其描述如下:

我从测试点沿着水平射线(增加 x,固定 y)运行一个半无限光线,并计算它穿过了多少条边缘。在每次交叉时,光线在内部和外部之间切换。这被称为乔丹曲线定理。

变量 c 每次水平射线穿过任何一条边缘时都会从 0 切换到 1 或从 1 切换到 0。因此,它基本上是跟踪穿过的边缘数量是偶数还是奇数。0 表示偶数,1 表示奇数。


10
它检查verty[i]verty[j]是否在testy两侧,因此它们永远不会相等。 - Peter Wood
3
确实有效。我创建了一个 github gist 来测试它(GLUT)。 - bobobobo
8
该代码不够健壮,我不建议使用它。这里有一篇文章提供了一些详细的分析:http://www-ma2.upc.es/geoc/Schirra-pointPolygon.pdf。 - Mikola
20
这种方法实际上确实有局限性(边缘情况):检查点(15,20)是否在多边形[(10,10),(10,20),(20,20),(20,10)]中会返回错误而非正确。同理,(10,20)或(20,15)也是如此。对于其他所有情况,算法都可以完美地工作,而边缘情况中的假负值对我的应用程序来说也没有问题。 - Alexander Pacha
18
@Alexander,实际上这是有意设计的:通过将左和下边界与顶部和右边界的处理方式相反,如果两个不同的多边形共享一条边,沿着这条边的任何点都只会位于一个多边形中。这是一个有用的特性。 - wardw
显示剩余18条评论

96

这是一个C#版本的nirg提出的解答,它来自于这位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;
    }

    // https://wrf.ecse.rpi.edu/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;
}

5
非常好,谢谢。我已将其转换为JavaScript。https://dev59.com/p3VC5IYBdhLWcg3wrDJd#17490923 - Philipp Lenssen
2
这比使用GraphicsPath.IsVisible快了1000倍以上!!边界框检查使函数变慢了约70%。 - James Brown
2
@GDavoli 这是一个效率问题。如果点不在多边形的边界框内,它就不能在多边形内。这基于一个假设,即查找多边形边界框的循环比第二个循环要快得多。在现代硬件上可能不是这样的。但在实际系统中,缓存每个多边形的边界框可能是有意义的,在这种情况下,边界框检查绝对是有意义的。 - M Katz
1
inside = !inside; 之后断开?不行。当它遇到各种多边形段时,算法会多次更改。 - M Katz
1
我想这是因为当我们遍历一个由点组成的多边形列表时,最后一个点与第一个点相连,存在“循环”逻辑。 - NotX
显示剩余5条评论

67

以下是基于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(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;
}

41

计算点p与多边形各个角之间的有向角之和。如果总有向角为360度,则点在内部。如果总角度为0,则该点在外面。

我更喜欢这种方法,因为它更加健壮,不太依赖数值精度。

计算交点数量的方法受限于您在计算交点数量时可能会“命中”一个角。

编辑:顺便说一下,这种方法适用于凸多边形和凹多边形。

编辑:我最近发现了整个关于此主题的维基百科文章


2
@DarenW:每个顶点只有一个acos;另一方面,这个算法应该是最容易在SIMD中实现的,因为它完全没有分支。 - Jasper Bekkers
2
首先使用边界框检查来解决“点远”的情况。对于三角函数,您可以使用预生成的表格。 - JOM
1
根据维基百科(http://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm)的说法,这个算法并不一定慢。你可以简单地检查象限而不是计算角度。 - kaka
2
我知道这是老帖子了,不确定是否有人会看到,但是... David,什么是“角度的定向和”?这只是我和所讨论点之间的角度,0..360吗?如果是这样,难道你不需要考虑多边形中的点数吗?它不仅适用于四边形吗? - Maury Markowitz
5
这是最佳解决方案,因为它的时间复杂度为O(n),需要最少的计算量。适用于所有多边形。我在30年前在我的第一份IBM工作中研究了这个解决方案。他们批准了它,并且今天仍然在他们的GIS技术中使用它。 - Dominic Cerisano
显示剩余6条评论

35

这个问题非常有趣。我有一个与其他回答不同的可行想法,可以使用角度和来决定目标点是在内部还是外部。更为人熟知的是绕数算法

设x为目标点。将区域的所有点[0、1、... n]连接到每个边界点的目标点,如果目标点位于该区域内,则所有角度的总和将为360度,否则角度总和将小于360。

参考以下图像以基本了解该想法:enter image description here

我的算法假定顺时针是正方向。
下面是一个潜在的输入:

[[-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

25

bobobobo提到的Eric Haines文章非常好。特别有趣的是比较算法性能的表格;与其他方法相比,角度汇总法真的很差。另一个有趣的事情是,优化措施,如使用查找网格将多边形进一步细分为“内部”和“外部”部分,即使在具有> 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

13

这个问题中的大多数答案都没有很好地处理所有的边角情况。以下是一些微妙的边角情况: ray casting corner cases 这是一个JavaScript版本,它能够处理所有的边角情况。

/** Get relationship between a point and a polygon using ray-casting algorithm
 * @param {{x:number, y:number}} P: point to check
 * @param {{x:number, y:number}[]} polygon: the polygon
 * @returns -1: outside, 0: on edge, 1: inside
 */
function relationPP(P, polygon) {
    const between = (p, a, b) => p >= a && p <= b || p <= a && p >= b
    let inside = false
    for (let i = polygon.length-1, j = 0; j < polygon.length; i = j, j++) {
        const A = polygon[i]
        const B = polygon[j]
        // corner cases
        if (P.x == A.x && P.y == A.y || P.x == B.x && P.y == B.y) return 0
        if (A.y == B.y && P.y == A.y && between(P.x, A.x, B.x)) return 0

        if (between(P.y, A.y, B.y)) { // if P inside the vertical range
            // filter out "ray pass vertex" problem by treating the line a little lower
            if (P.y == A.y && B.y >= A.y || P.y == B.y && A.y >= B.y) continue
            // calc cross product `PA X PB`, P lays on left side of AB if c > 0 
            const c = (A.x - P.x) * (B.y - P.y) - (B.x - P.x) * (A.y - P.y)
            if (c == 0) return 0
            if ((A.y < B.y) == (c > 0)) inside = !inside
        }
    }

    return inside? 1 : -1
}

这是最佳答案。所有其他答案都忽略了边角情况。 - Ahmad
最快且处理边缘情况! - AndrejH

12

我真的很喜欢Nirg发布并由bobobobo编辑的解决方案。我只是使它变得更适合我的使用,加入了JavaScript并让它更易读:

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;
}

11

Swift版本的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
    }
}

1
在计算b时,存在潜在的除零问题。只有当计算a显示可能相交时,才需要计算“b”和下一行的“c”。如果两个点都在上方或下方,则没有可能性-这由对“a”的计算描述。 - anorskdev

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接