median downsampling in python

Question

I would like to downsample an image using a median filter. The block_reduce method in skimage together with func=numpy.median seems to be appropriate at first sight. However, I got the impression that block_reduce applies the func one axis at a time, while I would expect the median filter to be applied to the whole block at once - for nontrivial inputs the outcome is not the same.

A minimal demonstration example would be this:

from skimage.measure import block_reduce
import numpy as np

image = np.array([[4, 6, 6, 2],
 [6, 7, 9, 9],
 [3, 0, 9, 0],
 [0, 6, 6, 4]])

expected = np.array([[np.median(image[0:2,0:2]), np.median(image[0:2,2:4])],
            [np.median(image[2:4,0:2]), np.median(image[2:4,2:4])]])

actual = block_reduce(image, (2,2), func=np.median)

assert np.array_equal(expected, actual)

The last assertion fails because the median filter is not applied in the way I expected.

Can you suggest other means of downsampling an image using a median filter?

Answer 1

Here is some hack which solves your specific issue.

def clever_func(block, axis):
    # axis unused on purpose
    if len(block.shape) == 4:
        return np.median(block, axis=[2, 3])
    else:
        return block

actual = block_reduce(image, (2,2), func=clever_func)

Essentially, the argument func to block_reduce will be called as func(out, axis=-1) in the code (source). Instead of calling numpy.median, I have hijacked this, so that median gets called with the right axis argument.

In your example, block.shape will be (2, 2, 2, 2) in the first call, and then (2, 2, 2) on the second call of clever_func. I am only using np.median on the first call (if len(block.shape) == 4), and nothing on the second call.

This is not a pretty solution.