我试图创建一个快速的2D点内多边形算法,用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
当使用qt (qt 4.3+)时,可以使用QPolygon的函数containsPoint
其他回答
没有什么比归纳定义问题更美好的了。为了完整起见,你在序言中有一个版本,它可能也澄清了光线投射背后的思想:
基于仿真的简化算法在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提供的算法的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的链接。
这似乎在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)}
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是一个很大的帮助。
我认为这是迄今为止所有答案中最简洁的一个。
例如,假设我们有一个多边形,它带有多边形凹,看起来像这样:
大多边形顶点的二维坐标为
[[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]
推荐文章
- 确定记录是否存在的最快方法
- 阅读GHC核心
- Python: List vs Dict用于查找表
- 为什么MATLAB的矩阵乘法运算这么快?
- for循环和for-each循环在性能上有区别吗?
- 就性能而言,使用std::memcpy()还是std::copy()更好?
- 什么时候我应该(不)想要在我的代码中使用熊猫apply() ?
- 如何加速gwt编译器?
- MySQL OR与IN性能
- 应该……接住环内还是环外?
- 哪个更快/最好?SELECT *或SELECT columnn1, colum2, column3等
- 加快R中的循环操作
- INT和VARCHAR主键之间有真正的性能差异吗?
- c++标准是否要求iostreams的性能很差,或者我只是在处理一个糟糕的实现?
- 大概的成本访问各种缓存和主存储器?