numpy数组中标记组件之间的最小边对边欧氏距离

Question

我在大型numpy数组中具有许多不同的形式，我想使用numpy和scipy来计算它们之间的边到边的欧氏距离。

注意：我做了搜索，这是从现在开始堆以前其它问题不同，因为我想获得一个阵列中，而不是点或单独阵列之间的标记斑块之间的最小距离为其他问题都问。

我目前的方法使用KDTree，但对于大型阵列来说效率非常低。本质上，我正在查找每个标记组件的坐标并计算所有其他组件之间的距离。最后，以平均最小距离为例进行计算。

我正在寻找一个更聪明的方法使用python，最好没有任何额外的模块。

import numpy 
from scipy import spatial 
from scipy import ndimage 

# Testing array 
a = numpy.zeros((8,8), dtype=numpy.int) 
a[2,2] = a[3,1] = a[3,2] = 1 
a[2,6] = a[2,7] = a[1,6] = 1 
a[5,5] = a[5,6] = a[6,5] = a[6,6] = a[7,5] = a[7,6] = 1  

# label it 
labeled_array,numpatches = ndimage.label(a) 

# For number of patches 
closest_points = [] 
for patch in [x+1 for x in range(numpatches)]: 
# Get coordinates of first patch 
    x,y = numpy.where(labeled_array==patch) 
    coords = numpy.vstack((x,y)).T # transform into array 
    # Built a KDtree of the coords of the first patch 
    mt = spatial.cKDTree(coords) 

    for patch2 in [i+1 for i in range(numpatches)]: 
     if patch == patch2: # If patch is the same as the first, skip 
      continue 
     # Get coordinates of second patch 
     x2,y2 = numpy.where(labeled_array==patch2) 
     coords2 = numpy.vstack((x2,y2)).T 

     # Now loop through points 
     min_res = [] 
     for pi in range(len(coords2)): 
      dist, indexes = mt.query(coords2[pi]) # query the distance and index 
      min_res.append([dist,pi]) 
     m = numpy.vstack(min_res) 
     # Find minimum as closed point and get index of coordinates 
     closest_points.append(coords2[m[numpy.argmin(m,axis=0)[0]][1]]) 


# The average euclidean distance can then be calculated like this: 
spatial.distance.pdist(closest_points,metric = "euclidean").mean() 

---

编辑 只是测试@morningsun提出的解决方案，这是一个巨大的速度提升。但是返回的值略有不同：

# Consider for instance the following array 
a = numpy.zeros((8,8), dtype=numpy.int) 
a[2,2] = a[2,6] = a[5,5] = 1 

labeled_array, numpatches = ndimage.label(cl_array,s) 

# Previous approach using KDtrees and pdist 
b = kd(labeled_array,numpatches) 
spatial.distance.pdist(b,metric = "euclidean").mean() 
#> 3.0413115592767102 

# New approach using the lower matrix and selecting only lower distances 
b = numpy.tril(feature_dist(labeled_array)) 
b[b == 0 ] = numpy.nan 
numpy.nanmean(b) 
#> 3.8016394490958878 

EDIT 2

啊，想通了。 spatial.distance.pdist不返回适当的距离矩阵，因此这些值是错误的。

Answer 1

这里找到距离矩阵用于标记对象完全矢量化方式：

import numpy as np 
from scipy.spatial.distance import cdist 

def feature_dist(input): 
    """ 
    Takes a labeled array as returned by scipy.ndimage.label and 
    returns an intra-feature distance matrix. 
    """ 
    I, J = np.nonzero(input) 
    labels = input[I,J] 
    coords = np.column_stack((I,J)) 

    sorter = np.argsort(labels) 
    labels = labels[sorter] 
    coords = coords[sorter] 

    sq_dists = cdist(coords, coords, 'sqeuclidean') 

    start_idx = np.flatnonzero(np.r_[1, np.diff(labels)]) 
    nonzero_vs_feat = np.minimum.reduceat(sq_dists, start_idx, axis=1) 
    feat_vs_feat = np.minimum.reduceat(nonzero_vs_feat, start_idx, axis=0) 

    return np.sqrt(feat_vs_feat) 

这种方法需要O（N 2 ）存储器，其中，N是非零的像素的数量。如果这太苛刻，可以沿着一个轴“去矢量化”（添加for循环）。