From b05688d104f18d7ec9c3e24d785784269b1631f6 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Fri, 7 Jul 2017 22:56:39 +0200 Subject: [PATCH 01/11] Add Fast Median Filter --- src/mapwindow.jl | 243 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 233 insertions(+), 10 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index 8a70985..6a2529b 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -90,24 +90,110 @@ function _mapwindow{T,N}(f, 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) + R2=last(Rinner).I[2] + L1=first(Rinner).I[1] + R1=first(Rinner).I[2] + L2=last(Rinner).I[1] Rwin = CartesianRange(map(+, window, first(Rinner).I)) copy!(buf, Rbuf, img, Rwin) - out = similar(img, typeof(f(bufrs))) + out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) + prev_mode= 0 # Handle the interior - for I in Rinner - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - out[I] = f(bufrs) + m_histogram=zeros(Int64,(256,)) + for i=L1:L2 + if (i-L1)%2==0 + RinnerN= CartesianRange(CartesianIndex((i,R1)), CartesianIndex((i,R2))) + for I in RinnerN + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + if (I[2]== R2 && prev_mode==1) + out[I] = f(bufrs,m_histogram,3,window) + prev_mode=3 + continue + elseif (I[2]==R1 && prev_mode==2) + out[I] = f(bufrs,m_histogram,4,window) + prev_mode=4 + continue + elseif prev_mode == 3 + out[I] = f(bufrs,m_histogram,5,window) + prev_mode = 2 + continue + elseif prev_mode == 4 + + out[I] = f(bufrs,m_histogram,6,window) + prev_mode = 1 + continue + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=1 + continue + elseif prev_mode == 1 + out[I] = f(bufrs,m_histogram,1,window) + prev_mode=1 + continue + elseif prev_mode ==2 + out[I] = f(bufrs,m_histogram,2,window) + prev_mode=2 + continue + end + end + + else + + for k= R2:-1:R1 + I=CartesianIndex((i,k)) + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + if (I[2]== R2 && prev_mode==1) + out[I] = f(bufrs,m_histogram,3,window) + prev_mode=3 + continue + elseif (I[2]==R1 && prev_mode==2) + out[I] = f(bufrs,m_histogram,4,window) + prev_mode=4 + continue + elseif prev_mode == 3 + out[I] = f(bufrs,m_histogram,5,window) + prev_mode = 2 + continue + elseif prev_mode == 4 + out[I] = f(bufrs,m_histogram,6,window) + prev_mode = 1 + continue + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=2 + continue + elseif prev_mode == 1 + out[I] = f(bufrs,m_histogram,1,window) + prev_mode=1 + continue + elseif prev_mode == 2 + out[I] = f(bufrs,m_histogram,2,window) + prev_mode=2 + continue + end + end + + end + end else copy_win!(buf, img, first(CartesianRange(inds)), border, offset) - out = similar(img, typeof(f(bufrs))) + 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 0 + mode =0 + m_histogram=zeros(Int64,(256,)) copy_win!(buf, img, I, border, offset) - out[I] = f(bufrs) + out[I] = f(bufrs,m_histogram,mode,window) end else # TODO: implement :view @@ -274,12 +360,149 @@ end @inline cyclecache(b, x) = b[1], (Base.tail(b)..., x) replace_function(f) = f -replace_function(::typeof(median!)) = function(v) +replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) + # reshape the value of the window so that they are horizontal major which allows to do horizontal traversal + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) + v= reshape(p,size(v)) + m_index=-1 inds = indices(v,1) - Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering())) + if mode == 0 + for i = first(inds):last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode == 1 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode == 2 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode == 3 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode ==4 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode ==5 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if m_histogram[i]<0 + println("stop:",mode) + end + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode ==6 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + end end default_shape(::Any) = identity default_shape(::typeof(median!)) = vec -end +end \ No newline at end of file From f0351ac1d7b41e2c412869f16a67e46ddfb97040 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Sat, 8 Jul 2017 13:18:27 +0200 Subject: [PATCH 02/11] Change median_filter to a seperate function --- src/mapwindow.jl | 70 ++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 68 insertions(+), 2 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index 6a2529b..9551b28 100644 --- a/src/mapwindow.jl +++ b/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(uses_histogram(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}, @@ -83,6 +90,49 @@ function _mapwindow{T,N}(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) + if !isempty(Rinner) + Rwin = CartesianRange(map(+, window, first(Rinner).I)) + copy!(buf, Rbuf, img, Rwin) + out = similar(img, typeof(f(bufrs))) + # Handle the interior + for I in Rinner + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + out[I] = f(bufrs) + end + else + copy_win!(buf, img, first(CartesianRange(inds)), border, offset) + out = similar(img, typeof(f(bufrs))) + end + # Now pick up the edge points we skipped over above + for I in EdgeIterator(inds, inner) + copy_win!(buf, img, I, border, offset) + out[I] = f(bufrs) + end + else + # TODO: implement :view + error("callmode $callmode not supported") + end + out +end + + +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) @@ -91,6 +141,7 @@ function _mapwindow{T,N}(f, # 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,)) @@ -196,12 +247,22 @@ function _mapwindow{T,N}(f, out[I] = f(bufrs,m_histogram,mode,window) end else - # TODO: implement :view + # TODO: implement :view error("callmode $callmode not supported") end + + out end +function uses_histogram(::typeof(median!)) + true +end + +function uses_histogram(::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) @@ -229,6 +290,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) @@ -360,8 +423,11 @@ end @inline cyclecache(b, x) = b[1], (Base.tail(b)..., x) replace_function(f) = f + replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) # reshape the value of the window so that they are horizontal major which allows to do horizontal traversal + + width=window[2].stop-window[2].start+1 height= window[1].stop-window[1].start+1 v_reshape= reshape(v,(height,width)) From 376cfc4f2ba5255f7e52975655d1b164dcb36e23 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Sun, 9 Jul 2017 00:07:03 +0200 Subject: [PATCH 03/11] Add nD support for median filter --- src/mapwindow2.jl | 1076 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1076 insertions(+) create mode 100644 src/mapwindow2.jl diff --git a/src/mapwindow2.jl b/src/mapwindow2.jl new file mode 100644 index 0000000..7fbd20d --- /dev/null +++ b/src/mapwindow2.jl @@ -0,0 +1,1076 @@ +module MapWindow + +using DataStructures, TiledIteration +using ..ImageFiltering: BorderSpecAny, Pad, Fill, borderinstance, _interior, padindex, imfilter +using Base: Indices, tail + +export mapwindow + +""" + mapwindow(f, img, window, [border="replicate"]) -> imgf + +Apply `f` to sliding windows of `img`, with window size or indices +specified by `window`. For example, `mapwindow(median!, img, window)` +returns an `Array` of values similar to `img` (median-filtered, of +course), whereas `mapwindow(extrema, img, window)` returns an `Array` +of `(min,max)` tuples over a window of size `window` centered on each +point of `img`. + +The function `f` receives a buffer `buf` for the window of data +surrounding the current point. If `window` is specified as a +Dims-tuple (tuple-of-integers), then all the integers must be odd and +the window is centered around the current image point. For example, if +`window=(3,3)`, then `f` will receive an Array `buf` corresponding to +offsets `(-1:1, -1:1)` from the `imgf[i,j]` for which this is +currently being computed. Alternatively, `window` can be a tuple of +AbstractUnitRanges, in which case the specified ranges are used for +`buf`; this allows you to use asymmetric windows if needed. + +`border` specifies how the edges of `img` should be handled; see +`imfilter` for details. + +For functions that can only take `AbstractVector` inputs, you might have to +first specialize `default_shape`: + +```julia +f = v->quantile(v, 0.75) +ImageFiltering.MapWindow.default_shape(::typeof(f)) = vec +``` + +and then `mapwindow(f, img, (m,n))` should filter at the 75th quantile. + +See also: [`imfilter`](@ref). +""" +function mapwindow(f, img::AbstractArray, window::Dims, args...; kwargs...) + all(isodd(w) for w in window) || error("entries in window must be odd, got $window") + halfsize = map(w->w>>1, window) + mapwindow(f, img, map(h->-h:h, halfsize), args...; kwargs...) +end +function mapwindow(f, img::AbstractVector, window::Integer, args...; kwargs...) + isodd(window) || error("window must be odd, got $window") + h = window>>1 + mapwindow(f, img, (-h:h,), args...; kwargs...) +end + +mapwindow(f, img::AbstractArray, window::Indices; kwargs...) = + mapwindow(f, img, window, "replicate"; kwargs...) +mapwindow(f, img::AbstractVector, window::AbstractUnitRange; kwargs...) = + mapwindow(f, img, (window,); kwargs...) + +function mapwindow(f, img::AbstractArray, window::Indices, border::AbstractString; + kwargs...) + mapwindow(f, img, window, borderinstance(border); kwargs...) +end +function mapwindow(f, img::AbstractVector, window::AbstractUnitRange, border::AbstractString; + kwargs...) + mapwindow(f, img, (window,), border; kwargs...) +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!) + if(uses_histogram(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}, + 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) + if !isempty(Rinner) + Rwin = CartesianRange(map(+, window, first(Rinner).I)) + copy!(buf, Rbuf, img, Rwin) + out = similar(img, typeof(f(bufrs))) + # Handle the interior + for I in Rinner + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + out[I] = f(bufrs) + end + else + copy_win!(buf, img, first(CartesianRange(inds)), border, offset) + out = similar(img, typeof(f(bufrs))) + end + # Now pick up the edge points we skipped over above + for I in EdgeIterator(inds, inner) + copy_win!(buf, img, I, border, offset) + out[I] = f(bufrs) + end + else + # TODO: implement :view + error("callmode $callmode not supported") + end + out +end + +function dept_iter{T,N}(f,img::AbstractArray{T,N},buf,bufrs,Rbuf,offset,Rwin,out,level,depth,indice,window::Indices{N}, + border::BorderSpecAny, + shape=default_shape(f);callmode=:copy!) + inds = indices(img) + + inner = _interior(inds, window) + Rinner=CartesianRange(inner) + buf = Array{T}(map(length, window)) + bufrs = shape(buf) + Rbuf = CartesianRange(size(buf)) + offset = CartesianIndex(map(w->first(w)-1, window)) + + if(level==depth) + +# println(img[:,:,indice]) +# println( CartesianRange(img[:,:,indice])) + + inds = indices(img) + inner = _interior(inds, window) + if callmode == :copy! + Rinner = CartesianRange(inner) + mode =0 + m_histogram=zeros(Int64,(256,)) + + if !isempty(Rinner) + R2=last(Rinner).I[2] + L1=first(Rinner).I[1] + R1=first(Rinner).I[2] + L2=last(Rinner).I[1] + + prev_mode= 0 + # Handle the interior + m_histogram=zeros(Int64,(256,)) + for i=L1:L2 + if (i-L1)%2==0 + temp1=(i,R1) + temp2=(i,R2) + id1=(temp1...,indice...) + id2=(temp2...,indice...) + RinnerN= CartesianRange(CartesianIndex(id1), CartesianIndex(id2)) + #println(RinnerN) + for I in RinnerN + #println(I) + Rwin = CartesianRange(map(+, window, I.I)) + #println(Rwin) + copy!(buf, Rbuf, img, Rwin) + if (I[2]== R2 && prev_mode==1) + out[I] = f(bufrs,m_histogram,3,window) + prev_mode=3 + continue + elseif (I[2]==R1 && prev_mode==2) + out[I] = f(bufrs,m_histogram,4,window) + prev_mode=4 + continue + elseif prev_mode == 3 + out[I] = f(bufrs,m_histogram,5,window) + prev_mode = 2 + continue + elseif prev_mode == 4 + + out[I] = f(bufrs,m_histogram,6,window) + prev_mode = 1 + continue + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=1 + continue + elseif prev_mode == 1 + out[I] = f(bufrs,m_histogram,1,window) + prev_mode=1 + continue + elseif prev_mode ==2 + out[I] = f(bufrs,m_histogram,2,window) + prev_mode=2 + continue + end + end + + else + + for k= R2:-1:R1 + + I=CartesianIndex(tuple((i,k)...,indice...)) + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + if (I[2]== R2 && prev_mode==1) + out[I] = f(bufrs,m_histogram,3,window) + prev_mode=3 + continue + elseif (I[2]==R1 && prev_mode==2) + out[I] = f(bufrs,m_histogram,4,window) + prev_mode=4 + continue + elseif prev_mode == 3 + out[I] = f(bufrs,m_histogram,5,window) + prev_mode = 2 + continue + elseif prev_mode == 4 + out[I] = f(bufrs,m_histogram,6,window) + prev_mode = 1 + continue + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=2 + continue + elseif prev_mode == 1 + out[I] = f(bufrs,m_histogram,1,window) + prev_mode=1 + continue + elseif prev_mode == 2 + out[I] = f(bufrs,m_histogram,2,window) + prev_mode=2 + continue + end + end + + 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 0 + #mode =0 + #m_histogram=zeros(Int64,(256,)) + #copy_win!(buf, img, I, border, offset) + #out[I] = f(bufrs,m_histogram,mode,window) + out[I]=0.0 + end + + else + # TODO: implement :view + error("callmode $callmode not supported") + end + + return + else + for i = 2:size(img,ndims(img)-level)-1 + #println(i) + indice[ndims(img)-level-2]=i + dept_iter(f,img,buf,bufrs,Rbuf,offset,Rwin,out,level+1,depth,indice,window,border,shape) + + end + end +end + + +function iterate_points{T,N}(f,img::AbstractArray{T,N},window::Indices{N}, + border::BorderSpecAny, + shape=default_shape(f);callmode=:copy!) + + depth= ndims(img)-2 + indice=zeros(Int64,(depth,)) + inds = indices(img) + mode = 0 + m_histogram=zeros(Int64,(256,)) + + + inner = _interior(inds, window) + Rinner=CartesianRange(inner) + 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 + + # Initialise the mode to zero and histogram consisting of 255 bins to zeros + + + Rwin = CartesianRange(map(+, window, first(Rinner).I)) + copy!(buf, Rbuf, img, Rwin) + #out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) + out=similar(img) + dept_iter(f,img,buf,bufrs,Rbuf,offset,Rwin,out,0,depth,indice,window,border,shape) + return out + +end + + + + +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! + Rinner = CartesianRange(inner) + 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 + R2=last(Rinner).I[2] + L1=first(Rinner).I[1] + R1=first(Rinner).I[2] + L2=last(Rinner).I[1] + skipper= (R2-R1+1)*(L2-L1+1) + # Initialise the mode to zero and histogram consisting of 255 bins to zeros + counter=0 + mode =0 + m_histogram=zeros(Int64,(256,)) + if !isempty(Rinner) + out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) + for I in Rinner + if(counter%skipper==0) + flag=0 + if(ndims(img)==2) + flag=1 + end + indice=I.I[3:end] + #println(indice) + prev_mode= 0 + # Handle the interior + m_histogram=zeros(Int64,(256,)) + for i=L1:L2 + if (i-L1)%2==0 + temp1=(i,R1) + temp2=(i,R2) + + if(flag==0) + id1=(temp1...,indice...) + id2=(temp2...,indice...) + else + id1=temp1 + id2=temp2 + end + RinnerN= CartesianRange(CartesianIndex(id1), CartesianIndex(id2)) + #println(RinnerN) + for I in RinnerN + #println(I) + Rwin = CartesianRange(map(+, window, I.I)) + + #println(Rwin) + copy!(buf, Rbuf, img, Rwin) + if (I[2]== R2 && prev_mode==1) + out[I] = f(bufrs,m_histogram,3,window) + prev_mode=3 + continue + elseif (I[2]==R1 && prev_mode==2) + out[I] = f(bufrs,m_histogram,4,window) + prev_mode=4 + continue + elseif prev_mode == 3 + out[I] = f(bufrs,m_histogram,5,window) + prev_mode = 2 + continue + elseif prev_mode == 4 + + out[I] = f(bufrs,m_histogram,6,window) + prev_mode = 1 + continue + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=1 + continue + elseif prev_mode == 1 + out[I] = f(bufrs,m_histogram,1,window) + prev_mode=1 + continue + elseif prev_mode ==2 + out[I] = f(bufrs,m_histogram,2,window) + prev_mode=2 + continue + end + end + + else + + for k= R2:-1:R1 + if flag==1 + I=CartesianIndex((i,k)) + else + I=CartesianIndex(tuple((i,k)...,indice...)) + end + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + if (I[2]== R2 && prev_mode==1) + out[I] = f(bufrs,m_histogram,3,window) + prev_mode=3 + continue + elseif (I[2]==R1 && prev_mode==2) + out[I] = f(bufrs,m_histogram,4,window) + prev_mode=4 + continue + elseif prev_mode == 3 + out[I] = f(bufrs,m_histogram,5,window) + prev_mode = 2 + continue + elseif prev_mode == 4 + out[I] = f(bufrs,m_histogram,6,window) + prev_mode = 1 + continue + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=2 + continue + elseif prev_mode == 1 + out[I] = f(bufrs,m_histogram,1,window) + prev_mode=1 + continue + elseif prev_mode == 2 + out[I] = f(bufrs,m_histogram,2,window) + prev_mode=2 + continue + end + end + + end + end + + end + counter+=1 + 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 0 + mode =0 + m_histogram=zeros(Int64,(256,)) + copy_win!(buf, img, I, border, offset) + out[I] = f(bufrs,m_histogram,mode,window) + #out[I]=0.0 + end + + else + # TODO: implement :view + error("callmode $callmode not supported") + end + out +end + +function uses_histogram(::typeof(median!)) + true +end + +function uses_histogram(::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) + win_img_inds = map(intersect, indices(img), win_inds) + padinds = map((inner,outer)->padindex(border, inner, outer), win_img_inds, win_inds) + docopy!(buf, img, padinds) + buf +end +docopy!(buf, img, padinds::NTuple{1}) = buf[:] = view(img, padinds[1]) +docopy!(buf, img, padinds::NTuple{2}) = buf[:,:] = view(img, padinds[1], padinds[2]) +docopy!(buf, img, padinds::NTuple{3}) = buf[:,:,:] = view(img, padinds[1], padinds[2], padinds[3]) +@inline function docopy!{N}(buf, img, padinds::NTuple{N}) + @show N + colons = ntuple(d->Colon(), Val{N}) + buf[colons...] = view(img, padinds...) +end + +function copy_win!{T,N}(buf::AbstractArray{T,N}, img, I, border::Fill, offset) + R = CartesianRange(indices(img)) + Ioff = I+offset + for J in CartesianRange(indices(buf)) + K = Ioff+J + buf[J] = K ∈ R ? img[K] : convert(eltype(img), border.value) + end + buf +end + + + +### Optimizations for particular window-functions + +mapwindow(::typeof(extrema), A::AbstractArray, window::Dims) = extrema_filter(A, window) +mapwindow(::typeof(extrema), A::AbstractVector, window::Integer) = extrema_filter(A, (window,)) + +# Max-min filter + +# This is an implementation of the Lemire max-min filter +# http://arxiv.org/abs/cs.DS/0610046 + +# Monotonic wedge +immutable Wedge{T} + L::CircularDeque{T} + U::CircularDeque{T} +end +(::Type{Wedge{T}}){T}(n::Integer) = Wedge(CircularDeque{T}(n), CircularDeque{T}(n)) + +function Base.push!(W::Wedge, i::Integer) + push!(W.L, i) + push!(W.U, i) + W +end + +function addtoback!(W::Wedge, A, i, J) + mn, mx = A[i, J] + @inbounds while !isempty(W.L) && mn < A[back(W.L), J][1] + pop!(W.L) + end + @inbounds while !isempty(W.U) && mx > A[back(W.U), J][2] + pop!(W.U) + end + push!(W.L, i) + push!(W.U, i) + W +end + +function Base.empty!(W::Wedge) + empty!(W.L) + empty!(W.U) + W +end + +@inline function getextrema(A, W::Wedge, J) + (A[front(W.L), J][1], A[front(W.U), J][2]) +end + +""" + extrema_filter(A, window) --> Array{(min,max)} + +Calculate the running min/max over a window of width `window[d]` along +dimension `d`, centered on the current point. The returned array has +the same indices as the input `A`. +""" +function extrema_filter{T,N}(A::AbstractArray{T,N}, window::NTuple{N,Integer}) + _extrema_filter!([(a,a) for a in A], window...) +end +extrema_filter(A::AbstractArray, window::AbstractArray) = extrema_filter(A, (window...,)) +extrema_filter(A::AbstractArray, window) = error("`window` must have the same number of entries as dimensions of `A`") + +extrema_filter{T,N}(A::AbstractArray{T,N}, window::Integer) = extrema_filter(A, ntuple(d->window, Val{N})) + +function _extrema_filter!(A::Array, w1, w...) + if w1 > 1 + a = first(A) + cache = ntuple(i->a, w1>>1) + _extrema_filter1!(A, w1, cache) + end + _extrema_filter!(permutedims(A, [2:ndims(A);1]), w...) +end +_extrema_filter!(A::Array) = A + +# Extrema-filtering along "columns" (dimension 1). This implements Lemire +# Algorithm 1, with the following modifications: +# - multidimensional array support by looping over trailing dimensions +# - working with min/max pairs rather than plain values, to +# facilitate multidimensional processing +# - output for all points of the array, handling the edges as max-min +# over halfwindow on either side +function _extrema_filter1!{T}(A::AbstractArray{Tuple{T,T}}, window::Int, cache) + # Initialise the internal wedges + # U[1], L[1] are the location of the global (within the window) maximum and minimum + # U[2], L[2] are the maximum and minimum over (U1, end] and (L1, end], respectively + W = Wedge{Int}(window+1) + tmp = Array{Tuple{T,T}}(window) + c = z = first(cache) + + inds = indices(A) + inds1 = inds[1] + halfwindow = window>>1 + iw = min(last(inds1), first(inds1)+window-1) + for J in CartesianRange(tail(inds)) + empty!(W) + # Leading edge. We can't overwrite any values yet in A because + # we'll need them again in later computations. + for i = first(inds1):iw + addtoback!(W, A, i, J) + c, cache = cyclecache(cache, getextrema(A, W, J)) + end + # Process the rest of the "column" + for i = iw+1:last(inds1) + A[i-window, J] = c + if i == window+front(W.U) + shift!(W.U) + end + if i == window+front(W.L) + shift!(W.L) + end + addtoback!(W, A, i, J) + c, cache = cyclecache(cache, getextrema(A, W, J)) + end + for i = last(inds1)-window+1:last(inds1)-1 + if i >= first(inds1) + A[i, J] = c + end + if i == front(W.U) + shift!(W.U) + end + if i == front(W.L) + shift!(W.L) + end + c, cache = cyclecache(cache, getextrema(A, W, J)) + end + A[last(inds1), J] = c + end + A +end + +# This is slightly faster than a circular buffer +@inline cyclecache(b, x) = b[1], (Base.tail(b)..., x) + +function median_find(m_histogram,window_size) + tempsum = 0 + m_index=-1 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= window_size/2 + m_index=i-1 + break + end + end + return convert(Float64,m_index)/255 +end + + +function update_median(v_blob,m_histogram,mode,window,recursion_depth,indice,level) + v = v_blob[:,:,indice] + + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + #v_reshape=v + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) + v= reshape(p,(size(v,1)*size(v,2),)) + m_index=-1 + inds = indices(v,1) + if mode == 0 + for i = first(inds):last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 1 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 2 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 3 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==4 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==5 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==6 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + end +end + + +function clear_histogram(v_blob,m_histogram,mode,window,recursion_depth,indice,level) + v = v_blob[:,:,indice] + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) +# v= reshape(p,size(v)) + v= reshape(p,(size(v,1)*size(v,2),)) + + inds = indices(v,1) + if mode == 0 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 1 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 2 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 3 + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==4 + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==5 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==6 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + end + +end + + + +function depth_looper_delete(v_blob,m_histogram,mode,window,recursion_depth,indice,level) + if level == recursion_depth + clear_histogram(v_blob,m_histogram,mode,window,recursion_depth,indice,level) + else + for i = 1:size(v_blob,ndims(v_blob)-level) + indice[ndims(v_blob)-level-2]=i + depth_looper_delete(v_blob,m_histogram,mode,window,recursion_depth,indice,level+1) + end + end +end + +function depth_looper_update(v_blob,m_histogram,mode,window,recursion_depth,indice,level) + if level == recursion_depth + update_median(v_blob,m_histogram,mode,window,recursion_depth,indice,level) + else + for i = 1:size(v_blob,ndims(v_blob)-level) + indice[ndims(v_blob)-level-2]=i + depth_looper_update(v_blob,m_histogram,mode,window,recursion_depth,indice,level+1) + end + end +end + +replace_function(f) = f + +replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) + # reshape the value of the window so that they are horizontal major which allows to do horizontal traversal + + window_size=size(v,1) + dims = map(x->x.stop-x.start+1,window) + startp= ones(Int64,(size(dims,1),)) + v_blob=reshape(v,dims) + blob_range=CartesianRange(CartesianIndex(tuple(startp...)),CartesianIndex(tuple(dims...))) +#= recursion_depth = size(window,1)-2 + indice = zeros(Int64,(recursion_depth,)) + #println(recursion_depth) +# println("v_blob: ",size(v_blob)) + depth_looper_update(v_blob,m_histogram,mode,window,recursion_depth,indice,0) + current_median=median_find(m_histogram,window_size) + depth_looper_delete(v_blob,m_histogram,mode,window,recursion_depth,indice,0) + println(current_median) + current_median +=# + skipper= size(v_blob,1)*size(v_blob,2) + counter=0 + for I in blob_range + if(counter%skipper==0) + indice=I.I[3:end] + if(ndims(v_blob)==2) + v=v_blob + else + v=v_blob[:,:,indice...] + end + m_index=-1 + inds = indices(v,1) + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + #v_reshape=v + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) + v= reshape(p,(size(v,1)*size(v,2),)) + m_index=-1 + inds = indices(v,1) + if mode == 0 + for i = first(inds):last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 1 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 2 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 3 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==4 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==5 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==6 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + end + end + end +# Find the median using histogram + tempsum = 0 + m_index=-1 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= window_size/2 + m_index=i-1 + break + end + end + +#Clear histogram +counter=0 + for I in blob_range + if(counter%skipper==0) + indice=I.I[3:end] + if(ndims(v_blob)==2) + v=v_blob + else + v=v_blob[:,:,indice...] + end + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) + # v= reshape(p,size(v)) + v= reshape(p,(size(v,1)*size(v,2),)) + + inds = indices(v,1) + if mode == 0 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 1 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 2 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 3 + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==4 + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==5 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==6 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + end + end + end + + + + return convert(Float64,m_index)/255 + + + +#= + + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) + v= reshape(p,size(v)) + width=3 + height=3 + m_index=-1 + inds = indices(v,1) + if mode == 0 + for i = first(inds):last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode == 1 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode == 2 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode == 3 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode ==4 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode ==5 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if m_histogram[i]<0 + println("stop:",mode) + end + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + elseif mode ==6 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + tempsum = 0 + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= last(inds)/2 + m_index=i-1 + break + end + end + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + return convert(Float64,m_index)/255 + end +=# +end + +default_shape(::Any) = identity +default_shape(::typeof(median!)) = vec + +end \ No newline at end of file From d73c56ed6c44874b72be12feb24653d17a6646d0 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Sun, 9 Jul 2017 14:21:33 +0200 Subject: [PATCH 04/11] Improved time for nD fast median filter --- src/mapwindow.jl | 417 ++++++++---------- src/mapwindow2.jl | 1076 --------------------------------------------- 2 files changed, 194 insertions(+), 1299 deletions(-) delete mode 100644 src/mapwindow2.jl diff --git a/src/mapwindow.jl b/src/mapwindow.jl index 9551b28..7d00d2a 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -124,6 +124,7 @@ function _mapwindow{T,N}(f, end + function median_filter{T,N}(f, img::AbstractArray{T,N}, window::Indices{N}, @@ -132,126 +133,107 @@ function median_filter{T,N}(f, callmode=:copy!) inds = indices(img) inner = _interior(inds, window) - if callmode == :copy! + Rinner = CartesianRange(inner) 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) + R2=last(Rinner).I[2] + L1=first(Rinner).I[1] + R1=first(Rinner).I[2] + L2=last(Rinner).I[1] # Initialise the mode to zero and histogram consisting of 255 bins to zeros - mode =0 m_histogram=zeros(Int64,(256,)) - + Rinner_new= CartesianRange(CartesianIndex((first(Rinner).I[3:end])),CartesianIndex((last(Rinner).I[3:end]))) if !isempty(Rinner) - R2=last(Rinner).I[2] - L1=first(Rinner).I[1] - R1=first(Rinner).I[2] - L2=last(Rinner).I[1] - Rwin = CartesianRange(map(+, window, first(Rinner).I)) - copy!(buf, Rbuf, img, Rwin) out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) - prev_mode= 0 - # Handle the interior - m_histogram=zeros(Int64,(256,)) - for i=L1:L2 - if (i-L1)%2==0 - RinnerN= CartesianRange(CartesianIndex((i,R1)), CartesianIndex((i,R2))) - for I in RinnerN - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - if (I[2]== R2 && prev_mode==1) - out[I] = f(bufrs,m_histogram,3,window) - prev_mode=3 - continue - elseif (I[2]==R1 && prev_mode==2) - out[I] = f(bufrs,m_histogram,4,window) - prev_mode=4 - continue - elseif prev_mode == 3 - out[I] = f(bufrs,m_histogram,5,window) - prev_mode = 2 - continue - elseif prev_mode == 4 - - out[I] = f(bufrs,m_histogram,6,window) - prev_mode = 1 - continue - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=1 - continue - elseif prev_mode == 1 - out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 - continue - elseif prev_mode ==2 - out[I] = f(bufrs,m_histogram,2,window) - prev_mode=2 - continue + for I in Rinner_new + indice=I.I + prev_mode= 0 + # Handle the interior + m_histogram=zeros(Int64,(256,)) + for i=L1:L2 + if (i-L1)%2==0 + temp1=(i,R1) + temp2=(i,R2) + id1=(temp1...,indice...) + id2=(temp2...,indice...) + RinnerN= CartesianRange(CartesianIndex(id1), CartesianIndex(id2)) + for I in RinnerN + + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + if (I[2]== R2 && prev_mode==1) + out[I] = f(bufrs,m_histogram,3,window) + prev_mode=3 + continue + elseif prev_mode == 4 + out[I] = f(bufrs,m_histogram,6,window) + prev_mode = 1 + continue + + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=1 + continue + elseif prev_mode == 1 + out[I] = f(bufrs,m_histogram,1,window) + prev_mode=1 + continue + end + end - end - else - - for k= R2:-1:R1 - I=CartesianIndex((i,k)) - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - if (I[2]== R2 && prev_mode==1) - out[I] = f(bufrs,m_histogram,3,window) - prev_mode=3 - continue - elseif (I[2]==R1 && prev_mode==2) - out[I] = f(bufrs,m_histogram,4,window) - prev_mode=4 - continue - elseif prev_mode == 3 - out[I] = f(bufrs,m_histogram,5,window) - prev_mode = 2 - continue - elseif prev_mode == 4 - out[I] = f(bufrs,m_histogram,6,window) - prev_mode = 1 - continue - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=2 - continue - elseif prev_mode == 1 - out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 - continue - elseif prev_mode == 2 - out[I] = f(bufrs,m_histogram,2,window) - prev_mode=2 - continue + else + + for k= R2:-1:R1 + I=CartesianIndex(tuple((i,k)...,indice...)) + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + if (I[2]==R1 && prev_mode==2) + out[I] = f(bufrs,m_histogram,4,window) + prev_mode=4 + continue + elseif prev_mode == 3 + out[I] = f(bufrs,m_histogram,5,window) + prev_mode = 2 + continue + elseif prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=2 + continue + elseif prev_mode == 2 + out[I] = f(bufrs,m_histogram,2,window) + prev_mode=2 + continue + end end - end + end end - - 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 0 - mode =0 - m_histogram=zeros(Int64,(256,)) + # Handle the edge points with mode -1 + mode =-1 copy_win!(buf, img, I, border, offset) out[I] = f(bufrs,m_histogram,mode,window) + #out[I]=0 end + else - # TODO: implement :view - error("callmode $callmode not supported") + # TODO: implement :view + error("callmode $callmode not supported") end - - out end @@ -422,150 +404,139 @@ 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,m_histogram,mode,window) - # reshape the value of the window so that they are horizontal major which allows to do horizontal traversal - - - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) - v= reshape(p,size(v)) - m_index=-1 - inds = indices(v,1) - if mode == 0 - for i = first(inds):last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode == 1 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode == 2 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode == 3 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode ==4 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode ==5 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if m_histogram[i]<0 - println("stop:",mode) - end - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode ==6 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 + # 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 + # reshape the value of the window so that they are horizontal major which allows to do horizontal traversal + window_size=size(v,1) + dims = map(x->x.stop-x.start+1,window) + startp= ones(Int64,(size(dims,1),)) + v_blob=reshape(v,dims) + blob_range=CartesianRange(CartesianIndex(tuple(startp...)),CartesianIndex(tuple(dims...))) + blob_range_new=CartesianRange(CartesianIndex(first(blob_range).I[3:end]),CartesianIndex(last(blob_range).I[3:end])) + # Update the histogram according to + for I in blob_range_new + indice=I.I + v=v_blob[:,:,indice...] + m_index=-1 + inds = indices(v,1) + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) + v= reshape(p,(size(v,1)*size(v,2),)) + m_index=-1 + inds = indices(v,1) + if mode == 0 + for i = first(inds):last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 1 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 2 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode == 3 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==4 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==5 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + elseif mode ==6 + for i = last(inds)-width+1:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]+=1 + end + end + end + # Find the median using histogram tempsum = 0 + m_index=-1 for i = 1:256 tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 + if tempsum>= window_size/2 m_index=i-1 break end end - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 + + # Clear histogram from useless past items + for I in blob_range_new + indice=I.I + v=v_blob[:,:,indice...] + + width=window[2].stop-window[2].start+1 + height= window[1].stop-window[1].start+1 + v_reshape= reshape(v,(height,width)) + p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) + transpose!(p,v_reshape) + v= reshape(p,(size(v,1)*size(v,2),)) + + inds = indices(v,1) + if mode == 0 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 1 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 2 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode == 3 + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==4 + for i = first(inds):first(inds)+width-1 + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==5 + for i = width:width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + elseif mode ==6 + for i = first(inds):width:last(inds) + id= trunc(Int64,(v[i]*255))+1 + m_histogram[id]-=1 + end + end + end return convert(Float64,m_index)/255 end + end default_shape(::Any) = identity diff --git a/src/mapwindow2.jl b/src/mapwindow2.jl deleted file mode 100644 index 7fbd20d..0000000 --- a/src/mapwindow2.jl +++ /dev/null @@ -1,1076 +0,0 @@ -module MapWindow - -using DataStructures, TiledIteration -using ..ImageFiltering: BorderSpecAny, Pad, Fill, borderinstance, _interior, padindex, imfilter -using Base: Indices, tail - -export mapwindow - -""" - mapwindow(f, img, window, [border="replicate"]) -> imgf - -Apply `f` to sliding windows of `img`, with window size or indices -specified by `window`. For example, `mapwindow(median!, img, window)` -returns an `Array` of values similar to `img` (median-filtered, of -course), whereas `mapwindow(extrema, img, window)` returns an `Array` -of `(min,max)` tuples over a window of size `window` centered on each -point of `img`. - -The function `f` receives a buffer `buf` for the window of data -surrounding the current point. If `window` is specified as a -Dims-tuple (tuple-of-integers), then all the integers must be odd and -the window is centered around the current image point. For example, if -`window=(3,3)`, then `f` will receive an Array `buf` corresponding to -offsets `(-1:1, -1:1)` from the `imgf[i,j]` for which this is -currently being computed. Alternatively, `window` can be a tuple of -AbstractUnitRanges, in which case the specified ranges are used for -`buf`; this allows you to use asymmetric windows if needed. - -`border` specifies how the edges of `img` should be handled; see -`imfilter` for details. - -For functions that can only take `AbstractVector` inputs, you might have to -first specialize `default_shape`: - -```julia -f = v->quantile(v, 0.75) -ImageFiltering.MapWindow.default_shape(::typeof(f)) = vec -``` - -and then `mapwindow(f, img, (m,n))` should filter at the 75th quantile. - -See also: [`imfilter`](@ref). -""" -function mapwindow(f, img::AbstractArray, window::Dims, args...; kwargs...) - all(isodd(w) for w in window) || error("entries in window must be odd, got $window") - halfsize = map(w->w>>1, window) - mapwindow(f, img, map(h->-h:h, halfsize), args...; kwargs...) -end -function mapwindow(f, img::AbstractVector, window::Integer, args...; kwargs...) - isodd(window) || error("window must be odd, got $window") - h = window>>1 - mapwindow(f, img, (-h:h,), args...; kwargs...) -end - -mapwindow(f, img::AbstractArray, window::Indices; kwargs...) = - mapwindow(f, img, window, "replicate"; kwargs...) -mapwindow(f, img::AbstractVector, window::AbstractUnitRange; kwargs...) = - mapwindow(f, img, (window,); kwargs...) - -function mapwindow(f, img::AbstractArray, window::Indices, border::AbstractString; - kwargs...) - mapwindow(f, img, window, borderinstance(border); kwargs...) -end -function mapwindow(f, img::AbstractVector, window::AbstractUnitRange, border::AbstractString; - kwargs...) - mapwindow(f, img, (window,), border; kwargs...) -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!) - if(uses_histogram(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}, - 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) - if !isempty(Rinner) - Rwin = CartesianRange(map(+, window, first(Rinner).I)) - copy!(buf, Rbuf, img, Rwin) - out = similar(img, typeof(f(bufrs))) - # Handle the interior - for I in Rinner - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - out[I] = f(bufrs) - end - else - copy_win!(buf, img, first(CartesianRange(inds)), border, offset) - out = similar(img, typeof(f(bufrs))) - end - # Now pick up the edge points we skipped over above - for I in EdgeIterator(inds, inner) - copy_win!(buf, img, I, border, offset) - out[I] = f(bufrs) - end - else - # TODO: implement :view - error("callmode $callmode not supported") - end - out -end - -function dept_iter{T,N}(f,img::AbstractArray{T,N},buf,bufrs,Rbuf,offset,Rwin,out,level,depth,indice,window::Indices{N}, - border::BorderSpecAny, - shape=default_shape(f);callmode=:copy!) - inds = indices(img) - - inner = _interior(inds, window) - Rinner=CartesianRange(inner) - buf = Array{T}(map(length, window)) - bufrs = shape(buf) - Rbuf = CartesianRange(size(buf)) - offset = CartesianIndex(map(w->first(w)-1, window)) - - if(level==depth) - -# println(img[:,:,indice]) -# println( CartesianRange(img[:,:,indice])) - - inds = indices(img) - inner = _interior(inds, window) - if callmode == :copy! - Rinner = CartesianRange(inner) - mode =0 - m_histogram=zeros(Int64,(256,)) - - if !isempty(Rinner) - R2=last(Rinner).I[2] - L1=first(Rinner).I[1] - R1=first(Rinner).I[2] - L2=last(Rinner).I[1] - - prev_mode= 0 - # Handle the interior - m_histogram=zeros(Int64,(256,)) - for i=L1:L2 - if (i-L1)%2==0 - temp1=(i,R1) - temp2=(i,R2) - id1=(temp1...,indice...) - id2=(temp2...,indice...) - RinnerN= CartesianRange(CartesianIndex(id1), CartesianIndex(id2)) - #println(RinnerN) - for I in RinnerN - #println(I) - Rwin = CartesianRange(map(+, window, I.I)) - #println(Rwin) - copy!(buf, Rbuf, img, Rwin) - if (I[2]== R2 && prev_mode==1) - out[I] = f(bufrs,m_histogram,3,window) - prev_mode=3 - continue - elseif (I[2]==R1 && prev_mode==2) - out[I] = f(bufrs,m_histogram,4,window) - prev_mode=4 - continue - elseif prev_mode == 3 - out[I] = f(bufrs,m_histogram,5,window) - prev_mode = 2 - continue - elseif prev_mode == 4 - - out[I] = f(bufrs,m_histogram,6,window) - prev_mode = 1 - continue - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=1 - continue - elseif prev_mode == 1 - out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 - continue - elseif prev_mode ==2 - out[I] = f(bufrs,m_histogram,2,window) - prev_mode=2 - continue - end - end - - else - - for k= R2:-1:R1 - - I=CartesianIndex(tuple((i,k)...,indice...)) - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - if (I[2]== R2 && prev_mode==1) - out[I] = f(bufrs,m_histogram,3,window) - prev_mode=3 - continue - elseif (I[2]==R1 && prev_mode==2) - out[I] = f(bufrs,m_histogram,4,window) - prev_mode=4 - continue - elseif prev_mode == 3 - out[I] = f(bufrs,m_histogram,5,window) - prev_mode = 2 - continue - elseif prev_mode == 4 - out[I] = f(bufrs,m_histogram,6,window) - prev_mode = 1 - continue - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=2 - continue - elseif prev_mode == 1 - out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 - continue - elseif prev_mode == 2 - out[I] = f(bufrs,m_histogram,2,window) - prev_mode=2 - continue - end - end - - 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 0 - #mode =0 - #m_histogram=zeros(Int64,(256,)) - #copy_win!(buf, img, I, border, offset) - #out[I] = f(bufrs,m_histogram,mode,window) - out[I]=0.0 - end - - else - # TODO: implement :view - error("callmode $callmode not supported") - end - - return - else - for i = 2:size(img,ndims(img)-level)-1 - #println(i) - indice[ndims(img)-level-2]=i - dept_iter(f,img,buf,bufrs,Rbuf,offset,Rwin,out,level+1,depth,indice,window,border,shape) - - end - end -end - - -function iterate_points{T,N}(f,img::AbstractArray{T,N},window::Indices{N}, - border::BorderSpecAny, - shape=default_shape(f);callmode=:copy!) - - depth= ndims(img)-2 - indice=zeros(Int64,(depth,)) - inds = indices(img) - mode = 0 - m_histogram=zeros(Int64,(256,)) - - - inner = _interior(inds, window) - Rinner=CartesianRange(inner) - 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 - - # Initialise the mode to zero and histogram consisting of 255 bins to zeros - - - Rwin = CartesianRange(map(+, window, first(Rinner).I)) - copy!(buf, Rbuf, img, Rwin) - #out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) - out=similar(img) - dept_iter(f,img,buf,bufrs,Rbuf,offset,Rwin,out,0,depth,indice,window,border,shape) - return out - -end - - - - -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! - Rinner = CartesianRange(inner) - 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 - R2=last(Rinner).I[2] - L1=first(Rinner).I[1] - R1=first(Rinner).I[2] - L2=last(Rinner).I[1] - skipper= (R2-R1+1)*(L2-L1+1) - # Initialise the mode to zero and histogram consisting of 255 bins to zeros - counter=0 - mode =0 - m_histogram=zeros(Int64,(256,)) - if !isempty(Rinner) - out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) - for I in Rinner - if(counter%skipper==0) - flag=0 - if(ndims(img)==2) - flag=1 - end - indice=I.I[3:end] - #println(indice) - prev_mode= 0 - # Handle the interior - m_histogram=zeros(Int64,(256,)) - for i=L1:L2 - if (i-L1)%2==0 - temp1=(i,R1) - temp2=(i,R2) - - if(flag==0) - id1=(temp1...,indice...) - id2=(temp2...,indice...) - else - id1=temp1 - id2=temp2 - end - RinnerN= CartesianRange(CartesianIndex(id1), CartesianIndex(id2)) - #println(RinnerN) - for I in RinnerN - #println(I) - Rwin = CartesianRange(map(+, window, I.I)) - - #println(Rwin) - copy!(buf, Rbuf, img, Rwin) - if (I[2]== R2 && prev_mode==1) - out[I] = f(bufrs,m_histogram,3,window) - prev_mode=3 - continue - elseif (I[2]==R1 && prev_mode==2) - out[I] = f(bufrs,m_histogram,4,window) - prev_mode=4 - continue - elseif prev_mode == 3 - out[I] = f(bufrs,m_histogram,5,window) - prev_mode = 2 - continue - elseif prev_mode == 4 - - out[I] = f(bufrs,m_histogram,6,window) - prev_mode = 1 - continue - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=1 - continue - elseif prev_mode == 1 - out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 - continue - elseif prev_mode ==2 - out[I] = f(bufrs,m_histogram,2,window) - prev_mode=2 - continue - end - end - - else - - for k= R2:-1:R1 - if flag==1 - I=CartesianIndex((i,k)) - else - I=CartesianIndex(tuple((i,k)...,indice...)) - end - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - if (I[2]== R2 && prev_mode==1) - out[I] = f(bufrs,m_histogram,3,window) - prev_mode=3 - continue - elseif (I[2]==R1 && prev_mode==2) - out[I] = f(bufrs,m_histogram,4,window) - prev_mode=4 - continue - elseif prev_mode == 3 - out[I] = f(bufrs,m_histogram,5,window) - prev_mode = 2 - continue - elseif prev_mode == 4 - out[I] = f(bufrs,m_histogram,6,window) - prev_mode = 1 - continue - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=2 - continue - elseif prev_mode == 1 - out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 - continue - elseif prev_mode == 2 - out[I] = f(bufrs,m_histogram,2,window) - prev_mode=2 - continue - end - end - - end - end - - end - counter+=1 - 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 0 - mode =0 - m_histogram=zeros(Int64,(256,)) - copy_win!(buf, img, I, border, offset) - out[I] = f(bufrs,m_histogram,mode,window) - #out[I]=0.0 - end - - else - # TODO: implement :view - error("callmode $callmode not supported") - end - out -end - -function uses_histogram(::typeof(median!)) - true -end - -function uses_histogram(::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) - win_img_inds = map(intersect, indices(img), win_inds) - padinds = map((inner,outer)->padindex(border, inner, outer), win_img_inds, win_inds) - docopy!(buf, img, padinds) - buf -end -docopy!(buf, img, padinds::NTuple{1}) = buf[:] = view(img, padinds[1]) -docopy!(buf, img, padinds::NTuple{2}) = buf[:,:] = view(img, padinds[1], padinds[2]) -docopy!(buf, img, padinds::NTuple{3}) = buf[:,:,:] = view(img, padinds[1], padinds[2], padinds[3]) -@inline function docopy!{N}(buf, img, padinds::NTuple{N}) - @show N - colons = ntuple(d->Colon(), Val{N}) - buf[colons...] = view(img, padinds...) -end - -function copy_win!{T,N}(buf::AbstractArray{T,N}, img, I, border::Fill, offset) - R = CartesianRange(indices(img)) - Ioff = I+offset - for J in CartesianRange(indices(buf)) - K = Ioff+J - buf[J] = K ∈ R ? img[K] : convert(eltype(img), border.value) - end - buf -end - - - -### Optimizations for particular window-functions - -mapwindow(::typeof(extrema), A::AbstractArray, window::Dims) = extrema_filter(A, window) -mapwindow(::typeof(extrema), A::AbstractVector, window::Integer) = extrema_filter(A, (window,)) - -# Max-min filter - -# This is an implementation of the Lemire max-min filter -# http://arxiv.org/abs/cs.DS/0610046 - -# Monotonic wedge -immutable Wedge{T} - L::CircularDeque{T} - U::CircularDeque{T} -end -(::Type{Wedge{T}}){T}(n::Integer) = Wedge(CircularDeque{T}(n), CircularDeque{T}(n)) - -function Base.push!(W::Wedge, i::Integer) - push!(W.L, i) - push!(W.U, i) - W -end - -function addtoback!(W::Wedge, A, i, J) - mn, mx = A[i, J] - @inbounds while !isempty(W.L) && mn < A[back(W.L), J][1] - pop!(W.L) - end - @inbounds while !isempty(W.U) && mx > A[back(W.U), J][2] - pop!(W.U) - end - push!(W.L, i) - push!(W.U, i) - W -end - -function Base.empty!(W::Wedge) - empty!(W.L) - empty!(W.U) - W -end - -@inline function getextrema(A, W::Wedge, J) - (A[front(W.L), J][1], A[front(W.U), J][2]) -end - -""" - extrema_filter(A, window) --> Array{(min,max)} - -Calculate the running min/max over a window of width `window[d]` along -dimension `d`, centered on the current point. The returned array has -the same indices as the input `A`. -""" -function extrema_filter{T,N}(A::AbstractArray{T,N}, window::NTuple{N,Integer}) - _extrema_filter!([(a,a) for a in A], window...) -end -extrema_filter(A::AbstractArray, window::AbstractArray) = extrema_filter(A, (window...,)) -extrema_filter(A::AbstractArray, window) = error("`window` must have the same number of entries as dimensions of `A`") - -extrema_filter{T,N}(A::AbstractArray{T,N}, window::Integer) = extrema_filter(A, ntuple(d->window, Val{N})) - -function _extrema_filter!(A::Array, w1, w...) - if w1 > 1 - a = first(A) - cache = ntuple(i->a, w1>>1) - _extrema_filter1!(A, w1, cache) - end - _extrema_filter!(permutedims(A, [2:ndims(A);1]), w...) -end -_extrema_filter!(A::Array) = A - -# Extrema-filtering along "columns" (dimension 1). This implements Lemire -# Algorithm 1, with the following modifications: -# - multidimensional array support by looping over trailing dimensions -# - working with min/max pairs rather than plain values, to -# facilitate multidimensional processing -# - output for all points of the array, handling the edges as max-min -# over halfwindow on either side -function _extrema_filter1!{T}(A::AbstractArray{Tuple{T,T}}, window::Int, cache) - # Initialise the internal wedges - # U[1], L[1] are the location of the global (within the window) maximum and minimum - # U[2], L[2] are the maximum and minimum over (U1, end] and (L1, end], respectively - W = Wedge{Int}(window+1) - tmp = Array{Tuple{T,T}}(window) - c = z = first(cache) - - inds = indices(A) - inds1 = inds[1] - halfwindow = window>>1 - iw = min(last(inds1), first(inds1)+window-1) - for J in CartesianRange(tail(inds)) - empty!(W) - # Leading edge. We can't overwrite any values yet in A because - # we'll need them again in later computations. - for i = first(inds1):iw - addtoback!(W, A, i, J) - c, cache = cyclecache(cache, getextrema(A, W, J)) - end - # Process the rest of the "column" - for i = iw+1:last(inds1) - A[i-window, J] = c - if i == window+front(W.U) - shift!(W.U) - end - if i == window+front(W.L) - shift!(W.L) - end - addtoback!(W, A, i, J) - c, cache = cyclecache(cache, getextrema(A, W, J)) - end - for i = last(inds1)-window+1:last(inds1)-1 - if i >= first(inds1) - A[i, J] = c - end - if i == front(W.U) - shift!(W.U) - end - if i == front(W.L) - shift!(W.L) - end - c, cache = cyclecache(cache, getextrema(A, W, J)) - end - A[last(inds1), J] = c - end - A -end - -# This is slightly faster than a circular buffer -@inline cyclecache(b, x) = b[1], (Base.tail(b)..., x) - -function median_find(m_histogram,window_size) - tempsum = 0 - m_index=-1 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= window_size/2 - m_index=i-1 - break - end - end - return convert(Float64,m_index)/255 -end - - -function update_median(v_blob,m_histogram,mode,window,recursion_depth,indice,level) - v = v_blob[:,:,indice] - - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - #v_reshape=v - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) - v= reshape(p,(size(v,1)*size(v,2),)) - m_index=-1 - inds = indices(v,1) - if mode == 0 - for i = first(inds):last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 1 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 2 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 3 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==4 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==5 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==6 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - end -end - - -function clear_histogram(v_blob,m_histogram,mode,window,recursion_depth,indice,level) - v = v_blob[:,:,indice] - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) -# v= reshape(p,size(v)) - v= reshape(p,(size(v,1)*size(v,2),)) - - inds = indices(v,1) - if mode == 0 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 1 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 2 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 3 - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==4 - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==5 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==6 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - end - -end - - - -function depth_looper_delete(v_blob,m_histogram,mode,window,recursion_depth,indice,level) - if level == recursion_depth - clear_histogram(v_blob,m_histogram,mode,window,recursion_depth,indice,level) - else - for i = 1:size(v_blob,ndims(v_blob)-level) - indice[ndims(v_blob)-level-2]=i - depth_looper_delete(v_blob,m_histogram,mode,window,recursion_depth,indice,level+1) - end - end -end - -function depth_looper_update(v_blob,m_histogram,mode,window,recursion_depth,indice,level) - if level == recursion_depth - update_median(v_blob,m_histogram,mode,window,recursion_depth,indice,level) - else - for i = 1:size(v_blob,ndims(v_blob)-level) - indice[ndims(v_blob)-level-2]=i - depth_looper_update(v_blob,m_histogram,mode,window,recursion_depth,indice,level+1) - end - end -end - -replace_function(f) = f - -replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) - # reshape the value of the window so that they are horizontal major which allows to do horizontal traversal - - window_size=size(v,1) - dims = map(x->x.stop-x.start+1,window) - startp= ones(Int64,(size(dims,1),)) - v_blob=reshape(v,dims) - blob_range=CartesianRange(CartesianIndex(tuple(startp...)),CartesianIndex(tuple(dims...))) -#= recursion_depth = size(window,1)-2 - indice = zeros(Int64,(recursion_depth,)) - #println(recursion_depth) -# println("v_blob: ",size(v_blob)) - depth_looper_update(v_blob,m_histogram,mode,window,recursion_depth,indice,0) - current_median=median_find(m_histogram,window_size) - depth_looper_delete(v_blob,m_histogram,mode,window,recursion_depth,indice,0) - println(current_median) - current_median -=# - skipper= size(v_blob,1)*size(v_blob,2) - counter=0 - for I in blob_range - if(counter%skipper==0) - indice=I.I[3:end] - if(ndims(v_blob)==2) - v=v_blob - else - v=v_blob[:,:,indice...] - end - m_index=-1 - inds = indices(v,1) - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - #v_reshape=v - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) - v= reshape(p,(size(v,1)*size(v,2),)) - m_index=-1 - inds = indices(v,1) - if mode == 0 - for i = first(inds):last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 1 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 2 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 3 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==4 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==5 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==6 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - end - end - end -# Find the median using histogram - tempsum = 0 - m_index=-1 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= window_size/2 - m_index=i-1 - break - end - end - -#Clear histogram -counter=0 - for I in blob_range - if(counter%skipper==0) - indice=I.I[3:end] - if(ndims(v_blob)==2) - v=v_blob - else - v=v_blob[:,:,indice...] - end - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) - # v= reshape(p,size(v)) - v= reshape(p,(size(v,1)*size(v,2),)) - - inds = indices(v,1) - if mode == 0 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 1 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 2 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 3 - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==4 - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==5 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==6 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - end - end - end - - - - return convert(Float64,m_index)/255 - - - -#= - - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) - v= reshape(p,size(v)) - width=3 - height=3 - m_index=-1 - inds = indices(v,1) - if mode == 0 - for i = first(inds):last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode == 1 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode == 2 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode == 3 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode ==4 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode ==5 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if m_histogram[i]<0 - println("stop:",mode) - end - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - elseif mode ==6 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - tempsum = 0 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= last(inds)/2 - m_index=i-1 - break - end - end - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - return convert(Float64,m_index)/255 - end -=# -end - -default_shape(::Any) = identity -default_shape(::typeof(median!)) = vec - -end \ No newline at end of file From 3bbafaf5d690a4d16bde41f186d68daf3a49c2f3 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Sun, 9 Jul 2017 21:37:11 +0200 Subject: [PATCH 05/11] Changed traversal to column major --- src/mapwindow.jl | 265 ++++++++++++++--------------------------------- 1 file changed, 78 insertions(+), 187 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index 7d00d2a..aca7951 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -134,88 +134,42 @@ function median_filter{T,N}(f, inds = indices(img) inner = _interior(inds, window) if callmode == :copy! - Rinner = CartesianRange(inner) 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 - R2=last(Rinner).I[2] - L1=first(Rinner).I[1] - R1=first(Rinner).I[2] - L2=last(Rinner).I[1] + Rinner = CartesianRange(inner) # Initialise the mode to zero and histogram consisting of 255 bins to zeros - mode =0 + mode = 0 m_histogram=zeros(Int64,(256,)) - Rinner_new= CartesianRange(CartesianIndex((first(Rinner).I[3:end])),CartesianIndex((last(Rinner).I[3:end]))) if !isempty(Rinner) out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) - for I in Rinner_new - indice=I.I - prev_mode= 0 - # Handle the interior - m_histogram=zeros(Int64,(256,)) - for i=L1:L2 - if (i-L1)%2==0 - temp1=(i,R1) - temp2=(i,R2) - id1=(temp1...,indice...) - id2=(temp2...,indice...) - RinnerN= CartesianRange(CartesianIndex(id1), CartesianIndex(id2)) - for I in RinnerN - - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - if (I[2]== R2 && prev_mode==1) - out[I] = f(bufrs,m_histogram,3,window) - prev_mode=3 - continue - elseif prev_mode == 4 - out[I] = f(bufrs,m_histogram,6,window) - prev_mode = 1 - continue - - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=1 - continue - elseif prev_mode == 1 - out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 - continue - end - - end - - else - - for k= R2:-1:R1 - I=CartesianIndex(tuple((i,k)...,indice...)) - Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) - if (I[2]==R1 && prev_mode==2) - out[I] = f(bufrs,m_histogram,4,window) - prev_mode=4 - continue - elseif prev_mode == 3 - out[I] = f(bufrs,m_histogram,5,window) - prev_mode = 2 - continue - elseif prev_mode == 0 - out[I] = f(bufrs,m_histogram,0,window) - prev_mode=2 - continue - elseif prev_mode == 2 - out[I] = f(bufrs,m_histogram,2,window) - prev_mode=2 - continue - end - end - - end + 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(curr_col - prev_col!=0) + m_histogram=zeros(Int64,(256,)) + prev_mode=0 end - end - + prev_col=curr_col + + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + if prev_mode == 0 + out[I] = f(bufrs,m_histogram,0,window) + prev_mode=1 + continue + 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) @@ -227,7 +181,6 @@ function median_filter{T,N}(f, mode =-1 copy_win!(buf, img, I, border, offset) out[I] = f(bufrs,m_histogram,mode,window) - #out[I]=0 end else @@ -413,129 +366,67 @@ replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) inds = indices(v,1) return Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering())) else - # reshape the value of the window so that they are horizontal major which allows to do horizontal traversal window_size=size(v,1) dims = map(x->x.stop-x.start+1,window) - startp= ones(Int64,(size(dims,1),)) - v_blob=reshape(v,dims) - blob_range=CartesianRange(CartesianIndex(tuple(startp...)),CartesianIndex(tuple(dims...))) - blob_range_new=CartesianRange(CartesianIndex(first(blob_range).I[3:end]),CartesianIndex(last(blob_range).I[3:end])) - # Update the histogram according to - for I in blob_range_new - indice=I.I - v=v_blob[:,:,indice...] + 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) + for j= i: dims[1]*dims[2]:last(inds) + id= trunc(Int64,(v[j]*255))+1 + m_histogram[id]+=1 + end + end + # Compute the median + tempsum = 0 m_index=-1 - inds = indices(v,1) - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) - v= reshape(p,(size(v,1)*size(v,2),)) + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= window_size/2 + m_index=i-1 + break + end + end + # Clear the histogram from previous value + for i = first(inds):width:last(inds) + for j= i: dims[1]*dims[2]:last(inds) + id= trunc(Int64,(v[j]*255))+1 + m_histogram[id]-=1 + end + end + return convert(Float64,m_index)/255 + + elseif mode == 1 + # Update the histogram according to new entries + for i = width:width:last(inds) + for j= i: dims[1]*dims[2]:last(inds) + id= trunc(Int64,(v[j]*255))+1 + m_histogram[id]+=1 + end + end + # Compute the median + tempsum = 0 m_index=-1 - inds = indices(v,1) - if mode == 0 - for i = first(inds):last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 1 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 2 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode == 3 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==4 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==5 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end - elseif mode ==6 - for i = last(inds)-width+1:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]+=1 - end + for i = 1:256 + tempsum+= m_histogram[i] + if tempsum>= window_size/2 + m_index=i-1 + break end - - end - # Find the median using histogram - tempsum = 0 - m_index=-1 - for i = 1:256 - tempsum+= m_histogram[i] - if tempsum>= window_size/2 - m_index=i-1 - break end - end - - # Clear histogram from useless past items - for I in blob_range_new - indice=I.I - v=v_blob[:,:,indice...] - - width=window[2].stop-window[2].start+1 - height= window[1].stop-window[1].start+1 - v_reshape= reshape(v,(height,width)) - p=zeros(eltype(v_reshape), (size(v_reshape,2),size(v_reshape,1))) - transpose!(p,v_reshape) - v= reshape(p,(size(v,1)*size(v,2),)) - - inds = indices(v,1) - if mode == 0 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 1 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 2 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode == 3 - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==4 - for i = first(inds):first(inds)+width-1 - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==5 - for i = width:width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end - elseif mode ==6 - for i = first(inds):width:last(inds) - id= trunc(Int64,(v[i]*255))+1 - m_histogram[id]-=1 - end + # Clear the histogram from previous value + for i = first(inds):width:last(inds) + for j= i: dims[1]*dims[2]:last(inds) + id= trunc(Int64,(v[j]*255))+1 + m_histogram[id]-=1 end - + end + return convert(Float64,m_index)/255 end - return convert(Float64,m_index)/255 + end + end From e1735ff22ccca56346b6e04a036974fbd02a3f42 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Mon, 10 Jul 2017 12:23:28 +0200 Subject: [PATCH 06/11] Clear Definitions --- src/mapwindow.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index aca7951..a62bb5e 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -75,7 +75,7 @@ function mapwindow{T,N}(f, window::Indices{N}, border::BorderSpecAny; callmode=:copy!) - if(uses_histogram(f)) + 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) @@ -125,6 +125,8 @@ 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}, @@ -151,7 +153,7 @@ function median_filter{T,N}(f, prev_col=1 for I in Rinner curr_col=I.I[2] - if(curr_col - prev_col!=0) + if(column_change(curr_col,prev_col)) m_histogram=zeros(Int64,(256,)) prev_mode=0 end @@ -159,10 +161,12 @@ function median_filter{T,N}(f, 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 @@ -190,11 +194,15 @@ function median_filter{T,N}(f, out end -function uses_histogram(::typeof(median!)) +function column_change(curr_col,prev_col) + return (curr_col - prev_col!=0) +end + +function is_medianfilter(::typeof(median!)) true end -function uses_histogram(::Any) +function is_medianfilter(::Any) false end From 5ecf8ef2158be6baf076bcfabdec5aae83ac3619 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Mon, 10 Jul 2017 13:07:29 +0200 Subject: [PATCH 07/11] changed median statistic --- src/mapwindow.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index a62bb5e..b8ddae7 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -391,7 +391,7 @@ replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) m_index=-1 for i = 1:256 tempsum+= m_histogram[i] - if tempsum>= window_size/2 + if tempsum>=trunc(Int64,window_size/2)+1 m_index=i-1 break end @@ -418,7 +418,7 @@ replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) m_index=-1 for i = 1:256 tempsum+= m_histogram[i] - if tempsum>= window_size/2 + if tempsum>= trunc(Int64,window_size/2)+1 m_index=i-1 break end From 99da897c424b9a855301f2a8a9daa86a2a97f531 Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Mon, 10 Jul 2017 13:35:28 +0200 Subject: [PATCH 08/11] changed median statistic --- src/mapwindow.jl | 36 ++++++++++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 8 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index b8ddae7..f080b78 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -374,6 +374,7 @@ replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) 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 @@ -381,10 +382,13 @@ replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) if mode == 0 # Update the histogram according to new entries for i = first(inds):last(inds) - for j= i: dims[1]*dims[2]:last(inds) - id= trunc(Int64,(v[j]*255))+1 - m_histogram[id]+=1 - end + 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 @@ -397,22 +401,34 @@ replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) end end # Clear the histogram from previous value - for i = first(inds):width:last(inds) - for j= i: dims[1]*dims[2]:last(inds) + 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:last(inds) + 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 @@ -424,10 +440,14 @@ replace_function(::typeof(median!)) = function(v,m_histogram,mode,window) end end # Clear the histogram from previous value - for i = first(inds):width:last(inds) + 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 From 56c4457123d7fc8a0da1aeedb62c8ec9f9b64ebe Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Mon, 24 Jul 2017 02:10:45 +0530 Subject: [PATCH 09/11] Update program structure --- src/ImageFiltering.jl | 4 +- src/mapwindow.jl | 317 ++++++++++++++++++++++-------------------- 2 files changed, 167 insertions(+), 154 deletions(-) diff --git a/src/ImageFiltering.jl b/src/ImageFiltering.jl index aac8f34..10c6572 100644 --- a/src/ImageFiltering.jl +++ b/src/ImageFiltering.jl @@ -7,7 +7,7 @@ using ColorVectorSpace # for filtering RGB arrays using Compat using Base: Indices, tail, fill_to_length, @pure, depwarn -export Kernel, KernelFactors, Pad, Fill, Inner, NA, NoPad, Algorithm, imfilter, imfilter!, mapwindow, imgradients, padarray, centered, kernelfactors, reflect +export Kernel, KernelFactors, Pad, Fill, Inner, NA, NoPad, Algorithm, imfilter, imfilter!, median_direct!, median_histogram!, mapwindow, imgradients, padarray, centered, kernelfactors, reflect @compat FixedColorant{T<:Normed} = Colorant{T} @compat StaticOffsetArray{T,N,A<:StaticArray} = OffsetArray{T,N,A} @@ -68,8 +68,10 @@ include("imfilter.jl") include("specialty.jl") include("mapwindow.jl") + using .MapWindow + function __init__() # See ComputationalResources README for explanation push!(LOAD_PATH, dirname(@__FILE__)) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index f080b78..11acd81 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -3,8 +3,11 @@ module MapWindow using DataStructures, TiledIteration using ..ImageFiltering: BorderSpecAny, Pad, Fill, borderinstance, _interior, padindex, imfilter using Base: Indices, tail +using FixedPointNumbers,Colors export mapwindow +export median_direct! +export median_histogram! """ mapwindow(f, img, window, [border="replicate"]) -> imgf @@ -41,6 +44,82 @@ and then `mapwindow(f, img, (m,n))` should filter at the 75th quantile. See also: [`imfilter`](@ref). """ +function median_direct!(v::AbstractVector) + inds = indices(v,1) + Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering())) +end + +function median_histogram!(v::AbstractVector,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->last(x)-first(x)+1,window) + width=last(window[1])-first(window[1])+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 + # Compute the median + tempsum = 0 + m_index=-1 + for i in linearindices(m_histogram) + 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 + 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 + + # Compute the median + tempsum = 0 + m_index=-1 + for i in linearindices(m_histogram) + 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 + end + end + return convert(Float64,m_index)/255 + end + + end + +end + + function mapwindow(f, img::AbstractArray, window::Dims, args...; kwargs...) all(isodd(w) for w in window) || error("entries in window must be odd, got $window") halfsize = map(w->w>>1, window) @@ -70,64 +149,105 @@ mapwindow(f, img, window::AbstractArray, args...; kwargs...) = mapwindow(f, img, -function mapwindow{T,N}(f, - img::AbstractArray{T,N}, - window::Indices{N}, - border::BorderSpecAny; - callmode=:copy!) - 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}, +# This is a implementation of Fast 2D median filter +# http://ieeexplore.ieee.org/document/1163188/ + +function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, + img::Union{AbstractArray{N0f8,N},AbstractArray{ColorTypes.Gray{N0f8},N}}, window::Indices{N}, border::BorderSpecAny, shape=default_shape(f); callmode=:copy!) + f = replace_function(f,window) + img= map(x->Float64(x),img) + # This ensures different method of image traversal for median_direct + if(typeof(f)==typeof(median_direct!)) + return _mapwindow(median_direct!, img, window, border, default_shape(f); callmode=callmode) + end inds = indices(img) inner = _interior(inds, window) if callmode == :copy! - buf = Array{T}(map(length, window)) + buf = Array{eltype(img)}(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) - out = similar(img, typeof(f(bufrs))) - # Handle the interior + prev_mode=0 + prev_col=1 for I in Rinner - Rwin = CartesianRange(map(+, window, I.I)) + 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) - out[I] = f(bufrs) + # Mode 0 corresponds to refilling the empty histogram with all the points in the window + if prev_mode == 0 + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + 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 + Rwin = CartesianRange(map(+, window, I.I)) + copy!(buf, Rbuf, img, Rwin) + 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))) + 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) + out[I] = f(bufrs,m_histogram,mode,window) end + else - # TODO: implement :view - error("callmode $callmode not supported") + # TODO: implement :view + error("callmode $callmode not supported") end out end +function mapwindow{T,N}(f::typeof(median!), + img::AbstractArray{T,N}, + window::Indices{N}, + border::BorderSpecAny; + callmode=:copy!) + _mapwindow(median_direct!, img, window, border, default_shape(f); callmode=callmode) + +end +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) + +end -# This is a implementation of Fast 2D median filter -# http://ieeexplore.ieee.org/document/1163188/ -function median_filter{T,N}(f, +function _mapwindow{T,N}(f, img::AbstractArray{T,N}, window::Indices{N}, border::BorderSpecAny, @@ -142,69 +262,39 @@ function median_filter{T,N}(f, 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 + out = similar(img, typeof(f(bufrs))) + # Handle the interior 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)) + 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 - + out[I] = f(bufrs) end - else copy_win!(buf, img, first(CartesianRange(inds)), border, offset) - out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) + out = similar(img, typeof(f(bufrs))) 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) + out[I] = f(bufrs) end - else - # TODO: implement :view - error("callmode $callmode not supported") + # 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) @@ -366,99 +456,20 @@ end @inline cyclecache(b, x) = b[1], (Base.tail(b)..., x) +replace_function(f,window) = replace_function(f) replace_function(f) = f -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 - +function replace_function(::typeof(median!), window) + # The threshold is defined for changing the method of finding median [Using sorting or Using Histogram] + if prod(map(length, window)) <= 81 + return median_direct! end - - + median_histogram! end + default_shape(::Any) = identity -default_shape(::typeof(median!)) = vec +default_shape(::Union{typeof(median!),typeof(median_direct!),typeof(median_histogram!)}) = vec + end \ No newline at end of file From 6c9bf747af784e3a7e21cb280699446e344c5c9f Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Fri, 28 Jul 2017 00:58:43 +0530 Subject: [PATCH 10/11] Separate dispatch for continuos and discrete input --- src/mapwindow.jl | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index 11acd81..a35553e 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -44,11 +44,16 @@ and then `mapwindow(f, img, (m,n))` should filter at the 75th quantile. See also: [`imfilter`](@ref). """ + + function median_direct!(v::AbstractVector) inds = indices(v,1) Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering())) end +# This is a implementation of Fast 2D median filter +# http://ieeexplore.ieee.org/document/1163188/ + function median_histogram!(v::AbstractVector,m_histogram,mode,window) # Handle the boundary points with mode -1 @@ -150,8 +155,8 @@ mapwindow(f, img, window::AbstractArray, args...; kwargs...) = mapwindow(f, img, -# This is a implementation of Fast 2D median filter -# http://ieeexplore.ieee.org/document/1163188/ + +# This dispatch takes input of type Union{AbstractArray{N0f8,N},AbstractArray{ColorTypes.Gray{N0f8},N}} and decides whether to use median_histogram or median_direct method function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, img::Union{AbstractArray{N0f8,N},AbstractArray{ColorTypes.Gray{N0f8},N}}, @@ -160,15 +165,15 @@ function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, shape=default_shape(f); callmode=:copy!) f = replace_function(f,window) - img= map(x->Float64(x),img) + img_f= map(x->Float64(x),img) # This ensures different method of image traversal for median_direct if(typeof(f)==typeof(median_direct!)) - return _mapwindow(median_direct!, img, window, border, default_shape(f); callmode=callmode) + return _mapwindow(median_direct!, img_f, window, border, default_shape(f); callmode=callmode) end - inds = indices(img) + inds = indices(img_f) inner = _interior(inds, window) if callmode == :copy! - buf = Array{eltype(img)}(map(length, window)) + buf = Array{Float64}(map(length, window)) bufrs = shape(buf) Rbuf = CartesianRange(size(buf)) offset = CartesianIndex(map(w->first(w)-1, window)) @@ -178,9 +183,9 @@ function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, mode = 0 m_histogram=zeros(Int64,(256,)) if !isempty(Rinner) - out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) + out = similar(img_f, typeof(f(bufrs,m_histogram,mode,window))) Rwin = CartesianRange(map(+, window, first(Rinner).I)) - copy!(buf, Rbuf, img, Rwin) + copy!(buf, Rbuf, img_f, Rwin) prev_mode=0 prev_col=1 for I in Rinner @@ -191,18 +196,18 @@ function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, end prev_col=curr_col Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) + copy!(buf, Rbuf, img_f, Rwin) # Mode 0 corresponds to refilling the empty histogram with all the points in the window if prev_mode == 0 Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) + copy!(buf, Rbuf, img_f, Rwin) 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 Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img, Rwin) + copy!(buf, Rbuf, img_f, Rwin) out[I] = f(bufrs,m_histogram,1,window) prev_mode=1 continue @@ -211,14 +216,14 @@ function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, end else - copy_win!(buf, img, first(CartesianRange(inds)), border, offset) - out = similar(img, typeof(f(bufrs,m_histogram,mode,window))) + copy_win!(buf, img_f, first(CartesianRange(inds)), border, offset) + out = similar(img_f, 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 + # Handle the edge points with mode -1 [finding median using simple sort] mode =-1 - copy_win!(buf, img, I, border, offset) + copy_win!(buf, img_f, I, border, offset) out[I] = f(bufrs,m_histogram,mode,window) end @@ -226,6 +231,8 @@ function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, # TODO: implement :view error("callmode $callmode not supported") end + + out = map(x->convert(eltype(img),x),out) out end From 7e80f2630becf3b43731662afc1ea25133c603ae Mon Sep 17 00:00:00 2001 From: hari-sikchi Date: Sat, 5 Aug 2017 20:25:18 +0530 Subject: [PATCH 11/11] removed type incompatibility --- src/mapwindow.jl | 86 +++++++++++++++++++++++------------------------- 1 file changed, 42 insertions(+), 44 deletions(-) diff --git a/src/mapwindow.jl b/src/mapwindow.jl index a35553e..57be3d1 100644 --- a/src/mapwindow.jl +++ b/src/mapwindow.jl @@ -4,6 +4,7 @@ using DataStructures, TiledIteration using ..ImageFiltering: BorderSpecAny, Pad, Fill, borderinstance, _interior, padindex, imfilter using Base: Indices, tail using FixedPointNumbers,Colors +using ImageCore:normedview export mapwindow export median_direct! @@ -51,30 +52,29 @@ function median_direct!(v::AbstractVector) Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering())) end -# This is a implementation of Fast 2D median filter +# This is a implementation of Fast 2D median filter extended for nD # http://ieeexplore.ieee.org/document/1163188/ -function median_histogram!(v::AbstractVector,m_histogram,mode,window) +function median_histogram!(v::AbstractVector{UInt8},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())) + return convert(eltype(v),Base.middle(Base.select!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering()))) else - window_size=size(v,1) dims = map(x->last(x)-first(x)+1,window) - width=last(window[1])-first(window[1])+1 + width = last(window[1])-first(window[1])+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 + id = v[i]+1; + m_histogram[id]+= 1 end # Compute the median tempsum = 0 - m_index=-1 + m_index = -1 for i in linearindices(m_histogram) tempsum+= m_histogram[i] if tempsum>=trunc(Int64,window_size/2)+1 @@ -84,42 +84,39 @@ function median_histogram!(v::AbstractVector,m_histogram,mode,window) 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 + for j = i: dims[1]*dims[2]:last(inds) + id = v[j]+1 + m_histogram[id]-= 1 end end - return convert(Float64,m_index)/255 - + return convert(UInt8,m_index) 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 + id = v[j]+1 + m_histogram[id]+= 1 end end - # Compute the median tempsum = 0 m_index=-1 for i in linearindices(m_histogram) tempsum+= m_histogram[i] if tempsum>= trunc(Int64,window_size/2)+1 - m_index=i-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 + for j = i: dims[1]*dims[2]:last(inds) + id = v[j]+1 + m_histogram[id]-= 1 end end - return convert(Float64,m_index)/255 + return convert(UInt8,m_index) end - end end @@ -154,8 +151,6 @@ mapwindow(f, img, window::AbstractArray, args...; kwargs...) = mapwindow(f, img, - - # This dispatch takes input of type Union{AbstractArray{N0f8,N},AbstractArray{ColorTypes.Gray{N0f8},N}} and decides whether to use median_histogram or median_direct method function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, @@ -165,15 +160,18 @@ function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, shape=default_shape(f); callmode=:copy!) f = replace_function(f,window) - img_f= map(x->Float64(x),img) + img_n=map(x->gray(x),img) + img_n=reinterpret.(img_n) # This ensures different method of image traversal for median_direct - if(typeof(f)==typeof(median_direct!)) - return _mapwindow(median_direct!, img_f, window, border, default_shape(f); callmode=callmode) + if(typeof(f) == typeof(median_direct!)) + out = _mapwindow(median_direct!, img_n, window, border, default_shape(f); callmode=callmode) + out = normedview(map(x->convert(UInt8,x),out)) + return map(x->convert(eltype(img),x),out) end - inds = indices(img_f) + inds = indices(img_n) inner = _interior(inds, window) if callmode == :copy! - buf = Array{Float64}(map(length, window)) + buf = Array{UInt8}(map(length, window)) bufrs = shape(buf) Rbuf = CartesianRange(size(buf)) offset = CartesianIndex(map(w->first(w)-1, window)) @@ -183,55 +181,55 @@ function mapwindow{N}(f::Union{typeof(median!),typeof(median_histogram!)}, mode = 0 m_histogram=zeros(Int64,(256,)) if !isempty(Rinner) - out = similar(img_f, typeof(f(bufrs,m_histogram,mode,window))) + out = similar(img_n, typeof(f(bufrs,m_histogram,mode,window))) Rwin = CartesianRange(map(+, window, first(Rinner).I)) - copy!(buf, Rbuf, img_f, Rwin) - prev_mode=0 - prev_col=1 + copy!(buf, Rbuf, img_n, 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 + prev_col = curr_col Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img_f, Rwin) + copy!(buf, Rbuf, img_n, Rwin) # Mode 0 corresponds to refilling the empty histogram with all the points in the window if prev_mode == 0 Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img_f, Rwin) + copy!(buf, Rbuf, img_n, Rwin) 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 Rwin = CartesianRange(map(+, window, I.I)) - copy!(buf, Rbuf, img_f, Rwin) + copy!(buf, Rbuf, img_n, Rwin) out[I] = f(bufrs,m_histogram,1,window) - prev_mode=1 + prev_mode=1 continue end end else - copy_win!(buf, img_f, first(CartesianRange(inds)), border, offset) - out = similar(img_f, typeof(f(bufrs,m_histogram,mode,window))) + copy_win!(buf, img_n, first(CartesianRange(inds)), border, offset) + out = similar(img_n, 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 [finding median using simple sort] mode =-1 - copy_win!(buf, img_f, I, border, offset) + copy_win!(buf, img_n, I, border, offset) out[I] = f(bufrs,m_histogram,mode,window) end else - # TODO: implement :view - error("callmode $callmode not supported") + # TODO: implement :view + error("callmode $callmode not supported") end - + out = normedview(out) out = map(x->convert(eltype(img),x),out) out end