下面是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
}

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

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

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

这似乎在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)}

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

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

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