通过在numpy中设置一些条件来检索元素的位置

通过在numpy中设置一些条件来检索元素的位置

使用scipy，您可以将这样的点描述为既是其邻域的最大值又是最小值的点：

import numpy as npimport scipy.ndimage.filters as filtersdef using_filters(data):    return np.where(np.logical_and.reduce(        [data == f(data, footprint=np.ones((3,3)), mode='constant', cval=np.inf)         for f in (filters.maximum_filter, filters.minimum_filter)]))using_filters(data)# (array([2, 3]), array([5, 9]))

仅使用numpy，您可以将

data

其自身的8个移位切片进行比较以找到相等的点：

def using_eight_shifts(data):    h, w = data.shape    data2 = np.empty((h+2, w+2))    data2[(0,-1),:] = np.nan    data2[:,(0,-1)] = np.nan    data2[1:1+h,1:1+w] = data    result = np.where(np.logical_and.reduce([        (data2[i:i+h,j:j+w] == data)        for i in range(3)        for j in range(3)        if not (i==1 and j==1)]))    return result

如您在上面所看到的，此策略创建了一个扩展数组，该数组在数据周围具有NaN边界。这允许将移位的切片表示为

data2[i:i+h,j:j+w]

。

如果您知道要与邻居进行比较，则可能应该

data

从一开始就用NaN的边界进行定义，这样就不必像上面那样做第二个数组了。

使用八移位（和比较）比循环遍历每个单元格

data

并将其与相邻单元格进行比较要快得多：

def using_quadratic_loop(data):    return np.array([[i,j] for i in range(1,np.shape(data)[0]-1) for j in range(1,np.shape(data)[1]-1) if np.all(data[i-1:i+2,j-1:j+2]==data[i,j])]).T

这是一个基准：

using_filters : 0.130using_eight_shifts       : 0.340using_quadratic_loop     : 18.794

---

这是用于产生基准的代码：

import timeitimport operatorimport numpy as npimport scipy.ndimage.filters as filtersimport matplotlib.pyplot as pltdata  = np.array([    [0,1,2,3,4,7,6,7,8,9,10],     [3,3,3,4,7,7,7,8,11,12,11],      [3,3,3,5,7,7,7,9,11,11,11],    [3,4,3,6,7,7,7,10,11,11,11],    [4,5,6,7,7,9,10,11,11,11,11]    ])data = np.tile(data, (50,50))def using_filters(data):    return np.where(np.logical_and.reduce(        [data == f(data, footprint=np.ones((3,3)), mode='constant', cval=np.inf)         for f in (filters.maximum_filter, filters.minimum_filter)]))def using_eight_shifts(data):    h, w = data.shape    data2 = np.empty((h+2, w+2))    data2[(0,-1),:] = np.nan    data2[:,(0,-1)] = np.nan    data2[1:1+h,1:1+w] = data    result = np.where(np.logical_and.reduce([        (data2[i:i+h,j:j+w] == data)        for i in range(3)        for j in range(3)        if not (i==1 and j==1)]))    return resultdef using_quadratic_loop(data):    return np.array([[i,j] for i in range(1,np.shape(data)[0]-1) for j in range(1,np.shape(data)[1]-1) if np.all(data[i-1:i+2,j-1:j+2]==data[i,j])]).Tnp.testing.assert_equal(using_quadratic_loop(data), using_filters(data))np.testing.assert_equal(using_eight_shifts(data), using_filters(data))timing = dict()for f in ('using_filters', 'using_eight_shifts', 'using_quadratic_loop'):    timing[f] = timeit.timeit('{f}(data)'.format(f=f),        'from __main__ import data, {f}'.format(f=f),        number=10)for f, t in sorted(timing.items(), key=operator.itemgetter(1)):    print('{f:25}: {t:.3f}'.format(f=f, t=t))