Code for median filter based on 'A Fast Two-Dimensional Median Filtering Algorithm' by Huang, Yang and Tang.

median_filter.jl

using Images, TestImages, Statistics, BenchmarkTools, Test

img = testimage("lena_gray")
high_resolution_img = convert(Array{Gray{Normed{UInt16,16}}, 2}, img)

function median_filter(img::Array{Gray{Normed{T,f}}, 2}, window_size::Tuple{Int64,Int64}) where {T,f}

    function update_median(median_val, n_pixels_below_median, hist, half_pixels)
        if n_pixels_below_median <= half_pixels
            for val in median_val:2^f
                if n_pixels_below_median + hist[val] > half_pixels
                    median_val = val
                    break
                else
                    n_pixels_below_median += hist[val]
                end
            end
        elseif n_pixels_below_median > half_pixels
            for val in median_val-1:-1:1
                n_pixels_below_median -= hist[val]
                if n_pixels_below_median <= half_pixels
                    median_val = val
                    break
                end
            end
        end
        return median_val, n_pixels_below_median
    end

    result = copy(img)

    if iseven(window_size[1]) || iseven(window_size[2])
        error("window width and height must be odd.")
    end

    half_window_size = window_size.>>1
    img_size = size(img)
    half_pixels = (window_size[1]*window_size[2])>>1

    for j in 1 + half_window_size[2]: img_size[2] - half_window_size[2]
        hist = zeros(Int, 2^f)
        for I in CartesianIndices((1:window_size[1], j-half_window_size[2]:j+half_window_size[2]))
            hist[img[I].val.i+1] += 1
        end

        median_val, n_pixels_below_median = update_median(1, 0, hist, half_pixels)
        result[1+half_window_size[2], j] = Gray{Normed{T,f}}((median_val-1)/(2^f-1))

        for i in 2 + half_window_size[1]: img_size[1] - half_window_size[1]
            for I in CartesianIndices((i-half_window_size[1]-1, j-half_window_size[2]:j+half_window_size[2]))
                hist[img[I].val.i+1] -= 1
                if img[I].val.i + 1 < median_val
                    n_pixels_below_median -= 1
                end
            end

            for I in CartesianIndices((i+half_window_size[1], j-half_window_size[2]:j+half_window_size[2]))
                hist[img[I].val.i+1] += 1
                if img[I].val.i + 1 < median_val
                    n_pixels_below_median += 1
                end
            end
            
            median_val, n_pixels_below_median = update_median(median_val, n_pixels_below_median, hist, half_pixels)
            result[i, j] = Gray{Normed{T,f}}((median_val-1)/(2^f-1))
        end
    end
    result[1+half_window_size[1]:img_size[1]-half_window_size[1], 1+half_window_size[2]:img_size[2]-half_window_size[2]]
end

#Testing

@test median_filter(high_resolution_img, (3,3)) == mapwindow(median, high_resolution_img, (3,3))[2:255, 2:255]
@test median_filter(high_resolution_img, (5, 5)) == mapwindow(median, high_resolution_img, (5,5))[3:254, 3:254]

@test median_filter(img, (3,3)) == mapwindow(median, img, (3,3))[2:255, 2:255]
@test median_filter(img, (5,5)) == mapwindow(median, img, (5,5))[3:254, 3:254]

#Benchmarking

@btime median_filter($img, (3,3));
@btime median_filter($img, (5,5));
@btime median_filter($img, (7,7));
@btime median_filter($img, (11,11));
@btime median_filter($img, (21,21));

@btime mapwindow($median, $img, (3,3));
@btime mapwindow($median, $img, (5,5));
@btime mapwindow($median, $img, (7,7));
@btime mapwindow($median, $img, (11,11));
@btime mapwindow($median, $img, (21,21));

@btime median_filter($high_resolution_img, (3,3));
@btime median_filter($high_resolution_img, (5,5));
@btime median_filter($high_resolution_img, (7,7));
@btime median_filter($high_resolution_img, (11,11));
@btime median_filter($high_resolution_img, (21,21));

@btime mapwindow($median, $high_resolution_img, (3,3));
@btime mapwindow($median, $high_resolution_img, (5,5));
@btime mapwindow($median, $high_resolution_img, (7,7));
@btime mapwindow($median, $high_resolution_img, (11,11));
@btime mapwindow($median, $high_resolution_img, (21,21));