以下是使用块的另一种(循环自由)方式的一些示例:
import numpy as np
from numpy.lib.stride_tricks import as_strided as ast
A= np.arange(36).reshape(6, 6)
print A
#[[ 0 1 2 3 4 5]
# [ 6 7 8 9 10 11]
# ...
# [30 31 32 33 34 35]]
# 2x2 block view
B= ast(A, shape= (3, 3, 2, 2), strides= (48, 8, 24, 4))
print B[1, 1]
#[[14 15]
# [20 21]]
# for preserving original shape
B[:, :]= np.dot(B[:, :], np.array([[0, 1], [1, 0]]))
print A
#[[ 1 0 3 2 5 4]
# [ 7 6 9 8 11 10]
# ...
# [31 30 33 32 35 34]]
print B[1, 1]
#[[15 14]
# [21 20]]
# for reducing shape, processing in 3D is enough
C= B.reshape(3, 3, -1)
print C.sum(-1)
#[[ 14 22 30]
# [ 62 70 78]
# [110 118 126]]
所以只是试图简单地将matlab功能复制到numpy并不是所有的方法进行最好的方法。有时需要“脱帽”思维。
警告:
一般来说,基于步幅技巧的实现可能(但不是必须)遭受一些性能惩罚。所以做好准备,以各种方式衡量你的表现。无论如何,首先检查所需的功能(或类似的功能,以便轻松适应)已经准备好以numpy或scipy的方式实现了。
更新:
请注意,这里并没有任何真正的魔力,因此我将提供一个简单的函数来获取任何适合的2D numpy数组的block_view。所以这里我们去:
from numpy.lib.stride_tricks import as_strided as ast
def block_view(A, block= (3, 3)):
"""Provide a 2D block view to 2D array. No error checking made.
Therefore meaningful (as implemented) only for blocks strictly
compatible with the shape of A."""
# simple shape and strides computations may seem at first strange
# unless one is able to recognize the 'tuple additions' involved ;-)
shape= (A.shape[0]/ block[0], A.shape[1]/ block[1])+ block
strides= (block[0]* A.strides[0], block[1]* A.strides[1])+ A.strides
return ast(A, shape= shape, strides= strides)
if __name__ == '__main__':
from numpy import arange
A= arange(144).reshape(12, 12)
print block_view(A)[0, 0]
#[[ 0 1 2]
# [12 13 14]
# [24 25 26]]
print block_view(A, (2, 6))[0, 0]
#[[ 0 1 2 3 4 5]
# [12 13 14 15 16 17]]
print block_view(A, (3, 12))[0, 0]
#[[ 0 1 2 3 4 5 6 7 8 9 10 11]
# [12 13 14 15 16 17 18 19 20 21 22 23]
# [24 25 26 27 28 29 30 31 32 33 34 35]]