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

我曾在Michael Stonebraker名下担任研究员时，对此进行了一些工作 - 你知道的，他是Ingres、PostgreSQL等项目的创始人。

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

如果你想要一个很棒的算法，可以查看开源项目PostgreSQL的源代码来处理地理信息...

我想指出的是，我们从未得到过关于左右手问题（也可表达为“内部”与“外部”问题）的任何见解...

---

更新

BKB的链接提供了许多合理的算法。我正在处理地球科学问题，因此需要一种在纬度/经度下工作的解决方案，并且它有一个特殊的问题 - 小区域内部还是大区域内部？答案是顶点的“方向”很重要 - 它可以是左手或右手，以此方式您可以将任何给定多边形的任一区域指示为“内部”。因此，我的工作使用了该页面上枚举的第三种解决方案。

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

...由于有人问：当顶点数超过某个数字时，我们发现边界框测试效果最好 - 如果需要，在进行更长时间的测试之前进行非常快速的测试...边界框通过简单地取最大x，最小x，最大y和最小y并将它们放在一起来创建一个框的四个点...

对于那些跟随的人的另一个提示：我们所有更复杂和“光暗”计算都在平面上的正点格网空间中完成，然后重新投影回“真实”的经度/纬度，从而避免了可能出现的错误，如越过经线180时的包裹错误以及处理极地区域时的错误。效果非常好！

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

}

一个简单的解决方案是将多边形分割成三角形，并像这里所解释的那样对其进行命中测试。

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

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

像许多优化一样，这些都是基于特定而不是通用情况的，并且根据摊销时间而不是单个使用效益。

在这个领域工作时，我发现Joseph O'Rourke的《C语言计算几何》（ISBN 0-521-44034-3）对我很有帮助。

我知道这已经是老话题了，但是这里提供了一个在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;
}

没有比归纳定义问题更美的了。为了完整起见，这里有一个用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]).

给出两点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答案的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");
    }
}

我已经用Python实现了nirg的C++代码（链接在此：nirg's code）：
输入：
- 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

再次强调，这个想法来自于这里。

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

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