如何确定2D点是否在多边形内？

Question

我正在尝试在多边形算法中创建一个快速 2D点,用于命中测试(例如Polygon.contains(p:Point)).对于有效技术的建议将不胜感激.

Answer 1

对于图形,我宁愿不喜欢整数.许多系统使用整数进行UI绘制(毕竟像素是整数),但macOS例如使用float来表示一切.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!
}

这是第一个针对任何一点运行的测试.正如您所看到的,此测试速度非常快,但也非常粗糙.要处理边界矩形内的点,我们需要更复杂的算法.有几种方法可以计算出来.如果多边形可以有孔或者总是是实心的,那么哪种方法也起作用.以下是实体(一个凸面,一个凹面)的示例:

而这里有一个洞:

绿色的中间有一个洞!

最简单的算法,可以处理上述所有三种情况,并且仍然非常快,称为光线投射.该算法的想法非常简单:从多边形外部的任何位置绘制虚拟光线到您的点,并计算它撞击多边形一侧的频率.如果命中数是偶数,则它在多边形之外,如果它是奇数,则它在内部.

该绕数的算法将是一种选择,它更准确的点非常接近多边形线,但它也慢得多.由于有限的浮点精度和圆角问题,射线投射可能因太靠近多边形的点而失败,但实际上这几乎不是问题,就好像一个点靠近一侧,它通常在视觉上甚至不可能观众识别它是否已经在里面或者仍在外面.

你还有上面的边框,还记得吗？只需在边界框外选取一个点,并将其用作光线的起点.例如,该点(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)
:

你需要测试所有侧面的光线.将光线视为矢量,将每一边视为矢量.射线必须完全击中每一侧或从不击中每一侧.它不能两次击中同一侧.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(和)和向量2(和)的两个端点.所以你有2个向量,4个点,8个坐标.而且很清楚.增加交叉点,什么都不做.v1x1/v1y1v1x2/v1y2v2x1/v2y1v2x2/v2y2YESNOYESNO

COLLINEAR怎么样？这意味着两个矢量位于相同的无限线上,取决于位置和长度,它们根本不相交或者它们在无数个点中相交.我不完全确定如何处理这种情况,我不认为它是交叉点.嗯,这种情况在实践中相当罕见,因为浮点舍入错误; 更好的代码可能不会测试,== 0.0f而是用于类似的东西< epsilon,其中epsilon是一个相当小的数字.

如果你需要测试更多的点,你可以通过将多边形边的线性方程式标准形式保存在内存中来加快整个过程,所以你不必每次都重新计算这些.这将在每次测试中为您节省两次浮点乘法和三次浮点减法,以换取在内存中存储每个多边形侧的三个浮点值.这是典型的内存与计算时间的权衡.

最后但并非最不重要:如果您可以使用3D硬件来解决问题,那么有一个有趣的选择.让GPU为您完成所有工作.创建一个屏幕外的绘画表面.用黑色完全填充.现在让OpenGL或Direct3D绘制多边形(或者甚至是所有多边形,如果你只想测试点是否在任何一个内,但你不关心哪一个)并用不同的方法填充多边形颜色,例如白色.要检查点是否在多边形内,请从绘图表面获取此点的颜色.这只是一个O(1)内存提取.

Of course this method is only usable if your drawing surface doesn't have to be huge. If it cannot fit into the GPU memory, this method is slower than doing it on the CPU. If it would have to be huge and your GPU supports modern shaders, you can still use the GPU by implementing the ray casting shown above as a GPU shader, which absolutely is possible. For a larger number of polygons or a large number of points to test, this will pay off, consider some GPUs will be able to test 64 to 256 points in parallel. Note however that transferring data from CPU to GPU and back is always expensive, so for just testing a couple of points against a couple of simple polygons, where either the points or the polygons are dynamic and will change frequently, a GPU approach will rarely pay off.

Answer 2

我认为以下代码是最好的解决方案(取自此处):

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表示奇数.

Answer 3

这是由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;
}

Answer 4

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

Answer 5

计算点p和每个多边形顶点之间的定向角度和.如果总取向角度是360度,则该点在内部.如果总数为0,则该点在外面.

我更喜欢这种方法,因为它更稳健,更少依赖于数值精度.

计算交叉点数均匀度的方法是有限的,因为您可以在计算交叉点数量时"击中"顶点.

编辑:通过The Way,这种方法适用于凹凸多边形.

编辑:我最近发现了一篇关于该主题的整篇维基百科文章.

Answer 6

bobobobo引用的Eric Haines文章非常出色.特别有趣的是比较算法性能的表格; 角度求和方法与其他方法相比非常糟糕.同样有趣的是,使用查找网格进一步将多边形细分为"in"和"out"扇区的优化可以使测试速度极快,即使在具有> 1000边的多边形上也是如此.

无论如何,这是早期的,但我的投票是针对"交叉"方法,这正是Mecki所描述的.然而,我发现它最由大卫伯克描述和编纂.我喜欢没有真正的三角函数,它适用于凸面和凹面,并且随着边数的增加,它表现得相当好.

顺便说一句,这是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

Answer 7

这个问题非常有趣.我有另一个可行的想法不同于这篇文章的其他答案.

我们的想法是使用角度之和来确定目标是在内部还是外部.如果目标位于区域内,则目标和每两个边界点形成的角度之和将为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

Answer 8

这个问题的大多数答案都没有很好地处理所有极端情况。一些微妙的极端情况如下： 这是一个 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
}

Answer 9

nirg的答案的 Swift版本:

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

Answer 10

真的很喜欢由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;
}

Answer 11

最简单的解决方案是将多边形划分为三角形并按此处所述对三角形进行命中测试

如果您的多边形是凸多边形，则可能有更好的方法。将多边形视为无限线的集合。每条线将空间一分为二。对于每个点，很容易判断它是在线的一侧还是另一侧。如果一个点位于所有线的同一侧，则它位于多边形内部。

Answer 12

当我是Michael Stonebraker的研究员时,我在这方面做了一些工作- 你知道,这位教授想出了Ingres,PostgreSQL等.

我们意识到最快的方法是首先做一个边界框,因为它的速度非常快.如果它在边界框之外,它就在外面.否则,你会做更艰苦的工作......

如果你想要一个很好的算法,请查看开源项目PostgreSQL的地理工作源代码......

我想指出,我们从来没有对左右手的任何见解(也可以表达为"内部"与"外部"问题......

UPDATE

BKB的链接提供了大量合理的算法.我正在研究地球科学问题,因此需要一个在纬度/经度上工作的解决方案,它有一个特殊的手性问题 - 是较小区域内还是较大区域内的区域？答案是顶点的"方向"很重要 - 它既可以是左手也可以是右手,这样你就可以将任一区域指定为任何给定多边形的"内部".因此,我的工作使用了该页面上列举的解决方案三.

另外,我的工作使用了"在线"测试的单独功能.

...自从有人问:我们发现当顶点数量超过一定数量时,边界框测试最好 - 如果有必要,在进行更长时间的测试之前做一个非常快速的测试...通过简单地获取边界框来创建边界框最大x,最小x,最大y和最小y并将它们组合在一起制作一个盒子的四个点......

对于接下来的人的另一个提示:我们在网格空间中进行了所有更复杂和"光线调暗"的计算,所有这些都在平面上的正点上,然后重新投射回"真实"经度/纬度,从而避免了可能的错误当经过一条经线180和处理极地区域时,环绕着.工作得很好!

Answer 13

David Segond的答案几乎是标准的一般答案,而Richard T是最常见的优化,尽管其他一些优化.其他强有力的优化基于不太通用的解决方案.例如,如果要检查具有大量点的相同多边形,对多边形进行三角测量可以大大加快速度,因为有许多非常快速的TIN搜索算法.另一个是如果多边形和点在低分辨率的有限平面上,比如屏幕显示,您可以将多边形绘制到给定颜色的内存映射显示缓冲区中,并检查给定像素的颜色以查看它是否位于在多边形中.

与许多优化一样,这些优化基于特定情况而非一般情况,并且基于摊销时间而非单次使用产生收益.

在这个领域工作,我发现Joeseph O'Rourkes的计算几何在C'ISBN 0-521-44034-3中是一个很好的帮助.

Answer 14

爪哇版：

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

}

Answer 15

我已经用Python实现了nirg的c++代码：

输入

• bounding_points：构成多边形的节点。

• bounding_box_positions：要过滤的候选点。（在我的实现中，从边界框创建。

  （输入是格式为的元组列表[(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

同样，这个想法来自这里