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

Question

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

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

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

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

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

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)

Answer 1

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

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 中实现是明智的。

更新:
请注意，这里 strides 并没有真正涉及到真正的 magic，所以我将提供一个简单的函数来获取任何合适的 block_view二维 numpy-数组。所以我们开始:

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