python - 跨 kd 树的双重递归以查找两组点之间最接近的方法

Question

我为两组点构建了 kd 树，以便找到两组点之间最接近的双色配对:

kd 树存储为 python 字典，可以在下面的代码中找到，并传递给一个函数('closest')，该函数旨在同时递归地分析两棵树以找到集合之间最接近的方法。这是为了防止必须暴力解决问题。

我的第一次尝试是基于this question的答案。通过这次尝试，我找不到强制函数在遇到叶子时“反弹”的条件，即旨在返回叶子与现有最小值之间的最小距离的 if 语句是从未达到。

第一次尝试 - 为上下文提供了完整的代码，这个问题仅与“最接近”的函数有关:

from operator import itemgetter
import math
import time
import pprint
import numpy as np


# builds the trees
def build_kd_tree(ar, depth=0, k=2):
    if len(ar) <= 0:
        return None
    axis = depth % k
    sorted_ar = sorted(ar, key=itemgetter(axis))
    idx = int(math.floor(len(ar)/2))
    return {
       'point': sorted_ar[idx],
       'left': build_kd_tree(sorted_ar[:idx], depth + 1),
       'right': build_kd_tree(sorted_ar[idx+1:], depth + 1)
    }


def min_dist(p1, p2):
    d1 = math.hypot(p1[0] - p2[0], p1[1] - p2[1])
    return d1


# function designed to simultaneously recurse two trees to find the closest approach
def closest(k1,k2,lim=float("inf")):

    cc1 = [k1[value] for value in k1 if k1[value] is not None and type(k1[value]) == dict]
    cc2 = [k2[value] for value in k2 if k2[value] is not None and type(k2[value]) == dict]

    if len(cc1) == 0 and len(cc2) == 0:
        return min(lim, min_dist(k1['point'], k2['point']))

    for md, c1, c2 in sorted((min_dist(c1['point'], c2['point']), c1, c2) for c1 in cc1 for c2 in cc2):
        if md >= lim: break
        lim = min(lim, closest(c1, c2, lim))
    return lim

# some example coordinates
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])

# some more example coordinates
#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])

# this is all just formatting the coordinates - only important thing to know is that p_midpoints and q_midpoints are two distinct sets of points, and are the targets in this question
px_edges = np.stack((px_coords, np.roll(px_coords, -1)),1)
px_midpoints = np.array(abs(px_coords + np.roll(px_coords, -1))/2)
py_edges = np.stack((py_coords, np.roll(py_coords, -1)),1)
py_midpoints = np.array(abs(py_coords + np.roll(py_coords, -1))/2)

p_edges = np.stack((px_edges, py_edges), axis=-1)[:-1]
p_midpoints = np.stack((px_midpoints, py_midpoints), axis=-1)[:-1]

qx_edges = np.stack((qx_coords, np.roll(qx_coords, -1)),1)
qx_midpoints = np.array(abs(qx_coords + np.roll(qx_coords, -1))/2)
qy_edges = np.stack((qy_coords, np.roll(qy_coords, -1)),1)
qy_midpoints = np.array(abs(qy_coords + np.roll(qy_coords, -1))/2)

q_edges = np.stack((qx_edges, qy_edges), axis=-1)[:-1]
q_midpoints = np.stack((qx_midpoints, qy_midpoints), axis=-1)[:-1]

# where the tree is actually built
p_tree = build_kd_tree(p_midpoints)
q_tree = build_kd_tree(q_midpoints)

# uncommect to see structure of tree
#pprint.pprint(p_tree)

near_distance = closest(p_tree, q_tree)

# brute force for testing
#distances = []
#for p_point in p_midpoints:
#    for q_point in q_midpoints:
#        distances.append(min_dist(p_point, q_point))
#
#m_dist = sorted(distances)[0]
#print(m_dist)

在我的第二次尝试中，我试图强制函数在遇到树叶时停止递归。这适用于两个样本坐标集中较小的一个，但不适用于两个样本坐标集中较大的一个，并因同样的问题而失败。

第二次尝试 - 只有“最接近”的函数，可以与上面代码中的同名函数进行类似的交换:

def closest(k1,k2,lim=float("inf")):
    cc1 = [k1]
    cc1 = cc1 + [k1[value] for value in k1 if k1[value] is not None and type(k1[value]) == dict]
    cc2 = [k2]
    cc2 = cc2 + [k2[value] for value in k2 if k2[value] is not None and type(k2[value]) == dict]

    if len(cc1) == 1 and len(cc2) == 1:
        return min(lim, min_dist(k1['point'], k2['point']))

    md = [[min_dist(cc1[i]['point'], cc2[j]['point']), i, j, (cc1[i]['point'], cc2[j]['point'])] for i in range(len(cc1) >> 1, len(cc1)) for j in range(len(cc1) >> 1, len(cc2))]
    md = sorted(md, key=itemgetter(0))
    for h in range(0, len(md)):
        lim = min(lim, closest(cc1[md[h][1]], cc2[md[h][2]],lim))
    return lim

我知道存在开箱即用的解决方案来解决这个问题，但这是一个我想通过从头开始构建自己的解决方案来更好地理解的领域。任何帮助表示赞赏。

Answer 1

kD 树的工作原理是，您可以快速找到查询点(假设它是红色的)到已知矩形(假设排列在一个矩形中)中包含的点的子集的最短和最长距离的边界。蓝树)。此外，矩形是通过连续划分获得的，这使得估计的计算更加简单。

如果你想适应双​​色情况，你可以处理红色树生成的矩形而不是单个红色点，并调整规则来估计到蓝色的最短距离(在重叠的情况下为 0)和最长距离矩形。

有不同的方式来组织两棵树的分割，例如

• 对于红树的每个分割级别，将蓝树分割到叶子，

• 相反，对于蓝色树的每个分割级别，将红色树分割到叶子，

• 或者在每个分割级别上，分割红色和蓝色并考虑所有组合。

我不知道如何在这些选项中进行选择(除了充分尝试它们)。