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

Question

给定一个二维 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 的每个元素定义的邻域进行排序，并仅返回非屏蔽元素的中值(如果数组是偶数)。

现在，目前我正在非 Python 循环中执行此操作，使用 bottleneck.nanmedian 通过将掩码映射到 NaN。这正是我所需要的，但我希望依靠二维数组操作例程。

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 循环。 IE。我们炸毁了原始数组，以便我们可以一次对所有子数组应用过滤器。这有点类似于 Numpy broadcasting.

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

在我的代码中，我使用了 NaN 进行屏蔽，但是通过一些额外的代码行，您还可以使用 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)