我试图创建一个快速的2D点内多边形算法,用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
如果你正在使用谷歌Map SDK,想要检查一个点是否在一个多边形内,你可以尝试使用GMSGeometryContainsLocation。效果很好!!它是这样运作的,
if GMSGeometryContainsLocation(point, polygon, true) {
print("Inside this polygon.")
} else {
print("outside this polygon")
}
这里是参考资料:https://developers.google.com/maps/documentation/ios-sdk/reference/group___geometry_utils#gaba958d3776d49213404af249419d0ffd
其他回答
没有什么比归纳定义问题更美好的了。为了完整起见,你在序言中有一个版本,它可能也澄清了光线投射背后的思想:
基于仿真的简化算法在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
nirg的c#版本的答案在这里:我只分享代码。这可能会节省一些时间。
public static bool IsPointInPolygon(IList<Point> polygon, Point testPoint) {
bool result = false;
int j = polygon.Count() - 1;
for (int i = 0; i < polygon.Count(); i++) {
if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) {
if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) {
result = !result;
}
}
j = i;
}
return result;
}
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
下面是Rust版本的@nirg答案(Philipp Lenssen javascript版本) 我给出这个答案是因为我从这个网站得到了很多帮助,我翻译javascript版本rust作为一个练习,希望可以帮助一些人,最后一个原因是,在我的工作中,我会把这段代码翻译成一个wasm,以提高我的画布的性能,这是一个开始。我的英语很差……,请原谅我 `
pub struct Point {
x: f32,
y: f32,
}
pub fn point_is_in_poly(pt: Point, polygon: &Vec<Point>) -> bool {
let mut is_inside = false;
let max_x = polygon.iter().map(|pt| pt.x).reduce(f32::max).unwrap();
let min_x = polygon.iter().map(|pt| pt.x).reduce(f32::min).unwrap();
let max_y = polygon.iter().map(|pt| pt.y).reduce(f32::max).unwrap();
let min_y = polygon.iter().map(|pt| pt.y).reduce(f32::min).unwrap();
if pt.x < min_x || pt.x > max_x || pt.y < min_y || pt.y > max_y {
return is_inside;
}
let len = polygon.len();
let mut j = len - 1;
for i in 0..len {
let y_i_value = polygon[i].y > pt.y;
let y_j_value = polygon[j].y > pt.y;
let last_check = (polygon[j].x - polygon[i].x) * (pt.y - polygon[i].y)
/ (polygon[j].y - polygon[i].y)
+ polygon[i].x;
if y_i_value != y_j_value && pt.x < last_check {
is_inside = !is_inside;
}
j = i;
}
is_inside
}
let pt = Point {
x: 1266.753,
y: 97.655,
};
let polygon = vec![
Point {
x: 725.278,
y: 203.586,
},
Point {
x: 486.831,
y: 441.931,
},
Point {
x: 905.77,
y: 445.241,
},
Point {
x: 1026.649,
y: 201.931,
},
];
let pt1 = Point {
x: 725.278,
y: 203.586,
};
let pt2 = Point {
x: 872.652,
y: 321.103,
};
println!("{}", point_is_in_poly(pt, &polygon));// false
println!("{}", point_is_in_poly(pt1, &polygon)); // true
println!("{}", point_is_in_poly(pt2, &polygon));// 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的性能很差,或者我只是在处理一个糟糕的实现?
- 大概的成本访问各种缓存和主存储器?