我试图创建一个快速的2D点内多边形算法,用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
当使用qt (qt 4.3+)时,可以使用QPolygon的函数containsPoint
其他回答
这似乎在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)}
计算点p与每个多边形顶点之间的有向角和。如果总倾斜角是360度,那么这个点在里面。如果总数为0,则点在外面。
我更喜欢这种方法,因为它更健壮,对数值精度的依赖更小。
计算交集数量的均匀性的方法是有限的,因为你可以在计算交集数量的过程中“击中”一个顶点。
编辑:顺便说一下,这种方法适用于凹凸多边形。
编辑:我最近在维基百科上找到了一篇关于这个话题的完整文章。
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是一个很大的帮助。
以下是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;
}
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