如何确定二维点是否在多边形内？

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

这是针对任何点运行的第一个测试。如您所见，此测试非常快，但也非常粗糙。为了处理边界矩形内的点，我们需要更复杂的算法。有几种计算方法。哪种方法有效还取决于多边形是否可以有孔或始终是实心的。以下是实心的示例（一凸一凹）：

https://i.stack.imgur.com/G76ta.jpg

这是一个有洞的：

https://i.stack.imgur.com/oDCCd.gif

绿色的中间有个洞！

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

https://i.stack.imgur.com/2EzSM.png

缠绕数算法将是一种替代方法，对于非常接近多边形线的点来说它更准确，但它也慢得多。由于有限的浮点精度和舍入问题，对于太靠近多边形边的点，光线投射可能会失败，但实际上这几乎不是问题，就好像一个点离边很近，通常在视觉上甚至不可能查看器来识别它是已经在里面还是仍然在外面。

你还有上面的边界框，记得吗？只需在边界框外选择一个点并将其用作射线的起点。例如，点 (Xmin - e/p.y) 肯定在多边形之外。

但什么是 e？嗯，e（实际上是 epsilon）为边界框提供了一些 填充。正如我所说，如果我们开始离多边形线太近，光线追踪就会失败。由于边界框可能等于多边形（如果多边形是轴对齐的矩形，则边界框等于多边形本身！），我们需要一些填充来确保安全，仅此而已。您应该选择多大的 e？不会太大。这取决于您用于绘图的坐标系比例。如果您的像素步长为 1.0，则只需选择 1.0（但 0.1 也可以）

现在我们有了射线及其起点和终点坐标，问题从“是多边形内的点”转移到“射线与多边形边相交的频率”。因此，我们不能像以前那样只使用多边形点，现在我们需要实际的边。边总是由两点定义。

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

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

// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
    // Test if current side intersects with ray.
    // If yes, intersections++;
}
if ((intersections & 1) == 1) {
    // Inside of polygon
} else {
    // Outside of polygon
}

到目前为止一切顺利，但是如何测试两个向量是否相交？这是一些 C 代码（未经测试），应该可以解决问题：

#define NO 0
#define YES 1
#define COLLINEAR 2

int areIntersecting(
    float v1x1, float v1y1, float v1x2, float v1y2,
    float v2x1, float v2y1, float v2x2, float v2y2
) {
    float d1, d2;
    float a1, a2, b1, b2, c1, c2;

    // Convert vector 1 to a line (line 1) of infinite length.
    // We want the line in linear equation standard form: A*x + B*y + C = 0
    // See: http://en.wikipedia.org/wiki/Linear_equation
    a1 = v1y2 - v1y1;
    b1 = v1x1 - v1x2;
    c1 = (v1x2 * v1y1) - (v1x1 * v1y2);

    // Every point (x,y), that solves the equation above, is on the line,
    // every point that does not solve it, is not. The equation will have a
    // positive result if it is on one side of the line and a negative one 
    // if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
    // 2 into the equation above.
    d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
    d2 = (a1 * v2x2) + (b1 * v2y2) + c1;

    // If d1 and d2 both have the same sign, they are both on the same side
    // of our line 1 and in that case no intersection is possible. Careful, 
    // 0 is a special case, that's why we don't test ">=" and "<=", 
    // but "<" and ">".
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // The fact that vector 2 intersected the infinite line 1 above doesn't 
    // mean it also intersects the vector 1. Vector 1 is only a subset of that
    // infinite line 1, so it may have intersected that line before the vector
    // started or after it ended. To know for sure, we have to repeat the
    // the same test the other way round. We start by calculating the 
    // infinite line 2 in linear equation standard form.
    a2 = v2y2 - v2y1;
    b2 = v2x1 - v2x2;
    c2 = (v2x2 * v2y1) - (v2x1 * v2y2);

    // Calculate d1 and d2 again, this time using points of vector 1.
    d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
    d2 = (a2 * v1x2) + (b2 * v1y2) + c2;

    // Again, if both have the same sign (and neither one is 0),
    // no intersection is possible.
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // If we get here, only two possibilities are left. Either the two
    // vectors intersect in exactly one point or they are collinear, which
    // means they intersect in any number of points from zero to infinite.
    if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;

    // If they are not collinear, they must intersect in exactly one point.
    return YES;
}

输入值是向量 1（v1x1/v1y1 和 v1x2/v1y2）和向量 2（v2x1/v2y1 和 v2x2/v2y2）的两个端点。所以你有 2 个向量，4 个点，8 个坐标。 YES 和 NO 是明确的。 YES 增加交叉点，NO 什么也不做。

COLLINEAR 怎么样？这意味着两个向量都位于同一条无限线上，取决于位置和长度，它们根本不相交或相交于无数个点。我不确定如何处理这种情况，无论哪种方式，我都不会将其视为交集。好吧，由于浮点舍入错误，这种情况在实践中相当罕见；更好的代码可能不会测试 == 0.0f，而是测试 < epsilon 之类的东西，其中 epsilon 是一个相当小的数字。

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

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

当然，这种方法仅在您的绘图表面不必很大时才可用。如果它无法适应 GPU 内存，则此方法比在 CPU 上执行要慢。如果它必须很大并且您的 GPU 支持现代着色器，您仍然可以通过将上面显示的光线投射实现为 GPU 着色器来使用 GPU，这绝对是可能的。对于要测试的大量多边形或大量点，这将获得回报，考虑一些 GPU 将能够并行测试 64 到 256 个点。但是请注意，将数据从 CPU 传输到 GPU 并返回总是很昂贵，因此仅针对几个简单的多边形测试几个点，其中点或多边形是动态的并且会经常变化，GPU 方法很少会支付离开。

---

我认为以下代码是最好的解决方案（取自 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;
}

论据

nvert：多边形中的顶点数。是否重复最后的第一个顶点已经在上面提到的文章中讨论过了。

vertx, verty：包含多边形顶点的 x 和 y 坐标的数组。

testx, testy：测试点的 X 和 y 坐标。

它既短又高效，适用于凸多边形和凹多边形。如前所述，您应该首先检查边界矩形并分别处理多边形孔。

这背后的想法很简单。作者是这样描述的：

我从测试点水平射出一条半无限光线（增加 x，固定 y），并计算它穿过多少条边。在每个交叉点，光线在内部和外部之间切换。这被称为乔丹曲线定理。

每次水平光线穿过任何边缘时，变量 c 都会从 0 切换到 1 和 1 到 0。所以基本上它会跟踪交叉的边数是偶数还是奇数。 0 表示偶数，1 表示奇数。

---

这是来自 this RPI professor 的 answer given by nirg 的 C# 版本。请注意，使用来自该 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;
}

---

这是 M. Katz 基于 Nirg 方法的答案的 JavaScript 变体：

function pointIsInPoly(p, polygon) {
    var isInside = false;
    var minX = polygon[0].x, maxX = polygon[0].x;
    var minY = polygon[0].y, maxY = polygon[0].y;
    for (var n = 1; n < polygon.length; n++) {
        var q = polygon[n];
        minX = Math.min(q.x, minX);
        maxX = Math.max(q.x, maxX);
        minY = Math.min(q.y, minY);
        maxY = Math.max(q.y, maxY);
    }

    if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
        return false;
    }

    var i = 0, j = polygon.length - 1;
    for (i, j; i < polygon.length; j = i++) {
        if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
                p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
            isInside = !isInside;
        }
    }

    return isInside;
}

---

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

我更喜欢这种方法，因为它更健壮且对数值精度的依赖更少。

计算交叉点数量均匀性的方法受到限制，因为您可以在计算交叉点数量期间“命中”顶点。

编辑：顺便说一句，这种方法适用于凹凸多边形。

编辑：我最近在该主题上找到了完整的 Wikipedia article。

---

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

设 x 为目标点。令数组 [0, 1, .... n] 为该区域的所有点。用一条线将目标点与每个边界点连接起来。如果目标点在该区域内。所有角度的总和将是 360 度。如果不是，角度将小于 360。

https://i.stack.imgur.com/Zbcec.jpg

我的算法假设顺时针方向是正方向。这是一个潜在的输入：

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

---

bobobobo 引用的Eric Haines article 非常出色。特别有趣的是比较算法性能的表格；与其他方法相比，角度求和方法确实很糟糕。同样有趣的是，使用查找网格将多边形进一步细分为“输入”和“输出”扇区等优化可以使测试速度非常快，即使在具有 > 的多边形上也是如此。 1000 面。

无论如何，现在还为时过早，但我的投票是“交叉”方法，这几乎就是 Mecki 所描述的我的想法。但是我发现它最简洁described and codified by 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

---

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

---

当我还是 Michael Stonebraker 的研究员时，我在这方面做了一些工作 - 你知道，提出 Ingres、PostgreSQL 等的教授。

我们意识到最快的方法是先做一个边界框，因为它非常快。如果它在边界框之外，它就在外面。否则，你会更努力地工作......

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

我想指出，我们从来没有深入了解过左撇子和右撇子（也可以表达为“内部”与“外部”问题......

更新

BKB 的链接提供了大量合理的算法。我正在研究地球科学问题，因此需要一个适用于纬度/经度的解决方案，它有一个特殊的惯用手问题——是较小区域内的区域还是较大区域内的区域？答案是顶点的“方向”很重要——它是左撇子或右撇子，通过这种方式，您可以将任一区域指示为任何给定多边形的“内部”。因此，我的工作使用了该页面上列举的解决方案三。

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

...既然有人问：我们发现当顶点的数量超过某个数字时，边界框测试是最好的——如果需要，在进行更长的测试之前做一个非常快速的测试......一个边界框是通过简单地获取最大的 x、最小的 x、最大的 y 和最小的 y 并将它们放在一起形成一个盒子的四个点...

给后面的人的另一个提示：我们在网格空间中进行了所有更复杂和“调光”的计算，所有这些计算都在平面上的正点上，然后重新投影回“真实”经度/纬度，从而避免可能的错误当一根越过经度线 180 和处理极地地区时环绕。工作得很好！

---

https://i.stack.imgur.com/20cli.png

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

---

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

像许多优化一样，这些优化基于特定而不是一般情况，并且基于摊销时间而不是单一使用来产生收益。

在这个领域工作时，我发现 Joeseph O'Rourkes 'Computation Geometry in C' ISBN 0-521-44034-3 对我很有帮助。

---

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

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

---

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

}

---

我意识到这是旧的，但这里有一个用 Cocoa 实现的光线投射算法，以防有人感兴趣。不确定这是最有效的做事方式，但它可能会帮助某人。

- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
    NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
    BOOL result;
    float aggregateX = 0; //I use these to calculate the centroid of the shape
    float aggregateY = 0;
    NSPoint firstPoint[1];
    [currentPath elementAtIndex:0 associatedPoints:firstPoint];
    float olderX = firstPoint[0].x;
    float olderY = firstPoint[0].y;
    NSPoint interPoint;
    int noOfIntersections = 0;

    for (int n = 0; n < [currentPath elementCount]; n++) {
        NSPoint points[1];
        [currentPath elementAtIndex:n associatedPoints:points];
        aggregateX += points[0].x;
        aggregateY += points[0].y;
    }

    for (int n = 0; n < [currentPath elementCount]; n++) {
        NSPoint points[1];

        [currentPath elementAtIndex:n associatedPoints:points];
        //line equations in Ax + By = C form
        float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;  
        float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
        float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);

        float _A_BAR = olderY - points[0].y;
        float _B_BAR = points[0].x - olderX;
        float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);

        float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
        if (det != 0) {
            //intersection points with the edges
            float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
            float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
            interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
            if (olderX <= points[0].x) {
                //doesn't matter in which direction the ray goes, so I send it right-ward.
                if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {  
                    noOfIntersections++;
                }
            } else {
                if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
                     noOfIntersections++;
                } 
            }
        }
        olderX = points[0].x;
        olderY = points[0].y;
    }
    if (noOfIntersections % 2 == 0) {
        result = FALSE;
    } else {
        result = TRUE;
    }
    return result;
}

---

nirg 答案的 Obj-C 版本，带有用于测试点的示例方法。 Nirg 的回答对我来说效果很好。

- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
    NSUInteger nvert = [vertices count];
    NSInteger i, j, c = 0;
    CGPoint verti, vertj;

    for (i = 0, j = nvert-1; i < nvert; j = i++) {
        verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
        vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
        if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
        ( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
            c = !c;
    }

    return (c ? YES : NO);
}

- (void)testPoint {

    NSArray *polygonVertices = [NSArray arrayWithObjects:
        [NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
        [NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
        [NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
        [NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
        [NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
        [NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
        nil
    ];

    CGPoint tappedPoint = CGPointMake(23.0, 70.0);

    if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
        NSLog(@"YES");
    } else {
        NSLog(@"NO");
    }
}

https://i.stack.imgur.com/vCFl0.png

---

没有什么比对问题的归纳定义更美好的了。为了完整起见，这里有一个 prolog 版本，它也可能阐明光线投射背后的想法：

基于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。由于这是一个更强的限制，它需要首先检查，所以我们首先只取那些满足这个要求的行，然后检查它的位置。根据 Jordan 曲线定理，任何投射到多边形的射线都必须在偶数条线上相交。所以我们完成了，我们将把光线向右抛，然后每次它与一条线相交时，切换它的状态。然而，在我们的实现中，我们将检查满足给定限制的解决方案包的长度并决定其内部结构。对于多边形中的每条线，都必须这样做。

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

---

我制作了 nirg's c++ code 的 Python 实现：

输入

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

同样，这个想法来自 here

---

C# 版本的 nirg 的答案在这里：我将分享代码。它可能会节省一些时间。

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

---

.Net 端口：

    static void Main(string[] args)
    {

        Console.Write("Hola");
        List<double> vertx = new List<double>();
        List<double> verty = new List<double>();

        int i, j, c = 0;

        vertx.Add(1);
        vertx.Add(2);
        vertx.Add(1);
        vertx.Add(4);
        vertx.Add(4);
        vertx.Add(1);

        verty.Add(1);
        verty.Add(2);
        verty.Add(4);
        verty.Add(4);
        verty.Add(1);
        verty.Add(1);

        int nvert = 6;  //Vértices del poligono

        double testx = 2;
        double testy = 5;


        for (i = 0, j = nvert - 1; i < nvert; j = i++)
        {
            if (((verty[i] > testy) != (verty[j] > testy)) &&
             (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
                c = 1;
        }
    }

---

VBA 版本：

注意：请记住，如果您的多边形是地图中的一个区域，则纬度/经度是 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

---

很惊讶之前没有人提出这个问题，但是对于需要数据库的实用主义者来说：MongoDB 对包括这个在内的地理查询有很好的支持。

您正在寻找的是：

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

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

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

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

---

这是 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 倍）。

---

如果您使用 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

---

这可能是来自 from this page 的 here 的 C 代码的优化程度稍差的版本。

我的 C++ 版本使用 std::vector<std::pair<double, double>> 和两个 double 作为 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;
}

---

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

如果 Point 在多边形区域内。 （可以通过光线投射算法完成）如果点在多边形边界上（可以通过与折线（线）上的点检测相同的算法完成）。

---

在 Ray casting algorithm 中处理以下特殊情况：

光线与多边形的一侧重叠。该点位于多边形内部，光线穿过多边形的一个顶点。该点位于多边形之外，光线刚好接触到多边形的一个角度。

检查Determining Whether A Point Is Inside A Complex Polygon。本文提供了一种简单的方法来解决它们，因此上述情况不需要特殊处理。

---

您可以通过检查通过将所需点连接到多边形的顶点形成的区域是否与多边形本身的区域匹配来做到这一点。

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

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

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

---

如果您正在寻找一个 java 脚本库，那么 Polygon 类有一个 javascript google maps v3 扩展来检测一个点是否位于其中。

var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);

Google Extention Github

---

使用 qt (Qt 4.3+) 时，可以使用 QPolygon 的函数 containsPoint

---