Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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.
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));
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment