我正在尝试创建一个快速的2D点在多边形内部算法,用于命中测试(例如Polygon.contains(p:Point)
)。欢迎提供有效技巧建议。
我正在尝试创建一个快速的2D点在多边形内部算法,用于命中测试(例如Polygon.contains(p:Point)
)。欢迎提供有效技巧建议。
(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!
}
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
}
#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;
}
我认为下面这段代码是最好的解决方案(来自这里):
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
两侧,因此它们永远不会相等。 - Peter Wood这是一个C#版本的
在顶部添加了包围框检查。然而,正如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;
}
inside = !inside;
之后断开?不行。当它遇到各种多边形段时,算法会多次更改。 - M Katz以下是基于Nirg方法的M. Katz答案的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;
}
计算点p与多边形各个角之间的有向角之和。如果总有向角为360度,则点在内部。如果总角度为0,则该点在外面。
我更喜欢这种方法,因为它更加健壮,不太依赖数值精度。
计算交点数量的方法受限于您在计算交点数量时可能会“命中”一个角。
编辑:顺便说一下,这种方法适用于凸多边形和凹多边形。
编辑:我最近发现了整个关于此主题的维基百科文章。
这个问题非常有趣。我有一个与其他回答不同的可行想法,可以使用角度和来决定目标点是在内部还是外部。更为人熟知的是绕数算法。
设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
bobobobo提到的Eric Haines文章非常好。特别有趣的是比较算法性能的表格;与其他方法相比,角度汇总法真的很差。另一个有趣的事情是,优化措施,如使用查找网格将多边形进一步细分为“内部”和“外部”部分,即使在具有> 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
这个问题中的大多数答案都没有很好地处理所有的边角情况。以下是一些微妙的边角情况:
这是一个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
}
我真的很喜欢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;
}
Swift版本的nirg的答案:
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
}
}
point.within(polygon)
,它会返回一个布尔值,表示点是否在多边形内部。 - J Agustin Barrachina