6k次浏览,没有带有代码示例的答案。
我曾经专门研究过这个问题,所以我觉得我应该把它从死亡中带回来,帮助下一个碰巧遇到这个问题的人;
Marco13的答案是正确的方法;使用四叉树。
(我的工作基于
https://github.com/Josephbakulikira/QuadTree-Visualization-with-python--pygame)
我在Python中实现了一个四叉树,这样我就可以在地图数据库上检查离特定地址最近的街道。我认为这基本上代表了OP想要做的事情。
我的四叉树接受线段和点。在我的情况下,它们是道路段和房屋位置。
然后,我可以在特定点上使用一定半径查询四叉树,以选择靠近该点的道路段,并获得最近的段。
以下是代码:
quadtree.py:
from shapes import Rectangle, Vector2, Line, Circle
class QuadTree:
"""
Python representation of a QuadTree
"""
def __init__(self, capacity : int, boundary : Rectangle) -> None:
self.capacity : int = capacity
self.boundary : Rectangle = boundary
self.elements = set([])
self.northWest : QuadTree = None
self.northEast : QuadTree = None
self.southWest : QuadTree = None
self.southEast : QuadTree = None
return
@property
def children(self) -> list | bool:
"""
Checker if tree has children.
If yes, returns them in list
"""
if self.northWest != None:
return [self.northWest, self.northEast, self.southWest, self.southEast]
else:
return False
def subdivide(self) -> None:
"""
Divide current tree into quadrants (children)
"""
parent = self.boundary
differencex = (parent.bottomRight.x - parent.topLeft.x) / 2
differencey = (parent.bottomRight.y - parent.topLeft.y) / 2
boundary_nw = Rectangle(
Vector2(parent.topLeft.x , parent.topLeft.y),
Vector2(parent.bottomRight.x - differencex, parent.bottomRight.y - differencey)
)
boundary_ne = Rectangle(
Vector2(parent.topLeft.x + differencex, parent.topLeft.y),
Vector2(parent.bottomRight.x, parent.bottomRight.y - differencey)
)
boundary_sw = Rectangle(
Vector2(parent.topLeft.x, parent.topLeft.y + differencey),
Vector2(parent.topLeft.x + differencex, parent.bottomRight.y)
)
boundary_se = Rectangle(
Vector2(parent.topLeft.x + differencex, parent.topLeft.y + differencey),
Vector2(parent.bottomRight.x, parent.bottomRight.y)
)
self.northWest = QuadTree(self.capacity, boundary_nw)
self.northEast = QuadTree(self.capacity, boundary_ne)
self.southWest = QuadTree(self.capacity, boundary_sw)
self.southEast = QuadTree(self.capacity, boundary_se)
for element in self.elements:
self.insertinchildren(element)
self.elements = set([])
return
def insert(self, element) -> bool:
"""
Insert element into tree.
If tree is at max capacity, subdivide it, and add elements to children
"""
if self.boundary.containsElement(element) == False:
return False
if len(self.elements) < self.capacity and not self.children:
self.elements.add(element)
return True
else:
if len(self.elements) == self.capacity and not self.children:
if type(element) == Line:
if element in self.elements:
return True
if not self.children:
self.subdivide()
return self.insertinchildren(element)
def insertinchildren(self, element) -> bool:
"""
Insert element into children
(Insert lines everywhere they intersect, while
inserting points(Vectors) only in one quadrant they are in)
"""
if isinstance(element, Line):
self.northWest.insert(element)
self.northEast.insert(element)
self.southWest.insert(element)
self.southEast.insert(element)
return True
else:
if not self.northWest.insert(element):
if not self.northEast.insert(element):
if not self.southWest.insert(element):
if not self.southEast.insert(element):
return False
return True
def query(self, _range : Rectangle | Circle) -> set:
"""
Use a circle or rectangle to query the tree for
elements in range
"""
elements_in_range = set([])
if _range.intersects(self.boundary):
if self.children:
elements_in_range.update(self.northWest.query(_range))
elements_in_range.update(self.northEast.query(_range))
elements_in_range.update(self.southWest.query(_range))
elements_in_range.update(self.southEast.query(_range))
else:
for element in self.elements:
if _range.containsElement(element):
elements_in_range.add(element)
return elements_in_range
shapes.py
import numpy as np
class Vector2(np.ndarray):
"""
Numpy array subclass to represent a 2D vector
"""
def __new__(cls, x, y) -> np.ndarray:
return np.array([x, y]).view(cls)
@property
def x(self) -> float:
return self[0]
@property
def y(self) -> float:
return self[1]
def __hash__(self) -> int:
return id(self)
def __eq__(self, other) -> bool:
if not isinstance(other, Vector2):
return False
return np.array_equal(self, other)
class Line:
"""
Representation of a line segment
"""
def __init__(self, point1 : Vector2, point2 : Vector2) -> None:
self.a = point1
self.b = point2
return
def __str__(self) -> str:
return f"Line({self.a}, {self.b})"
def distancetoPoint(self, p : Vector2) -> float:
"""
Get the shortest distance between point p and the line segment
"""
return np.linalg.norm(self.closestPoint(p) - p)
def closestPoint(self, p : Vector2) -> float:
"""
Get the closest point on the line segment from point p
"""
a = self.a
b = self.b
cp = None
ab = b - a
ap = p - a
proj = np.dot(ap, ab)
abLenSq = ab[0]**2 + ab[1]**2
normalizedProj = proj / abLenSq
if normalizedProj <= 0:
cp = a
elif normalizedProj >= 1:
cp = b
else:
cp = a + normalizedProj * ab
return cp
def onSegment(self, p, q, r) -> bool:
"""
Given three collinear points p, q, r, the function checks if
point q lies on line segment 'pr'
"""
if ( (q.x <= max(p.x, r.x)) and (q.x >= min(p.x, r.x)) and
(q.y <= max(p.y, r.y)) and (q.y >= min(p.y, r.y))):
return True
return False
def orientation(self, p, q, r) -> int:
"""
To find the orientation of an ordered triplet (p,q,r)
function returns the following values:
0 : Collinear points
1 : Clockwise points
2 : Counterclockwise
See https://www.geeksforgeeks.org/orientation-3-ordered-points/amp/
for details of below formula.
"""
val = (float(q.y - p.y) * (r.x - q.x)) - (float(q.x - p.x) * (r.y - q.y))
if (val > 0):
return 1
elif (val < 0):
return 2
else:
return 0
def intersects(self, p1 : Vector2, p2 : Vector2) -> bool:
"""
Check if current line intersects with another line
made of point1 and point2
"""
o1 = self.orientation(self.a, self.b, p1)
o2 = self.orientation(self.a, self.b, p2)
o3 = self.orientation(p1, p2, self.a)
o4 = self.orientation(p1, p2, self.b)
if ((o1 != o2) and (o3 != o4)):
return True
if ((o1 == 0) and self.onSegment(self.a, p1, self.b)):
return True
if ((o2 == 0) and self.onSegment(self.a, p2, self.b)):
return True
if ((o3 == 0) and self.onSegment(p1, self.a, p2)):
return True
if ((o4 == 0) and self.onSegment(p1, self.b, p2)):
return True
return False
class Rectangle:
"""
Rectangle used as boundaries for the QuadTree
or as a query range for the QuadTree
"""
def __init__(self, topLeft : Vector2, bottomRight : Vector2) -> None:
self.topLeft = topLeft
self.bottomRight = bottomRight
self.topRight = Vector2(self.bottomRight.x, self.topLeft.y)
self.bottomLeft = Vector2(self.topLeft.x, self.bottomRight.y)
self.center = (self.topLeft + self.bottomRight) / 2
self.top = Line(self.topLeft, Vector2(self.bottomRight.x, self.topLeft.y))
self.bottom = Line(Vector2(self.topLeft.x, self.bottomRight.y), self.bottomRight)
self.left = Line(self.topLeft, Vector2(self.topLeft.x, self.bottomRight.y))
self.right = Line(Vector2(self.bottomRight.x, self.topLeft.y), self.bottomRight)
return
def __str__(self) -> str:
return f"Rectangle({self.topLeft}, {self.bottomRight})"
def containsElement(self, element : Vector2 | Line) -> bool:
"""
Checks if the element is contained by this rectangle
"""
if isinstance(element, Line):
return self.containsLine(element)
x, y = element
bx , by = self.topLeft
cx, cy = self.bottomRight
if x >= bx and x <= cx and y >= by and y <= cy:
return True
else:
return False
def containsLine(self, line : Line) -> bool:
"""
Checks if the line goes through this rectangle
"""
if self.containsElement(line.a) or self.containsElement(line.b):
return True
else:
if not line.intersects(self.top.a, self.top.b):
if not line.intersects(self.bottom.a, self.bottom.b):
if not line.intersects(self.left.a, self.left.b):
if not line.intersects(self.right.a, self.right.b):
return False
return True
def intersects(self, other : 'Rectangle') -> bool:
"""
Check if this rectangle intersects with another rectangle
ie. if any of the sides of the rectangle crosses another rectangle
"""
return not (self.bottomRight.x < other.topLeft.x
or self.topLeft.x > other.bottomRight.x
or self.bottomRight.y < other.topLeft.y
or self.topLeft.y > other.bottomRight.y)
class Circle:
"""
Circle used as range to query the QuadTree
"""
def __init__(self, center : Vector2, radius : float) -> None:
self.center = center
self.radius = radius
return
def __str__(self) -> str:
return f"Circle({self.center}, r{self.radius})"
def containsElement(self, element : Vector2 | Line) -> bool:
"""
Checks if the element is contained by this circle
"""
if isinstance(element, Line):
return self.containsLine(element)
dist = np.linalg.norm(self.center - element)
if dist <= self.radius:
return True
else:
return False
def containsLine(self, line : Line) -> bool:
"""
Checks if the line goes through this circle
"""
if line.distancetoPoint(self.center) <= self.radius:
return True
else:
return False
def intersects(self, rectangle : Rectangle) -> bool:
"""
Check if this circle intersects with a rectangle
"""
if not rectangle.containsElement(self.center):
if not self.containsLine(rectangle.top):
if not self.containsLine(rectangle.bottom):
if not self.containsLine(rectangle.left):
if not self.containsLine(rectangle.right):
return False
return True
main.py
from quadtree import QuadTree
from shapes import Rectangle, Circle, Vector2, Line
node0 = Vector2(107.1360616, 46.1032011)
node1 = Vector2(107.13563750000003, 46.102939)
node2 = Vector2(107.13538800000003, 46.1028447)
node3 = Vector2(107.13493060000002, 46.1027272)
node4 = Vector2(107.13318809999998, 46.1021657)
node5 = Vector2(107.13255700000002, 46.1019075)
node6 = Vector2(107.1316817, 46.1014684)
node7 = Vector2(107.13088449999998, 46.1011467)
node8 = Vector2(107.1304002, 46.1010027)
node9 = Vector2(107.1299055, 46.1009219)
node10 = Vector2(107.12915320000002, 46.1008729)
node11 = Vector2(107.12877689999999, 46.1008901)
node12 = Vector2(107.12828999999999, 46.10097)
node13 = Vector2(107.1274209, 46.1012472)
node14 = Vector2(107.1270212, 46.1014094)
node15 = Vector2(107.1265209, 46.101732)
node16 = Vector2(107.12579249999999, 46.1022694)
node17 = Vector2(107.1256219, 46.1024435)
node18 = Vector2(107.12546250000003, 46.1026809)
node19 = Vector2(107.12500690000002, 46.1034946)
nodes = [node0, node1, node2, node3, node4, node5, node6, node7, node8, node9,
node10, node11, node12, node13, node14, node15, node16, node17, node18, node19]
lines = []
for i in range(len(nodes)-1):
lines.append(Line(nodes[i], nodes[i+1]))
house0 = Vector2(107.13595197486893, 46.10285045021604)
house1 = Vector2(107.1300706177878, 46.101115083412374)
house2 = Vector2(107.13043003379954, 46.100917943095524)
house3 = Vector2(107.13046758472615, 46.101234111186926)
house4 = Vector2(107.13071434795808, 46.10100721426972)
house5 = Vector2(107.13108449280605, 46.10108904605243)
house6 = Vector2(107.12480471447458, 46.10391049893959)
house7 = Vector2(107.12635081112013, 46.10232085692584)
houses = [house0, house1, house2, house3, house4, house5, house6, house7]
tree = QuadTree(2, Rectangle(Vector2(107.12383, 46.09690), Vector2(107.13709, 46.10700)))
for l in lines:
tree.insert(l)
for h in houses:
tree.insert(h)
queryAT = Vector2(107.12739, 46.10358)
radQuery = Circle(queryAT, 0.0022)
found_elements = tree.query(radQuery)
print("Found elements:")
for element in found_elements:
print(element)
shortest_dist = 999999
closest_line = None
for element in found_elements:
if isinstance(element, Line):
dist = element.distancetoPoint(queryAT)
if dist < shortest_dist:
shortest_dist = dist
closest_line = element
print("\nClosest line: " + str(closest_line))
import matplotlib.pyplot as plt
import matplotlib.patches as patches
fig, ax = plt.subplots()
ax.set_xlim(tree.boundary.topLeft.x, tree.boundary.bottomRight.x)
ax.set_ylim(tree.boundary.topLeft.y, tree.boundary.bottomRight.y)
def plot_tree(ax, quadtree):
"""
Traverse the quadtree recursively and
plot the boundaries and elements of each node
"""
if quadtree.children:
for child in quadtree.children:
plot_tree(ax, child)
else:
ax.add_patch(patches.Rectangle((quadtree.boundary.topLeft.x, quadtree.boundary.topLeft.y), quadtree.boundary.bottomRight.x-quadtree.boundary.topLeft.x, quadtree.boundary.bottomRight.y - quadtree.boundary.topLeft.y, fill=False))
for element in quadtree.elements:
clr = 'black'
if element in found_elements:
clr = 'magenta'
if element is closest_line:
clr = 'red'
if isinstance(element, Line):
ax.plot([element.a.x, element.b.x], [element.a.y, element.b.y], color=clr)
ax.plot(element.a.x, element.a.y, 'o', color=clr)
ax.plot(element.b.x, element.b.y, 'o', color=clr)
elif isinstance(element, Vector2):
ax.plot(element.x, element.y, 'go')
return
plot_tree(ax, tree)
ax.add_patch(patches.Circle((queryAT.x, queryAT.y), radQuery.radius, fill=False, color='orange'))
ax.plot(queryAT.x, queryAT.y, 'o', color='orange')
ax.set_box_aspect(1)
plt.show()
在main.py中,我使用坐标点并设置了一个面积约为1平方公里的四叉树,其中只有一条由多个线段组成的街道。
然后,我查询了一个特定的点(107.12739,46.10358)经度,纬度,距离为0.0022,约为170米。
这给我带来了以下结果:
{{link1:
![enter image description here](https://istack.dev59.com/ha3Hy.webp)
}}
在matplotlib中的可视化:
橙色是查询圆
洋红色是查询找到的线段
红线是最近的线段。
![enter image description here](https://istack.dev59.com/0GB9Z.webp)
![enter image description here](https://istack.dev59.com/uU80x.webp)
真的对于1000个片段来说更有效吗?我不知道。不管怎样,给你了。