我试图创建一个快速的2D点内多边形算法,用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
计算点p与每个多边形顶点之间的有向角和。如果总倾斜角是360度,那么这个点在里面。如果总数为0,则点在外面。
我更喜欢这种方法,因为它更健壮,对数值精度的依赖更小。
计算交集数量的均匀性的方法是有限的,因为你可以在计算交集数量的过程中“击中”一个顶点。
编辑:顺便说一下,这种方法适用于凹凸多边形。
编辑:我最近在维基百科上找到了一篇关于这个话题的完整文章。
其他回答
当我还是Michael Stonebraker手下的一名研究员时,我做了一些关于这方面的工作——你知道,就是那位提出了Ingres、PostgreSQL等的教授。
我们意识到最快的方法是首先做一个边界框,因为它非常快。如果它在边界框之外,它就在外面。否则,你就得做更辛苦的工作……
如果你想要一个伟大的算法,看看开源项目PostgreSQL的源代码的地理工作…
我想指出的是,我们从来没有深入了解过左撇子和右撇子(也可以表达为“内”和“外”的问题……
更新
BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.
此外,我的工作使用单独的函数进行“在线”测试。
...因为有人问:我们发现当垂直的数量超过某个数字时,边界盒测试是最好的——如果有必要,在做更长的测试之前做一个非常快速的测试……边界框是通过简单地将最大的x,最小的x,最大的y和最小的y放在一起,组成一个框的四个点来创建的……
另一个提示是:我们在网格空间中进行了所有更复杂的“调光”计算,都是在平面上的正点上进行的,然后重新投影到“真实”的经度/纬度上,从而避免了在经度180线交叉时和处理极地时可能出现的环绕错误。工作好了!
为了完整性,这里是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的链接。
计算点p与每个多边形顶点之间的有向角和。如果总倾斜角是360度,那么这个点在里面。如果总数为0,则点在外面。
我更喜欢这种方法,因为它更健壮,对数值精度的依赖更小。
计算交集数量的均匀性的方法是有限的,因为你可以在计算交集数量的过程中“击中”一个顶点。
编辑:顺便说一下,这种方法适用于凹凸多边形。
编辑:我最近在维基百科上找到了一篇关于这个话题的完整文章。
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
对于检测多边形上的命中,我们需要测试两件事:
如果点在多边形区域内。(可通过Ray-Casting算法实现) 如果点在多边形边界上(可以用与在折线(线)上检测点相同的算法来完成)。
推荐文章
- 如何加速gwt编译器?
- MySQL OR与IN性能
- 应该……接住环内还是环外?
- 哪个更快/最好?SELECT *或SELECT columnn1, colum2, column3等
- 加快R中的循环操作
- INT和VARCHAR主键之间有真正的性能差异吗?
- c++标准是否要求iostreams的性能很差,或者我只是在处理一个糟糕的实现?
- 大概的成本访问各种缓存和主存储器?
- 模拟慢速互联网连接
- 如何检查表上持有哪些锁
- 检查字符串是否包含字符串列表中的元素
- Scala vs Python的Spark性能
- 现代c++能让你免费获得性能吗?
- 对于PostgreSQL表来说,多大才算太大?
- 即使从未抛出异常,使用try-catch块的代价是否昂贵?