python/numpy最快的方法,用于对掩码数组进行2d内核排名过滤(和/或选择性排名)

Question

给定2D numpy数组

MyArray = np.array([[ 8.02,  9.54,  0.82,  7.56,  2.26,  9.47],
           [ 2.68,  7.3 ,  2.74,  3.03,  2.25,  8.84],
           [ 2.21,  3.62,  0.55,  2.94,  5.77,  0.21],
           [ 5.78,  5.72,  8.85,  0.24,  5.37,  9.9 ],
           [ 9.1 ,  7.21,  4.14,  9.95,  6.73,  6.08],
           [ 1.8 ,  5.14,  5.02,  6.52,  0.3 ,  6.11]])

和一个掩码阵列

MyMask =  np.array([[ 0.,  0.,  1.,  1.,  0.,  1.],
            [ 1.,  0.,  0.,  0.,  0.,  1.],
            [ 0.,  0.,  0.,  1.,  0.,  0.],
            [ 0.,  1.,  1.,  1.,  1.,  0.],
            [ 0.,  1.,  0.,  1.,  0.,  0.],
            [ 0.,  1.,  0.,  0.,  1.,  1.]])

我希望运行一个"漏洞"中值过滤器,忽略蒙面元素.

例如,带有内核的秩过滤器

k = np.array([[ 1, 1, 1],
              [ 1, 0, 1],
              [ 1, 1, 1]]); 

将运行MyArray:对每个元素的内核定义的邻域进行排序,MyArray并返回仅非掩码元素的中值(如果数组是偶数则求平均值).

现在,我现在在unpythonic循环中这样做,使用bottleneck.nanmedian将掩码映射到NaNs.这正是我所需要的,但我希望依赖于2D数组操作例程.

scipy.signal.order_filter并且scipy.ndimage.filters.rank_filter都可用(rank_filter看起来要快得多),但看起来它们在返回排名和偏置结果之前排序NaN并Inf位于数组的顶部.似乎这些方法都不支持numpy.ma数组(屏蔽),也不接受选择性排名数组(然后我可以用0填充所有掩码并抵消我的排名),也没有明显的方法来改变内核每个地点.

我想知道我是否错过了组合和/或python功能,或者我是否应该在Cython中实现新的例程.

忽略边界处理,上述问题的内部要点将是

[[ 0.     0.     0.     0.     0.     0.   ]
 [ 0.     3.18   3.62   2.26   2.645  0.   ]
 [ 0.     2.74   3.325  2.74   2.64   0.   ]
 [ 0.     3.88   3.62   4.955  6.08   0.   ]
 [ 0.     5.02   5.77   5.77   6.52   0.   ]
 [ 0.     0.     0.     0.     0.     0.   ]]

Answer 1

一种方法是牺牲 RAM 使用量来放弃 Python 循环。即我们炸毁原始数组，以便我们可以立即将过滤器应用于所有子数组。这有点类似于Numpy 广播。

在我的测试中，对于 1000x1000 的数组，矢量化函数的执行速度大约快 100 倍。

在我的代码中，我使用NaN's 进行屏蔽，但通过一些额外的代码行，您也可以使用numpy.ma数组。而且我没有nanmedian函数，所以我使用了nanmean，但性能应该是可比的。

import numpy as np
from numpy.lib.stride_tricks import as_strided

# test data
N = 1000
A = np.random.rand(N, N)*10
mask = np.random.choice([True, False], size=(N, N))

def filter_loop(A, mask):
    kernel = np.array([[1,1,1],[1,0,1],[1,1,1]], bool)
    A = A.copy()
    A[mask] = np.nan
    N = A.shape[0] - 2  # assuming square matrix
    out = np.empty((N, N))
    for i in xrange(N):
        for j in xrange(N):
            out[i,j] = np.nanmean(A[i:i+3, j:j+3][kernel])
    return out    

def filter_broadcast(A, mask):
    A = A.copy()
    A[mask] = np.nan
    N = A.shape[0] - 2
    B = as_strided(A, (N, N, 3, 3), A.strides+A.strides)
    B = B.copy().reshape((N, N, 3*3))
    B[:,:,4] = np.nan
    return np.nanmean(B, axis=2)