Median filter for image Python3

Question

I wanted to implement a radial median filter. I have the following picture (size = (Nx,Ny))

I want to derive radius for each pixels. For each radius compute median value and put it to a new matrix in the place of all pixels with the same radius. I found Image Smoothing Using Median Filter, but it isn't fast enough. And I created my own script, unfortunately, it isn't fast too. I tested it on some generatic data:

import cv2
from PIL import Image
from scipy import stats, ndimage, misc
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.image as mpimg
from scipy import stats


a = np.array([[0.,0.,0.,0.,0.],[0.,5.,1.,9.,0.],[0.,10.,2.,10.,0.],[0.,9.,1.,5.,0.],[0.,0.,0.,0.,0.]])

b = a.copy().flatten()

y,x = np.indices((a.shape))
center = [len(x)//2, len(y)//2]
r = np.hypot(x-center[0],y-center[1])

r = r.astype(np.int) # integer part of radii (bin size = 1)

set_r = set(r.flatten()) # get the list of r without duplication
max_r = max(set_r) # determine the maximum r

median_r = np.array([0.]*len(r.flatten())) # array of median I for each r


for j in set_r:
    result = np.where(r.flatten() == j) 
    median_r[result[0]] = np.median(b[result[0]])



a_med = median_r.reshape(a.shape)

am_med = ndimage.median_filter(a, 3)

plt.figure(figsize=(16, 5))

plt.subplot(141)
plt.imshow(a, interpolation='nearest')
plt.axis('off')
plt.title('Original image', fontsize=20)
plt.subplot(142)
plt.imshow(am_med, interpolation='nearest', vmin=0, vmax=5)
plt.axis('off')
plt.title('Median filter', fontsize=20)
plt.subplot(143)
plt.imshow(a_med, interpolation='nearest')
plt.axis('off')
plt.title('Own median', fontsize=20)


plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0, left=0,
                    right=1)

plt.show()

I'd like to find a convenient way for solving this issue

user_2836621 · Answer

I think you want to replace all pixels around the radius of each circle of the image with the mean of the pixels on that same radius in the input image.

I propose to warp the image to cartesian coordinates, calculate the mean and then warp back to polar coordinates.

I generated some test data of a decent size like this:

#!/usr/bin/env python3

import cv2
from PIL import Image
from scipy import stats, ndimage, misc
import matplotlib.image as mpimg
from scipy import stats
import numpy as np

w, h = 600, 600
a = np.zeros((h,w),np.uint8)

# Generate some arcs
for s in range(1,6):
    radius = int(s*w/14)
    centre = (int(w/2), int(w/2))
    axes = (radius, radius)
    angle = 360
    startAngle = 0
    endAngle = 72*s

    cv2.ellipse(a, centre, axes, angle, startAngle, endAngle, 255, 2)

That gives this:

Image.fromarray(a.astype(np.uint8)).save('start.png')

def orig(a):
    b = a.copy().flatten()
    y,x = np.indices((a.shape))
    center = [len(x)//2, len(y)//2]
    r = np.hypot(x-center[0],y-center[1])
    r = r.astype(np.int) # integer part of radii (bin size = 1)
    set_r = set(r.flatten()) # get the list of r without duplication
    max_r = max(set_r) # determine the maximum r
    median_r = np.array([0.]*len(r.flatten())) # array of median I for each r
    for j in set_r:
        result = np.where(r.flatten() == j) 
        median_r[result[0]] = np.median(b[result[0]])
    return median_r

def me(a):
    h, w = a.shape
    centre = (int(h/2), int(w/2))
    maxRad = np.sqrt(((h/2.0)**2.0)+((w/2.0)**2.0))
    pol = cv2.warpPolar(a.astype(np.float), a.shape, centre, maxRad, flags=cv2.WARP_POLAR_LINEAR+cv2.WARP_FILL_OUTLIERS)
    polmed = np.median(pol,axis=0,keepdims=True)
    polmed = np.broadcast_to(polmed,a.shape)
    res = cv2.warpPolar(polmed, a.shape, centre,  maxRad, cv2.WARP_INVERSE_MAP)
    return res.astype(np.uint8)

a_med = orig(a).reshape(a.shape)

Image.fromarray(a_med.astype(np.uint8)).save('result.png')

r = me(a)
Image.fromarray(r).save('result-me.png')

The result is the same as yours, i.e. it removes all arcs less than 180 degrees and fills all arcs over 180 degrees:

But the timing for mine is 10x faster:

In [58]: %timeit a_med = orig(a).reshape(a.shape)                                                                               
287 ms ± 17.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [59]: %timeit r = me(a)                                                                                                      
29.9 ms ± 107 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In case you are having difficulty imagining what I get after warpPolar(), it looks like this. Then I use np.mean() to take the mean down the columns, i.e. axis=0:

Keywords: Python, radial mean, radial median, cartesian coordinates, polar coordinates, rectangular, warpPolar, linearPolar, OpenCV, image, image processing

user_2221270 · Answer

Most of the answers here seem to center on performance optimizations of the naive median filtering algorithm. It's worth noting that the median filters you would find in imaging packages like OpenCV/scikit-image/MATLAB/etc. implement faster algorithms.

http://nomis80.org/ctmf.pdf

If you are median filtering uint8 data, there are a lot of clever tricks to be played with reusing histograms as you move from neighborhood to neighborhood.

I would use the median filter in an imaging package rather than trying to roll one yourself if you care about speed.

Median filter for image Python3

2 Answers