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

发布于 2024-07-08 06:18:56 字数 106 浏览 16 评论 0 原文

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

I'm trying to create a fast 2D point inside polygon algorithm, for use in hit-testing (e.g. Polygon.contains(p:Point)). Suggestions for effective techniques would be appreciated.

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(30

飘然心甜 2024-07-15 06:18:57

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

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

I think the following piece of code is the best solution (taken from here):

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

Arguments

  • nvert: Number of vertices in the polygon. Whether to repeat the first vertex at the end has been discussed in the article referred above.
  • vertx, verty: Arrays containing the x- and y-coordinates of the polygon's vertices.
  • testx, testy: X- and y-coordinate of the test point.

It's both short and efficient and works both for convex and concave polygons. As suggested before, you should check the bounding rectangle first and treat polygon holes separately.

The idea behind this is pretty simple. The author describes it as follows:

I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point, and count how many edges it crosses. At each crossing, the ray switches between inside and outside. This is called the Jordan curve theorem.

The variable c is switching from 0 to 1 and 1 to 0 each time the horizontal ray crosses any edge. So basically it's keeping track of whether the number of edges crossed are even or odd. 0 means even and 1 means odd.

罪歌 2024-07-15 06:18:57

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

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

Here is a C# version of the answer given by nirg, which comes from this RPI professor. Note that use of the code from that RPI source requires attribution.

A bounding box check has been added at the top. However, as James Brown points out, the main code is almost as fast as the bounding box check itself, so the bounding box check can actually slow the overall operation, in the case that most of the points you are checking are inside the bounding box. So you could leave the bounding box check out, or an alternative would be to precompute the bounding boxes of your polygons if they don't change shape too often.

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;
}
东北女汉子 2024-07-15 06:18:57

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

Here is a JavaScript variant of the answer by M. Katz based on Nirg's approach:

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;
}
情归归情 2024-07-15 06:18:57

计算点 p 与每个多边形顶点之间的角度的定向和。 如果总定向角度为 360 度,则该点位于内部。 如果总数为 0,则该点位于外部。

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

计算交叉点数量均匀度的方法受到限制,因为在计算交叉点数量期间您可能会“击中”顶点。

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

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

Compute the oriented sum of angles between the point p and each of the polygon apices. If the total oriented angle is 360 degrees, the point is inside. If the total is 0, the point is outside.

I like this method better because it is more robust and less dependent on numerical precision.

Methods that compute evenness of number of intersections are limited because you can 'hit' an apex during the computation of the number of intersections.

EDIT: By The Way, this method works with concave and convex polygons.

EDIT: I recently found a whole Wikipedia article on the topic.

铁轨上的流浪者 2024-07-15 06:18:57

这个问题太有趣了。 我有另一个可行的想法,不同于这篇文章的其他答案。 这个想法是使用角度之和来确定目标是在内部还是在外部。 更广为人知的名称是绕数

设 x 为目标点。 令数组 [0, 1, .... n] 为该区域的所有点。 用一条线将目标点与每个边界点连接起来。 如果目标点在该区域内。 所有角度的总和将为 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

This question is so interesting. I have another workable idea different from other answers to this post. The idea is to use the sum of angles to decide whether the target is inside or outside. Better known as winding number.

Let x be the target point. Let array [0, 1, .... n] be the all the points of the area. Connect the target point with every border point with a line. If the target point is inside of this area. The sum of all angles will be 360 degrees. If not the angles will be less than 360.

Refer to this image to get a basic understanding of the idea:
enter image description here

My algorithm assumes the clockwise is the positive direction. Here is a potential input:

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

The following is the python code that implements the idea:

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
音栖息无 2024-07-15 06:18:57

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

The Eric Haines article cited by bobobobo is really excellent. Particularly interesting are the tables comparing performance of the algorithms; the angle summation method is really bad compared to the others. Also interesting is that optimisations like using a lookup grid to further subdivide the polygon into "in" and "out" sectors can make the test incredibly fast even on polygons with > 1000 sides.

Anyway, it's early days but my vote goes to the "crossings" method, which is pretty much what Mecki describes I think. However I found it most succintly described and codified by David Bourke. I love that there is no real trigonometry required, and it works for convex and concave, and it performs reasonably well as the number of sides increases.

By the way, here's one of the performance tables from the Eric Haines' article for interest, testing on random polygons.

                       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
咽泪装欢 2024-07-15 06:18:57

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

Most of the answers in this question are not handling all corner cases well. Some subtle corner cases like below:
ray casting corner cases
This is a javascript version with all corner cases well handled.

/** 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
}
屋檐 2024-07-15 06:18:57

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

Really like the solution posted by Nirg and edited by bobobobo. I just made it javascript friendly and a little more legible for my use:

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;
}
燕归巢 2024-07-15 06:18:57

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

Swift version of the answer by 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
    }
}
耳根太软 2024-07-15 06:18:57

当我还是 Michael Stonebraker 的研究员时,我在这方面做了一些工作 - 你知道,教授谁想出了安格尔PostgreSQL 等。

我们意识到最快的方法是首先做一个边界框,因为它超级快。 如果它在边界框之外,则它在边界框之外。 否则,你会做更困难的工作...

如果你想要一个出色的算法,请查看开源项目 PostgreSQL 源代码以进行地理工作...

我想指出,我们从未深入了解右手与左手习惯(也可以表示为“内部”与“外部”问题...


更新

BKB 的链接提供了大量合理的算法,我正在研究地球科学问题,因此需要一个适用于纬度/经度的解决方案,并且它具有独特的功能。旋向问题 - 较小区域内的区域还是较大区域内的区域?答案是顶点的“方向”很重要 - 它是左旋或右旋,这样您就可以将任一区域指示为“内部”因此,我的工作使用了该页面上列举的解决方案三。

此外,我的工作使用了单独的函数进行“在线”测试

……因为有人问:我们发现边界框测试在什么情况下是最好的。顶点的数量超出了某个数字 - 如果需要的话,在进行更长的测试之前先进行一次非常快速的测试...通过简单地取最大的 x、最小的 x、最大的 y 和最小的 y 并将它们放在一起来创建边界框一个盒子的四个点...

给后面的人的另一个提示:我们在网格空间中进行了所有更复杂和“调光”的计算,所有计算都在平面上的正点上,然后重新投影回“真实”经度/纬度,从而避免了当跨越经度 180 度线以及处理极地区域时可能出现的绕行错误。 效果很好!

I did some work on this back when I was a researcher under Michael Stonebraker - you know, the professor who came up with Ingres, PostgreSQL, etc.

We realized that the fastest way was to first do a bounding box because it's SUPER fast. If it's outside the bounding box, it's outside. Otherwise, you do the harder work...

If you want a great algorithm, look to the open source project PostgreSQL source code for the geo work...

I want to point out, we never got any insight into right vs left handedness (also expressible as an "inside" vs "outside" problem...


UPDATE

BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.

In addition, my work used separate functions for "on the line" tests.

...Since someone asked: we figured out that bounding box tests were best when the number of verticies went beyond some number - do a very quick test before doing the longer test if necessary... A bounding box is created by simply taking the largest x, smallest x, largest y and smallest y and putting them together to make four points of a box...

Another tip for those that follow: we did all our more sophisticated and "light-dimming" computing in a grid space all in positive points on a plane and then re-projected back into "real" longitude/latitude, thus avoiding possible errors of wrapping around when one crossed line 180 of longitude and when handling polar regions. Worked great!

一身骄傲 2024-07-15 06:18:57

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

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

The trivial solution would be to divide the polygon to triangles and hit test the triangles as explained here

If your polygon is CONVEX there might be a better approach though. Look at the polygon as a collection of infinite lines. Each line dividing space into two. for every point it's easy to say if its on the one side or the other side of the line. If a point is on the same side of all lines then it is inside the polygon.

以为你会在 2024-07-15 06:18:57

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

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

在这个领域工作时,我发现 Joeseph O'Rourkes 的《C 语言计算几何》ISBN 0-521-44034-3 给了我很大的帮助。

David Segond's answer is pretty much the standard general answer, and Richard T's is the most common optimization, though therre are some others. Other strong optimizations are based on less general solutions. For example if you are going to check the same polygon with lots of points, triangulating the polygon can speed things up hugely as there are a number of very fast TIN searching algorithms. Another is if the polygon and points are on a limited plane at low resolution, say a screen display, you can paint the polygon onto a memory mapped display buffer in a given colour, and check the color of a given pixel to see if it lies in the polygons.

Like many optimizations, these are based on specific rather than general cases, and yield beneifits based on amortized time rather than single usage.

Working in this field, i found Joeseph O'Rourkes 'Computation Geometry in C' ISBN 0-521-44034-3 to be a great help.

古镇旧梦 2024-07-15 06:18:57

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

}

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

}
泪意 2024-07-15 06:18:57

我已经用 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

同样,这个想法取自此处

I've made a Python implementation of nirg's c++ code:

Inputs

  • bounding_points: nodes that make up the polygon.
  • bounding_box_positions: candidate points to filter. (In my implementation created from the bounding box.

    (The inputs are lists of tuples in the format: [(xcord, ycord), ...])

Returns

  • All the points that are inside the polygon.
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

Again, the idea is taken from here

小霸王臭丫头 2024-07-15 06:18:57

我意识到这已经很旧了,但这里有一个在 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;
}

I realize this is old, but here is a ray casting algorithm implemented in Cocoa, in case anyone is interested. Not sure it is the most efficient way to do things, but it may help someone out.

- (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;
}
烂人 2024-07-15 06:18:57

nilg 答案的 Obj-C 版本,带有测试点的示例方法。 尼尔格的回答对我来说效果很好。

- (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");
    }
}

样本多边形

Obj-C version of nirg's answer with sample method for testing points. Nirg's answer worked well for me.

- (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");
    }
}

sample polygon

守不住的情 2024-07-15 06:18:57

没有什么比问题的归纳定义更美妙的了。 为了完整起见,您在序言中提供了一个版本,它也可能澄清光线投射背后的想法:

基于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 (Line(A,B)) 的直线方程为:

                    (YB-YA)
           Y - YA = ------- * (X - XA) 
                    (XB-YB) 

重要的是,直线的旋转方向为
对于边界设置为顺时针,对于孔设置为逆时针。
我们要检查点(X,Y),即测试点是否在左边
我们生产线的半平面(这是一个品味问题,也可以是
右侧,而且边界线的方向也必须改变
这种情况),这是将光线从该点投影到右侧(或左侧)
并确认与线的交点。 我们选择了项目
水平方向的光线(这又是一个品味问题,
它也可以在垂直方向上完成,但有类似的限制),所以我们有:

               (XB-XA)
           X < ------- * (Y - YA) + XA
               (YB-YA) 

现在我们需要知道该点是否位于
只是线段,而不是整个平面,所以我们需要
将搜索限制为仅此段,但这很容易,因为
在线段内只有一个点可以更高
比垂直轴上的 Y 更重要。 由于这是一个更强的限制
需要第一个检查,所以我们只首先检查这些行
满足此要求,然后检查其位置。 约旦河边
曲线定理 投射到多边形上的任何射线必须相交于
偶数行。 这样我们就完成了,我们将把光线扔到
向右,然后每次它与一条线相交时,切换其状态。
然而,在我们的实现中,我们要检查的长度
满足给定限制的一袋解决方案并决定
其内在。 对于多边形中的每条线都必须完成此操作。

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).

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), i.e 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).
禾厶谷欠 2024-07-15 06:18:57

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

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;
        }
酒中人 2024-07-15 06:18:57

VBA 版本:

注意:请记住,如果您的多边形是地图内的一个区域,则纬度/经度是 Y/X 值,而不是 X/Y(纬度 = Y,经度 = X),因为我所理解的是很久以前经度还不是一种测量方法时的历史意义。

类别模块: 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

模块:

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

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(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
尤怨 2024-07-15 06:18:57

我认为这是迄今为止所有答案中最简洁的另一个 numpyic 实现。

例如,假设我们有一个带有多边形空心的多边形,如下所示:
输入图片这里的描述

大多边形顶点的 2D 坐标为

[[139, 483], [227, 792], [482, 849], [523, 670], [352, 330]]

空心方形顶点的坐标为

[[248, 518], [336, 510], [341, 614], [250, 620]]

空心三角形顶点的坐标为

[[416, 531], [505, 517], [495, 616]]

假设我们要测试两个点 [ 296, 557][422, 730] 如果它们位于红色区域内(不包括边缘)。 如果我们找到这两个点,它将如下所示:

[296, 557] 不在读取区域内,而 [422, 730] 则在读取区域内。

我的解决方案基于绕数算法。 下面是我仅使用 numpy 的 4 行 python 代码:

def detect(points, *polygons):
    import numpy as np
    endpoint1 = np.r_[tuple(np.roll(p, 1, 0) for p in polygons)][:, None] - points
    endpoint2 = np.r_[polygons][:, None] - points
    p1, p2 = np.cross(endpoint1, endpoint2), np.einsum('...i,...i', endpoint1, endpoint2)
    return ~((p1.sum(0) < 0) ^ (abs(np.arctan2(p1, p2).sum(0)) > np.pi) | ((p1 == 0) & (p2 <= 0)).any(0))

测试实现:

points = [[296, 557], [422, 730]]
polygon1 = [[139, 483], [227, 792], [482, 849], [523, 670], [352, 330]]
polygon2 = [[248, 518], [336, 510], [341, 614], [250, 620]]
polygon3 = [[416, 531], [505, 517], [495, 616]]

print(detect(points, polygon1, polygon2, polygon3))

输出:

[False  True]

Yet another numpyic implementation which I believe is the most concise one out of all the answers so far.

For example, let's say we have a polygon with polygon hollows that looks like this:
enter image description here

The 2D coordinates for the vertices of the large polygon are

[[139, 483], [227, 792], [482, 849], [523, 670], [352, 330]]

The coordinates for the vertices of the square hollow are

[[248, 518], [336, 510], [341, 614], [250, 620]]

The coordinates for the vertices of the triangle hollow are

[[416, 531], [505, 517], [495, 616]]

Say we want to test two points [296, 557] and [422, 730] if they are within the red area (excluding the edges). If we locate the two points, it will look like this:
enter image description here

Obviously, [296, 557] is not inside the read area, whereas [422, 730] is.

My solution is based on the winding number algorithm. Below is my 4-line python code using only numpy:

def detect(points, *polygons):
    import numpy as np
    endpoint1 = np.r_[tuple(np.roll(p, 1, 0) for p in polygons)][:, None] - points
    endpoint2 = np.r_[polygons][:, None] - points
    p1, p2 = np.cross(endpoint1, endpoint2), np.einsum('...i,...i', endpoint1, endpoint2)
    return ~((p1.sum(0) < 0) ^ (abs(np.arctan2(p1, p2).sum(0)) > np.pi) | ((p1 == 0) & (p2 <= 0)).any(0))

To test the implementation:

points = [[296, 557], [422, 730]]
polygon1 = [[139, 483], [227, 792], [482, 849], [523, 670], [352, 330]]
polygon2 = [[248, 518], [336, 510], [341, 614], [250, 620]]
polygon3 = [[416, 531], [505, 517], [495, 616]]

print(detect(points, polygon1, polygon2, polygon3))

Output:

[False  True]
飘过的浮云 2024-07-15 06:18:57

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

.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;
        }
    }
满地尘埃落定 2024-07-15 06:18:57

令人惊讶的是,之前没有人提出这一点,但对于需要数据库的实用主义者来说:MongoDB 对地理查询(包括这个查询)具有出色的支持。

您正在寻找的是:

db.neighborhoods.findOne({ 几何: { $geoIntersects: { $geometry: {
类型:“点”,坐标:[“经度”,“纬度”] } } }
})

Neighborhoods 是一个以标准 GeoJson 格式存储一个或多个多边形的集合。 如果查询返回 null,则表示不相交,否则表示相交。

这里有很好的记录:
https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/< /a>

在 330 个不规则多边形网格中分类的 6,000 多个点的性能不到一分钟,根本没有任何优化,包括用各自的多边形更新文档的时间。

Surprised nobody brought this up earlier, but for the pragmatists requiring a database: MongoDB has excellent support for Geo queries including this one.

What you are looking for is:

db.neighborhoods.findOne({ geometry: { $geoIntersects: { $geometry: {
type: "Point", coordinates: [ "longitude", "latitude" ] } } }
})

Neighborhoods is the collection that stores one or more polygons in standard GeoJson format. If the query returns null it is not intersected otherwise it is.

Very well documented here:
https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

The performance for more than 6,000 points classified in a 330 irregular polygon grid was less than one minute with no optimization at all and including the time to update documents with their respective polygon.

黑色毁心梦 2024-07-15 06:18:57

这是 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 倍)交集法)。

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).

最终幸福 2024-07-15 06:18:57

如果您使用 Google Map SDK 并且想要检查某个点是否位于多边形内,可以尝试使用 GMSGeometryContainsLocation。 效果很好! 这是它的工作原理,

if GMSGeometryContainsLocation(point, polygon, true) {
    print("Inside this polygon.")
} else {
    print("outside this polygon")
}

这是参考:https:// /developers.google.com/maps/documentation/ios-sdk/reference/group___geometry_utils#gaba958d3776d49213404af249419d0ffd

If you're using Google Map SDK and want to check if a point is inside a polygon, you can try to use GMSGeometryContainsLocation. It works great!! Here is how that works,

if GMSGeometryContainsLocation(point, polygon, true) {
    print("Inside this polygon.")
} else {
    print("outside this polygon")
}

Here is the reference: https://developers.google.com/maps/documentation/ios-sdk/reference/group___geometry_utils#gaba958d3776d49213404af249419d0ffd

半寸时光 2024-07-15 06:18:57

这可能是来自此处的 C 代码的优化版本,其来源为从此页面

我的 C++ 版本使用 std::vector> 和两个双精度数作为 x 和 y。 逻辑应该与原始 C 代码完全相同,但我发现我的代码更容易阅读。 对于表演,我无话可说。

bool point_in_poly(std::vector<std::pair<double, double>>& verts, double point_x, double point_y)
{
    bool in_poly = false;
    auto num_verts = verts.size();
    for (int i = 0, j = num_verts - 1; i < num_verts; j = i++) {
        double x1 = verts[i].first;
        double y1 = verts[i].second;
        double x2 = verts[j].first;
        double y2 = verts[j].second;

        if (((y1 > point_y) != (y2 > point_y)) &&
            (point_x < (x2 - x1) * (point_y - y1) / (y2 - y1) + x1))
            in_poly = !in_poly;
    }
    return in_poly;
}

原来的C代码是

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

This is a presumably slightly less optimized version of the C code from here which was sourced from this page.

My C++ version uses a std::vector<std::pair<double, double>> and two doubles as an x and y. The logic should be exactly the same as the original C code, but I find mine easier to read. I can't speak for the performance.

bool point_in_poly(std::vector<std::pair<double, double>>& verts, double point_x, double point_y)
{
    bool in_poly = false;
    auto num_verts = verts.size();
    for (int i = 0, j = num_verts - 1; i < num_verts; j = i++) {
        double x1 = verts[i].first;
        double y1 = verts[i].second;
        double x2 = verts[j].first;
        double y2 = verts[j].second;

        if (((y1 > point_y) != (y2 > point_y)) &&
            (point_x < (x2 - x1) * (point_y - y1) / (y2 - y1) + x1))
            in_poly = !in_poly;
    }
    return in_poly;
}

The original C code is

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;
}
感受沵的脚步 2024-07-15 06:18:57

在大多数情况下,这是一个比其他人的算法更快的算法。

它新颖而优雅。 我们花费O(n * log(n))时间构建一个表,该表将允许我们在O(log(n) + k)中测试多边形内的点时间。

使用扫描光束表对同一多边形进行多次检查时,您可以明显更快地获得结果,而不是光线追踪或角度。 我们必须预先构建一个扫描光束活动边缘表,这是大多数代码正在做的事情。

scanbeam 图

我们计算 y 方向上该位置的扫描光束和活动边缘。 我们创建一个按 y 分量排序的点列表,并针对两个事件迭代该列表。 Start-Y 和 End-Y,我们在处理列表时跟踪活动边缘。 我们按顺序处理事件,对于每个扫描光束,我们记录事件的 y 值以及每个事件的活动边缘(事件为 start-y 和 end-y),但我们仅在事件 y 与事件 y 不同时记录这些事件上次(因此在我们将事件点标记到表中之前,会处理事件点的所有内容)。

所以我们得到了我们的表:

  1. []
  2. p6p5, p6p7
  3. p6p5, p6p7, p2p3, p2p1
  4. p6p7, p5p4, p2p3, p3p1
  5. p7p8, p5p4, p2p3, p2p1
  6. p7p8, p5p4, p3p4, p2p1
  7. p7p8, p2p1
  8. p7p8, p1p0
  9. p8p0, p1p0
  10. []

构建此表后,执行工作的实际代码只需几行。

  • 注意:这里的代码使用复数值作为点。 因此,.real.x.imag.y
def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

我们对事件进行二分搜索以查找特定值的 actives_at_y。 对于该 y 处的所有活动,我们计算该点的特定 y 处的段 x 值。 每当这个 x 截距大于我们点的 x 分量时,我们就加 1。 然后我们将该总数除以 2。(这是奇偶填充规则,您可以轻松地将其适应任何其他填充规则)。

完整代码:


from bisect import bisect

def build_edge_list(polygon):
    edge_list = []
    for i in range(1, len(polygon)):
        if (polygon[i].imag, polygon[i].real) < (polygon[i - 1].imag, polygon[i - 1].real):
            edge_list.append((polygon[i], i))
            edge_list.append((polygon[i - 1], ~i))
        else:
            edge_list.append((polygon[i], ~i))
            edge_list.append((polygon[i - 1], i))

    def sort_key(e):
        return e[0].imag, e[0].real, ~e[1]

    edge_list.sort(key=sort_key)
    return edge_list


def build_scanbeam(edge_list):
    actives_table = []
    events = []
    y = -float("inf")
    actives = []
    for pt, index in edge_list:
        if y != pt.imag:
            actives_table.append(list(actives))
            events.append(y)
        if index >= 0:
            actives.append(index)
        else:
            actives.remove(~index)
        y = pt.imag
    return actives_table, events

def point_in_polygon(polygon, point):
    def x_intercept(e, y):
        pt0 = polygon[e-1]
        pt1 = polygon[e]
        if pt1.real - pt0.real == 0:
            return pt0.real
        m = (pt1.imag - pt0.imag) / (pt1.real - pt0.real)
        b = pt0.imag - (m * pt0.real)
        return (y - b) / m

    edge_list = build_edge_list(polygon)
    actives_table, events = build_scanbeam(edge_list)
    try:
        if len(point):
            return [point_in_scantable(actives_table, events, x_intercept, p) for p in point]
    except TypeError:
        return point_in_scantable(actives_table, events, x_intercept, point)

def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

如果我们忽略扫描表的 O(n*log(n)) 构建时间。 我们实际上是在 O(log(n) + k) 时间内查找内容。 其中,n 是多边形中线段数量的大小,k 是所述多边形中活动边的典型数量。 其他光线跟踪方法实际上需要 O(n) 时间。 每次我们检查一个点时,它都会迭代整个多边形。 因此,即使这种实现方式明显不理想,它也轻而易举地击败了其他所有方法。


有一些性能技巧可以实现,例如,我们可以将时间复杂度降低到 O(log(n) + log(k)) 时间。 为此,我们将在扫描线中实施 Bentley-Ottmann,而不是将交叉点作为不同的事件进行处理,而是在交叉点处分割线。 然后,我们还按 x 轴截距对活动边进行排序。 然后我们知道在扫描光束期间不会发生相交,并且因为我们对片段进行了排序(请注意在扫描光束内正确排序它们,即使它们从相同的初始点开始(您需要查看斜率,或者只比较片段的中点)然后,我们有一个排序的无交集的活动列表扫描束表,这意味着我们也可以对活动边缘列表进行二分搜索,尽管这听起来像是对 k 值进行大量工作。通常为 2 或可能 4。

此外,由于这基本上变成了一个查找表和一些 x 截距的最小计算,因此更可以使用 GPU 来完成,因此您不再需要循环遍历多边形。实际上,您可以使用 numpy 之类的工具批量计算这些点,同时执行所有计算可以提高性能。

Here's an algorithm faster than everybody else's algorithm for most cases.

It's new and elegant. We spend O(n * log(n)) time building a table that will allow us to test point-in-polygon in O(log(n) + k) time.

Rather than ray-tracing or angles, you can get significantly faster results for multiple checks of the same polygon using a scanbeam table. We have to prebuild a scanbeam active edge table which is what most of the code is doing.

scanbeam diagram

We calculate the scanbeam and the active edges for that position in the y-direction. We make a list of points sorted by their y-component and we iterate through this list, for two events. Start-Y and End-Y, we track the active edges as we process the list. We process the events in order and for each scanbeam we record the y-value of the event and the active edges at each event (events being start-y and end-y) but we only record these when our event-y is different than last time (so everything at the event point is processed before we mark it in our table).

So we get our table:

  1. []
  2. p6p5, p6p7
  3. p6p5, p6p7, p2p3, p2p1
  4. p6p7, p5p4, p2p3, p3p1
  5. p7p8, p5p4, p2p3, p2p1
  6. p7p8, p5p4, p3p4, p2p1
  7. p7p8, p2p1
  8. p7p8, p1p0
  9. p8p0, p1p0
  10. []

The actual code doing the work is just a couple lines after this table is built.

  • Note: the code here is using complex values as points. So .real is .x and .imag is .y.
def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

We binary search into the events to find the actives_at_y for the specific value. For all the actives at that y we calculate the segments x value at the specific y of our point. We add 1 each time this x-intercept is bigger than our point's x-component. We then mod that total by 2. (This is the even-odd fill rule, you can easily adapt this to any other fill rule).

Complete Code:


from bisect import bisect

def build_edge_list(polygon):
    edge_list = []
    for i in range(1, len(polygon)):
        if (polygon[i].imag, polygon[i].real) < (polygon[i - 1].imag, polygon[i - 1].real):
            edge_list.append((polygon[i], i))
            edge_list.append((polygon[i - 1], ~i))
        else:
            edge_list.append((polygon[i], ~i))
            edge_list.append((polygon[i - 1], i))

    def sort_key(e):
        return e[0].imag, e[0].real, ~e[1]

    edge_list.sort(key=sort_key)
    return edge_list


def build_scanbeam(edge_list):
    actives_table = []
    events = []
    y = -float("inf")
    actives = []
    for pt, index in edge_list:
        if y != pt.imag:
            actives_table.append(list(actives))
            events.append(y)
        if index >= 0:
            actives.append(index)
        else:
            actives.remove(~index)
        y = pt.imag
    return actives_table, events

def point_in_polygon(polygon, point):
    def x_intercept(e, y):
        pt0 = polygon[e-1]
        pt1 = polygon[e]
        if pt1.real - pt0.real == 0:
            return pt0.real
        m = (pt1.imag - pt0.imag) / (pt1.real - pt0.real)
        b = pt0.imag - (m * pt0.real)
        return (y - b) / m

    edge_list = build_edge_list(polygon)
    actives_table, events = build_scanbeam(edge_list)
    try:
        if len(point):
            return [point_in_scantable(actives_table, events, x_intercept, p) for p in point]
    except TypeError:
        return point_in_scantable(actives_table, events, x_intercept, point)

def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

If we ignore then O(n*log(n)) buildtime for the scantable. We're actually looking things up in O(log(n) + k) time. Where n is the size of the number of segments in the polygon and k is the typical number of active edges in said polygon. The other methods for raytracing actually needs O(n) time. Each time we check a point it iterates the entire polygon. So even with this notably suboptimal implementation of this, it beats all the rest of them hands down.


There's a few performance tricks that could be done, for example, we can lower the time complexity to O(log(n) + log(k)) time. To do this we would implement Bentley-Ottmann into the sweep line, and rather than processing the intersections as different events, we split the lines at the intersections. We then also sort the active edges by their x-intercepts. We then know that no intersections occur during a scanbeam and since we sorted our segments (taking care to order them correctly within the scanbeam even if they start at the same initial point (you need to look at the slopes, or just compare midpoints of the segments). We then have a sorted intersection-less actives lists scanbeam table which means we can binary search into active edge list as well. Though that sounds like a lot of work for a value of k which is going to be typically 2 or maybe 4.

Also, since this basically becomes a lookup table and some minimal calculations for x-intercept it's much more able to be done with a GPU. You don't need to loop over the polygon anymore. So you can actually mass calculate these points with things like numpy where you get performance boosts for doing all the calculations at once.

半边脸i 2024-07-15 06:18:57

为了检测多边形上的命中,我们需要测试两件事:

  1. 点是否在多边形区域内。 (可以通过光线投射算法来完成)
  2. 如果点位于多边形边界上(可以通过用于折线(线)上的点检测的相同算法来完成)。

For Detecting hit on Polygon we need to test two things:

  1. If Point is inside polygon area. (can be accomplished by Ray-Casting Algorithm)
  2. If Point is on the polygon border(can be accomplished by same algorithm which is used for point detection on polyline(line)).
泪痕残 2024-07-15 06:18:57

光线投射算法

  1. 射线与多边形的一侧重叠。
  2. 该点位于多边形内部,并且射线穿过多边形的顶点。
  3. 该点位于多边形之外,并且射线仅触及多边形的一个角度。

检查确定点是否在复杂多边形内。 本文提供了解决这些问题的简单方法,因此上述情况无需特殊处理。

To deal with the following special cases in Ray casting algorithm:

  1. The ray overlaps one of the polygon's side.
  2. The point is inside of the polygon and the ray passes through a vertex of the polygon.
  3. 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.

想挽留 2024-07-15 06:18:57

您可以通过检查将所需点连接到多边形顶点形成的面积是否与多边形本身的面积相匹配来完成此操作。

或者您可以检查从您的点到每对两个连续多边形顶点到您的检查点的内角总和是否为 360,但我感觉第一个选项更快,因为它不涉及除法或计算三角函数的反函数。

我不知道如果你的多边形里面有一个洞会发生什么,但在我看来,主要思想可以适应这种情况

你也可以在数学社区中发布问题。 我打赌他们有一百万种方法可以做到这一点

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

尬尬 2024-07-15 06:18:56

对于图形,我宁愿不喜欢整数。 许多系统使用整数进行 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 withouthole

这是一个有孔的:

带孔的多边形

绿色的中间有一个孔!

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

演示光线如何穿过多边形

绕数算法将是一种替代方案,它是对于非常接近折线的点来说更准确,但速度也慢得多。 由于有限的浮点精度和舍入问题,对于太靠近多边形边的点,光线投射可能会失败,但实际上这几乎不是问题,就好像一个点距离边很近,通常在视觉上甚至不可能观察者可以识别它是否已经在内部或仍在外部。

你还记得上面的边界框吗? 只需在边界框外选择一个点并将其用作光线的起点即可。 例如,点 (Xmin - e/py) 肯定位于多边形之外。

但是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/v1y1v1x2/v1y2 )和向量 2(v2x1/v2y1v2x2/v2y2)。 所以你有 2 个向量、4 个点、8 个坐标。 YESNO 很清楚。 YES 增加交叉点,NO 不执行任何操作。

共线怎么样? 这意味着两个向量位于同一条无限直线上,具体取决于位置和长度,它们根本不相交,或者相交于无数个点。 我不太确定如何处理这种情况,无论如何我都不会将其算作交叉点。 嗯,无论如何,由于浮点舍入错误,这种情况在实践中相当罕见; 更好的代码可能不会测试 == 0.0f ,而是测试类似 的内容。 epsilon,其中 epsilon 是一个相当小的数字。

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

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

当然,只有当您的绘图表面不必很大时,此方法才可用。 如果 GPU 内存无法容纳,则此方法比在 CPU 上执行速度慢。 如果它必须很大并且您的 GPU 支持现代着色器,您仍然可以通过将上面所示的光线投射实现为 GPU 着色器来使用 GPU,这绝对是可能的。 对于大量的多边形或大量的点进行测试,这将会带来回报,考虑到某些 GPU 将能够并行测试 64 到 256 个点。 但请注意,将数据从 CPU 传输到 GPU 并返回的成本总是昂贵的,因此,如果仅针对几个简单的多边形测试几个点(其中点或多边形是动态的并且会频繁变化),GPU 方法很少会带来回报离开。

For graphics, I'd rather not prefer integers. Many systems use integers for UI painting (pixels are ints after all), but macOS, for example, uses float for everything. macOS only knows points and a point can translate to one pixel, but depending on monitor resolution, it might translate to something else. On retina screens half a point (0.5/0.5) is pixel. Still, I never noticed that macOS UIs are significantly slower than other UIs. After all, 3D APIs (OpenGL or Direct3D) also work with floats and modern graphics libraries very often take advantage of GPU acceleration.

Now you said speed is your main concern, okay, let's go for speed. Before you run any sophisticated algorithm, first do a simple test. Create an axis aligned bounding box around your polygon. This is very easy, fast and can already save you a lot of calculations. How does that work? Iterate over all points of the polygon and find the min/max values of X and Y.

E.g. you have the points (9/1), (4/3), (2/7), (8/2), (3/6). This means Xmin is 2, Xmax is 9, Ymin is 1 and Ymax is 7. A point outside of the rectangle with the two edges (2/1) and (9/7) cannot be within the polygon.

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

This is the first test to run for any point. As you can see, this test is ultra fast but it's also very coarse. To handle points that are within the bounding rectangle, we need a more sophisticated algorithm. There are a couple of ways how this can be calculated. Which method works also depends on whether the polygon can have holes or will always be solid. Here are examples of solid ones (one convex, one concave):

Polygon without hole

And here's one with a hole:

Polygon with hole

The green one has a hole in the middle!

The easiest algorithm, that can handle all three cases above and is still pretty fast is named ray casting. The idea of the algorithm is pretty simple: Draw a virtual ray from anywhere outside the polygon to your point and count how often it hits a side of the polygon. If the number of hits is even, it's outside of the polygon, if it's odd, it's inside.

Demonstrating how the ray cuts through a polygon

The winding number algorithm would be an alternative, it is more accurate for points being very close to a polygon line but it's also much slower. Ray casting may fail for points too close to a polygon side because of limited floating point precision and rounding issues, but in reality that is hardly a problem, as if a point lies that close to a side, it's often visually not even possible for a viewer to recognize if it is already inside or still outside.

You still have the bounding box of above, remember? Just pick a point outside the bounding box and use it as starting point for your ray. E.g. the point (Xmin - e/p.y) is outside the polygon for sure.

But what is e? Well, e (actually epsilon) gives the bounding box some padding. As I said, ray tracing fails if we start too close to a polygon line. Since the bounding box might equal the polygon (if the polygon is an axis aligned rectangle, the bounding box is equal to the polygon itself!), we need some padding to make this safe, that's all. How big should you choose e? Not too big. It depends on the coordinate system scale you use for drawing. If your pixel step width is 1.0, then just choose 1.0 (yet 0.1 would have worked as well)

Now that we have the ray with its start and end coordinates, the problem shifts from "is the point within the polygon" to "how often does the ray intersects a polygon side". Therefore we can't just work with the polygon points as before, now we need the actual sides. A side is always defined by two points.

side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:

You need to test the ray against all sides. Consider the ray to be a vector and every side to be a vector. The ray has to hit each side exactly once or never at all. It can't hit the same side twice. Two lines in 2D space will always intersect exactly once, unless they are parallel, in which case they never intersect. However since vectors have a limited length, two vectors might not be parallel and still never intersect because they are too short to ever meet each other.

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

So far so well, but how do you test if two vectors intersect? Here's some C code (not tested), that should do the trick:

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

The input values are the two endpoints of vector 1 (v1x1/v1y1 and v1x2/v1y2) and vector 2 (v2x1/v2y1 and v2x2/v2y2). So you have 2 vectors, 4 points, 8 coordinates. YES and NO are clear. YES increases intersections, NO does nothing.

What about COLLINEAR? It means both vectors lie on the same infinite line, depending on position and length, they don't intersect at all or they intersect in an endless number of points. I'm not absolutely sure how to handle this case, I would not count it as intersection either way. Well, this case is rather rare in practice anyway because of floating point rounding errors; better code would probably not test for == 0.0f but instead for something like < epsilon, where epsilon is a rather small number.

If you need to test a larger number of points, you can certainly speed up the whole thing a bit by keeping the linear equation standard forms of the polygon sides in memory, so you don't have to recalculate these every time. This will save you two floating point multiplications and three floating point subtractions on every test in exchange for storing three floating point values per polygon side in memory. It's a typical memory vs computation time trade off.

Last but not least: If you may use 3D hardware to solve the problem, there is an interesting alternative. Just let the GPU do all the work for you. Create a painting surface that is off screen. Fill it completely with the color black. Now let OpenGL or Direct3D paint your polygon (or even all of your polygons if you just want to test if the point is within any of them, but you don't care for which one) and fill the polygon(s) with a different color, e.g. white. To check if a point is within the polygon, get the color of this point from the drawing surface. This is just a O(1) memory fetch.

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.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文