在numpy中查找两个多边形之间的距离

问题描述 投票:5回答:2

我有两个多边形P和Q,其中多边形的外部线性环由两个封闭的点集(以numpy数组存储)定义,它们以逆时针方向连接。 P和Q的格式如下:

P['x_coords'] = [299398.56 299402.16 299410.25 299419.7  299434.97 299443.75 299454.1 299465.3  299477.   299488.25 299496.8  299499.5  299501.28 299504. 299511.62 299520.62 299527.8  299530.06 299530.06 299525.12 299520.2 299513.88 299508.5  299500.84 299487.34 299474.78 299458.6  299444.66 299429.8  299415.4  299404.84 299399.47 299398.56 299398.56] 
P['y_coords'] = [822975.2  822989.56 823001.25 823005.3  823006.7  823005.06 823001.06 822993.4  822977.2  822961.   822943.94 822933.6  822925.06 822919.7 822916.94 822912.94 822906.6  822897.6  822886.8  822869.75 822860.75 822855.8  822855.4  822857.2  822863.44 822866.6  822870.6  822876.94 822886.8  822903.   822920.3  822937.44 822954.94 822975.2]

Q['x_coords'] = [292316.94 292317.94 292319.44 292322.47 292327.47 292337.72 292345.75 292350.   292352.75 292353.5  292352.25 292348.75 292345.75 292342.5 292338.97 292335.97 292333.22 292331.22 292329.72 292324.72 292319.44 292317.2  292316.2  292316.94]
Q['y_coords'] = [663781.   663788.25 663794.   663798.06 663800.06 663799.3  663796.56 663792.75 663788.5  663782.   663773.25 663766.   663762.   663758.25 663756.5  663756.25 663757.5  663761.   663763.75 663767.5  663769.5 663772.25 663777.5  663781.  ]

## SIMPLIFIED AND FORMATTED FOR EASY TESTING:
import numpy as np

px_coords = np.array([299398,299402,299410.25,299419.7,299398])
py_coords = np.array([822975.2,822920.3,822937.44,822954.94,822975.2])

qx_coords = np.array([292316,292331.22,292329.72,292324.72,292319.44,292317.2,292316])
qy_coords = np.array([663781,663788.25,663794,663798.06,663800.06,663799.3,663781])

P的外环是通过连接P['x_coords'][0], P['y_coords'][0] -> P['x_coords'][1], P['y_coords'][1]等形成的。每个数组的最后一个坐标与第一个相同,表示形状在拓扑上是封闭的。

是否可以使用numpy计算出P和Q外圈之间的简单最小距离?我一直在搜索SO上下两步,而没有找到任何明确的内容,因此我怀疑这可能是一个非常复杂的问题的过度简化。我知道可以使用现成的空间库(例如GDAL或Shapely)来完成距离计算,但是我很想通过在numpy中从头开始构建一些东西来了解它们的工作原理。

我已经考虑或尝试过的一些事情:

  • 计算两个数组中每个点之间的距离。这不起作用,因为P和Q之间的最接近点可以是边顶点对。使用scipy.spatial计算的每个形状的凸包都存在相同的问题。
  • 一种低效的蛮力方法,计算每对点之间以及边缘点对的每种组合之间的距离

是否有解决此问题的更好方法?

python numpy geometry
2个回答
1
投票

[k-d树上有many variations,用于存储具有范围的对象,例如多边形的[[edges。我最熟悉(但没有链接)的方法涉及将轴对齐的边界框与每个节点相关联。叶子与对象相对应,内部节点的框是包围两个子节点的框(通常重叠)的最小框。通常的中位数切割方法适用于对象框的中点(对于线段,这是它们的中点)。

已经为每个多边形构建了这些,下面的[[双重递归找到了最接近的方法:def closest(k1,k2,true_dist): return _closest(k1,0,k2,0,true_dist,float("inf")) def _closest(k1,i1,k2,i2,true_dist,lim): b1=k1.bbox[i1] b2=k2.bbox[i2] cc1=k1.child[i1] or (i1,) cc2=k2.child[i2] or (i2,) if len(cc1)==1 and len(cc2)==1: return min(lim,true_dist(i1,i2)) # Consider 2 or 4 pairs of children, possibly-closest first: for md,c1,c2 in sorted((min_dist(k1.bbox[c1],k2.bbox[c2]),c1,c2) for c1 in cc1 for c2 in cc2): if md>=lim: break lim=min(lim,_closest(k1,c1,k2,c2,true_dist,lim) return lim

注意:

两个不相交的线段之间的最接近方法必须涉及至少一个端点。

点与线段之间的距离可以大于点与包含线段的线之间的距离。

    不需要点对点检查:可以通过相邻的边缘找到这样的一对(四次)。
  • min_dist的参数可能重叠,在这种情况下,它必须返回0。

0
投票
pure python shape a edges: 33 shape b edges: 15 total loops: 1000 total time = 6.889256715774536 average time per loop = 0.006896152868643179 max time per loop = 0.022176027297973633 min time per loop = 0.0 cython loop shape a edges: 33 shape b edges: 15 total loops: 1000 total time = 0.18484902381896973 average time per loop = 0.00018503405787684656 max time per loop = 0.015654802322387695 min time per loop = 0.0

为了清楚起见,我在下面附加了纯Python版本的代码,如果需要,可以提供cython。

import numpy as np import time import math def segments_distance(x11, y11, x12, y12, x21, y21, x22, y22): if segments_intersect(x11, y11, x12, y12, x21, y21, x22, y22): return 0 distances = [] distances.append(point_segment_distance(x11, y11, x21, y21, x22, y22)) distances.append(point_segment_distance(x12, y12, x21, y21, x22, y22)) distances.append(point_segment_distance(x21, y21, x11, y11, x12, y12)) distances.append(point_segment_distance(x22, y22, x11, y11, x12, y12)) return min(distances) def segments_intersect(x11, y11, x12, y12, x21, y21, x22, y22): dx1 = x12 - x11 dy1 = y12 - y11 dx2 = x22 - x21 dy2 = y22 - y21 delta = dx2 * dy1 - dy2 * dx1 if delta == 0: return False # parallel segments s = (dx1 * (y21 - y11) + dy1 * (x11 - x21)) / delta t = (dx2 * (y11 - y21) + dy2 * (x21 - x11)) / (-delta) return (0 <= s <= 1) and (0 <= t <= 1) def point_segment_distance(px, py, x1, y1, x2, y2): dx = x2 - x1 dy = y2 - y1 if dx == dy == 0: # the segment's just a point return math.hypot(px - x1, py - y1) # Calculate the t that minimizes the distance. t = ((px - x1) * dx + (py - y1) * dy) / (dx * dx + dy * dy) # See if this represents one of the segment's # end points or a point in the middle. if t < 0: dx = px - x1 dy = py - y1 elif t > 1: dx = px - x2 dy = py - y2 else: near_x = x1 + t * dx near_y = y1 + t * dy dx = px - near_x dy = py - near_y return math.hypot(dx, dy) px_coords=np.array([299398.56,299402.16,299410.25,299419.7,299434.97,299443.75,299454.1,299465.3,299477.,299488.25,299496.8,299499.5,299501.28,299504.,299511.62,299520.62,299527.8,299530.06,299530.06,299525.12,299520.2,299513.88,299508.5,299500.84,299487.34,299474.78,299458.6,299444.66,299429.8,299415.4,299404.84,299399.47,299398.56,299398.56]) py_coords=np.array([822975.2,822989.56,823001.25,823005.3,823006.7,823005.06,823001.06,822993.4,822977.2,822961.,822943.94,822933.6,822925.06,822919.7,822916.94,822912.94,822906.6,822897.6,822886.8,822869.75,822860.75,822855.8,822855.4,822857.2,822863.44,822866.6,822870.6,822876.94,822886.8,822903.,822920.3,822937.44,822954.94,822975.2]) qx_coords=np.array([384072.1,384073.2,384078.9,384085.7,384092.47,384095.3,384097.12,384097.12,384093.9,384088.9,384082.47,384078.9,384076.03,384074.97,384073.53,384072.1]) qy_coords=np.array([780996.8,781001.1,781003.6,781003.6,780998.25,780993.25,780987.9,780981.8,780977.5,780974.7,780974.7,780977.2,780982.2,780988.25,780992.5,780996.8]) p_coords = np.dstack([px_coords.ravel(),py_coords.ravel()])[0][:-1] q_coords = np.dstack([qx_coords.ravel(),qy_coords.ravel()])[0][:-1] px_edges = np.stack((px_coords, np.roll(px_coords, -1)),1) py_edges = np.stack((py_coords, np.roll(py_coords, -1)),1) p_edges = np.stack((px_edges, py_edges), axis=-1)[:-1] qx_edges = np.stack((qx_coords, np.roll(qx_coords, -1)),1) qy_edges = np.stack((qy_coords, np.roll(qy_coords, -1)),1) q_edges = np.stack((qx_edges, qy_edges), axis=-1)[:-1] timings = [] for i in range(1,1000): start = time.time() edge_distances = [segments_distance(p_edges[n][0][0],p_edges[n][0][1],p_edges[n][1][0],p_edges[n][1][1],q_edges[m][0][0],q_edges[m][0][1],q_edges[m][1][0],q_edges[m][1][1]) for m in range(0,len(q_edges)) for n in range(0,len(p_edges))] end = time.time() - start timings.append(end) print(f'shape a edges: {len(px_coords)}') print(f'shape b edges: {len(qy_coords)}') print(f'total loops: {i+1}') print(f'total time = {sum(timings)}') print(f'average time per loop = {sum(timings)/len(timings)}') print(f'max time per loop = {max(timings)}') print(f'min time per loop = {min(timings)}')

© www.soinside.com 2019 - 2024. All rights reserved.