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

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

当我还是Michael Stonebraker手下的一名研究员时，我做了一些关于这方面的工作——你知道，就是那位提出了Ingres、PostgreSQL等的教授。

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

如果你想要一个伟大的算法，看看开源项目PostgreSQL的源代码的地理工作…

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

更新

BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.

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

.．.因为有人问:我们发现当垂直的数量超过某个数字时，边界盒测试是最好的——如果有必要，在做更长的测试之前做一个非常快速的测试……边界框是通过简单地将最大的x，最小的x，最大的y和最小的y放在一起，组成一个框的四个点来创建的……

另一个提示是:我们在网格空间中进行了所有更复杂的“调光”计算，都是在平面上的正点上进行的，然后重新投影到“真实”的经度/纬度上，从而避免了在经度180线交叉时和处理极地时可能出现的环绕错误。工作好了!

计算点p与每个多边形顶点之间的有向角和。如果总倾斜角是360度，那么这个点在里面。如果总数为0，则点在外面。

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

计算交集数量的均匀性的方法是有限的，因为你可以在计算交集数量的过程中“击中”一个顶点。

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

编辑:我最近在维基百科上找到了一篇关于这个话题的完整文章。

David Segond's answer is pretty much the standard general answer, and Richard T's is the most common optimization, though therre are some others. Other strong optimizations are based on less general solutions. For example if you are going to check the same polygon with lots of points, triangulating the polygon can speed things up hugely as there are a number of very fast TIN searching algorithms. Another is if the polygon and points are on a limited plane at low resolution, say a screen display, you can paint the polygon onto a memory mapped display buffer in a given colour, and check the color of a given pixel to see if it lies in the polygons.

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

在这个领域工作，我发现约瑟夫·奥鲁克斯的《计算几何》在C' ISBN 0-521-44034-3是一个很大的帮助。

For graphics, I'd rather not prefer integers. Many systems use integers for UI painting (pixels are ints after all), but macOS, for example, uses float for everything. macOS only knows points and a point can translate to one pixel, but depending on monitor resolution, it might translate to something else. On retina screens half a point (0.5/0.5) is pixel. Still, I never noticed that macOS UIs are significantly slower than other UIs. After all, 3D APIs (OpenGL or Direct3D) also work with floats and modern graphics libraries very often take advantage of GPU acceleration.

现在你说速度是你最关心的，好吧，让我们追求速度。在运行任何复杂的算法之前，先做一个简单的测试。在多边形周围创建一个轴对齐的包围框。这是非常简单，快速的，已经可以节省你很多计算。这是怎么做到的呢?遍历多边形的所有点，找到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!
}

这是对任意点运行的第一个测试。正如你所看到的，这个测试非常快，但也非常粗糙。要处理边界矩形内的点，我们需要更复杂的算法。有几种计算方法。哪种方法有效还取决于多边形是否有孔或始终是固体。以下是实体的例子(一个凸面，一个凹面):

这里有一个洞:

绿色的中间有个洞!

最简单的算法，可以处理上述三种情况，并且仍然非常快，叫做射线投射。该算法的思想非常简单:从多边形外的任何地方绘制一条虚拟光线到你的点，并计算它击中多边形一侧的频率。如果命中次数是偶数，则在多边形外，如果是奇数，则在多边形内。

圈数算法是另一种选择，它对非常接近多边形线的点更准确，但也慢得多。由于有限的浮点精度和舍入问题，光线投射可能会因为太靠近多边形一侧的点而失败，但在现实中这几乎不是问题，因为如果一个点靠近一侧，在视觉上甚至不可能让观看者识别它是否已经在内部或仍然在外部。

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

But what is e? Well, e (actually epsilon) gives the bounding box some padding. As I said, ray tracing fails if we start too close to a polygon line. Since the bounding box might equal the polygon (if the polygon is an axis aligned rectangle, the bounding box is equal to the polygon itself!), we need some padding to make this safe, that's all. How big should you choose e? Not too big. It depends on the coordinate system scale you use for drawing. If your pixel step width is 1.0, then just choose 1.0 (yet 0.1 would have worked as well)

现在我们有了光线的起始坐标和结束坐标，问题从“是多边形内的点”转移到“光线与多边形边相交的频率”。因此，我们不能像以前那样只处理多边形点，现在我们需要实际的边。一条边总是由两点来定义。

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很清楚。“是”增加了交叉路口，“否”什么都不做。

What about COLLINEAR? It means both vectors lie on the same infinite line, depending on position and length, they don't intersect at all or they intersect in an endless number of points. I'm not absolutely sure how to handle this case, I would not count it as intersection either way. Well, this case is rather rare in practice anyway because of floating point rounding errors; better code would probably not test for == 0.0f but instead for something like < epsilon, where epsilon is a rather small number.

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

Last but not least: If you may use 3D hardware to solve the problem, there is an interesting alternative. Just let the GPU do all the work for you. Create a painting surface that is off screen. Fill it completely with the color black. Now let OpenGL or Direct3D paint your polygon (or even all of your polygons if you just want to test if the point is within any of them, but you don't care for which one) and fill the polygon(s) with a different color, e.g. white. To check if a point is within the polygon, get the color of this point from the drawing surface. This is just a O(1) memory fetch.

Of course this method is only usable if your drawing surface doesn't have to be huge. If it cannot fit into the GPU memory, this method is slower than doing it on the CPU. If it would have to be huge and your GPU supports modern shaders, you can still use the GPU by implementing the ray casting shown above as a GPU shader, which absolutely is possible. For a larger number of polygons or a large number of points to test, this will pay off, consider some GPUs will be able to test 64 to 256 points in parallel. Note however that transferring data from CPU to GPU and back is always expensive, so for just testing a couple of points against a couple of simple polygons, where either the points or the polygons are dynamic and will change frequently, a GPU approach will rarely pay off.

bobobobo引用的Eric Haines的文章真的很棒。特别有趣的是比较算法性能的表格;角度求和法和其他方法比起来真的很差。同样有趣的是，使用查找网格将多边形进一步细分为“in”和“out”扇区的优化可以使测试非常快，即使是在> 1000条边的多边形上。

不管怎样，现在还为时过早，但我的投票倾向于“交叉”方法，我认为这几乎就是Mecki所描述的。然而，我发现大卫·伯克(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

我认为下面这段代码是最好的解决方案(从这里开始):

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表示奇数。

我知道这是旧的，但这里是一个在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;
}

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

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

Obj-C版本nirg的答案与样本方法测试点。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");
    }
}

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

基于仿真的简化算法在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).

下面是nirg给出的答案的c#版本，它来自RPI教授。请注意，使用来自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;
}

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

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

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

}

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

这只适用于凸形状，但是Minkowski Portal Refinement和GJK也是测试一个点是否在多边形中的很好的选择。您使用闵可夫斯基减法从多边形中减去点，然后运行这些算法来查看多边形是否包含原点。

另外，有趣的是，你可以用支持函数更隐式地描述你的形状，它以一个方向向量作为输入，并输出沿该向量的最远点。这可以让你描述任何凸形状..弯曲的，由多边形制成的，或混合的您还可以执行一些操作，将简单支持函数的结果组合起来，以生成更复杂的形状。

更多信息: http://xenocollide.snethen.com/mpr2d.html

此外，game programming gems 7讨论了如何在3d中做到这一点(:

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

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

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

在Ray casting算法中处理以下特殊情况:

射线与多边形的一条边重叠。 点在多边形的内部，光线穿过多边形的顶点。 该点在多边形的外部，光线只接触到多边形的一个角。

检查确定一个点是否在一个复杂多边形内。本文提供了一种简单的解决方法，因此对于上述情况不需要特殊处理。

您可以通过检查将所需点连接到多边形顶点所形成的面积是否与多边形本身的面积相匹配来实现这一点。

或者你可以检查从你的点到每一对连续的多边形顶点到你的检查点的内角之和是否为360，但我有一种感觉，第一种选择更快，因为它不涉及除法，也不计算三角函数的反函数。

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

你也可以把问题贴在数学社区里。我打赌他们有一百万种方法

如果你正在寻找一个java脚本库，有一个javascript谷歌maps v3扩展的Polygon类，以检测是否有一个点驻留在它里面。

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

谷歌扩展Github

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

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实现这两种算法的优秀文章。我试过了，效果不错。

http://geomalgorithms.com/a03-_inclusion.html

这个问题很有趣。我有另一个可行的想法，不同于这篇文章的其他答案。其原理是利用角度之和来判断目标是在内部还是外部。也就是圈数。

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

参考这张图来对这个概念有一个基本的了解:

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

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

Scala版本的解决方案由nirg(假设边界矩形预检查是单独完成的):

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

我已经做了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

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

下面是golang版本的@nirg答案(灵感来自于@@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
}

令人惊讶的是之前没有人提出这个问题，但是对于需要数据库的实用主义者来说:MongoDB对Geo查询提供了出色的支持，包括这个查询。

你需要的是:

db.neighborhoods。findOne({geometry: {$geoIntersects: {$geometry: { type: "Point"，坐标:["经度"，"纬度"]}}} })

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

这里有详细的记录: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

在330个不规则多边形网格中，超过6000个点分类的性能不到一分钟，没有任何优化，包括用各自的多边形更新文档的时间。

这似乎在R中工作(为丑陋道歉，希望看到更好的版本!)。

pnpoly <- function(nvert,vertx,verty,testx,testy){
          c <- FALSE
          j <- nvert 
          for (i in 1:nvert){
              if( ((verty[i]>testy) != (verty[j]>testy)) && 
   (testx < (vertx[j]-vertx[i])*(testy-verty[i])/(verty[j]-verty[i])+vertx[i]))
            {c <- !c}
             j <- i}
   return(c)}

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

为了完整性，这里是nirg提供的算法的lua实现，由Mecki讨论:

function pnpoly(area, test)
    local inside = false
    local tx, ty = table.unpack(test)
    local j = #area
    for i=1, #area do
        local vxi, vyi = table.unpack(area[i])
        local vxj, vyj = table.unpack(area[j])
        if (vyi > ty) ~= (vyj > ty)
        and tx < (vxj - vxi)*(ty - vyi)/(vyj - vyi) + vxi
        then
            inside = not inside
        end
        j = i
    end
    return inside
end

变量区域是一个点的表，这些点依次存储为2D表。例子:

> A = {{2, 1}, {1, 2}, {15, 3}, {3, 4}, {5, 3}, {4, 1.5}}
> T = {2, 1.1}
> pnpoly(A, T)
true

GitHub Gist的链接。

这大概是一个稍微不那么优化的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;
}

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

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

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

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

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

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

Like David Segonds' answer suggests I use an approach of angle summation derived from my concave polygon drawing algorithm. It relies of adding up the approximate angles of subtriangles around the point to obtain a weight. A weight around 1.0 means the point is inside the triangle, a weight around 0.0 means outside, a weight around -1.0 is what happens when inside the polygon but in reverse order (like with one of the halves of a bowtie-shaped tetragon) and a weight of NAN if exactly on an edge. The reason it's not slow is that angles don't need to be estimated accurately at all. Holes can be handled by treating them as separate polygons and subtracting the weights.

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

因此，对于多边形的每一段，都形成一个子三角形，并计算点，然后旋转每个子三角形以计算其近似角度并添加到权重。

调用subtriangle_angle_approx(y, x)可以替换为atan2(y, x) / (2.*pi)，但是一个非常粗略的近似值就足够精确了:

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

下面是Rust版本的@nirg答案(Philipp Lenssen javascript版本) 我给出这个答案是因为我从这个网站得到了很多帮助，我翻译javascript版本rust作为一个练习，希望可以帮助一些人，最后一个原因是，在我的工作中，我会把这段代码翻译成一个wasm，以提高我的画布的性能，这是一个开始。我的英语很差……,请原谅我 ｀

pub struct Point {
    x: f32,
    y: f32,
}
pub fn point_is_in_poly(pt: Point, polygon: &Vec<Point>) -> bool {
    let mut is_inside = false;

    let max_x = polygon.iter().map(|pt| pt.x).reduce(f32::max).unwrap();
    let min_x = polygon.iter().map(|pt| pt.x).reduce(f32::min).unwrap();
    let max_y = polygon.iter().map(|pt| pt.y).reduce(f32::max).unwrap();
    let min_y = polygon.iter().map(|pt| pt.y).reduce(f32::min).unwrap();

    if pt.x < min_x || pt.x > max_x || pt.y < min_y || pt.y > max_y {
        return is_inside;
    }

    let len = polygon.len();
    let mut j = len - 1;

    for i in 0..len {
        let y_i_value = polygon[i].y > pt.y;
        let y_j_value = polygon[j].y > pt.y;
        let last_check = (polygon[j].x - polygon[i].x) * (pt.y - polygon[i].y)
            / (polygon[j].y - polygon[i].y)
            + polygon[i].x;
        if y_i_value != y_j_value && pt.x < last_check {
            is_inside = !is_inside;
        }
        j = i;
    }
    is_inside
}


let pt = Point {
    x: 1266.753,
    y: 97.655,
};
let polygon = vec![
    Point {
        x: 725.278,
        y: 203.586,
    },
    Point {
        x: 486.831,
        y: 441.931,
    },
    Point {
        x: 905.77,
        y: 445.241,
    },
    Point {
        x: 1026.649,
        y: 201.931,
    },
];
let pt1 = Point {
    x: 725.278,
    y: 203.586,
};
let pt2 = Point {
    x: 872.652,
    y: 321.103,
};
println!("{}", point_is_in_poly(pt, &polygon));// false
println!("{}", point_is_in_poly(pt1, &polygon)); // true
println!("{}", point_is_in_poly(pt2, &polygon));// true

`

在大多数情况下，这是一个比其他算法都快的算法。

它又新又雅致。我们花费O(n * log(n))时间构建一个表，允许我们在O(log(n) + k)时间内测试多边形中的点。

与光线跟踪或角度不同，使用扫描光束表可以更快地对同一多边形进行多次检查。我们必须预先构建一个扫描束活动边表，这是大多数代码正在做的事情。

We calculate the scanbeam and the active edges for that position in the y-direction. We make a list of points sorted by their y-component and we iterate through this list, for two events. Start-Y and End-Y, we track the active edges as we process the list. We process the events in order and for each scanbeam we record the y-value of the event and the active edges at each event (events being start-y and end-y) but we only record these when our event-y is different than last time (so everything at the event point is processed before we mark it in our table).

我们得到我们的表格:

[] p6p5、p6p7 p6p5, p6p7, p2p3, p2p1 p6p7, p5p4, p2p3, p3p1 p7p8, p5p4, p2p3, p2p1 p7p8, p5p4, p3p4, p2p1 p7p8 p2p1、 p7p8、p1p0 p8p0、p1p0 []

在构建该表之后，实际执行工作的代码只有几行。

注意:这里的代码使用复数值作为点。所以。real是。x。imag是。y。

def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

我们对事件进行二进制搜索，以找到特定值的actives_at_y。对于在y点的所有活动，我们计算在我们点的特定y点的x段值。每次x截距大于点的x分量时加1。然后对总数乘以2。(这是偶数-奇数填充规则，你可以很容易地适应任何其他填充规则)。

完整的代码:

from bisect import bisect

def build_edge_list(polygon):
    edge_list = []
    for i in range(1, len(polygon)):
        if (polygon[i].imag, polygon[i].real) < (polygon[i - 1].imag, polygon[i - 1].real):
            edge_list.append((polygon[i], i))
            edge_list.append((polygon[i - 1], ~i))
        else:
            edge_list.append((polygon[i], ~i))
            edge_list.append((polygon[i - 1], i))

    def sort_key(e):
        return e[0].imag, e[0].real, ~e[1]

    edge_list.sort(key=sort_key)
    return edge_list


def build_scanbeam(edge_list):
    actives_table = []
    events = []
    y = -float("inf")
    actives = []
    for pt, index in edge_list:
        if y != pt.imag:
            actives_table.append(list(actives))
            events.append(y)
        if index >= 0:
            actives.append(index)
        else:
            actives.remove(~index)
        y = pt.imag
    return actives_table, events

def point_in_polygon(polygon, point):
    def x_intercept(e, y):
        pt0 = polygon[e-1]
        pt1 = polygon[e]
        if pt1.real - pt0.real == 0:
            return pt0.real
        m = (pt1.imag - pt0.imag) / (pt1.real - pt0.real)
        b = pt0.imag - (m * pt0.real)
        return (y - b) / m

    edge_list = build_edge_list(polygon)
    actives_table, events = build_scanbeam(edge_list)
    try:
        if len(point):
            return [point_in_scantable(actives_table, events, x_intercept, p) for p in point]
    except TypeError:
        return point_in_scantable(actives_table, events, x_intercept, point)

def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

如果忽略，则扫描表的构建时间为O(n*log(n))。我们实际上是在O(log(n) + k)时间内查到的。其中n是多边形中段数的大小，k是该多边形中典型的活动边数。其他的光线追踪方法实际上需要O(n)时间。每次我们检查一个点，它迭代整个多边形。所以即使有这个明显的次优实现，它也轻而易举地打败了其他所有的。

There's a few performance tricks that could be done, for example, we can lower the time complexity to O(log(n) + log(k)) time. To do this we would implement Bentley-Ottmann into the sweep line, and rather than processing the intersections as different events, we split the lines at the intersections. We then also sort the active edges by their x-intercepts. We then know that no intersections occur during a scanbeam and since we sorted our segments (taking care to order them correctly within the scanbeam even if they start at the same initial point (you need to look at the slopes, or just compare midpoints of the segments). We then have a sorted intersection-less actives lists scanbeam table which means we can binary search into active edge list as well. Though that sounds like a lot of work for a value of k which is going to be typically 2 or maybe 4.

此外，由于这基本上变成了一个查找表和一些x截距的最小计算，它更能用GPU完成。你不再需要在多边形上循环了。所以你可以用numpy这样的东西来大量计算这些点，这样你就可以一次做所有的计算，从而提高性能。