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循環）。