Vectorized extraction of submatrices in numpy-array

Question

My goal is to implement a median-filter, which is a function, that replaces each pixel in a (mostly) 2d-array with the median of its surrounding pixels. It can be used to denoise images.

My implementation extracts submatrices from the original matrix, that contain the pixel itself and its neighbors. This extraction is currently done with a for-loop and as you can imagine that for-loop takes away approximately 95% of execution-time.

Here is my current implementation:

def median_blur(img, fsize=3):
    img_intermediate = np.zeros((img.shape[0] + 2, img.shape[1] + 2), np.uint8)  # First intermediate img, used for padding original image
    img_intermediate[1:img.shape[0]+1, 1:img.shape[1]+1] = img
    img_result = np.empty((*img.shape, fsize, fsize), np.uint8)  # Will contain result, first receives kernel-submatrices

    # Extract submatrices from image
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            img_result[i, j] = img_intermediate[i:i+fsize, j:j+fsize]

    img_result = img_result.reshape(*img.shape, fsize**2)  # Reshape Result-Matrix
    img_result = np.median(img_result, axis=2)  # Calculate median in one go
    return img_result.astype('uint8')

How can i extract these submatrices using a vectorized operation?

As a bonus, if someone with experience in computer-vision is reading this: Is there an better approach of implementing a median-filter other than applying it to an intermediate zero-padded matrix?

Thank you very much.

@edoput it is covered by several libraries (scipy, scikit-image, opencv...), but my goal is to implement it myself in the most performant way. — Tim Hilt
– Tim Hilt, Commented May 2, 2020 at 10:09
Use view-as_windows : stackoverflow.com/questions/42831022 and then use np.median? — Divakar
– Divakar, Commented May 2, 2020 at 10:16
@Divakar thanks for mentioning, i didn't know that function! It wouldn't be an only-numpy-solution but certainly would do what i want. — Tim Hilt
– Tim Hilt, Commented May 2, 2020 at 10:25
view_as_windows is a standalone function codebase. So, you can just pick it up from source. — Divakar
– Divakar, Commented May 2, 2020 at 10:27

Ehsan · Accepted Answer · 2020-05-02 10:45:06Z

3

Here is a vectorized solution. However, you can come up with a faster solution by paying attention to memory order of the image array:

from numpy.lib.stride_tricks import as_strided

img_padded = np.pad(img, 1, mode='constant')
sub_shape = (fsize, fsize)
view_shape = tuple(np.subtract(img_padded.shape, sub_shape) + 1) + sub_shape
sub_img = as_strided(img_padded, view_shape, img_padded.strides * 2)
sub_img = sub_img.reshape(img.shape + (fsize**2,))
result = np.median(sub_img, axis=2).astype(int)

Use pad function to pad first and then as_strided to get the sub_matrices and finally apply median to your strides.

UPDATE: using view_as_windows suggested by @Divakar in comments:

from skimage.util.shape import view_as_windows

img_padded = np.pad(img, 1, mode='constant')
sub_shape = (fsize, fsize)
view_shape = view_as_windows(img_padded, sub_shape, 1)
sub_img = view_shape.reshape(img.shape + (fsize**2,))
result = np.median(sub_img, axis=2).astype(int)

view_as_windows also provides sub-matrices similar to strides.

sample image and output:

img: 
[[ 1  2  3  4  5  6]
 [ 7  8  9 10 11 12]
 [13 14 15 16 17 18]
 [19 20 21 22 23 24]]

median_filtered: 
[[ 0  2  3  4  5  0]
 [ 2  8  9 10 11  6]
 [ 8 14 15 16 17 12]
 [ 0 14 15 16 17  0]]

edited May 2, 2020 at 10:45

answered May 2, 2020 at 10:37

Ehsan

12.5k2 gold badges24 silver badges36 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Tim Hilt Over a year ago

Beautiful! This made my existing solution approx. 7x faster! Can you elaborate on the possible memory-order gains?

Ehsan Over a year ago

Glad it helped. There might be a faster way of implementing this using numba or strided loops if you know your datatype. You would probably benefit from it if your array is saved in a C-contiguous manner (meaning rows first). view_as_windows will create a new array if your array is not C-contiguous. I would expect as_strides to do the same, but I am not sure. Generally speaking anything faster than this probably needs to get into C-contiguous/F-contiguous memory allocation and knowing your datatypes.

Collectives™ on Stack Overflow

Vectorized extraction of submatrices in numpy-array

1 Answer 1

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related