Code for median filter based on 'A coarse-to-fine algorithm for fast median filtering of image data with a huge number of levels' by Alparone et al.

median_filter_multilevel.jl

using Images, TestImages, Statistics, BenchmarkTools, Test

img = testimage("lena_gray")

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

    nlevels = 2^f
    course_nlevels = 2^(f>>1)
    bin_bits = f>>1

    function update_median(median_val, lt_median, course_median_val, lt_course_median, hist, course_hist, half_pixels)
        if lt_median <= half_pixels
            for val in course_median_val:course_nlevels
                if lt_course_median + course_hist[val] > half_pixels
                    course_median_val = val
                    break
                else
                    lt_course_median += course_hist[val]
                end
            end
            
            if (median_val - 1) >> bin_bits + 1 != course_median_val
                median_val = (course_median_val - 1) << bin_bits + 1
                lt_median = lt_course_median
            end
            
            for val in median_val:nlevels
                if lt_median + hist[val] > half_pixels
                    median_val = val
                    break
                else
                    lt_median += hist[val]
                end
            end
            
        elseif lt_median > half_pixels
            if lt_course_median <= half_pixels
                for val in median_val-1:-1:1
                    lt_median -= hist[val]
                    if lt_median <= half_pixels
                        median_val = val
                        break
                    end
                end
            else
                for val in course_median_val-1:-1:1
                    lt_course_median -= course_hist[val]
                    if lt_course_median <= half_pixels
                        course_median_val = val
                        break
                    end
                end
                median_val = (course_median_val - 1) << bin_bits + 1
                lt_median = lt_course_median
                for val in median_val:nlevels
                    if lt_median + hist[val] > half_pixels
                        median_val = val
                        break
                    else
                        lt_median += hist[val]
                    end
                end
            end
        end
        return median_val, lt_median, course_median_val, lt_course_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, nlevels)
        course_hist = zeros(Int, course_nlevels+1)

        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
            course_hist[(img[I].val.i)>>bin_bits+1] += 1
        end

        median_val, lt_median, course_median_val, lt_course_median = update_median(1, 0, 1, 0, hist, course_hist, half_pixels)
        result[1+half_window_size[2], j] = Gray{Normed{T,f}}((median_val-1)/(nlevels-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
                course_hist[(img[I].val.i)>>bin_bits+1] -= 1
                if img[I].val.i + 1 < median_val
                    lt_median -= 1
                end

                if (img[I].val.i)>>bin_bits+1 < course_median_val
                    lt_course_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
                course_hist[(img[I].val.i)>>bin_bits+1] += 1
                if img[I].val.i + 1 < median_val
                    lt_median += 1
                end

                if (img[I].val.i)>>bin_bits+1 < course_median_val
                    lt_course_median += 1
                end
            end
            
            median_val, lt_median, course_median_val, lt_course_median = update_median(median_val, lt_median, course_median_val, lt_course_median, hist, course_hist, half_pixels)
            result[i, j] = Gray{Normed{T,f}}((median_val-1)/(nlevels-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_multilevel(img, (3,3)) == mapwindow(median, img, (3,3))[2:255, 2:255]
@test median_filter_multilevel(img, (5,5)) == mapwindow(median, img, (5,5))[3:254, 3:254]

#Benchmarking

@btime median_filter_multilevel($img, (3,3));
@btime median_filter_multilevel($img, (5,5));
@btime median_filter_multilevel($img, (7,7));
@btime median_filter_multilevel($img, (11,11));
@btime median_filter_multilevel($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));

high_resolution_img = convert(Array{Gray{Normed{UInt32,32}}, 2}, img)

@btime median_filter_multilevel($high_resolution_img, (3,3));
@btime median_filter_multilevel($high_resolution_img, (5,5));
@btime median_filter_multilevel($high_resolution_img, (7,7));
@btime median_filter_multilevel($high_resolution_img, (11,11));
@btime median_filter_multilevel($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));