你正在操作原始数组的滑动窗口。关于滑动窗口、numpy和python,SO上有很多问题和答案。通过操纵数组的步幅,这个过程可以大大加快。这是一个通用函数,它将返回具有或不具有重叠的数组的(x,y)窗口。使用这个步幅技巧似乎只比@mskimm的解决方案慢一点点。这是一个很好的工具,值得拥有。这个函数不是我的 - 它是在
Efficient Overlapping Windows with Numpy中发现的。
import numpy as np
from numpy.lib.stride_tricks import as_strided as ast
from itertools import product
def norm_shape(shape):
'''
Normalize numpy array shapes so they're always expressed as a tuple,
even for one-dimensional shapes.
Parameters
shape - an int, or a tuple of ints
Returns
a shape tuple
from http://www.johnvinyard.com/blog/?p=268
'''
try:
i = int(shape)
return (i,)
except TypeError:
pass
try:
t = tuple(shape)
return t
except TypeError:
pass
raise TypeError('shape must be an int, or a tuple of ints')
def sliding_window(a,ws,ss = None,flatten = True):
'''
Return a sliding window over a in any number of dimensions
Parameters:
a - an n-dimensional numpy array
ws - an int (a is 1D) or tuple (a is 2D or greater) representing the size
of each dimension of the window
ss - an int (a is 1D) or tuple (a is 2D or greater) representing the
amount to slide the window in each dimension. If not specified, it
defaults to ws.
flatten - if True, all slices are flattened, otherwise, there is an
extra dimension for each dimension of the input.
Returns
an array containing each n-dimensional window from a
from http://www.johnvinyard.com/blog/?p=268
'''
if None is ss:
ss = ws
ws = norm_shape(ws)
ss = norm_shape(ss)
ws = np.array(ws)
ss = np.array(ss)
shape = np.array(a.shape)
ls = [len(shape),len(ws),len(ss)]
if 1 != len(set(ls)):
error_string = 'a.shape, ws and ss must all have the same length. They were{}'
raise ValueError(error_string.format(str(ls)))
if np.any(ws > shape):
error_string = 'ws cannot be larger than a in any dimension. a.shape was {} and ws was {}'
raise ValueError(error_string.format(str(a.shape),str(ws)))
newshape = norm_shape(((shape - ws) // ss) + 1)
newshape += norm_shape(ws)
newstrides = norm_shape(np.array(a.strides) * ss) + a.strides
strided = ast(a,shape = newshape,strides = newstrides)
if not flatten:
return strided
meat = len(ws) if ws.shape else 0
firstdim = (np.product(newshape[:-meat]),) if ws.shape else ()
dim = firstdim + (newshape[-meat:])
dim = filter(lambda i : i != 1,dim)
return strided.reshape(dim)
使用方法:
b = sliding_window(arr, (2,2), flatten = False)
c = b.sum((1,2))
使用numpy.einsum
可以获得约24%的性能提升。
c = np.einsum('ijkl -> ij', b)
一个SO问答示例 如何高效地像Matlab的blkproc(blockproc)函数一样分块处理numpy数组,选择的答案会对你有帮助。
np.einsum('ijkl -> ik', a.reshape(2160, 2, 4320, 2))
获得更快的速度。 - wwiinp.einsum
的优势通常来自于求和和乘积操作的向量化(SIMD)。我认为最新版本的numpy在np.sum
中也具有相同的特性。 - Jaime