我正在尝试在多边形算法中创建一个快速 2D 点，用于命中测试（例如Polygon.contains(p:Point)）。对有效技术的建议将不胜感激。
【问题讨论】：
标签： performance graphics collision-detection polygon point-in-polygon
对于图形，我宁愿不喜欢整数。许多系统使用整数进行 UI 绘制（像素毕竟是整数），但例如 macOS，对所有内容都使用浮点数。 macOS 只知道点，一个点可以转换为一个像素，但根据显示器分辨率，它可能会转换为其他东西。在视网膜屏幕上，半个点 (0.5/0.5) 是像素。尽管如此，我从未注意到 macOS UI 比其他 UI 慢得多。毕竟，3D API（OpenGL 或 Direct3D）也适用于浮点数，现代图形库经常利用 GPU 加速。
现在你说速度是你最关心的问题，好吧，让我们追求速度。在你运行任何复杂的算法之前，首先做一个简单的测试。在多边形周围创建一个轴对齐边界框。这非常简单、快速，并且已经可以为您节省大量计算。这是如何运作的？遍历多边形的所有点并找到 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) 点肯定在多边形之外。
但是e 是什么？好吧，e（实际上是 epsilon）给边界框一些 padding。正如我所说，如果我们开始离多边形线太近，光线追踪就会失败。由于边界框可能等于多边形（如果多边形是轴对齐的矩形，则边界框等于多边形本身！），我们需要一些填充来确保安全，仅此而已。你应该选择多大的e？不会太大。这取决于您用于绘图的坐标系比例。如果您的像素步长为 1.0，则只需选择 1.0（但 0.1 也可以）
现在我们有了射线及其起点和终点坐标，问题从“是多边形内的点”转移到“射线与多边形边相交的频率”。因此，我们不能像以前那样只使用多边形点，现在我们需要实际的边。边总是由两个点定义的。
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 很清楚。 YES 增加交叉点，NO 什么都不做。
COLLINEAR 呢？这意味着两个向量都位于同一条无限线上，取决于位置和长度，它们根本不相交或相交于无数个点。我不确定如何处理这种情况，无论哪种方式，我都不会将其视为交集。好吧，由于浮点舍入错误，这种情况在实践中相当罕见；更好的代码可能不会测试== 0.0f，而是测试< epsilon 之类的东西，其中epsilon 是一个相当小的数字。
如果你需要测试更多的点，你当然可以通过将多边形边的线性方程标准形式保存在内存中来加快整个过程，这样你就不必每次都重新计算这些。这将在每次测试中为您节省两个浮点乘法和三个浮点减法，以换取在内存中存储每个多边形边的三个浮点值。这是典型的内存与计算时间的权衡。
最后但同样重要的是：如果您可以使用 3D 硬件来解决问题，还有一个有趣的替代方案。只需让 GPU 为您完成所有工作。创建一个屏幕外的绘画表面。用黑色完全填充它。现在让 OpenGL 或 Direct3D 绘制您的多边形（如果您只想测试该点是否在其中任何一个内，甚至可以绘制所有多边形，但您不关心哪个）并用不同的填充多边形颜色，例如白色的。要检查一个点是否在多边形内，请从绘图表面获取该点的颜色。这只是一个 O(1) 内存提取。
当然，这种方法仅在您的绘图表面不必很大时才可用。如果它无法适应 GPU 内存，则此方法比在 CPU 上执行要慢。如果它必须很大并且您的 GPU 支持现代着色器，您仍然可以通过将上面显示的光线投射实现为 GPU 着色器来使用 GPU，这绝对是可能的。对于要测试的大量多边形或大量点，这将获得回报，考虑一些 GPU 将能够并行测试 64 到 256 个点。但是请注意，将数据从 CPU 传输到 GPU 并返回总是很昂贵，因此仅针对几个简单的多边形测试几个点，其中点或多边形是动态的并且会经常变化，GPU 方法很少会支付关闭。
【讨论】：
我认为这是迄今为止所有答案中最简洁的另一个 numpyic 实现。
大多边形顶点的二维坐标是
[[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] 在。
我的解决方案基于winding number algorithm。以下是我仅使用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]
【讨论】：
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 版本是通用的。您可以为真/假结果输入单个 x 和单个 y 值，也可以使用 range 为 x 和 y 遍历整个点网格。如果使用范围，则返回所有 True 点的 x/y 对的 list。 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
【讨论】：
这个问题中的大多数答案都不能很好地处理所有极端情况。一些微妙的极端情况，如下所示： 这是一个 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
}
【讨论】：
这可能是来自here 的 C 代码的优化程度稍差的版本，来自from this page。
我的 C++ 版本使用 std::vector<std::pair<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;
}
【讨论】：
这是answer given by nirg 的C# 版本，它来自this RPI professor。请注意，使用来自该 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;
}
【讨论】：
for 和if 有什么意义？最后一个 for 适用于所有情况。
为了完整起见，这里是nirg提供并由Mecki讨论的算法的lua实现：
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
变量area 是一个点表，这些点又存储为二维表。示例：
> A = {{2, 1}, {1, 2}, {15, 3}, {3, 4}, {5, 3}, {4, 1.5}}
> T = {2, 1.1}
> pnpoly(A, T)
true
link 到 GitHub Gist。
【讨论】：
如果您使用 Google Map SDK 并想检查一个点是否在多边形内，您可以尝试使用GMSGeometryContainsLocation。效果很好！！这是它的工作原理，
if GMSGeometryContainsLocation(point, polygon, true) {
print("Inside this polygon.")
} else {
print("outside this polygon")
}
【讨论】：
这似乎在 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)}
【讨论】：
这个问题很有趣。我有另一个与这篇文章的其他答案不同的可行想法。这个想法是使用角度的总和来确定目标是在内部还是外部。更为人所知的是winding number。
设 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
【讨论】：
真的很喜欢 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;
}
【讨论】：
我认为下面这段代码是最好的解决方案（取自here）：
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;
}
它既短又高效，适用于凸多边形和凹多边形。如前所述，您应该首先检查边界矩形并分别处理多边形孔。
这背后的想法很简单。作者是这样描述的：
我从测试点水平射出一条半无限射线（增加 x，固定 y），并计算它穿过了多少条边。在每个交叉点，光线在内部和外部之间切换。这就是乔丹曲线定理。
每次水平光线穿过任何边缘时，变量 c 都会从 0 切换到 1 和 1 到 0。所以基本上它会跟踪交叉的边数是偶数还是奇数。 0 表示偶数，1 表示奇数。
【讨论】：
verty[i] 和verty[j] 是testy 的任一侧，所以它们永远不相等。
很惊讶之前没有人提出这个问题，但是对于需要数据库的实用主义者来说：MongoDB 对包括这个在内的地理查询有很好的支持。
你要找的是：
db.neighborhoods.findOne({ 几何: { $geoIntersects: { $geometry: { 类型：“点”，坐标：[“经度”，“纬度”] } } } })
Neighborhoods 是以标准 GeoJson 格式存储一个或多个多边形的集合。如果查询返回 null 则不相交，否则相交。
这里有很好的记录： https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
在 330 个不规则多边形网格中分类的 6,000 多个点的性能不到一分钟，完全没有优化，包括用各自的多边形更新文档的时间。
【讨论】：
这是@nirg 答案的golang 版本（灵感来自@@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
}
【讨论】：
我已经做了一个nirg's c++code的Python实现：
输入
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
同样，这个想法来自here
【讨论】：
nirg 解决方案的 Scala 版本（假设边界矩形预检查是单独完成的）：
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)
}
【讨论】：
答案取决于您是简单多边形还是复杂多边形。简单多边形不得有任何线段交点。所以它们可以有孔，但线不能相互交叉。复杂区域可以有线交点 - 因此它们可以有重叠区域，或者仅通过一个点相互接触的区域。
对于简单的多边形，最好的算法是光线投射（交叉数）算法。对于复杂的多边形，此算法不会检测重叠区域内的点。所以对于复杂的多边形，你必须使用绕组数算法。
这是一篇出色的文章，其中包含两种算法的 C 实现。我试过了，效果很好。
【讨论】：
answer by 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
}
}
【讨论】：
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
【讨论】：
当使用qt（Qt 4.3+）时，可以使用QPolygon的函数containsPoint
【讨论】：
如果您正在寻找一个 java-script 库，那么 Polygon 类有一个 javascript google maps v3 扩展来检测一个点是否驻留在其中。
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
【讨论】：
您可以通过检查将所需点连接到多边形顶点所形成的区域是否与多边形本身的区域相匹配来做到这一点。
或者您可以检查从您的点到每对两个连续多边形顶点到您的检查点的内角总和是否为 360，但我觉得第一个选项更快，因为它不涉及三角函数倒数的除法和计算。
我不知道如果你的多边形里面有一个洞会发生什么，但在我看来，主要思想可以适应这种情况
您也可以在数学社区中发布问题。我敢打赌他们有一百万种方法来做到这一点
【讨论】：
处理Ray casting algorithm中的以下特殊情况：
检查Determining Whether A Point Is Inside A Complex Polygon。本文提供了一种简单的解决方法，因此上述情况无需特殊处理。
【讨论】：
为了检测多边形上的命中，我们需要测试两件事：
【讨论】：
bobobobo 引用的Eric Haines article 真的很棒。特别有趣的是比较算法性能的表格；与其他方法相比，角度求和方法确实很糟糕。同样有趣的是，使用查找网格将多边形进一步细分为“内”和“外”扇区等优化可以使测试速度非常快，即使在边数超过 1000 的多边形上也是如此。
无论如何，现在还为时过早，但我的投票是“交叉”方法，我认为这与 Mecki 所描述的差不多。但是我发现它最简洁described and codified by 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
【讨论】：
这仅适用于凸形，但 Minkowski Portal Refinement 和 GJK 也是测试点是否在多边形中的绝佳选择。您使用 minkowski 减法从多边形中减去点，然后运行这些算法来查看多边形是否包含原点。
另外，有趣的是，您可以使用支持函数更隐含地描述您的形状，该支持函数将方向向量作为输入并沿该向量吐出最远的点。这使您可以描述任何凸面形状.. 弯曲的、由多边形制成的或混合的。您还可以进行操作，将简单支持函数的结果组合成更复杂的形状。
更多信息： http://xenocollide.snethen.com/mpr2d.html
此外，游戏编程 gems 7 讨论了如何在 3d 中执行此操作（：
【讨论】：
计算点 p 和每个多边形顶点之间的方向角和。如果总定向角度为 360 度，则该点在内部。如果总和为 0，则该点在外面。
我更喜欢这种方法，因为它更健壮且对数值精度的依赖性更小。
计算交叉点数量均匀度的方法受到限制，因为您可以在计算交叉点数量期间“命中”顶点。
编辑：顺便说一下，这种方法适用于凹凸多边形。
编辑：我最近在该主题上找到了整个Wikipedia article。
【讨论】：
这是 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 倍)。
【讨论】：
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;
}
}
【讨论】：