我怎样才能有效地处理类似于Matlab的blkproc(blockproc)函数块的numpy数组

我正在寻找一种有效的方法,将图像有效地分割成小区域,分别处理每个区域,然后将每个进程的结果重新组合为单个处理后的图像。 Matlab有一个叫做blkproc的工具(在Matlab的新版本中由blockproc取代)。

在一个理想的世界中,函数或类也会支持inputmatrix中各个分区之间的重叠。 在Matlab帮助中,blkproc被定义为:

B = blkproc(A,[mn],[mborder nborder],好玩,…)

  • A是你的inputmatrix,
  • [mn]是块大小
  • [mborder,nborder]是你的边界区域的大小(可选)
  • 乐趣是适用于每个块的function

我曾经在一起讨论过一个方法,但是它让我觉得很笨拙,我敢打赌还有更好的办法。 冒着我自己的尴尬风险,这是我的代码:

 import numpy as np def segmented_process(M, blk_size=(16,16), overlap=(0,0), fun=None): rows = [] for i in range(0, M.shape[0], blk_size[0]): cols = [] for j in range(0, M.shape[1], blk_size[1]): cols.append(fun(M[i:i+blk_size[0], j:j+blk_size[1]])) rows.append(np.concatenate(cols, axis=1)) return np.concatenate(rows, axis=0) R = np.random.rand(128,128) passthrough = lambda(x):x Rprime = segmented_process(R, blk_size=(16,16), overlap=(0,0), fun=passthrough) np.all(R==Rprime) 

6 Solutions collect form web for “我怎样才能有效地处理类似于Matlab的blkproc(blockproc)函数块的numpy数组”

这里有一些不同的(无循环)方式来处理块的一些例子:

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

所以只是简单的将matlabfunction复制到numpy并不是所有的方法都是最好的方法。 有时需要“脱帽”的想法。

警告
一般来说,基于步幅技巧的实现可能 (但不一定需要)遭受一些性能损失。 所以准备好各种方式来衡量你的performance。 在任何情况下,首先检查所需function(或足够类似的,以便于适应)已经准备好在numpyscipy实现是scipy

更新
请注意,这里没有真正的magic ,所以我会提供一个简单的函数来获得任何合适的二维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]] 

按切片/视图进行处理。 连接是非常昂贵的。

 for x in xrange(0, 160, 16): for y in xrange(0, 160, 16): view = A[x:x+16, y:y+16] view[:,:] = fun(view) 

我采取了两个input,以及我原来的方法,并比较结果。 正如@eat正确指出的那样,结果取决于input数据的性质。 令人惊讶的是,在less数情况下连接节拍视图处理。 每种方法都有一个甜蜜点。 这是我的基准代码:

 import numpy as np from itertools import product def segment_and_concatenate(M, fun=None, blk_size=(16,16), overlap=(0,0)): # truncate M to a multiple of blk_size M = M[:M.shape[0]-M.shape[0]%blk_size[0], :M.shape[1]-M.shape[1]%blk_size[1]] rows = [] for i in range(0, M.shape[0], blk_size[0]): cols = [] for j in range(0, M.shape[1], blk_size[1]): max_ndx = (min(i+blk_size[0], M.shape[0]), min(j+blk_size[1], M.shape[1])) cols.append(fun(M[i:max_ndx[0], j:max_ndx[1]])) rows.append(np.concatenate(cols, axis=1)) return np.concatenate(rows, axis=0) from numpy.lib.stride_tricks import as_strided 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 as_strided(A, shape= shape, strides= strides) def segmented_stride(M, fun, blk_size=(3,3), overlap=(0,0)): # This is some complex function of blk_size and M.shape stride = blk_size output = np.zeros(M.shape) B = block_view(M, block=blk_size) O = block_view(output, block=blk_size) for b,o in zip(B, O): o[:,:] = fun(b); return output def view_process(M, fun=None, blk_size=(16,16), overlap=None): # truncate M to a multiple of blk_size from itertools import product output = np.zeros(M.shape) dz = np.asarray(blk_size) shape = M.shape - (np.mod(np.asarray(M.shape), blk_size)) for indices in product(*[range(0, stop, step) for stop,step in zip(shape, blk_size)]): # Don't overrun the end of the array. #max_ndx = np.min((np.asarray(indices) + dz, M.shape), axis=0) #slices = [slice(s, s + f, None) for s,f in zip(indices, dz)] output[indices[0]:indices[0]+dz[0], indices[1]:indices[1]+dz[1]][:,:] = fun(M[indices[0]:indices[0]+dz[0], indices[1]:indices[1]+dz[1]]) return output if __name__ == "__main__": R = np.random.rand(128,128) squareit = lambda(x):x*2 from timeit import timeit t ={} kn = np.array(list(product((8,16,64,128), (128, 512, 2048, 4096)) ) ) methods = ("segment_and_concatenate", "view_process", "segmented_stride") t = np.zeros((kn.shape[0], len(methods))) for i, (k, N) in enumerate(kn): for j, method in enumerate(methods): t[i,j] = timeit("""Rprime = %s(R, blk_size=(%d,%d), overlap = (0,0), fun = squareit)""" % (method, k, k), setup=""" from segmented_processing import %s import numpy as np R = np.random.rand(%d,%d) squareit = lambda(x):x**2""" % (method, N, N), number=5 ) print "k =", k, "N =", N #, "time:", t[i] print (" Speed up (view vs. concat, stride vs. concat): %0.4f, %0.4f" % ( t[i][0]/t[i][1], t[i][0]/t[i][2])) 

结果如下:

numpy中以块为单位处理大型矩阵的三种方法的比较 请注意,对于小块大小,分段步幅方法赢得3-4倍。 只有在较大的块大小(128 x 128)和非常大的matrix(2048 x 2048或更大)时,视图处理方法才会赢,然后只有一小部分。 基于烘烤,它看起来像@eat得到的复选标记! 感谢你们俩的好例子!

甚至在游戏的后期。 有一个叫做Bob的瑞士image processing软件包可以在以下urlfind: https ://www.idiap.ch/software/bob/docs/releases/last/sphinx/html/index.html它有一个python命令bob.ip.block在: https : //www.idiap.ch/software/bob/docs/releases/last/sphinx/html/ip/generated/bob.ip.block.html#bob.ip.block似乎做的一切Matlab命令“blockproc”呢。 我没有testing过。

还有一个有趣的命令bob.ip.DCTFeatures,它包含了“块”function来提取或修改图像的DCT系数。

比赛迟到了,但这会做重叠的块。 我没有在这里做,但你可以很容易地适应步骤大小的窗口,我想:

 from numpy.lib.stride_tricks import as_strided def rolling_block(A, block=(3, 3)): shape = (A.shape[0] - block[0] + 1, A.shape[1] - block[1] + 1) + block strides = (A.strides[0], A.strides[1]) + A.strides return as_strided(A, shape=shape, strides=strides) 

我发现这个教程 – 最终的源代码提供了完全所需的function! 它甚至应该为任何维度工作(我没有testing它) http://www.johnvinyard.com/blog/?p=268

尽pipe源代码最后的“扁平”选项似乎有点小问题。 尽pipe如此,非常好的软件!

  • 如何规范一个NumPy数组在一定的范围内?
  • numpy,scipy,matplotlib和pylab之间的混淆
  • 我如何使用scipy进行二维插值?
  • 我应该使用scipy.pi,numpy.pi还是math.pi?
  • 如何在matplotlib中创build密度图?
  • Python的networking托pipe:Numpy,Matplotlib,科学计算
  • 无法通过pip安装Scipy
  • Python / SciPy的峰值searchalgorithm
  • 计算两个multidimensional array之间的相关系数
  • 直方图Matplotlib
  • 浓缩距离matrix如何工作? (pdist)