如何有效地处理类似于 Matlab 的 blkproc (blockproc) 函数的块中的 numpy 数组

Posted

技术标签:

【中文标题】如何有效地处理类似于 Matlab 的 blkproc (blockproc) 函数的块中的 numpy 数组【英文标题】:How can I efficiently process a numpy array in blocks similar to Matlab's blkproc (blockproc) function 【发布时间】:2011-07-01 17:17:00 【问题描述】:

我正在寻找一种很好的方法来有效地将图像划分为小区域,分别处理每个区域,然后将每个处理的结果重新组合成单​​个处理后的图像。 Matlab 有一个名为 blkproc 的工具(在较新版本的 Matlab 中由 blockproc 代替)。

在理想情况下,函数或类也将支持输入矩阵中的划分之间的重叠。在 Matlab 帮助中,blkproc 定义为:

B = blkproc(A,[m n],[mborder nborder],fun,...)

A 是你的输入矩阵, [m n] 是块大小 [mborder, nborder] 是边框区域的大小(可选) fun 是应用于每个块的函数

我已经拼凑了一种方法,但它让我觉得很笨拙,我敢打赌会有更好的方法。冒着让我尴尬的风险,这是我的代码:


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)

【问题讨论】:

【参考方案1】:

按切片/视图处理。串联非常昂贵。

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)

【讨论】:

这肯定有助于解决我的性能问题。我将不得不稍微调整它以适应一般大小的矩阵。另外,我怀疑我可以使用一些 itertools 使其成为单个 for 循环。那么也许它会支持并行化。 这种方法绝对是最容易理解的并且提供了实质性的改进(2-4 倍),但是它被 shape 和 stride 方法所迷惑。另一种方法对小图像区域展示了一个数量级的改进!【参考方案2】:

以下是使用块的不同(无循环)方式的一些示例:

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 并不是最好的方法。有时需要“现成”的想法。

警告: 一般来说,基于步幅技巧的实现可能(但不一定需要)遭受一些性能损失。所以准备好用各种方式衡量你的表现。无论如何,明智的做法是首先检查所需的功能(或足够类似的功能,以便轻松适应)是否已在numpyscipy 中实现。

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

【讨论】:

这是一个巧妙的技巧,我学到了一些关于 numpy internal structure 的知识,但它看起来很难推广到任意大小的矩阵。 @Carl F.:arbitrary sized matrices 是什么意思?如果您正在使用 2D 数组,并且需要将它们划分为任何兼容的块,那么 stridesshapes 的计算是微不足道的。你可能会问更多与这方面相关的问题。无论如何,我建议你给stride tricks 一个机会,从长远来看它可能会有所回报。谢谢 我想要一个可以将任何块大小作为输入的函数,但我不知道如何计算步幅值。你能提供一种基于块大小计算它们的算法吗?计算步幅(16,16)与计算步幅(2,2)?我通读了文档,但发现它们很难理解。老实说,这些天我更接近物理而不是计算机科学。 @Zoran Pavlovic:请详细说明“如何将各个块连接回单个数组”的意思!请注意,“块”只是原始数组的视图,因此无需加入任何块,单个数组只是原始数组;-) 再看这个例子,请注意步幅是为 'int32' 类型计算的。对于“int64”,它将不起作用【参考方案3】:

我采用了两种输入以及我的原始方法并比较了结果。正如@eat 正确指出的那样,结果取决于输入数据的性质。令人惊讶的是,在少数情况下,concatenate 会击败视图处理。每种方法都有一个甜蜜点。这是我的基准代码:

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

结果如下:

请注意,对于小块大小,分段步幅方法的优势是 3-4 倍。只有在大块大小(128 x 128)和非常大的矩阵(2048 x 2048 和更大)下,视图处理方法才能获胜,然后仅占一小部分。根据烘烤结果,看起来@eat 得到了复选标记!感谢你们两位的好例子!

【讨论】:

不错的基准!我认为它很好地展示了尝试制作通用代码时将面临的一些困难。实际上,生成始终以最佳方式执行的代码是非常重要的问题。谢谢【参考方案4】:

游戏有点晚了,但这会造成重叠块。我还没有在这里做过,但我认为你可以轻松地适应移动窗口的步长:

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)

【讨论】:

【参考方案5】:

甚至在游戏后期。 有一个名为 Bob 的瑞士图像处理包,可在以下位置获得: https://www.idiap.ch/software/bob/ 它有一些用于块的 python 命令,例如bob.ip.base.block 这似乎完成了 Matlab 命令“blockproc”所做的一切。 我没有测试过。

还有一些有趣的命令 bob.ip.base.DCTFeatures 结合了“块”功能来提取或修改图像的 DCT 系数。

【讨论】:

看起来很棒,但是链接坏了! @RobBW【参考方案6】:

我找到了本教程 - 最终的源代码提供了所需的功能! 它甚至应该适用于任何维度(我没有测试它) http://www.johnvinyard.com/blog/?p=268

虽然源代码末尾的“扁平化”选项似乎有点错误。不过,非常不错的软件!

【讨论】:

re: flatten,如果在 Python 3 中应该是dim = list(filter(lambda i : i != 1,dim)) 死链接,对于任何仍然好奇的人:web.archive.org/web/20161118125327/http://www.johnvinyard.com/…

以上是关于如何有效地处理类似于 Matlab 的 blkproc (blockproc) 函数的块中的 numpy 数组的主要内容,如果未能解决你的问题,请参考以下文章

图像提取基于matlab鱼眼图像有效区域提取含Matlab源码 2186期

信号处理高斯滤波器matlab源码

如何在 Python 中读取类似于 Matlab audioread 的音频文件?

如何在 for 循环中添加不同大小的 numpy 数组条目(类似于 Matlab 的单元格数组)?

在Matlab中有效地计算成对平方欧几里德距离

如何在 Fortran 中有效地计算矩阵内积?