Skip to content

Instantly share code, notes, and snippets.

@tejus-gupta
Last active October 4, 2018 13:33
Show Gist options
  • Save tejus-gupta/2f29e101e18a58dcf2477991151b3a42 to your computer and use it in GitHub Desktop.
Save tejus-gupta/2f29e101e18a58dcf2477991151b3a42 to your computer and use it in GitHub Desktop.
Code for median filter based on 'A Fast Two-Dimensional Median Filtering Algorithm' by Huang, Yang and Tang.
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));
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment