How can I efficiently process a numpy array in blocks similar to Matlab's blkproc (blockproc) function

I'm looking for a good approach for efficiently dividing an image into small regions, processing each region separately, and then re-assembling the results from each process into a single processed image. Matlab had a tool for this called blkproc (replaced by blockproc in newer versions of Matlab).

In an ideal world, the function or class would support overlap between the divisions in the input matrix too. In the Matlab help, blkproc is defined as:

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

• A is your input matrix,
• [m n] is the block size
• [mborder, nborder] is the size of your border region (optional)
• fun is a function to apply to each block

I have kluged together an approach, but it strikes me as clumsy and I bet there's a much better way. At the risk of my own embarrassment, here's my code:

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)