Adds optimised algorithm for Median Filter

[hari-sikchi]

Fast Median filter based on algorithm as explained in paper: http://ieeexplore.ieee.org/document/1163188/?reload=true

Algorithm in short: Algorithm : Huang’s O(n) median filtering algorithm.
Input: Image X of size m × n, kernel radius r
Output: Image Y of the same size as X
Initialize kernel histogram H for i = 1 to m do
for j = 1 to n do
for k = −r to r do
Remove X[i+k,j−r−1] from H
Add X[i+k,j+r] to H
end for
Y[i,j] ← median(H)
end for
end for

Improves the time complexity from O((ImageDims) x (kernelSize x kernelSize x log(kernelSize))) to constant(k) x O((ImageDims) x (kernelSize)

Currently, the filter is running slow for small filters than the original one and I have observed for 512x512 images, the new filter runs faster for kernels greater than 7x7 size. I think the code if written in a more optimized way can help make it run faster also for smaller kernels.

I am fairly new to Julia and would be grateful for help to make it faster.

[timholy]

Cool! I'd be excited to have a more efficient algorithm for large kernels.

Before I do a detailed review, maybe a couple of high-level comments. First, it looks like you've changed the generic implementation to use a histogram- and mode-based approach. This is clearly useful for median-filtering, but I doubt this will be helpful for, e.g, mapwindow(std, img, window). I'm wondering if it makes sense to introduce a trait-function

uses_histogram(::typeof(median!)) = true
uses_histogram(::Any) = false

and then branch to a different implementation whenever uses_histogram(f) evaluates to true. (It will be important to make sure that out is allocated as having the same type for either algorithm.)

Second, it sure would be nice to have this work in more than 2d, but I think this implementation is quite specific to 2d, right? If so we should also limit this to just that case (or generalize it to higher dimensions).

Third, perhaps @tejus-gupta might chime in? See https://github.com/JuliaImages/Images.jl/pull/639#issuecomment-308203538

[hari-sikchi]

Having a trait function would be definitely better. I'll make those changes and update the pull request.

Yes, for now it is specific to 2D. I'll see if I can generalise it to nD.

[timholy]

I'll see if I can generalise it to nD.

If you haven't read this blog post, it might help. Seems that if you just use an n-1-dimensional CartesianIndex for the trailing dimensions, it might do much of it for you?

[hari-sikchi]

I have updated the commit with support for N-Dimensional inputs also added the uses_histogram function to maintain the generic implementation for mapwindow. I had some doubts about the code:-

1. Can I use CartesianRange to traverse backward.
2. Allowing this algorithm to take n-Dimensional input is making it run much slower(maybe some problem in code) for 2-D input than it did before. Should algorithms be different for different dimensions?
3. This algorithm uses bucket sort for sorting and this makes up a large hidden constant(256). So it runs slower on smaller kernels than simple sorting and finding the median. Again, should there be different algorithms for kernels of different size?
4. This algorithm uses a zig zag path for traversal.eg.
   -------------------------------------> <------------------------------------ ---------------------------------------> Is there any way for traversing these type of way efficiently?
5. Also for ndims(input)>=4 .The original mapwindow function starts printining N=ndims(Input) in a loop. Is there a reason for that?

[tejus-gupta]

@hari-sikchi You should also read this paper. I think it would be best to write separate functions for cases where the color values are discrete and when they are continuous.

For discrete color levels, the paper you cited works in O(n) per pixel while this paper works in O(1) per pixel. A 2D implementation is simple but there are some trade-offs involved with n-D extensions. With constant time median filtering algorithm, the memory requirement grows exponentially with number of dimensions. For both the algorithms, we can't use histograms with continuous values images and need to store the values in a BST. This will increase the time complexity to O(nlogn) and O(logn). Any of these versions should have far better performance than the mapwindow version.

As far as I understand, you are primarily facing difficulty in iterating in zig-zag manner in higher dimensions. You may want to benchmark if there is any benefit of doing this rather than simply re-initializing the histogram for the next row.

[hari-sikchi]

I agree, I am converting the code to column major for now. Let's see whether reinitializing the histogram can make it faster.

[tejus-gupta]

@timholy Are there functions to compute union and difference for CartesianRange? I think it would be easier to keep track of pixels to add/subtract if we could simply compute union/difference of CartesianRange as we move to next pixel.

Wouldn't it be cleaner(and easier to debug) to write a separate median_filter function?

[hari-sikchi]

@tejus-gupta I have made changes in the code to be column major and obtained significant time reduction. Also, zig-zag traversal increases the time complexity. In the paper, for constant time median filter per pixel we don't need to have BST or union data structures. Only requirement is more memory which increases with dimension largely as:dim[2] x dim[3] x dim[4] x .... x dim[n] .

The point of inflection is around 9 x 9 filter size for 512 x 512 input, after that the new implementation runs faster.

[hari-sikchi]

I am on a flight, will fix this as soon as I reach back.

[timholy]

Let me know if you're still working on this or ready for a re-review.

[timholy]

Sorry, I hadn't noticed your responses to my review comments since they got folded. Since the file structure has been reorganized GitHub doesn't let me add on to the existing threads, so I'll just reply by quoting:

• Two new functions are added...Do suggest, if there can be a better dispatch method.

  Workable in principle, but currently your median_histogram! function makes too many unchecked assumptions (image values range between 0 and 1, m_histogram has 256 "slots") to be exportable on its own. If we really want to create something that works for any numeric/grayscale type, I think there probably have to be some pretty fundamental revisions. For example, rather than using discrete histograms, we could come up with some type that e.g. maintains a cyclic buffer of sorted columns and then allows you to compute the median by figuring out which column contains the true median (ideally not "median-of-medians," but a more sophisticated algorithm that will check non-median values in individual columns, if necessary, to compute the true median). I don't know of a reference for this but it seems possible (if complex) in principle.

  But the alternatives are (1) to support only N0f8, or (2) check the assumptions explicitly and throw informative errors when they are violated.

• Sadly, when I am trying to break Rwin ...is leading to more running time.

  Perhaps check the results of @code_warntype. I bet you have a type-instability somewhere.

• If the median_histogram method is restricted to N0f8, that won't be problem right?

  Correct. An even better way to phrase that line might be id = reinterpret(gray(v[i])) + 1, which exploits the "raw" integer representation underlying N0f8 directly.

[hari-sikchi]

I have changed the mapwindow dispatch function to work for types img::Union{AbstractArray{N0f8,N},AbstractArray{ColorTypes.Gray{N0f8},N}} and removed type incompatibilities. Still, when I am trying to make it run in O(n) complexity by breaking the copy method, It's running type is increasing significantly. Progress is slow as I am dealing with intern examinations right now. I have attached the @code_warntype result after breaking the copy function.

The code for complete O(n) complexity[after breaking the copy function] looks somewhat like this : (https://gist.github.com/hari-sikchi/b965eec72d45c43392512b2b9d970b82)