Minimaler euklidischer Abstand zwischen markierten Komponenten in einem Array

Question

Ich habe eine Reihe unterschiedlicher Formen in großen numpy Arrays, und ich möchte den euklidischen Abstand zwischen den Kanten unter Verwendung von numpy und scipy berechnen.

Hinweis:: Ich habe eine Suche durchgeführt und dies unterscheidet sich von vorherigen Fragen hier auf Stack, da ich den kleinsten Abstand zwischen markierten Patches innerhalb eines Arrays und nicht zwischen Punkten oder separaten Arrays erhalten möchte.

Mein aktueller Ansatz funktioniert mit einem KDTree, ist aber für große Arrays schrecklich ineffizient. Im Wesentlichen lese ich die Koordinaten jeder markierten Komponente nach und berechne den Abstand zwischen allen anderen Komponenten. Abschließend wird der durchschnittliche minimale Abstand als ein Beispiel berechnet.

Ich suche einen intelligenteren Ansatz mit Python und vorzugsweise ohne zusätzliche Module.

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 Gerade @morningsun vorgeschlagene Lösung getestet und es ist eine enorme Verbesserung der Geschwindigkeit. Allerdings kehrten die Werte sind etwas anders:

# 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

Ah, es herausgefunden. spatial.distance.pdist gibt keine korrekte Abstandsmatrix zurück und die Werte waren daher falsch.

Answer 1

Hier ist ein vollständig vektorisiert Weg, um die Distanzmatrix für die markierten Objekte zu finden:

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) 

Dieser Ansatz O (N) Speicher benötigt, wobei N die Anzahl von Nicht-Null-Pixel ist. Wenn dies zu anspruchsvoll ist, können Sie es entlang einer Achse "de-vektorisieren" (fügen Sie eine For-Schleife hinzu).