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


当前回答

我已经做了nirg的c++代码的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

同样,这个想法也是从这里得来的

其他回答

from typing import Iterable

def pnpoly(verts, x, y):
    #check if x and/or y is iterable
    xit, yit = isinstance(x, Iterable), isinstance(y, Iterable)
    #if not iterable, make an iterable of length 1
    X = x if xit else (x, )
    Y = y if yit else (y, )
    #store verts length as a range to juggle j
    r = range(len(verts))
    #final results if x or y is iterable
    results = []
    #traverse x and y coordinates
    for xp in X:
        for yp in Y:
            c = 0 #reset c at every new position
            for i in r:
                j = r[i-1] #set j to position before i
                #store a few arguments to shorten the if statement
                yneq       = (verts[i][1] > yp) != (verts[j][1] > yp)
                xofs, yofs = (verts[j][0] - verts[i][0]), (verts[j][1] - verts[i][1])
                #if we have crossed a line, increment c
                if (yneq and (xp < xofs * (yp - verts[i][1]) / yofs + verts[i][0])):
                    c += 1
            #if c is odd store the coordinates        
            if c%2:
                results.append((xp, yp))
    #return either coordinates or a bool, depending if x or y was an iterable
    return results if (xit or yit) else bool(c%2)

这个python版本是通用的。您可以为True/False结果输入单个x和单个y值,也可以使用x和y的范围来遍历整个点网格。如果使用范围,则返回所有True点的x/y对列表。vertices参数需要一个由x/y对组成的二维Iterable,例如:[(x1,y1), (x2,y2),…]

使用示例:

vertices = [(25,25), (75,25), (75,75), (25,75)]
pnpoly(vertices, 50, 50) #True
pnpoly(vertices, range(100), range(100)) #[(25,25), (25,26), (25,27), ...]

实际上,这些都可以。

pnpoly(vertices, 50, range(100)) #check 0 to 99 y at x of 50
pnpoly(vertices, range(100), 50) #check 0 to 99 x at y of 50

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

在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倍)。

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

如果点在多边形区域内。(可通过Ray-Casting算法实现) 如果点在多边形边界上(可以用与在折线(线)上检测点相同的算法来完成)。