我试图创建一个快速的2D点内多边形算法,用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
这个问题很有趣。我有另一个可行的想法,不同于这篇文章的其他答案。其原理是利用角度之和来判断目标是在内部还是外部。也就是圈数。
设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
其他回答
我知道这是旧的,但这里是一个在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;
}
简单的解决方案是将多边形划分为三角形,并按这里解释的那样对三角形进行测试
如果你的多边形是凸多边形,可能有更好的方法。把这个多边形看作是无限条线的集合。每一行将空间一分为二。对于每一个点,很容易判断它是在直线的一边还是另一边。如果一个点在所有直线的同一侧,那么它在多边形内。
在大多数情况下,这是一个比其他算法都快的算法。
它又新又雅致。我们花费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这样的东西来大量计算这些点,这样你就可以一次做所有的计算,从而提高性能。
以下是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;
}
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)
}
推荐文章
- 如何加速gwt编译器?
- MySQL OR与IN性能
- 应该……接住环内还是环外?
- 哪个更快/最好?SELECT *或SELECT columnn1, colum2, column3等
- 加快R中的循环操作
- INT和VARCHAR主键之间有真正的性能差异吗?
- c++标准是否要求iostreams的性能很差,或者我只是在处理一个糟糕的实现?
- 大概的成本访问各种缓存和主存储器?
- 模拟慢速互联网连接
- 如何检查表上持有哪些锁
- 检查字符串是否包含字符串列表中的元素
- Scala vs Python的Spark性能
- 现代c++能让你免费获得性能吗?
- 对于PostgreSQL表来说,多大才算太大?
- 即使从未抛出异常,使用try-catch块的代价是否昂贵?