Adds optimised algorithm for Median Filter

src/mapwindow.jl

@@ -68,12 +68,19 @@ end
 
 mapwindow(f, img, window::AbstractArray, args...; kwargs...) = mapwindow(f, img, (window...,), args...; kwargs...)
 
+
+
 function mapwindow{T,N}(f,
                         img::AbstractArray{T,N},
                         window::Indices{N},
                         border::BorderSpecAny;
                         callmode=:copy!)
-    _mapwindow(replace_function(f), img, window, border, default_shape(f); callmode=callmode)
+    if(is_medianfilter(f))
+        median_filter(replace_function(f), img, window, border, default_shape(f); callmode=callmode)
+    else
+        _mapwindow(replace_function(f), img, window, border, default_shape(f); callmode=callmode)
+    end 
+
 end
 function _mapwindow{T,N}(f,
                          img::AbstractArray{T,N},
@@ -116,6 +123,89 @@ function _mapwindow{T,N}(f,
     out
 end
 
+
+
+# This is a implementation of Fast 2D median filter
+# http://ieeexplore.ieee.org/document/1163188/
+function median_filter{T,N}(f,
+                         img::AbstractArray{T,N},
+                         window::Indices{N},
+                         border::BorderSpecAny,
+                         shape=default_shape(f);
+                         callmode=:copy!)
+    inds = indices(img)
+    inner = _interior(inds, window)
+    if callmode == :copy!
+        buf = Array{T}(map(length, window))
+        bufrs = shape(buf)
+        Rbuf = CartesianRange(size(buf))
+        offset = CartesianIndex(map(w->first(w)-1, window))
+        # To allocate the output, we have to evaluate f once
+        Rinner = CartesianRange(inner)
+        # Initialise the mode to zero and histogram consisting of 255 bins to zeros
+        mode = 0
+        m_histogram=zeros(Int64,(256,))
+        if !isempty(Rinner)
+            out = similar(img, typeof(f(bufrs,m_histogram,mode,window)))
+            Rwin = CartesianRange(map(+, window, first(Rinner).I))
+            copy!(buf, Rbuf, img, Rwin)
+            prev_mode=0
+            prev_col=1            
+            for I in Rinner
+                curr_col=I.I[2]
+                if(column_change(curr_col,prev_col))
+                    m_histogram=zeros(Int64,(256,))
+                    prev_mode=0
+                end
+                prev_col=curr_col
+
+                Rwin = CartesianRange(map(+, window, I.I))                                
+                copy!(buf, Rbuf, img, Rwin)
+                # Mode 0 corresponds to refilling the empty histogram with all the points in the window
+                if prev_mode == 0
+                    out[I] = f(bufrs,m_histogram,0,window)
+                    prev_mode=1
+                    continue
+                # Mode 1 corresponds to adding only the new points to the histogram and removing the old ones
+                elseif prev_mode == 1
+                    out[I] = f(bufrs,m_histogram,1,window)
+                    prev_mode=1
+                    continue
+                end
+
+            end
+
+        else
+            copy_win!(buf, img, first(CartesianRange(inds)), border, offset)
+            out = similar(img, typeof(f(bufrs,m_histogram,mode,window)))
+        end
+        # Now pick up the edge points we skipped over above
+        for I in EdgeIterator(inds, inner)
+            # Handle the edge points with mode -1
+            mode =-1
+            copy_win!(buf, img, I, border, offset)
+            out[I] = f(bufrs,m_histogram,mode,window)
+        end
+
+    else
+                # TODO: implement :view
+            error("callmode $callmode not supported")
+    end
+    out
+end
+
+function column_change(curr_col,prev_col)
+    return (curr_col - prev_col!=0)
+end
+
+function is_medianfilter(::typeof(median!))
+    true
+end
+
+function is_medianfilter(::Any)
+    false
+end
+
 # For copying along the edge of the image
 function copy_win!{T,N}(buf::AbstractArray{T,N}, img, I, border::Pad, offset)
     win_inds = map(+, indices(buf), (I+offset).I)
@@ -143,6 +233,8 @@ function copy_win!{T,N}(buf::AbstractArray{T,N}, img, I, border::Fill, offset)
     buf
 end
 
+
+
 ### Optimizations for particular window-functions
 
 mapwindow(::typeof(extrema), A::AbstractArray, window::Dims) = extrema_filter(A, window)
@@ -273,13 +365,100 @@ end
 # This is slightly faster than a circular buffer
 @inline cyclecache(b, x) = b[1], (Base.tail(b)..., x)
 
+
 replace_function(f) = f
-replace_function(::typeof(median!)) = function(v)
-    inds = indices(v,1)
-    Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering()))
+
+replace_function(::typeof(median!)) = function(v,m_histogram,mode,window)
+   # Handle the boundary points with mode -1
+    if(mode==-1)
+        inds = indices(v,1)
+        return Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering()))
+    else
+
+        window_size=size(v,1)
+        dims = map(x->x.stop-x.start+1,window)
+        width=window[1].stop-window[1].start+1
+        inds = indices(v,1)
+        if mode == 0
+        # Update the histogram according to new entries
+            for i = first(inds):last(inds)
+                id= trunc(Int64,(v[i]*255))+1
+                m_histogram[id]+=1
+            end
+            counter=0
+
+            for i =1:256
+                counter+=m_histogram[i]
+            end
+        # Compute the median
+            tempsum = 0
+            m_index=-1
+            for i = 1:256
+                tempsum+= m_histogram[i]
+                if tempsum>=trunc(Int64,window_size/2)+1
+                    m_index=i-1
+                    break
+                end
+            end
+        # Clear the histogram from previous value
+            for i = first(inds):width:dims[1]*dims[2]
+                for j= i: dims[1]*dims[2]:last(inds)
+                    println(j)    
+                    id= trunc(Int64,(v[j]*255))+1
+                    m_histogram[id]-=1
+                    if(m_histogram[id]<0)
+                        println("stop")
+                    end
+                end
+            end
+            return convert(Float64,m_index)/255
+
+        elseif mode == 1
+        # Update the histogram according to new entries
+            for i =  width:width:dims[1]*dims[2]
+                for j= i: dims[1]*dims[2]:last(inds)
+                    id= trunc(Int64,(v[j]*255))+1
+                    m_histogram[id]+=1
+                end
+            end
+            counter=0
+
+            for i =1:256
+                counter+=m_histogram[i]
+            end
+            println("mode:",mode)
+            println("counter:",counter)
+
+        # Compute the median        
+            tempsum = 0
+            m_index=-1
+            for i = 1:256
+                tempsum+= m_histogram[i]
+                if tempsum>= trunc(Int64,window_size/2)+1
+                    m_index=i-1
+                    break
+                end
+            end
+        # Clear the histogram from previous value
+            for i = first(inds):width:dims[1]*dims[2]
+                for j= i: dims[1]*dims[2]:last(inds)
+                    id= trunc(Int64,(v[j]*255))+1
+                    m_histogram[id]-=1
+                    if(m_histogram[id]<0)
+                        println("stop")
+                    end
+
+                end
+            end
+            return convert(Float64,m_index)/255
+        end
+
+    end
+
+
 end
 
 default_shape(::Any) = identity
 default_shape(::typeof(median!)) = vec
 
-end
+end