Блог

Я использую scikit.morphology для выполнения эрозии двумерного массива. Мне также нужно определить расстояние от каждой ячейки до минимального значения, указанного в эрозии.

Пример:

 np.reshape(np.arange(1,126,step=5),[5,5])
array([[  1,   6,  11,  16,  21],
       [ 26,  31,  36,  41,  46],
       [ 51,  56,  61,  66,  71],
       [ 76,  81,  86,  91,  96],
       [101, 106, 111, 116, 121]])

erosion(np.reshape(np.arange(1,126,step=5),[5,5]),selem=disk(3))
array([[ 1,  1,  1,  1,  6],
       [ 1,  1,  1,  6, 11],
       [ 1,  1,  1,  6, 11],
       [ 1,  6, 11, 16, 21],
       [26, 31, 36, 41, 46]])
  

Теперь я хочу также вернуть массив, который дает мне расстояние до минимума следующим образом:

 array([[ 0,  1,  2,  3,  3],
       [ 1,  1,  2,  3,  3],
       [ 2,  2,  3,  3,  3],
       [ 3,  3,  3,  3,  3],
       [ 3,  3,  3,  3,  3]])
  

Есть ли инструмент scikit, который может это сделать? Если нет, какие-либо советы о том, как эффективно достичь этого результата?

Вы можете найти расстояния от центра вашего следа, используя scipy.ndimage.distance_transform_cdt , а затем использовать SciPy ndimage.generic_filter для возврата этих значений:

 import numpy as np
from skimage.morphology import erosion, disk
from scipy import ndimage as ndi


input_arr = np.reshape(np.arange(1,126,step=5),[5,5])
footprint = disk(3)

def distance_from_min(values, distance_values):
    d = np.inf
    min_val = np.inf
    for i in range(len(values)):
        if values[i] <= min_val:
            min_val = values[i]
            d = distance_values[i]
    return d

full_footprint = np.ones_like(footprint, dtype=float)
full_footprint[tuple(i//2 for i in footprint.shape)] = 0
# use `ndi.distance_transform_edt` instead for the euclidean distance
distance_footprint = ndi.distance_transform_cdt(
    full_footprint, metric='taxicab'
)

# set values outside footprint to 0 for pretty-printing
distance_footprint[~footprint.astype(bool)] = 0
# then, extract it into values matching the values in generic_filter
distance_values = distance_footprint[footprint.astype(bool)]


output = ndi.generic_filter(
    input_arr.astype(float),
    distance_from_min,
    footprint=footprint,
    mode='constant',
    cval=np.inf,
    extra_arguments=(distance_values,),
)

print('input:n', input_arr)
print('footprint:n', footprint)
print('distance_footprint:n', distance_footprint)
print('output:n', output)
  

Что дает:

 input:
 [[  1   6  11  16  21]
 [ 26  31  36  41  46]
 [ 51  56  61  66  71]
 [ 76  81  86  91  96]
 [101 106 111 116 121]]
footprint:
 [[0 0 0 1 0 0 0]
 [0 1 1 1 1 1 0]
 [0 1 1 1 1 1 0]
 [1 1 1 1 1 1 1]
 [0 1 1 1 1 1 0]
 [0 1 1 1 1 1 0]
 [0 0 0 1 0 0 0]]
distance_footprint:
 [[0 0 0 3 0 0 0]
 [0 4 3 2 3 4 0]
 [0 3 2 1 2 3 0]
 [3 2 1 0 1 2 3]
 [0 3 2 1 2 3 0]
 [0 4 3 2 3 4 0]
 [0 0 0 3 0 0 0]]
output:
 [[0. 1. 2. 3. 3.]
 [1. 2. 3. 3. 3.]
 [2. 3. 4. 4. 4.]
 [3. 3. 3. 3. 3.]
 [3. 3. 3. 3. 3.]]
  

Однако эта функция будет очень медленной. Если вы хотите сделать это быстрее, вам понадобится (а) решение типа Numba или Cython для функции filter в сочетании с SciPy LowLevelCallables и (б) жестко закодировать массив расстояний в функцию distance, потому что для LowLevelCallables сложнее передать дополнительные аргументы. Вот полный пример, с llc-tools которым вы можете установить pip install numba llc-tools .

 import numpy as np
from scipy import ndimage as ndi
from skimage.morphology import erosion, disk
import llc


def filter_func_from_footprint(footprint):
    # first, create a footprint where the values are the distance from the
    # center
    full_footprint = np.ones_like(footprint, dtype=float)
    full_footprint[tuple(i//2 for i in footprint.shape)] = 0
    # use `ndi.distance_transform_edt` instead for the euclidean distance
    distance_footprint = ndi.distance_transform_cdt(
        full_footprint, metric='taxicab'
    )
    # then, extract it into values matching the values in generic_filter
    distance_footprint[~footprint.astype(bool)] = 0
    distance_values = distance_footprint[footprint.astype(bool)]

    # finally, create a filter function with the values hardcoded
    @llc.jit_filter_function
    def distance_from_min(values):
        d = np.inf
        min_val = np.inf
        for i in range(len(values)):
            if values[i] <= min_val:
                min_val = values[i]
                d = distance_values[i]
        return d

    return distance_from_min


if __name__ == '__main__':
    input_arr = np.reshape(np.arange(1,126,step=5),[5,5])
    footprint = disk(3)
    eroded = erosion(input_arr, selem=footprint)
    filter_func = filter_func_from_footprint(footprint)
    result = ndi.generic_filter(
        # use input_arr.astype(float) when using euclidean dist
        input_arr,
        filter_func,
        footprint=disk(3),
        mode='constant',
        cval=np.inf,
    )

    print('input:n', input_arr)
    print('output:n', result)
  

Что дает:

 input:
 [[  1   6  11  16  21]
 [ 26  31  36  41  46]
 [ 51  56  61  66  71]
 [ 76  81  86  91  96]
 [101 106 111 116 121]]
output:
 [[0 1 2 3 3]
 [1 2 3 3 3]
 [2 3 4 4 4]
 [3 3 3 3 3]
 [3 3 3 3 3]]
  

Для получения дополнительной информации о низкоуровневых вызываемых и llc-инструментах, в дополнение к документации LowLevelCallable на сайте SciPy (ссылка выше, плюс ссылки в нем), вы можете прочитать эти два сообщения в блоге, которые я написал несколько лет назад:

• Новый LowLevelCallable от SciPy меняет правила игры
• Более красивые вызовы низкого уровня с помощью Numba JIT и декораторов