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


问题内容

对于给定的2d数据数组,如何以粗体显示7和11的位置(索引)。因为只有它们是邻居中被相同值包围的元素

  import numpy as np
    data  = 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]
        ])

    print data

问题答案:

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

import numpy as np
import scipy.ndimage.filters as filters

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)]))

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.130
using_eight_shifts       : 0.340
using_quadratic_loop     : 18.794

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

import timeit
import operator
import numpy as np
import scipy.ndimage.filters as filters
import matplotlib.pyplot as plt

data  = 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 result


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

np.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))