Blockproc like function for Python image processing

Question

edit: it's an image so the suggested (How can I efficiently process a numpy array in blocks similar to Matlab's blkproc (blockproc) function) isn't really working for me

I have the following matlab code

fun = @(block_struct) ...
std2(block_struct.data) * ones(size(block_struct.data));

B=blockproc(im2double(Icorrected), [4 4], fun);

I want to remake my code, but this time in Python. I have installed Scikit and i'm trying to work around it like this

b = np.std(a, axis = 2)

The problem of course it's that i'm not applying the std for a number of blocks, just like above.

How can i do something like this? Start a loop and try to call the function for each X*X blocks? Then i wouldn't keep the size the it was.

Is there another more efficient way?

duplicate of http://stackoverflow.com/questions/5073767/how-can-i-efficiently-process-a-numpy-array-in-blocks-similar-to-matlabs-blkpro ? (Does that help?) — Andy Hayden, Dec 20 '14 at 09:33
@AndyHayden It's something i have checked, but considering the time asked and how libraries like Scikit are introduced,i hoped there was something more clear. — Rentonie, Dec 20 '14 at 12:26
The ```as_strided``` solution @AndyHayden pointed to is your best bet. Here is another using ```as_strided```: http://stackoverflow.com/a/26848312/2823755 — wwii, Dec 20 '14 at 16:26
@wwii My list was already in the format as shown before the B[:, :] part. Although i could not wrap my head around the rest of the solution, so went with my own mediocre thing. — Rentonie, Dec 20 '14 at 16:52

wwii · Answer 1 · 2016-11-10T16:56:53.990

If there is no overlap in the windows you can reshape the data to suit your needs:

Find the mean of 3x3 windows of a 9x9 array.

import numpy as np

>>> a
array([[ 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],
       [36, 37, 38, 39, 40, 41, 42, 43, 44],
       [45, 46, 47, 48, 49, 50, 51, 52, 53],
       [54, 55, 56, 57, 58, 59, 60, 61, 62],
       [63, 64, 65, 66, 67, 68, 69, 70, 71],
       [72, 73, 74, 75, 76, 77, 78, 79, 80]])

Find the new shape

>>> window_size = (3,3)
>>> tuple(np.array(a.shape) / window_size) + window_size
(3, 3, 3, 3)
>>> b = a.reshape(3,3,3,3)

Find the mean along the first and third axes.

>>> b.mean(axis = (1,3))
array([[ 10.,  13.,  16.],
       [ 37.,  40.,  43.],
       [ 64.,  67.,  70.]])
>>>

2x2 windows of a 4x4 array:

>>> a = np.arange(16).reshape((4,4))
>>> window_size = (2,2)
>>> tuple(np.array(a.shape) / window_size) + window_size
(2, 2, 2, 2)
>>> b = a.reshape(2,2,2,2)
>>> b.mean(axis = (1,3))
array([[  2.5,   4.5],
       [ 10.5,  12.5]])
>>>

It won't work if the window size doesn't divide into the array size evenly. In that case you need some overlap in the windows or if you just want overlap numpy.lib.stride_tricks.as_strided is the way to go - a generic N-D function can be found at Efficient Overlapping Windows with Numpy

Another option for 2d arrays is sklearn.feature_extraction.image.extract_patches_2d and for ndarray's - sklearn.feature_extraction.image.extract_patches. Each manipulate the array's strides to produce the patches/windows.

Rentonie · Accepted Answer · 2014-12-20T14:24:44.877

I did the following

io.use_plugin('pil', 'imread')
a = io.imread('C:\Users\Dimitrios\Desktop\polimesa\\arizona.jpg')

B = np.zeros((len(a)/2 +1, len(a[0])/2 +1))


for i in xrange(0, len(a), 2):
    for j in xrange(0, len(a[0]), 2):
        x.append(a[i][j])
        if i+1 < len(a):
            x.append(a[i+1][j])
        if j+1 < len(a[0]):
           x.append(a[i][j+1])
        if i+1 < len(a) and j+1 < len(a[0]):
           x.append(a[i+1][j+1])
        B[i/2][j/2] = np.std(x)
        x[:] = []

and i think it's correct. Iterating over the image by 2 and taking each neighbour node, adding them to a list and calculating std.

edit* later edited for 4x4 blocks.

Sandipan Dey · Answer 3 · 2022-06-18T16:28:12.147

We can implement blockproc() in python the following way:

def blockproc(im, block_sz, func):
    h, w = im.shape
    m, n = block_sz
    for x in range(0, h, m):
        for y in range(0, w, n):
            block = im[x:x+m, y:y+n]
            block[:,:] = func(block)
    return im

Now, let's apply it to implement contrast enhancement with local histogram equalization, with the low-contrast moon image (of size 512x512) as input and choosing 32x32 blocks:

from skimage import data, exposure
img = data.moon()
img = img / img.max()
m, n = 64, 64
img_eq = blockproc(img.copy(), (m, n), exposure.equalize_hist)

Display the input and output images:

Note that the function does in-place modification to the image, hence a copy of the input image is passed instead.

Blockproc like function for Python image processing

3 Answers3