如何确定一个二维点是否在多边形内？（涉及IT技术）

Question

如何确定一个二维点是否在多边形内？（涉及IT技术）

performancegraphicscollision-detectionpolygonpoint-in-polygon

605

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

- Scott Evernden

您忘了告诉我们关于左右利手的问题您的看法——这也可以被解释为“内部”与“外部”。 - Richard T

17

是的，我现在意识到这个问题留下了许多未指定的细节，但目前我有点想看看各种回答的多样性。 - Scott Evernden

5

一个拥有90个面的多边形被称为九十边形（或称九角形），一个拥有10,000个面的多边形被称为万边形。 - user263678

1

“最优雅”已经不是目标了，因为我一直在寻找一个“能够工作”的算法。我必须自己想出来：https://dev59.com/GmUq5IYBdhLWcg3wEcVZ#18190354 - davidkonrad

这个库已经实现了它，所以你只需要在 Python 中使用 point.within(polygon)，它会返回一个布尔值，表示点是否在多边形内部。 - J Agustin Barrachina

42个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Aladar · Answer 1

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

- Santiago M. Quintero · Answer 2

惊讶的是之前没有人提到，但对于需要数据库的实用主义者来说：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个点的性能，甚至没有任何优化并且包括更新其各自多边形所需的时间，少于一分钟。

- Yuan Fu · Answer 3

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

- ideasman42 · Answer 4

这是一个不使用射线法的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倍）。

- Shaun Han · Answer 5

这是另一种基于NumPy的实现，我认为它是目前所有答案中最简洁的一个。

例如，假设我们有一个带空洞的多边形，看起来像这样：

大多边形顶点的二维坐标为

[[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] 在其中。

我的解決方案基於 winding number algorithm。以下是我使用僅 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]

- TankorSmash · Answer 6

这是一个较不优化的C代码版本，此处提供了源代码，此页面提供了来源。

我的C ++版本使用std :: vector<std :: pair <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;
}

- Shanaka Rathnayaka · Answer 7

如果你正在寻找一个JavaScript库，那么有一个JavaScript Google Maps V3扩展可以用于Polygon类，以便检测一个点是否在其中。

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

Google Extension Github

- Michael-7 · Answer 8

以下是 nirg 提供的 Scala 版本解决方案（假设边界矩形预检已单独完成）：

def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {

  val length = polygon.length

  @tailrec
  def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
    if (i == length)
      tracker
    else {
      val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
      oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
    }
  }

  oddIntersections(0, length - 1, tracker = false)
}

- SamTech · Answer 9

这是@girg答案的golang版本（受@m-katz的C#代码启发）

func isPointInPolygon(polygon []point, testp point) bool {
    minX := polygon[0].X
    maxX := polygon[0].X
    minY := polygon[0].Y
    maxY := polygon[0].Y

    for _, p := range polygon {
        minX = min(p.X, minX)
        maxX = max(p.X, maxX)
        minY = min(p.Y, minY)
        maxY = max(p.Y, maxY)
    }

    if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
        return false
    }

    inside := false
    j := len(polygon) - 1
    for i := 0; i < len(polygon); i++ {
        if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
            inside = !inside
        }
        j = i
    }

    return inside
}

- Michel Rouzic · Answer 10

就像David Segonds' answer所建议的那样，我使用了一种角度求和方法，该方法源自我的凹多边形绘制算法。它依赖于将围绕点的子三角形的近似角度相加以获得权重。权重约为1.0表示该点在三角形内部，约为0.0表示在外部，约为-1.0的权重是发生在多边形内但顺序相反的情况（例如菱形的两半之一），如果恰好位于边缘，则为NAN。它之所以不慢，是因为角度根本不需要精确估计。通过将孔视为单独的多边形并减去权重来处理它们。

typedef struct { double x, y; } xy_t;

xy_t sub_xy(xy_t a, xy_t b)
{
    a.x -= b.x;
    a.y -= b.y;
    return a;
}

double calc_sharp_subtriangle_pixel_weight(xy_t p0, xy_t p1)
{
    xy_t rot, r0, r1;
    double weight;

    // Rotate points (unnormalised)
    rot = sub_xy(p1, p0);
    r0.x = rot.x*p0.y - rot.y*p0.x;
    r0.y = rot.x*p0.x + rot.y*p0.y;
    r1.y = rot.x*p1.x + rot.y*p1.y;

    // Calc weight
    weight = subtriangle_angle_approx(r1.y, r0.x) - subtriangle_angle_approx(r0.y, r0.x);

    return weight;
}

double calc_sharp_polygon_pixel_weight(xy_t p, xy_t *corner, int corner_count)
{
    int i;
    xy_t p0, p1;
    double weight = 0.;

    p0 = sub_xy(corner[corner_count-1], p);
    for (i=0; i < corner_count; i++)
    {
        // Transform corner coordinates
        p1 = sub_xy(corner[i], p);

        // Calculate weight for each subtriangle
        weight += calc_sharp_subtriangle_pixel_weight(p0, p1);
        p0 = p1;
    }

    return weight;
}

因此，对于多边形的每个部分，将形成一个与待评估点相关的子三角形，然后将旋转每个子三角形以评估其近似角度并加入权重。

可以用 atan2(y, x) / (2.*pi) 替换对 subtriangle_angle_approx(y, x) 的调用，但是粗略的近似就足够准确了：

double subtriangle_angle_approx(double y, double x)
{
    double angle, d;
    int obtuse;

    if (x == 0.)
        return NAN;

    obtuse = fabs(y) > fabs(x);
    if (obtuse)
        swap_double(&y, &x);

    // Core of the approximation, a very loosely approximate atan(y/x) / (2.*pi) over ]-1 , 1[
    d = y / x;
    angle = 0.13185 * d;

    if (obtuse)
        angle = sign(d)*0.25 - angle;

    return angle;
}