跨kdtrees的对偶递归以找到两组点之间的最接近点

Question

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

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

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

第一次尝试-为上下文提供完整代码，此问题仅适用于函数“closest”：

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)

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

第二次尝试-只有“最近的”函数，可以用上面代码中的namesake替换掉：

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）和最长距离

有不同的方法来组织这两棵树的细分，例如

• 对于红色树的每个细分级别，将蓝色树细分为树叶，

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

• 或者在每个细分级别上，同时细分红色和蓝色，并考虑所有的组合

我不知道如何在这些选项中进行选择（除了完全尝试之外）