Минимальное евклидово расстояние от края до края между помеченными компонентами в массиве numpy

У меня есть несколько различных форм в больших массивах 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

ИЗМЕНИТЬ 2

А, разобрался. space.distance.pdist не возвращает правильную матрицу расстояний, поэтому значения были неверными.


numpy scipy image-processing euclidean-distance connected-components
person Curlew    schedule 14.05.2016    source источник

Ответы (1)


arrow_upward
5
arrow_downward

Вот полностью векторизованный способ найти матрицу расстояний для помеченных объектов:

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(N2) памяти, где N — количество ненулевых пикселей. Если это слишком сложно, вы можете «девекторизовать» его вдоль одной оси (добавить цикл for).

person Community    schedule 14.05.2016
comment
Спасибо вам за это! Я только что проверил его на одном из своих наборов данных, и он работает почти на 89% быстрее. Сила векторизации. Хотя я не совсем понимаю, почему вычислялся «кваклид». Он также возвращает разные значения, если попытаться вычислить, например, среднее значение всех различий (см. редактирование в вопросе). - person Curlew; 14.05.2016
comment
Ааа, разобрался (см выше). Pdist не возвращает правильную матрицу расстояний, поэтому мои предыдущие значения были неправильными... Еще раз спасибо за ваше решение! - person Curlew; 14.05.2016
comment
@Curlew - евклидов в квадрате вычисляется быстрее. Обратите внимание, что я использовал его только для промежуточных результатов; квадратный корень берется в операторе return. - person ; 14.05.2016