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


当前回答

如果你正在使用谷歌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

其他回答

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

简单的解决方案是将多边形划分为三角形,并按这里解释的那样对三角形进行测试

如果你的多边形是凸多边形,可能有更好的方法。把这个多边形看作是无限条线的集合。每一行将空间一分为二。对于每一个点,很容易判断它是在直线的一边还是另一边。如果一个点在所有直线的同一侧,那么它在多边形内。

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

例如,假设我们有一个多边形,它带有多边形凹,看起来像这样:

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

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

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

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

基于仿真的简化算法在http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

一些helper谓词:

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

给定两点a和B的直线(直线(a,B))方程为:

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