JuliaImages / ImageFiltering.jl

Julia implementations of multidimensional array convolution and nonlinear stencil operations
Other
99 stars 49 forks source link

Optimized Median Filter #222

Open tejus-gupta opened 5 years ago

tejus-gupta commented 5 years ago

I have an implementation of median filter based on A Fast Two-Dimensional Median Filtering Algorithm. The algorithm scales very well with window size. Though the paper desribes a method specifically for 2-D images and median filter, it can easily be extended to work with N-D images and any percentile filter e.g, minima.

Comparison with mapwindow for 256x256 grayscale N0f8 image -

window size time for proposed method (ms) time for mapwindow (ms)
(3,3) 3.81 14.27
(5,5) 4.30 36.02
(7,7) 4.86 71.21
(11,11) 6.07 183.70
(21,21) 9.07 736.39

However, the algorithm only works for image types with discrete pixel values. It also doesn't scale well with increased precision of pixel values as the histogram size increases exponentially.

Comparison with mapwindow for 256x256 grayscale N0f16 image -

window size time for proposed method (ms) time for mapwindow (ms)
(3,3) 97.60 14.20
(5,5) 75.05 36.25
(7,7) 67.52 71.70
(11,11) 59.60 188.72
(21,21) 52.41 743.58

Can we use this median filter in Images.jl specifically for N0f8 images? I would be happy to submit a pull request. My code is available here.

timholy commented 5 years ago

Really nice to see you again, @tejus-gupta!

This would be great to fix. I am remembering https://github.com/JuliaImages/ImageFiltering.jl/pull/34 also. I haven't yet looked at your code, I confess.

I think it's reasonable to have a specialized algorithm for N0f8. We will also need a generic algorithm, and that's part of what has held me back. If I dig into some distant memories, last time I thought about this I seem to remember contemplating "slice accumulation" strategies: create a circular buffer of sorted values from each slice along dimension 1 (the "shift dimension"), and compute the median within the window by clever aggregation among these slices. That's a nice speedup in 2d but gets less exciting in 3d and 4d---if you discount the overhead of the aggregation, it effectively saves you one dimension.

tejus-gupta commented 5 years ago

@timholy I implemented the algorithm described in A coarse-to-fine algorithm for fast median filtering of image data with a huge number of levels by Alparone et al. and got some nice speedup for N0f16. This algorihm is faster than mapwindow on window sizes greater than 3x3.

Comparison with mapwindow for 256x256 grayscale N0f16 image -

window size time for proposed method (ms) time for mapwindow (ms)
(3,3) 30.31 14.05
(5,5) 31.16 36.15
(7,7) 31.48 72.16
(11,11) 34.67 189.71
(21,21) 37.85 754.36

This algorithm uses two-tiers of histograms. The course histogram stores frequency for higher 8 bits and the fine histogram stores frequency for full 16 bits data. Now, we can reach the median faster by first searching for median in course histogram and then using course median to define a small search region in fine histogram. This means that the worst-case number of search steps during median update is 2^8 + 2^8 instead of 2^16.

Right now, I am trying to improve performance further by using more histogram levels. I will send a pull request after this. The code for above algorithm is available here.

I also remember trying to do something similar to what you describe for continous pixel types last year. However, the overhead for maintaining the sorted buffer was quite big for my implementation.

timholy commented 5 years ago

Interesting how close to size-invariant the coarse-to-fine algorithm is. We can easily check the size and use the direct method for sufficiently small sizes.

All that sounds fantastic. Definitely interested in this!

Oblynx commented 5 years ago

@tejus-gupta you mention that the first algorithm can be extended to compute any percentiles. Have you implemented this extension?