Минимальное эвклидовое расстояние между мечеными компонентами в массиве 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() 

---

EDIT Просто испытания @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

Ах, понял это. spaces.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).