Custom median filters

With median filter as an example, this demo shows how you could construct a custom kernel and pass it to mapwindow for common stencil operations.

Unlike other image processing toolboxes, ImageFiltering.jl does not provide you an out-of-box function for median filters. It instead provides the more general mapwindow function, which allows you to apply any function to the window around each pixel.

using ImageFiltering, ImageCore, ImageShow # or you could just `using Images`
using TestImages
using Statistics

img = Float64[isodd(i) - isodd(j) for i = 1:5, j = 1:5]
img[3, 3] = 1000

5×5 Matrix{Float64}:
  0.0  1.0     0.0  1.0   0.0
 -1.0  0.0    -1.0  0.0  -1.0
  0.0  1.0  1000.0  1.0   0.0
 -1.0  0.0    -1.0  0.0  -1.0
  0.0  1.0     0.0  1.0   0.0

patch_size = (3, 3)
imgm = mapwindow(median, img, patch_size)

5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

mapwindow provides a high-level interface to loop over the image and call a function (median in this demo) on each patch/window. It's performing an operation that is nearly equivalent to the loop below, albeit more efficiently:

# For simplicity, border condition is not included here.
imgm = zeros(axes(img))
R = CartesianIndices(img)
I_first, I_last = first(R), last(R)
Δ = CartesianIndex(patch_size .÷ 2)
for I in R
    patch = max(I_first, I-Δ):min(I_last, I+Δ)
    imgm[I] = median(img[patch])
end
imgm

5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

Compared to this hand-written loop, mapwindow offers additional flexibility in handling the borders and allows you to process only a subset of the windows. For some functions (e.g., min, max, and extrema), mapwindow uses a custom implementation that is far more efficient than naively computing the function "freshly" on each window.

When the input array has NaN or some unwanted pixel values, a pre-filtering process is needed to exclude them first. This can be done quite easily by compositing a given kernel operation and a filter operation. Let's still take median filter as an example.

img[2,2] = NaN
imgm = mapwindow(median, img, (3,3))

5×5 Matrix{Float64}:
 NaN    NaN    NaN    0.0  0.0
 NaN    NaN    NaN    0.0  0.0
 NaN    NaN    NaN    0.0  0.0
   0.0    0.0    1.0  0.0  0.0
   0.0    0.0    0.0  0.0  0.0

NaN has polluted the output, which is usually not wanted. This can be fixed quite easily by compositing a new filter kernel. This way, we only compute the median value on the non-NaN subset of each patch/window.

_median(patch) = median(filter(x->!isnan(x), patch))
imgm = mapwindow(_median, img, patch_size)

5×5 Matrix{Float64}:
  0.0   0.0  0.5  0.0  0.0
  0.0   0.0  1.0  0.0  0.0
 -0.5  -0.5  0.0  0.0  0.0
  0.0   0.0  1.0  0.0  0.0
  0.0   0.0  0.0  0.0  0.0

Impulse (Salt and Pepper) noise removal

Median filters are quite robust to outliers (unusual values), this property can be used to remove salt&pepper noise. An image with n-level salt&papper noise is defined as: each pixel p has n/2 probabilty to be filled by 0(the minimal gamut value), n/2 probabilty to be filled by 1(the maximal gamut value), and 1-n probablity to be undamaged.

img = testimage("cameraman")

n = 0.2 # noise level
noisy_img = map(img) do p
    sp = rand()
    if sp < n/2
        eltype(img)(gamutmin(eltype(img))...)
    elseif sp < n
        eltype(img)(gamutmax(eltype(img))...)
    else
        p
    end
end

denoised_img = mapwindow(median!, noisy_img, (5, 5))

mosaicview(img, noisy_img, denoised_img; nrow=1)

This page was generated using DemoCards.jl and Literate.jl.