我试图创建一个快速的2D点内多边形算法，用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
在大多数情况下，这是一个比其他算法都快的算法。
它又新又雅致。我们花费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这样的东西来大量计算这些点，这样你就可以一次做所有的计算，从而提高性能。
其他回答
答案取决于你用的是简单多边形还是复杂多边形。简单多边形不能有任何线段交点。所以它们可以有洞，但线不能交叉。复杂区域可以有直线交点，所以它们可以有重叠的区域，或者只有一点相交的区域。
对于简单多边形，最好的算法是光线投射(交叉数)算法。对于复杂多边形，该算法不检测重叠区域内的点。所以对于复杂多边形你必须使用圈数算法。
下面是一篇用C实现这两种算法的优秀文章。我试过了，效果不错。
http://geomalgorithms.com/a03-_inclusion.html
下面是nirg给出的答案的c#版本，它来自RPI教授。请注意，使用来自RPI源代码的代码需要归属。
在顶部添加了一个边界框复选。然而，正如James Brown所指出的，主代码几乎和边界框检查本身一样快，所以边界框检查实际上会减慢整体操作，因为您正在检查的大多数点都在边界框内。所以你可以让边界框签出，或者另一种选择是预先计算多边形的边界框，如果它们不经常改变形状的话。
public bool IsPointInPolygon( Point p, Point[] polygon )
{
double minX = polygon[ 0 ].X;
double maxX = polygon[ 0 ].X;
double minY = polygon[ 0 ].Y;
double maxY = polygon[ 0 ].Y;
for ( int i = 1 ; i < polygon.Length ; i++ )
{
Point q = polygon[ i ];
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;
}
// https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
bool inside = false;
for ( int i = 0, j = polygon.Length - 1 ; 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 )
{
inside = !inside;
}
}
return inside;
}
在大多数情况下，这是一个比其他算法都快的算法。
它又新又雅致。我们花费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这样的东西来大量计算这些点，这样你就可以一次做所有的计算，从而提高性能。
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
Java版本:
public class Geocode {
private float latitude;
private float longitude;
public Geocode() {
}
public Geocode(float latitude, float longitude) {
this.latitude = latitude;
this.longitude = longitude;
}
public float getLatitude() {
return latitude;
}
public void setLatitude(float latitude) {
this.latitude = latitude;
}
public float getLongitude() {
return longitude;
}
public void setLongitude(float longitude) {
this.longitude = longitude;
}
}
public class GeoPolygon {
private ArrayList<Geocode> points;
public GeoPolygon() {
this.points = new ArrayList<Geocode>();
}
public GeoPolygon(ArrayList<Geocode> points) {
this.points = points;
}
public GeoPolygon add(Geocode geo) {
points.add(geo);
return this;
}
public boolean inside(Geocode geo) {
int i, j;
boolean c = false;
for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
(geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
c = !c;
}
return c;
}
}
推荐文章
- 确定记录是否存在的最快方法
- 阅读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的性能很差，或者我只是在处理一个糟糕的实现?
- 大概的成本访问各种缓存和主存储器?
