From 3f832dbb918c0e6f99574370577670c9b38ec80c Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 20 Aug 2016 07:48:16 -0500 Subject: [PATCH 01/21] A few fixes and performance enhancements --- src/ImagesFiltering.jl | 2 +- src/imfilter.jl | 21 +++++++++++++++------ src/specialty.jl | 15 ++++++--------- test/specialty.jl | 2 ++ 4 files changed, 24 insertions(+), 16 deletions(-) diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index 316283d..824501a 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -2,7 +2,7 @@ module ImagesFiltering using Colors, FixedPointNumbers, ImagesCore, MappedArrays, FFTViews, OffsetArrays, StaticArrays, ComputationalResources using ColorVectorSpace # in case someone filters RGB arrays -using Base: Indices, tail, fill_to_length +using Base: Indices, tail, fill_to_length, @pure export Kernel, KernelFactors, Pad, Fill, Inner, Algorithm, imfilter, imfilter!, padarray, centered diff --git a/src/imfilter.jl b/src/imfilter.jl index a430142..2922ded 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -120,7 +120,7 @@ imfilter # have to be a little more cautious to make sure that later methods # don't inadvertently call back to these: in methods that take an # AbstractResource argument, exclude `NoPad()` as a border option. -function imfilter!(out::AbstractArray, img::AbstractArray, kernel::AbstractArray, args...) +function imfilter!(out::AbstractArray, img::AbstractArray, kernel::Union{AbstractArray,Laplacian}, args...) imfilter!(out, img, factorkernel(kernel), args...) end @@ -321,11 +321,11 @@ function _imfilter_inbounds!(out, A::AbstractArray, kern::OffsetArray, R) ΔI = CartesianIndex(kern.offsets) _imfilter_inbounds!(out, (A, ΔI), parent(kern), R) end -function _imfilter_inbounds!(out, A::OffsetArray, kern::AbstractArray, R) - ΔI = CartesianIndex(kern.offsets) +function _imfilter_inbounds!(out, A::OffsetArray, kern, R) + ΔI = -CartesianIndex(A.offsets) _imfilter_inbounds!(out, (parent(A), ΔI), kern, R) end -function _imfilter_inbounds!(out, A::AbstractArray, kern::AbstractArray, R) +function _imfilter_inbounds!(out, A::AbstractArray, kern, R) indsk = indices(kern) p = first(A) * first(kern) TT = typeof(p+p) @@ -640,14 +640,23 @@ end ### Utilities filter_type{S,T}(img::AbstractArray{S}, kernel::AbstractArray{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) -filter_type{S,T}(img::AbstractArray{S}, kernel::Tuple{AbstractArray{T},Vararg{AbstractArray{T}}}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) -filter_type{S,T}(img::AbstractArray{S}, kernel::Tuple{IIRFilter{T},Vararg{IIRFilter{T}}}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) +filter_type{S,T}(img::AbstractArray{S}, kernel::IIRFilter{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) filter_type{S<:Union{UFixed,FixedColorant}}(img::AbstractArray{S}, ::Laplacian) = float32(S) filter_type{S<:Colorant}(img::AbstractArray{S}, ::Laplacian) = S filter_type{S<:AbstractFloat}(img::AbstractArray{S}, ::Laplacian) = S filter_type{S<:Signed}(img::AbstractArray{S}, ::Laplacian) = S filter_type{S<:Unsigned}(img::AbstractArray{S}, ::Laplacian) = signed(S) +@inline function filter_type(img, kernel::Tuple{Any,Vararg{Any}}) + T = filter_type(img, kernel[1]) + filter_type(T, img, tail(kernel)) +end +@inline function filter_type{T}(::Type{T}, img, kernel::Tuple{Any,Vararg{Any}}) + S = promote_type(T, filter_type(img, kernel[1])) + filter_type(S, img, tail(kernel)) +end +filter_type{T}(::Type{T}, img, kernel::Tuple{}) = T + factorkernel(kernel::AbstractArray) = (kernelshift(indices(kernel), kernel),) factorkernel(L::Laplacian) = (L,) diff --git a/src/specialty.jl b/src/specialty.jl index 96d6772..6377067 100644 --- a/src/specialty.jl +++ b/src/specialty.jl @@ -1,12 +1,9 @@ ## Laplacian -# function !{S,T,N}(r::AbstractResource, -# out::AbstractArray{S,N}, -# A::AbstractArray{T,N}, -# L::Tuple{Laplacian{N}}) -# _imfilter_padded!(r, out, A, L[1]) -# end - +function _imfilter_inbounds!{N}(out, A::OffsetArray, kern::Laplacian{N}, R) + ΔI = -CartesianIndex(A.offsets) + _imfilter_inbounds!(out, (parent(A), ΔI), kern, R) +end function _imfilter_inbounds!(out, A::AbstractArray, L::Laplacian, R) TT = eltype(out) # accumtype(eltype(out), eltype(A)) n = 2*length(L.offsets) @@ -22,11 +19,11 @@ function _imfilter_inbounds!(out, A::AbstractArray, L::Laplacian, R) end function _imfilter_inbounds!(out, Ashift::Tuple{AbstractArray,CartesianIndex}, L::Laplacian, R) A, ΔI = Ashift - TT = accumtype(eltype(out), eltype(A)) + TT = eltype(out) # accumtype(eltype(out), eltype(A)) n = 2*length(L.offsets) for I in R Ishift = I + ΔI - @inbounds tmp = convert(TT, - n * A[I]) + @inbounds tmp = convert(TT, - n * A[Ishift]) @inbounds for J in L.offsets tmp += A[Ishift+J] tmp += A[Ishift-J] diff --git a/test/specialty.jl b/test/specialty.jl index 00b1bc6..7998c3c 100644 --- a/test/specialty.jl +++ b/test/specialty.jl @@ -17,6 +17,8 @@ using Base.Test af = imfilter(a, kern) T = eltype(a) @test af == [zero(T),a[3],-2a[3],a[3],zero(T)] + af = imfilter(a, (kern,)) + @test af == [zero(T),a[3],-2a[3],a[3],zero(T)] end for a in (makeimpulse(Float64, (5,), 1), makeimpulse(Gray{U8}, (5,), 1), From 932fc27c82f2dadb41afe22b19a92c3db5d56086 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 20 Aug 2016 07:48:43 -0500 Subject: [PATCH 02/21] Add sobel and prewitt to Kernel, and improve docs --- src/kernel.jl | 77 +++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 63 insertions(+), 14 deletions(-) diff --git a/src/kernel.jl b/src/kernel.jl index 037c8ea..0304e6e 100644 --- a/src/kernel.jl +++ b/src/kernel.jl @@ -4,22 +4,63 @@ using StaticArrays, OffsetArrays using ..ImagesFiltering: centered, KernelFactors import ..ImagesFiltering: _reshape +function product2d(kf) + k1, k2 = kf + k1[1].*k1[2], k2[1].*k2[2] +end + +""" + diff1, diff2 = sobel() + +Return kernels for two-dimensional gradient compution using the Sobel +operator. `diff1` computes the gradient along the first (y) dimension, +and `diff2` computes the gradient along the second (x) dimension. + +See also: KernelFactors.sobel, Kernel.prewitt, Kernel.ando. +""" +sobel() = product2d(KernelFactors.sobel()) + +""" + diff1, diff2 = prewitt() + +Return kernels for two-dimensional gradient compution using the +Prewitt operator. `diff1` computes the gradient along the first (y) +dimension, and `diff2` computes the gradient along the second (x) +dimension. + +See also: KernelFactors.prewitt, Kernel.sobel, Kernel.ando. """ -`kern1, kern2 = ando3()` returns optimal 3x3 gradient filters for dimensions 1 and 2 of your image, as defined in -Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March 2000. +prewitt() = product2d(KernelFactors.prewitt()) + +""" + diff1, diff2 = ando3() + +Return 3x3 kernels for two-dimensional gradient compution using the +optimal "Ando" filters. `diff1` computes the gradient along the +y-axis (first dimension), and `diff2` computes the gradient along the +x-axis (second dimension). + +# Citation +Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March +2000 See also: KernelFactors.ando3, Kernel.ando4, Kernel.ando5. """ -function ando3() - k1, k2 = KernelFactors.ando3() - k1[1].*k1[2], k2[1].*k2[2] -end +ando3() = product2d(KernelFactors.ando3()) """ -`kern1, kern2 = ando4()` returns optimal 4x4 gradient filters for dimensions 1 and 2 of your image, as defined in -Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March 2000. + diff1, diff2 = ando4() + +Return 4x4 kernels for two-dimensional gradient compution using the +optimal "Ando" filters. `diff1` computes the gradient along the +y-axis (first dimension), and `diff2` computes the gradient along the +x-axis (second dimension). + +# Citation +Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March +2000 -See also: `KernelFactors.ando4`, `Kernel.ando3`, `Kernel.ando5`. +See also: KernelFactors.ando4, Kernel.ando3, Kernel.ando5. """ function ando4() f = centered(@SMatrix [ -0.022116 -0.025526 0.025526 0.022116 @@ -30,8 +71,16 @@ function ando4() end """ -`kern1, kern2 = ando5()` returns optimal 5x5 gradient filters for dimensions 1 and 2 of your image, as defined in -Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March 2000. + diff1, diff2 = ando5() + +Return 5x5 kernels for two-dimensional gradient compution using the +optimal "Ando" filters. `diff1` computes the gradient along the +y-axis (first dimension), and `diff2` computes the gradient along the +x-axis (second dimension). + +# Citation +Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March +2000 See also: `KernelFactors.ando5`, `Kernel.ando3`, `Kernel.ando4`. """ @@ -69,9 +118,9 @@ gaussian(σ::Tuple{}) = reshape([1]) gaussian(σ::Real) = gaussian((σ, σ)) """ - DoG((σp1, σp2, ...), (σm1, σm2, ...), [l]) -> k - DoG((σ1, σ2, ...)) -> k - DoG(σ::Real) -> k + DoG((σp1, σp2, ...), (σm1, σm2, ...), [l1, l2, ...]) -> k + DoG((σ1, σ2, ...)) -> k + DoG(σ::Real) -> k Construct a multidimensional difference-of-gaussian kernel `k`, equal to `gaussian(σp, l)-gaussian(σm, l)`. When only a single `σ` is From 0ad7e32467980ad3017aac508ad47d0a53a5c181 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 20 Aug 2016 14:10:07 -0500 Subject: [PATCH 03/21] More bugfixes and performance improvements --- src/imfilter.jl | 11 +++++++++-- src/utils.jl | 2 +- test/2d.jl | 5 +++-- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/imfilter.jl b/src/imfilter.jl index 2922ded..494ddf6 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -251,7 +251,13 @@ function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, ::NoPad) copy!(out, R, A1, R) end -function _imfilter!(r, out, A1, A2, kernel::Tuple, ::NoPad) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, ::NoPad) + kern = kernel[1] + iscopy(kern) && return imfilter!(r, out, A, (), NoPad()) + imfilter!(r, out, A1, samedims(out, kern), NoPad()) +end + +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, ::NoPad) kern = kernel[1] iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), NoPad()) kernN = samedims(A2, kern) @@ -343,10 +349,11 @@ function _imfilter_inbounds!(out, Ashift::Tuple{AbstractArray,CartesianIndex}, k indsk = indices(kern) p = first(A) * first(kern) TT = typeof(p+p) + Rk = CartesianRange(indsk) for I in R Ishift = I + ΔI tmp = zero(TT) - @inbounds for J in CartesianRange(indsk) + @inbounds for J in Rk tmp += A[Ishift+J]*kern[J] end @inbounds out[I] = tmp diff --git a/src/utils.jl b/src/utils.jl index 78ed1a6..647d0c7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -35,7 +35,7 @@ function _vec(a::OffsetArray) if i == 0 i = 1 end - OffsetArray(vec(parent(a)), inds[1]) + OffsetArray(vec(parent(a)), inds[i]) end samedims{N}(::Type{Val{N}}, kernel) = _reshape(kernel, Val{N}) diff --git a/test/2d.jl b/test/2d.jl index a6f6eaf..b57f87d 100644 --- a/test/2d.jl +++ b/test/2d.jl @@ -28,7 +28,7 @@ using Base.Test @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32(targetimg) end end - targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:3) + targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:5) @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner @test @inferred(imfilter(f32type(img), img, kernel, Inner())) ≈ float32(targetimg_inner) for alg in (Algorithm.FIR(), Algorithm.FFT()) @@ -52,7 +52,8 @@ using Base.Test @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32(targetimg) end end - targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:3) + targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:5) + @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner @test @inferred(imfilter(f32type(img), img, kernel, Inner())) ≈ float32(targetimg_inner) for alg in (Algorithm.FIR(), Algorithm.FFT()) From 1920ced26f6eb75db7562aa6eeec80736fd82839 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 20 Aug 2016 16:20:58 -0500 Subject: [PATCH 04/21] Test kernels with Rational coefficients A good test for the eltype generality of our operations --- src/imfilter.jl | 7 ++++--- test/2d.jl | 33 ++++++++++++++++++++++++++++++++- 2 files changed, 36 insertions(+), 4 deletions(-) diff --git a/src/imfilter.jl b/src/imfilter.jl index 494ddf6..c9aed28 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -333,7 +333,8 @@ function _imfilter_inbounds!(out, A::OffsetArray, kern, R) end function _imfilter_inbounds!(out, A::AbstractArray, kern, R) indsk = indices(kern) - p = first(A) * first(kern) + Rk = CartesianRange(indsk) + p = A[first(R)+first(Rk)] * first(kern) TT = typeof(p+p) for I in R tmp = zero(TT) @@ -347,9 +348,9 @@ end function _imfilter_inbounds!(out, Ashift::Tuple{AbstractArray,CartesianIndex}, kern, R) A, ΔI = Ashift indsk = indices(kern) - p = first(A) * first(kern) - TT = typeof(p+p) Rk = CartesianRange(indsk) + p = A[first(R)+first(Rk)+ΔI] * first(kern) + TT = typeof(p+p) for I in R Ishift = I + ΔI tmp = zero(TT) diff --git a/test/2d.jl b/test/2d.jl index b57f87d..2833ea8 100644 --- a/test/2d.jl +++ b/test/2d.jl @@ -54,13 +54,44 @@ using Base.Test end targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:5) @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner - @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner @test @inferred(imfilter(f32type(img), img, kernel, Inner())) ≈ float32(targetimg_inner) for alg in (Algorithm.FIR(), Algorithm.FFT()) @test @inferred(imfilter(img, kernel, Inner(), alg)) ≈ targetimg_inner @test @inferred(imfilter(f32type(img), img, kernel, Inner(), alg)) ≈ float32(targetimg_inner) end end + # Rational filter coefficients + kernel = (centered([1//3, 1//3, 1//3]), centered([1//3, 1//3, 1//3]')) + kern = kernel[1].*kernel[2] + for img in (imgf, imgi, imgg, imgc) + targetimg = zeros(typeof(img[1]*kern[1]), size(img)) + targetimg[2:4,3:5] = img[3,4]*(1//9) + @test @inferred(imfilter(img, kernel)) ≈ targetimg + @test @inferred(imfilter(f32type(img), img, kernel)) ≈ float32(targetimg) + for border in (Pad{:replicate}(), Pad{:circular}(), Pad{:symmetric}(), Pad{:reflect}(), Fill(zero(eltype(img)))) + @test @inferred(imfilter(img, kernel, border)) ≈ targetimg + @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32(targetimg) + for alg in (Algorithm.FIR(), Algorithm.FFT()) + if alg == Algorithm.FFT() && eltype(img) == Int + @test @inferred(imfilter(Float64, img, kernel, border, alg)) ≈ targetimg + else + @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg + end + @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32(targetimg) + end + end + targetimg_inner = OffsetArray(targetimg[2:end-1, 2:end-1], 2:4, 2:6) + @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner + @test @inferred(imfilter(f32type(img), img, kernel, Inner())) ≈ float32(targetimg_inner) + for alg in (Algorithm.FIR(), Algorithm.FFT()) + if alg == Algorithm.FFT() && eltype(img) == Int + @test @inferred(imfilter(Float64, img, kernel, Inner(), alg)) ≈ targetimg_inner + else + @test @inferred(imfilter(img, kernel, Inner(), alg)) ≈ targetimg_inner + end + @test @inferred(imfilter(f32type(img), img, kernel, Inner(), alg)) ≈ float32(targetimg_inner) + end + end ## Images for which the boundary conditions matter imgf = zeros(5, 7); imgf[1,2] = 1 imgi = zeros(Int, 5, 7); imgi[1,2] = 1 From 72030d0fee1b19853430c5b56789cd28375c2871 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 21 Aug 2016 08:36:04 -0500 Subject: [PATCH 05/21] Greatly improve test coverage --- src/ImagesFiltering.jl | 2 +- src/border.jl | 10 ++++--- src/imfilter.jl | 61 +++++++++++++++++++++------------------- src/kernel.jl | 3 ++ src/kernelfactors.jl | 2 ++ src/utils.jl | 18 +++++++----- test/2d.jl | 37 ++++++++++++++++++++++--- test/basic.jl | 32 +++++++++++++++++++++ test/border.jl | 60 +++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 1 + test/specialty.jl | 63 ++++++++++++++++++++++++++++++++++++++++-- 11 files changed, 241 insertions(+), 48 deletions(-) create mode 100644 test/basic.jl diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index 824501a..69827b1 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -4,7 +4,7 @@ using Colors, FixedPointNumbers, ImagesCore, MappedArrays, FFTViews, OffsetArray using ColorVectorSpace # in case someone filters RGB arrays using Base: Indices, tail, fill_to_length, @pure -export Kernel, KernelFactors, Pad, Fill, Inner, Algorithm, imfilter, imfilter!, padarray, centered +export Kernel, KernelFactors, Pad, Fill, Inner, NoPad, Algorithm, imfilter, imfilter!, padarray, centered typealias FixedColorant{T<:UFixed} Colorant{T} diff --git a/src/border.jl b/src/border.jl index ab106b6..e5427ea 100644 --- a/src/border.jl +++ b/src/border.jl @@ -47,6 +47,7 @@ Pad the input image symmetrically, `m` pixels at the lower and upper edge of dim Pad the input image symmetrically, `m` pixels at the lower and upper edge of dimension 1, `n` pixels for dimension 2. """ +(::Type{Pad{Style}}){Style}(::Tuple{}) = Pad{Style}() (::Type{Pad{Style}}){Style,N}(both::Dims{N}) = Pad{Style,N}(both, both) (::Type{Pad{Style}}){Style }(lo::Tuple{}, hi::Tuple{}) = Pad{Style,0}(lo, hi) @@ -178,10 +179,10 @@ Fill{T,N}(value::T, lo::Dims{N}, hi::Dims{N}) = Fill{T,N}(value, lo, hi) Fill(value, lo::AbstractVector, hi::AbstractVector) = Fill(value, (lo...,), (hi...,)) Fill{T,N}(value::T, inds::Base.Indices{N}) = Fill{T,N}(value, map(lo,inds), map(hi,inds)) Fill(value, kernel::AbstractArray) = Fill(value, indices(kernel)) +Fill(value, kernel::Laplacian) = Fill(value, indices(kernel)) Fill(value, factkernel::Tuple) = Fill(value, accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) -(p::Fill)(kernel::AbstractArray, img, ::FIR) = Fill(p.value, kernel) -(p::Fill)(factkernel::Tuple, img, ::FIR) = Fill(p.value, factkernel) +(p::Fill)(kernel, img, ::FIR) = Fill(p.value, kernel) function (p::Fill)(factkernel::Tuple, img, ::FFT) inds = accumulate_padding(indices(factkernel[1]), tail(factkernel)...) newinds = map(padfft, inds, map(length, indices(img))) @@ -248,8 +249,10 @@ end # @inline _flatten(t1, t...) = (t1, flatten(t)...) # _flatten() = () -accumulate_padding(inds::Indices, kernel1::AbstractArray, kernels...) = +accumulate_padding(inds::Indices, kernel1, kernels...) = accumulate_padding(_accumulate_padding(inds, indices(kernel1)), kernels...) +accumulate_padding(inds::Indices, kernel1::TriggsSdika, kernels...) = + accumulate_padding(inds, kernels...) accumulate_padding(inds) = inds _accumulate_padding(inds1, inds2) = (__accumulate_padding(inds1[1], inds2[1]), _accumulate_padding(tail(inds1), tail(inds2))...) _accumulate_padding(::Tuple{}, ::Tuple{}) = () @@ -263,7 +266,6 @@ modrange(A::AbstractArray, r::AbstractUnitRange) = map(x->modrange(x, r), A) arraytype{T}(A::AbstractArray, ::Type{T}) = Array{T} # fallback arraytype(A::BitArray, ::Type{Bool}) = BitArray -interior(r::AbstractResource{FFT}, A::AbstractArray, kernel) = indices(A) # periodic boundary conditions interior(r::AbstractResource, A::AbstractArray, kernel) = interior(A, kernel) interior(A::AbstractArray, kernel::Union{AbstractArray,Laplacian}) = _interior(indices(A), indices(kernel)) diff --git a/src/imfilter.jl b/src/imfilter.jl index c9aed28..d077351 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -222,39 +222,29 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke copy!(out, R, A, R) end -iscopy(kernel::AbstractArray) = all(x->x==0:0, indices(kernel)) && first(kernel) == 1 -iscopy(kernel::Laplacian) = false - -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{AbstractArray}, ::NoPad) +function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any}, ::NoPad) kern = kernel[1] iscopy(kern) && return imfilter!(r, out, A, (), NoPad()) imfilter!(r, out, A, samedims(out, kern), NoPad()) end -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Laplacian}, ::NoPad) - imfilter!(r, out, A, samedims(out, kernel[1]), NoPad()) -end - -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple, ::NoPad) +function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, ::NoPad) kern = kernel[1] iscopy(kern) && return imfilter!(r, out, A, tail(kernel), NoPad()) # For multiple stages of filtering, we introduce a second buffer # and swap them at each stage. The first of the two is the one # that holds the most recent result. - A2 = similar(A) # let's hope it's really the *same* type... + A2 = similar(A) # for type-stability, let's hope it's really the *same* type... _imfilter!(r, out, A, A2, kernel, NoPad()) return out end function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, ::NoPad) - R = CartesianRange(indices(out)) - copy!(out, R, A1, R) + imfilter!(r, out, A1, kernel, NoPad()) end function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, ::NoPad) - kern = kernel[1] - iscopy(kern) && return imfilter!(r, out, A, (), NoPad()) - imfilter!(r, out, A1, samedims(out, kern), NoPad()) + imfilter!(r, out, A1, kernel, NoPad()) end function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, ::NoPad) @@ -525,7 +515,7 @@ _imfilter_inplace_tuple!(r, out, img, ::Tuple{}, Rbegin, inds, Rend, border) = o out, img, kernel::TriggsSdika{T,k,l}, Rbegin::CartesianRange, ind::AbstractUnitRange, Rend::CartesianRange, border::BorderSpec) - if all(x->x==0, kernel.a) && all(x->x==0, kernel.b) && kernel.scale == 1 + if iscopy(kernel) if !(out === img) copy!(out, img) end @@ -647,23 +637,30 @@ end ### Utilities -filter_type{S,T}(img::AbstractArray{S}, kernel::AbstractArray{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) -filter_type{S,T}(img::AbstractArray{S}, kernel::IIRFilter{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) -filter_type{S<:Union{UFixed,FixedColorant}}(img::AbstractArray{S}, ::Laplacian) = float32(S) -filter_type{S<:Colorant}(img::AbstractArray{S}, ::Laplacian) = S -filter_type{S<:AbstractFloat}(img::AbstractArray{S}, ::Laplacian) = S -filter_type{S<:Signed}(img::AbstractArray{S}, ::Laplacian) = S -filter_type{S<:Unsigned}(img::AbstractArray{S}, ::Laplacian) = signed(S) +filter_type{S}(img::AbstractArray{S}, kernel) = filter_type(S, kernel) + +filter_type{S,T}(::Type{S}, kernel::AbstractArray{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) +filter_type{S,T}(::Type{S}, kernel::IIRFilter{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) +filter_type{S<:Union{UFixed,FixedColorant}}(::Type{S}, ::Laplacian) = float32(S) +filter_type{S<:Colorant}(::Type{S}, kernel::Laplacian) = S +filter_type{S<:AbstractFloat}(::Type{S}, ::Laplacian) = S +filter_type{S<:Signed}(::Type{S}, ::Laplacian) = S +filter_type{S<:Unsigned}(::Type{S}, ::Laplacian) = signed_type(S) +filter_type(::Type{Bool}, ::Laplacian) = Int8 + +signed_type(::Type{UInt8}) = Int16 +signed_type(::Type{UInt16}) = Int32 +signed_type(::Type{UInt32}) = Int64 -@inline function filter_type(img, kernel::Tuple{Any,Vararg{Any}}) +@inline function filter_type(img::AbstractArray, kernel::Tuple{Any,Vararg{Any}}) T = filter_type(img, kernel[1]) filter_type(T, img, tail(kernel)) end -@inline function filter_type{T}(::Type{T}, img, kernel::Tuple{Any,Vararg{Any}}) +@inline function filter_type{T}(::Type{T}, img::AbstractArray, kernel::Tuple{Any,Vararg{Any}}) S = promote_type(T, filter_type(img, kernel[1])) filter_type(S, img, tail(kernel)) end -filter_type{T}(::Type{T}, img, kernel::Tuple{}) = T +filter_type{T}(::Type{T}, img::AbstractArray, kernel::Tuple{}) = T factorkernel(kernel::AbstractArray) = (kernelshift(indices(kernel), kernel),) factorkernel(L::Laplacian) = (L,) @@ -705,14 +702,16 @@ function prod_kernel{N}(::Type{Val{N}}, kernel::AbstractArray) reshape(kernel, newinds) end +kernelshift(::Tuple{}, A::StridedArray) = A +kernelshift(::Tuple{}, A) = A kernelshift{N}(inds::NTuple{N,Base.OneTo}, A::StridedArray) = _kernelshift(inds, A) kernelshift{N}(inds::NTuple{N,Base.OneTo}, A) = _kernelshift(inds, A) function _kernelshift(inds, A) Base.depwarn("assuming that the origin is at the center of the kernel; to avoid this warning, call `centered(kernel)` or use an OffsetArray", :_kernelshift) # this may be necessary long-term? centered(A) end -kernelshift(inds::Any, A::StridedArray) = OffsetArray(A, inds...) -function kernelshift(inds::Any, A) +kernelshift(inds::Indices, A::StridedArray) = OffsetArray(A, inds...) +function kernelshift(inds::Indices, A) @assert indices(A) == inds A end @@ -721,7 +720,7 @@ filter_algorithm(out, img, kernel) = FIR() filter_algorithm(out, img, kernel::Tuple{IIRFilter,Vararg{IIRFilter}}) = IIR() isseparable(kernels::Tuple{Vararg{TriggsSdika}}) = true -isseparable(kernels::Tuple) = all(x->extendeddims(x)==1, kernels) +isseparable(kernels::Tuple) = all(x->nextendeddims(x)==1, kernels) function imfilter_na_inseparable!{T}(r, out::AbstractArray{T}, img, nanflag, kernel::Tuple{Vararg{TriggsSdika}}) fc, fn = Fill(zero(T)), Fill(zero(eltype(T))) # color, numeric @@ -783,3 +782,7 @@ function normalize_dims!{T,N}(A::AbstractArray{T,N}, factors::NTuple{N}) end A end + +iscopy(kernel::AbstractArray) = all(x->x==0:0, indices(kernel)) && first(kernel) == 1 +iscopy(kernel::Laplacian) = false +iscopy(kernel::TriggsSdika) = all(x->x==0, kernel.a) && all(x->x==0, kernel.b) && kernel.scale == 1 diff --git a/src/kernel.jl b/src/kernel.jl index 0304e6e..50a1028 100644 --- a/src/kernel.jl +++ b/src/kernel.jl @@ -110,8 +110,11 @@ See also: KernelFactors.gaussian. broadcast(.*, KernelFactors.gaussian(σs, ls)...) gaussian(σ::Tuple{Real}, l::Tuple{Integer}) = KernelFactors.gaussian(σ[1], l[1]) gaussian(σ::Tuple{}, l::Tuple{}) = reshape([1]) # 0d +gaussian{T<:Real,I<:Integer}(σs::AbstractVector{T}, ls::AbstractVector{I}) = + gaussian((σs...,), (ls...,)) @inline gaussian{N}(σs::NTuple{N,Real}) = broadcast(.*, KernelFactors.gaussian(σs)...) +gaussian{T<:Real}(σs::AbstractVector{T}) = gaussian((σs...,)) gaussian(σ::Tuple{Real}) = KernelFactors.gaussian(σ[1]) gaussian(σ::Tuple{}) = reshape([1]) diff --git a/src/kernelfactors.jl b/src/kernelfactors.jl index c640a5e..a95a2fd 100644 --- a/src/kernelfactors.jl +++ b/src/kernelfactors.jl @@ -6,6 +6,8 @@ using Base: tail abstract IIRFilter{T} +Base.eltype{T}(kernel::IIRFilter{T}) = T + #### FIR filters ## gradients diff --git a/src/utils.jl b/src/utils.jl index 647d0c7..9d3b317 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -18,23 +18,27 @@ dummyind(::AbstractUnitRange) = 0:0 dummykernel{N}(inds::Indices{N}) = similar(dims->ones(ntuple(d->1,Val{N})), map(dummyind, inds)) -extendeddims(k::AbstractArray) = sum(ind->length(ind)>1, indices(k)) +nextendeddims(inds::Indices) = sum(ind->length(ind)>1, inds) +nextendeddims(a::AbstractArray) = nextendeddims(indices(a)) + +function checkextended(inds::Indices, n) + dimstr = n == 1 ? "dimension" : "dimensions" + nextendeddims(inds) != n && throw(ArgumentError("need $n extended $dimstr, got indices $inds")) + nothing +end +checkextended(a::AbstractArray, n) = checkextended(indices(a), n) _reshape{_,N}(A::OffsetArray{_,N}, ::Type{Val{N}}) = A _reshape{N}(A::OffsetArray, ::Type{Val{N}}) = OffsetArray(reshape(parent(A), Val{N}), fill_to_length(A.offsets, -1, Val{N})) _reshape{N}(A::AbstractArray, ::Type{Val{N}}) = reshape(A, Val{N}) _vec(a::AbstractVector) = a -_vec(a::AbstractArray) = vec(a) -_vec{_}(a::OffsetArray{_,0}) = a +_vec(a::AbstractArray) = (checkextended(a, 1); a) _vec{_}(a::OffsetArray{_,1}) = a function _vec(a::OffsetArray) inds = indices(a) - extendeddims(a) < 2 || throw(ArgumentError("_vec only works on arrays with just 1 extended dimension, got indices $inds")) + checkextended(inds, 1) i = find(ind->length(ind)>1, inds) - if i == 0 - i = 1 - end OffsetArray(vec(parent(a)), inds[i]) end diff --git a/test/2d.jl b/test/2d.jl index 2833ea8..7e7adac 100644 --- a/test/2d.jl +++ b/test/2d.jl @@ -1,4 +1,4 @@ -using ImagesFiltering, ImagesCore, OffsetArrays, Colors, FFTViews, ColorVectorSpace +using ImagesFiltering, ImagesCore, OffsetArrays, Colors, FFTViews, ColorVectorSpace, ComputationalResources using Base.Test @testset "FIR/FFT" begin @@ -18,15 +18,26 @@ using Base.Test for img in (imgf, imgi, imgg, imgc) targetimg = zeros(typeof(img[1]*kern[1]), size(img)) targetimg[3:4,2:3] = rot180(kern)*img[3,4] - @test @inferred(imfilter(img, kernel)) ≈ targetimg + @test (ret = @inferred(imfilter(img, kernel))) ≈ targetimg + @test @inferred(imfilter(img, (kernel,))) ≈ targetimg @test @inferred(imfilter(f32type(img), img, kernel)) ≈ float32(targetimg) + fill!(ret, zero(eltype(ret))) + @test @inferred(imfilter!(ret, img, kernel)) ≈ targetimg + fill!(ret, zero(eltype(ret))) + @test @inferred(imfilter!(CPU1(Algorithm.FIR()), ret, img, kernel)) ≈ targetimg for border in (Pad{:replicate}(), Pad{:circular}(), Pad{:symmetric}(), Pad{:reflect}(), Fill(zero(eltype(img)))) - @test @inferred(imfilter(img, kernel, border)) ≈ targetimg + @test (ret = @inferred(imfilter(img, kernel, border))) ≈ targetimg @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32(targetimg) + fill!(ret, zero(eltype(ret))) + @test @inferred(imfilter!(ret, img, kernel, border)) ≈ targetimg for alg in (Algorithm.FIR(), Algorithm.FFT()) @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg + @test @inferred(imfilter(img, (kernel,), border, alg)) ≈ targetimg @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32(targetimg) + fill!(ret, zero(eltype(ret))) + @test @inferred(imfilter!(CPU1(alg), ret, img, kernel, border)) ≈ targetimg end + @test_throws MethodError imfilter!(CPU1(Algorithm.FIR()), ret, img, kernel, border, Algorithm.FFT()) end targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:5) @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner @@ -34,6 +45,7 @@ using Base.Test for alg in (Algorithm.FIR(), Algorithm.FFT()) @test @inferred(imfilter(img, kernel, Inner(), alg)) ≈ targetimg_inner @test @inferred(imfilter(f32type(img), img, kernel, Inner(), alg)) ≈ float32(targetimg_inner) + @test @inferred(imfilter(CPU1(alg), img, kernel, Inner())) ≈ targetimg_inner end end # Factored kernel @@ -51,6 +63,7 @@ using Base.Test @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32(targetimg) end + @test_throws MethodError imfilter!(CPU1(Algorithm.FIR()), ret, img, kernel, border, Algorithm.FFT()) end targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:5) @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner @@ -173,6 +186,14 @@ using Base.Test @test @inferred(zerona!(imfilter(f32type(img), img, kernel, border, alg), nanflag)) ≈ float32(targetimg) end end + # filtering with a 0d kernel + a = rand(3,3) + # @test_approx_eq imfilter(a, reshape([2])) 2a + # imfilter! as a copy + ret = zeros(3, 3) + b = @inferred(imfilter!(CPU1(Algorithm.FIR()), ret, a, (), NoPad())) + @test b == a + @test !(b === a) end @testset "TriggsSdika 1d" begin @@ -209,12 +230,20 @@ end x = -2:2 y = (-3:3)' kernel = KernelFactors.IIRGaussian((σ, σ)) - for img in (copy(imgf), copy(imgg), copy(imgc)) + for img in (imgf, imgg, imgc) imgcmp = img[3,4]*exp(-(x.^2 .+ y.^2)/(2*σ^2))/(σ^2*2*pi) border = Fill(zero(eltype(img))) + img0 = copy(img) imgf = @inferred(imfilter(img, kernel, border)) @test sumabs2(imgcmp - imgf) < 0.2^2*sumabs2(imgcmp) + @test img == img0 + @test imgf != img end + ret = imfilter(imgf, KernelFactors.IIRGaussian((0,σ)), Pad{:replicate}()) + @test ret != imgf + out = similar(ret) + imfilter!(CPU1(Algorithm.IIR()), out, imgf, KernelFactors.IIRGaussian(σ), 2, Pad{:replicate}()) + @test out == ret end diff --git a/test/basic.jl b/test/basic.jl new file mode 100644 index 0000000..43b5d3b --- /dev/null +++ b/test/basic.jl @@ -0,0 +1,32 @@ +using ImagesFiltering, OffsetArrays, Base.Test + +@testset "basic" begin + @test eltype(KernelFactors.IIRGaussian(3)) == Float64 + @test eltype(KernelFactors.IIRGaussian(Float32, 3)) == Float32 + @test isa(KernelFactors.IIRGaussian([1,2.0f0]), Tuple{KernelFactors.TriggsSdika,KernelFactors.TriggsSdika}) + + @test KernelFactors.kernelfactors(([0,3], [1,7])) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2)) + @test KernelFactors.kernelfactors(([0,3], [1,7]')) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2)) + + # Warnings + const OLDERR = STDERR + fname = tempname() + open(fname, "w") do f + redirect_stderr(f) + try + KernelFactors.IIRGaussian(0.5, emit_warning=false) + finally + redirect_stderr(OLDERR) + end + end + @test isempty(chomp(readstring(fname))) + open(fname, "w") do f + redirect_stderr(f) + try + KernelFactors.IIRGaussian(0.5) + finally + redirect_stderr(OLDERR) + end + end + @test contains(readstring(fname), "too small for accuracy") +end diff --git a/test/border.jl b/test/border.jl index 513a090..06cfcb0 100644 --- a/test/border.jl +++ b/test/border.jl @@ -1,4 +1,4 @@ -using ImagesFiltering, OffsetArrays +using ImagesFiltering, OffsetArrays, Colors using Base.Test @testset "padarray" begin @@ -122,6 +122,64 @@ using Base.Test B = view(A,1:8,1:8,1:8) @test isa(parent(padarray(A, Pad((1,1,1)))), BitArray{3}) @test isa(parent(padarray(B, Pad((1,1,1)))), BitArray{3}) + A = reshape(1:15, 3, 5) + B = @inferred(padarray(A, Inner((1,1)))) + @test B == A + # test that it's a copy + B[1,1] = 0 + @test B != A + A = rand(RGB{U8}, 3, 5) + ret = @test_throws ErrorException padarray(A, Fill(0, (0,0), (0,0))) + @test contains(ret.value.msg, "element type ColorTypes.RGB") + A = bitrand(3, 5) + ret = @test_throws ErrorException padarray(A, Fill(7, (0,0), (0,0))) + @test contains(ret.value.msg, "element type Bool") + @test isa(parent(padarray(A, Fill(false, (1,1), (1,1)))), BitArray) +end + +@testset "Pad" begin + @test Pad{:replicate}([1,2], [5,3]) == Pad{:replicate}((1,2), (5,3)) + @test @inferred(Pad{:replicate,2}([1,2], [5,3])) == Pad{:replicate}((1,2), (5,3)) + @test_throws TypeError Pad{:replicate,3}([1,2], [5,3]) + @test @inferred(Pad{:circular}()(rand(3,5))) == Pad{:circular,2}((0,0),(3,5)) + @test @inferred(Pad{:circular}()(centered(rand(3,5)))) == Pad{:circular,2}((1,2),(1,2)) + @test @inferred(Pad{:symmetric}()(Kernel.Laplacian())) == Pad{:symmetric,2}((1,1),(1,1)) + @test @inferred(Pad{:symmetric}()(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Pad{:symmetric,2}((1,1),(1,1)) + + a = reshape(1:15, 3, 5) + targetinds = (OffsetArray([1,1,2,3,3], 0:4), OffsetArray([1:5;], 1:5)) + ERROLD = STDERR + tmpfile = tempname() + try + open(tmpfile, "w") do f + redirect_stderr(f) + @test @inferred(ImagesFiltering.padindices(rand(3,5), Pad{:replicate}((1,)))) == targetinds + end + finally + redirect_stderr(ERROLD) + end + @test contains(readstring(tmpfile), "lacks the proper padding sizes") +end + +@testset "Inner" begin + @test @inferred(Inner((1,2))) == Inner((1,2), (1,2)) + @test @inferred(Inner((1,2), ())) == Inner((1,2), (0,0)) + @test @inferred(Inner((), (1,2))) == Inner((0,0), (1,2)) + @test @inferred(Inner{2}([1,2],[5,6])) == Inner((1,2), (5,6)) + @test Inner([1,2],[5,6]) == Inner((1,2), (5,6)) + @test @inferred(Inner()(rand(3,5))) == Inner((0,0),(3,5)) +end + +@testset "Fill" begin + @test Fill(1,[1,2],[5,6]) == Fill(1, (1,2), (5,6)) + @test @inferred(Fill(-1, rand(3,5))) == Fill(-1, (0,0), (3,5)) + @test @inferred(Fill(2)(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Fill(2, (1,1),(1,1)) +end + +@testset "misc" begin + a0 = reshape([1]) # 0-dimensional + @test ImagesFiltering.accumulate_padding((), a0) == () + @test ImagesFiltering.accumulate_padding((0:1, -1:1), a0) == (0:1, -1:1) end nothing diff --git a/test/runtests.jl b/test/runtests.jl index 87a0002..23981e9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ include("border.jl") include("2d.jl") include("cascade.jl") include("specialty.jl") +include("basic.jl") include("deprecated.jl") nothing diff --git a/test/specialty.jl b/test/specialty.jl index 7998c3c..947d302 100644 --- a/test/specialty.jl +++ b/test/specialty.jl @@ -3,6 +3,21 @@ using Base.Test @testset "specialty" begin @testset "Laplacian" begin + L1 = OffsetArray([1,-2,1],-1:1) + L2 = OffsetArray([0 1 0; 1 -4 1; 0 1 0], -1:1, -1:1) + kern = Kernel.Laplacian((true,)) + @test !isempty(kern) + @test convert(AbstractArray, kern) == L1 + kern = Kernel.Laplacian() + @test convert(AbstractArray, kern) == L2 + kern = Kernel.Laplacian((true,false)) + @test convert(AbstractArray, kern) == reshape(L1, -1:1, 0:0) + kern = Kernel.Laplacian((false,true)) + @test convert(AbstractArray, kern) == reshape(L1, 0:0, -1:1) + kern = Kernel.Laplacian([1], 2) + @test convert(AbstractArray, kern) == reshape(L1, -1:1, 0:0) + kern = Kernel.Laplacian([2], 2) + @test convert(AbstractArray, kern) == reshape(L1, 0:0, -1:1) function makeimpulse(T, sz, x) A = zeros(T, sz) A[x] = one(T) @@ -10,10 +25,15 @@ using Base.Test end # 1d kern = Kernel.Laplacian((true,)) - @test convert(AbstractArray, kern) == OffsetArray([1,-2,1],-1:1) for a in (makeimpulse(Float64, (5,), 3), + makeimpulse(Int, (5,), 3), + makeimpulse(UInt32, (5,), 3), + makeimpulse(UInt16, (5,), 3), + makeimpulse(UInt8, (5,), 3), + makeimpulse(Bool, (5,), 3), makeimpulse(Gray{U8}, (5,), 3), - makeimpulse(RGB{Float32}, (5,), 3)) + makeimpulse(RGB{Float32}, (5,), 3), + makeimpulse(RGB{U8}, (5,), 3)) af = imfilter(a, kern) T = eltype(a) @test af == [zero(T),a[3],-2a[3],a[3],zero(T)] @@ -127,6 +147,45 @@ using Base.Test end end end + + @testset "gaussian" begin + function gaussiancmp(σ, xr) + cmp = [exp(-x^2/(2σ^2)) for x in xr] + OffsetArray(cmp/sum(cmp), xr) + end + @test_approx_eq KernelFactors.gaussian(2, 9) gaussiancmp(2, -4:4) + k = KernelFactors.gaussian((2,3), (9,7)) + @test_approx_eq k[1] gaussiancmp(2, -4:4) + @test_approx_eq k[2] gaussiancmp(3, -3:3) + @test_approx_eq sum(KernelFactors.gaussian(5)) 1 + for k = (KernelFactors.gaussian((2,3)), KernelFactors.gaussian([2,3]), KernelFactors.gaussian([2,3], [9,7])) + @test_approx_eq sum(k[1]) 1 + @test_approx_eq sum(k[2]) 1 + end + @test_approx_eq Kernel.gaussian((2,), (9,)) gaussiancmp(2, -4:4) + @test_approx_eq Kernel.gaussian((2,3), (9,7)) gaussiancmp(2, -4:4).*gaussiancmp(3, -3:3)' + @test_approx_eq sum(Kernel.gaussian(5)) 1 + for k = (Kernel.gaussian((2,3)), Kernel.gaussian([2,3]), Kernel.gaussian([2,3], [9,7])) + @test_approx_eq sum(k) 1 + end + end + + @testset "DoG" begin + function gaussiancmp(σ, xr) + cmp = [exp(-x^2/(2σ^2)) for x in xr] + OffsetArray(cmp/sum(cmp), xr) + end + @test_approx_eq Kernel.DoG((2,), (3,), (9,)) gaussiancmp(2, -4:4) - gaussiancmp(3, -4:4) + k = Kernel.DoG((2,3), (4,3.5), (9,7)) + k1 = gaussiancmp(2, -4:4) .* gaussiancmp(3, -3:3)' + k2 = gaussiancmp(4, -4:4) .* gaussiancmp(3.5, -3:3)' + @test_approx_eq k k1-k2 + @test abs(sum(Kernel.DoG(5))) < 1e-8 + @test Kernel.DoG(5) == Kernel.DoG((5,5)) + @test abs(sum(Kernel.DoG((5,), (7,), (21,)))) < 1e-8 + @test indices(Kernel.DoG((5,), (7,), (21,))) == (-10:10,) + end + end nothing From 4ad88aa441c85b97ad87ef6bd58f7d27346dc110 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 23 Aug 2016 09:23:09 -0500 Subject: [PATCH 06/21] Introduce ReshapedVector to unify TriggsSdika with array kernels Now it's possible to construct kernel-tuples with a mix of TS and AbstractArray kernels, and things Just Work. --- src/ImagesFiltering.jl | 5 +- src/border.jl | 20 +++++-- src/imfilter.jl | 131 ++++++++++++++++++----------------------- src/kernelfactors.jl | 111 +++++++++++++++++++++++++++++----- src/specialty.jl | 29 ++------- src/utils.jl | 2 + test/basic.jl | 2 +- 7 files changed, 179 insertions(+), 121 deletions(-) diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index 69827b1..2b6bac0 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -22,10 +22,13 @@ Alg{A<:Alg}(r::AbstractResource{A}) = r.settings include("utils.jl") include("kernelfactors.jl") -using .KernelFactors: TriggsSdika, IIRFilter +using .KernelFactors: TriggsSdika, IIRFilter, ReshapedVector, iterdims include("kernel.jl") using .Kernel using .Kernel: Laplacian + +typealias ArrayLike{T} Union{AbstractArray{T}, IIRFilter{T}, ReshapedVector{T}} + include("border.jl") include("deprecated.jl") include("imfilter.jl") diff --git a/src/border.jl b/src/border.jl index e5427ea..962fd60 100644 --- a/src/border.jl +++ b/src/border.jl @@ -4,15 +4,21 @@ using OffsetArrays, CatIndices abstract AbstractBorder -immutable NoPad <: AbstractBorder end +immutable NoPad{T} <: AbstractBorder + border::T +end +NoPad() = NoPad(nothing) """ NoPad() + NoPad(border) -Indicates that no padding should be applied to the input array. +Indicates that no padding should be applied to the input array. Passing a `border` object allows you to preserve "memory" of a border choice; it can be retrieved by indexing with `[]`. """ NoPad +Base.getindex(np::NoPad) = np.border + """ `Pad{Style,N}` is a type that stores choices about padding. `Style` is a Symbol specifying the boundary conditions of the image, one of: @@ -125,8 +131,12 @@ padarray(img::AbstractArray, border::Pad) = padarray(eltype(img), img, border) function padarray{T}(::Type{T}, img::AbstractArray, border::Pad) inds = padindices(img, border) # like img[inds...] except that we can control the element type - dest = similar(img, T, map(Base.indices1, inds)) - Base._unsafe_getindex!(dest, img, inds...) + newinds = map(Base.indices1, inds) + dest = similar(img, T, newinds) + @unsafe for I in CartesianRange(newinds) + J = CartesianIndex(map((i,x)->x[i], I.I, inds)) + dest[I] = img[J] + end dest end padarray{P}(img, ::Type{P}) = img[padindices(img, P)...] # just to throw the nice error @@ -268,7 +278,7 @@ arraytype(A::BitArray, ::Type{Bool}) = BitArray interior(r::AbstractResource, A::AbstractArray, kernel) = interior(A, kernel) -interior(A::AbstractArray, kernel::Union{AbstractArray,Laplacian}) = _interior(indices(A), indices(kernel)) +interior(A::AbstractArray, kernel::Union{ArrayLike,Laplacian}) = _interior(indices(A), indices(kernel)) interior(A, factkernel::Tuple) = _interior(indices(A), accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) function _interior{N}(indsA::NTuple{N}, indsk) indskN = fill_to_length(indsk, 0:0, Val{N}) diff --git a/src/imfilter.jl b/src/imfilter.jl index d077351..5c0d6f3 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -2,7 +2,8 @@ typealias BorderSpec{Style,T} Union{Pad{Style,0}, Fill{T,0}, Inner{0}} typealias BorderSpecNoNa{T} Union{Pad{:replicate,0}, Pad{:circular,0}, Pad{:symmetric,0}, Pad{:reflect,0}, Fill{T,0}, Inner{0}} typealias BorderSpecAny Union{BorderSpec,NoPad} -typealias NDimKernel{N,K} Union{AbstractArray{K,N},Laplacian{N}} +typealias NDimKernel{N,K} Union{AbstractArray{K,N},ReshapedVector{K,N},Laplacian{N}} +typealias ReshapedTriggsSdika{T,N,Npre,V<:TriggsSdika} ReshapedVector{T,N,Npre,V} typealias ProcessedKernel Tuple @@ -14,7 +15,7 @@ typealias ProcessedKernel Tuple end # Step 2: if necessary, put the kernel into cannonical (factored) form -@inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::Union{AbstractArray,Laplacian}, args...) +@inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}, args...) imfilter(T, img, factorkernel(kernel), args...) end @@ -38,7 +39,7 @@ end imfilter(r, filter_type(img, kernel), img, kernel, args...) end -@inline function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::AbstractArray, args...) +@inline function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ArrayLike, args...) imfilter(r, T, img, factorkernel(kernel), args...) end @@ -120,15 +121,15 @@ imfilter # have to be a little more cautious to make sure that later methods # don't inadvertently call back to these: in methods that take an # AbstractResource argument, exclude `NoPad()` as a border option. -function imfilter!(out::AbstractArray, img::AbstractArray, kernel::Union{AbstractArray,Laplacian}, args...) +function imfilter!(out::AbstractArray, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}, args...) imfilter!(out, img, factorkernel(kernel), args...) end -function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::Union{AbstractArray,Laplacian}) +function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}) imfilter!(r, out, img, factorkernel(kernel)) end -function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::Union{AbstractArray,Laplacian}, border::BorderSpec) +function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}, border::BorderSpec) imfilter!(r, out, img, factorkernel(kernel), border) end @@ -214,7 +215,7 @@ function imfilter!{S,T,N}(r::AbstractResource, # By specifying NoPad(), we ensure that dispatch will never # accidentally "go back" to an earlier routine and apply more # padding - imfilter!(r, out, A, kernel, NoPad()) + imfilter!(r, out, A, kernel, NoPad(border)) end function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{}, ::NoPad) @@ -222,37 +223,37 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke copy!(out, R, A, R) end -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any}, ::NoPad) +function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any}, border::NoPad) kern = kernel[1] - iscopy(kern) && return imfilter!(r, out, A, (), NoPad()) - imfilter!(r, out, A, samedims(out, kern), NoPad()) + iscopy(kern) && return imfilter!(r, out, A, (), border) + imfilter!(r, out, A, samedims(out, kern), border) end -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, ::NoPad) +function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) kern = kernel[1] - iscopy(kern) && return imfilter!(r, out, A, tail(kernel), NoPad()) + iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border) # For multiple stages of filtering, we introduce a second buffer # and swap them at each stage. The first of the two is the one # that holds the most recent result. A2 = similar(A) # for type-stability, let's hope it's really the *same* type... - _imfilter!(r, out, A, A2, kernel, NoPad()) + _imfilter!(r, out, A, A2, kernel, border) return out end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, ::NoPad) - imfilter!(r, out, A1, kernel, NoPad()) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, border::NoPad) + imfilter!(r, out, A1, kernel, border) end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, ::NoPad) - imfilter!(r, out, A1, kernel, NoPad()) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, border::NoPad) + imfilter!(r, out, A1, kernel, border) end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, ::NoPad) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) kern = kernel[1] - iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), NoPad()) + iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), border) kernN = samedims(A2, kern) - imfilter!(r, A2, A1, kernN, NoPad(), interior(r, A2, kernN)) # store result in A2 - _imfilter!(r, out, A2, A1, tail(kernel), NoPad()) # swap the buffers + imfilter!(r, A2, A1, kernN, border, interior(r, A2, kernN)) # store result in A2 + _imfilter!(r, out, A2, A1, tail(kernel), border) # swap the buffers end ## FIR filtering @@ -281,11 +282,11 @@ each stage. See also: imfilter. """ -function imfilter!{S,T,N}(::CPU1{FIR}, +function imfilter!{S,T,N}(r::CPU1{FIR}, out::AbstractArray{S,N}, A::AbstractArray{T,N}, kern::NDimKernel{N}, - ::NoPad, + border::NoPad, R::CartesianRange=CartesianRange(indices(out))) (isempty(A) || isempty(kern)) && return out indso, indsA, indsk = indices(out), indices(A), indices(kern) @@ -303,56 +304,32 @@ function imfilter!{S,T,N}(::CPU1{FIR}, throw(DimensionMismatch("requested range $R and kernel indices $indsk do not agree with indices of padded input, $indsA")) end end - # For performance reasons, we now dispatch to a method that strips + # For performance reasons, we now dispatch to a method that may strip # off index-shift containers (in Julia 0.5 we shouldn't remove # boundschecks from OffsetArray). - _imfilter_inbounds!(out, A, kern, R) + _imfilter_inbounds!(r, out, A, kern, border, R) end -function _imfilter_inbounds!(out, A::OffsetArray, kern::OffsetArray, R) - ΔI = CartesianIndex(kern.offsets) - CartesianIndex(A.offsets) - _imfilter_inbounds!(out, (parent(A), ΔI), parent(kern), R) +function _imfilter_inbounds!{T,N,Npre,V<:TriggsSdika}(r::AbstractResource, out, A::AbstractArray, kern::ReshapedVector{T,N,Npre,V}, border::NoPad, R) + indspre, ind, indspost = iterdims(ranges(R), kern) + _imfilter_dim!(r, out, A, kern.data, CartesianRange(indspre), ind, CartesianRange(indspost), border[]) end -function _imfilter_inbounds!(out, A::AbstractArray, kern::OffsetArray, R) - ΔI = CartesianIndex(kern.offsets) - _imfilter_inbounds!(out, (A, ΔI), parent(kern), R) -end -function _imfilter_inbounds!(out, A::OffsetArray, kern, R) - ΔI = -CartesianIndex(A.offsets) - _imfilter_inbounds!(out, (parent(A), ΔI), kern, R) -end -function _imfilter_inbounds!(out, A::AbstractArray, kern, R) + +function _imfilter_inbounds!(r, out, A::AbstractArray, kern, border::NoPad, R) indsk = indices(kern) Rk = CartesianRange(indsk) p = A[first(R)+first(Rk)] * first(kern) TT = typeof(p+p) for I in R tmp = zero(TT) - @inbounds for J in CartesianRange(indsk) + @unsafe for J in CartesianRange(indsk) tmp += A[I+J]*kern[J] end - @inbounds out[I] = tmp - end - out -end -function _imfilter_inbounds!(out, Ashift::Tuple{AbstractArray,CartesianIndex}, kern, R) - A, ΔI = Ashift - indsk = indices(kern) - Rk = CartesianRange(indsk) - p = A[first(R)+first(Rk)+ΔI] * first(kern) - TT = typeof(p+p) - for I in R - Ishift = I + ΔI - tmp = zero(TT) - @inbounds for J in Rk - tmp += A[Ishift+J]*kern[J] - end - @inbounds out[I] = tmp + @unsafe out[I] = tmp end out end - ### FFT filtering """ @@ -377,30 +354,30 @@ function imfilter!{S,T,K,N}(r::AbstractCPU{FFT}, out::AbstractArray{S,N}, img::AbstractArray{T,N}, kernel::AbstractArray{K,N}, - ::NoPad) - imfilter!(r, out, img, (kernel,), NoPad()) + border::NoPad) + imfilter!(r, out, img, (kernel,), border) end function imfilter!{S,T,N}(r::AbstractCPU{FFT}, out::AbstractArray{S,N}, A::AbstractArray{T,N}, kernel::Tuple{AbstractArray}, - ::NoPad) - _imfilter_fft!(r, out, A, kernel, NoPad()) # ambiguity resolution + border::NoPad) + _imfilter_fft!(r, out, A, kernel, border) # ambiguity resolution end function imfilter!{S,T,N}(r::AbstractCPU{FFT}, out::AbstractArray{S,N}, A::AbstractArray{T,N}, kernel::Tuple{AbstractArray,Vararg{AbstractArray}}, - ::NoPad) - _imfilter_fft!(r, out, A, kernel, NoPad()) + border::NoPad) + _imfilter_fft!(r, out, A, kernel, border) end function _imfilter_fft!{S,T,N}(r::AbstractCPU{FFT}, out::AbstractArray{S,N}, A::AbstractArray{T,N}, kernel::Tuple{AbstractArray,Vararg{AbstractArray}}, - ::NoPad) + border::NoPad) kern = prod_kernel(Val{N}, kernel...) krn = FFTView(zeros(eltype(kern), map(length, indices(A)))) for I in CartesianRange(indices(kern)) @@ -534,7 +511,7 @@ _imfilter_inplace_tuple!(r, out, img, ::Tuple{}, Rbegin, inds, Rend, border) = o # is Pad{:replicate}, so that the original value is still # available. rightborder! will handle that point. for i = ind[k]+1:ind[end-1] - @inbounds for Ibegin in Rbegin + @unsafe for Ibegin in Rbegin tmp = one(T)*img[Ibegin, i, Iend] for j = 1:k tmp += kernel.a[j]*out[Ibegin, i-j, Iend] @@ -549,7 +526,7 @@ _imfilter_inplace_tuple!(r, out, img, ::Tuple{}, Rbegin, inds, Rend, border) = o end # Propagate backwards for i = ind[end-l]:-1:ind[1] - @inbounds for Ibegin in Rbegin + @unsafe for Ibegin in Rbegin tmp = one(T)*out[Ibegin, i, Iend] for j = 1:l tmp += kernel.b[j]*out[Ibegin, i+j, Iend] @@ -559,7 +536,7 @@ _imfilter_inplace_tuple!(r, out, img, ::Tuple{}, Rbegin, inds, Rend, border) = o end # Final scaling for i in ind - @inbounds for Ibegin in Rbegin + @unsafe for Ibegin in Rbegin out[Ibegin, i, Iend] *= kernel.scale end end @@ -639,8 +616,7 @@ end filter_type{S}(img::AbstractArray{S}, kernel) = filter_type(S, kernel) -filter_type{S,T}(::Type{S}, kernel::AbstractArray{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) -filter_type{S,T}(::Type{S}, kernel::IIRFilter{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) +filter_type{S,T}(::Type{S}, kernel::ArrayLike{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T)) filter_type{S<:Union{UFixed,FixedColorant}}(::Type{S}, ::Laplacian) = float32(S) filter_type{S<:Colorant}(::Type{S}, kernel::Laplacian) = S filter_type{S<:AbstractFloat}(::Type{S}, ::Laplacian) = S @@ -652,15 +628,15 @@ signed_type(::Type{UInt8}) = Int16 signed_type(::Type{UInt16}) = Int32 signed_type(::Type{UInt32}) = Int64 -@inline function filter_type(img::AbstractArray, kernel::Tuple{Any,Vararg{Any}}) - T = filter_type(img, kernel[1]) - filter_type(T, img, tail(kernel)) +@inline function filter_type{S}(::Type{S}, kernel::Tuple{Any,Vararg{Any}}) + T = filter_type(S, kernel[1]) + filter_type(T, S, tail(kernel)) end -@inline function filter_type{T}(::Type{T}, img::AbstractArray, kernel::Tuple{Any,Vararg{Any}}) - S = promote_type(T, filter_type(img, kernel[1])) - filter_type(S, img, tail(kernel)) +@inline function filter_type{T,S}(::Type{T}, ::Type{S}, kernel::Tuple{Any,Vararg{Any}}) + Tnew = promote_type(T, filter_type(S, kernel[1])) + filter_type(Tnew, S, tail(kernel)) end -filter_type{T}(::Type{T}, img::AbstractArray, kernel::Tuple{}) = T +filter_type{T,S}(::Type{T}, ::Type{S}, kernel::Tuple{}) = T factorkernel(kernel::AbstractArray) = (kernelshift(indices(kernel), kernel),) factorkernel(L::Laplacian) = (L,) @@ -757,6 +733,7 @@ function imfilter_na_separable!{T}(r, out::AbstractArray{T}, img, kernel::Tuple) normalize_separable!(r, out, kernel, fn) end +normalize_separable!(r::AbstractResource, A, ::Tuple{}, border) = error("this shouldn't happen") function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,TriggsSdika}, border) inds = indices(A) function imfilter_inplace!(r, a, kern, border) @@ -765,6 +742,9 @@ function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,Trigg filtdims = ntuple(d->imfilter_inplace!(r, similar(dims->ones(dims), inds[d]), kernels[d], border), Val{N}) normalize_dims!(A, filtdims) end +function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,ReshapedTriggsSdika}, border) + normalize_separable!(r, A, map(_vec, kernels), border) +end function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N}, border) inds = indices(A) @@ -786,3 +766,4 @@ end iscopy(kernel::AbstractArray) = all(x->x==0:0, indices(kernel)) && first(kernel) == 1 iscopy(kernel::Laplacian) = false iscopy(kernel::TriggsSdika) = all(x->x==0, kernel.a) && all(x->x==0, kernel.b) && kernel.scale == 1 +iscopy(kernel::ReshapedVector) = iscopy(kernel.data) diff --git a/src/kernelfactors.jl b/src/kernelfactors.jl index a95a2fd..251f32f 100644 --- a/src/kernelfactors.jl +++ b/src/kernelfactors.jl @@ -1,13 +1,82 @@ module KernelFactors using StaticArrays, OffsetArrays -using ..ImagesFiltering: centered, dummyind, _reshape -using Base: tail +using ..ImagesFiltering: centered, dummyind +import ..ImagesFiltering: _reshape, _vec, nextendeddims +using Base: tail, Indices abstract IIRFilter{T} Base.eltype{T}(kernel::IIRFilter{T}) = T +""" + ReshapedVector{N,Npre}(data) + +Return an object of dimensionality `N`, where `data` must have +dimensionality 1. The indices are `0:0` for the first `Npre` +dimensions, have the indices of `data` for dimension `Npre+1`, and are +`0:0` for the remaining dimensions. + +`data` must support `eltype` and `ndims`, but does not have to be an +AbstractArray. + +ReshapedVectors allow one to specify a "filtering dimension" for a +1-dimensional filter. +""" +immutable ReshapedVector{T,N,Npre,V} # note not <: AbstractArray{T,N} (more general) + data::V + + function ReshapedVector(data::V) + ndims(V) == 1 || throw(DimensionMismatch("must be one dimensional, got $(ndims(V))")) + new(data) + end +end + +(::Type{ReshapedVector{N,Npre}}){N,Npre,V}(data::V) = ReshapedVector{eltype(data),N,Npre,V}(data) +# Convenient loop constructor that uses dummy NTuple{N,Bool} to keep +# track of dimensions for type-stability +@inline function ReshapedVector{Npre,Npost}(pre::NTuple{Npre,Bool}, data, post::NTuple{Npost,Bool}) + total = (pre..., true, post...) + _reshapedvector(total, pre, data) +end +_reshapedvector{N,Npre}(::NTuple{N,Bool}, ::NTuple{Npre,Bool}, data) = ReshapedVector{eltype(data),N,Npre,typeof(data)}(data) + +Base.eltype{T}(A::ReshapedVector{T}) = T +Base.ndims{_,N}(A::ReshapedVector{_,N}) = N +Base.isempty(A::ReshapedVector) = isempty(A.data) + +@inline Base.indices{_,N,Npre}(A::ReshapedVector{_,N,Npre}) = Base.fill_to_length((Base.ntuple(d->0:0, Val{Npre})..., Base.indices1(A.data)), 0:0, Val{N}) + +_reshape{T,N}(A::ReshapedVector{T,N}, ::Type{Val{N}}) = A +_vec(A::ReshapedVector) = A.data +nextendeddims(a::ReshapedVector) = 1 + +""" + iterdims(inds, v::ReshapedVector{T,N,Npre}) -> Rpre, ind, Rpost + +Return a pair `Rpre`, `Rpost` of CartesianRanges that correspond to +pre- and post- dimensions for iterating over an array with indices +`inds`. `Rpre` corresponds to the `Npre` "pre" dimensions, and `Rpost` +to the trailing dimensions (not including the vector object wrapped in +`v`). Concretely, + + Rpre = CartesianRange(inds[1:Npre]) + ind = inds[Npre+1] + Rpost = CartesianRange(inds[Npre+2:end]) + +although the implementation differs for reason of type-stability. +""" +function iterdims{_,N,Npre}(inds::Indices{N}, v::ReshapedVector{_,N,Npre}) + indspre, ind, indspost = _iterdims((), (), inds, v) + CartesianRange(indspre), ind, CartesianRange(indspost) +end +@inline function _iterdims(indspre, ::Tuple{}, inds, v) + _iterdims((indspre..., inds[1]), (), tail(inds), v) # consume inds and push to indspre +end +@inline function _iterdims{_,N,Npre}(indspre::NTuple{Npre}, ::Tuple{}, inds, v::ReshapedVector{_,N,Npre}) + indspre, inds[1], tail(inds) # return the "central" and trailing dimensions +end + #### FIR filters ## gradients @@ -150,6 +219,11 @@ function TriggsSdika{T}(a::SVector{3,T}, scale) TriggsSdika(a, a, scale, M/Mdenom) end Base.vec(kernel::TriggsSdika) = kernel +Base.ndims(kernel::TriggsSdika) = 1 +Base.ndims{T<:TriggsSdika}(::Type{T}) = 1 +Base.indices1(kernel::TriggsSdika) = 0:0 +Base.indices(kernel::TriggsSdika) = (indices1(kernel),) +Base.isempty(kernel::TriggsSdika) = false # Note that there's a sign reversal between Young & Triggs. """ @@ -169,14 +243,14 @@ which case `Float64` will be chosen). I. T. Young, L. J. van Vliet, and M. van Ginkel, "Recursive Gabor Filtering". IEEE Trans. Sig. Proc., 50: 2798-2805 (2002). """ -function IIRGaussian{T}(::Type{T}, sigma::Real; emit_warning::Bool = true) - if emit_warning && sigma < 1 && sigma != 0 - warn("sigma is too small for accuracy") +function IIRGaussian{T}(::Type{T}, σ::Real; emit_warning::Bool = true) + if emit_warning && σ < 1 && σ != 0 + warn("σ is too small for accuracy") end m0 = convert(T,1.16680) m1 = convert(T,1.10783) m2 = convert(T,1.40586) - q = convert(T,1.31564*(sqrt(1+0.490811*sigma*sigma) - 1)) + q = convert(T,1.31564*(sqrt(1+0.490811*σ*σ) - 1)) ascale = (m0+q)*(m1*m1 + m2*m2 + 2m1*q + q*q) B = (m0*(m1*m1 + m2*m2)/ascale)^2 # This is what Young et al call -b, but in filt() notation would be called a @@ -186,19 +260,26 @@ function IIRGaussian{T}(::Type{T}, sigma::Real; emit_warning::Bool = true) a = SVector(a1,a2,a3) TriggsSdika(a, B) end -IIRGaussian(sigma::Real; emit_warning::Bool = true) = IIRGaussian(iirgt(sigma), sigma; emit_warning=emit_warning) +IIRGaussian(σ::Real; emit_warning::Bool = true) = IIRGaussian(iirgt(σ), σ; emit_warning=emit_warning) -function IIRGaussian{T}(::Type{T}, sigma::Tuple; emit_warning::Bool = true) - map(s->IIRGaussian(T, s; emit_warning=emit_warning), sigma) +function IIRGaussian{T,N}(::Type{T}, σs::NTuple{N,Real}; emit_warning::Bool = true) + iirg(T, (), σs, tail(ntuple(d->true, Val{N})), emit_warning) end -IIRGaussian(sigma::Tuple; emit_warning::Bool = true) = IIRGaussian(iirgt(sigma), sigma; emit_warning=emit_warning) +IIRGaussian(σs::Tuple; emit_warning::Bool = true) = IIRGaussian(iirgt(σs), σs; emit_warning=emit_warning) + +IIRGaussian(σs::AbstractVector; kwargs...) = IIRGaussian((σs...,); kwargs...) +IIRGaussian{T}(::Type{T}, σs::AbstractVector; kwargs...) = IIRGaussian(T, (σs...,); kwargs...) -IIRGaussian(sigma::AbstractVector; kwargs...) = IIRGaussian((sigma...,); kwargs...) -IIRGaussian{T}(::Type{T}, sigma::AbstractVector; kwargs...) = IIRGaussian(T, (sigma...,); kwargs...) +iirgt(σ::AbstractFloat) = typeof(σ) +iirgt(σ::Real) = Float64 +iirgt(σs::Tuple) = promote_type(map(iirgt, σs)...) -iirgt(sigma::AbstractFloat) = typeof(sigma) -iirgt(sigma::Real) = Float64 -iirgt(sigma::Tuple) = promote_type(map(iirgt, sigma)...) +@inline function iirg{T}(::Type{T}, pre, σs, post, emit_warning) + kern = ReshapedVector(pre, IIRGaussian(T, σs[1]; emit_warning=emit_warning), post) + (kern, iirg(T, (pre..., post[1]), tail(σs), tail(post), emit_warning)...) +end +iirg{T}(::Type{T}, pre, σs::Tuple{Real}, ::Tuple{}, emit_warning) = + (ReshapedVector(pre, IIRGaussian(T, σs[1]; emit_warning=emit_warning), ()),) ###### Utilities diff --git a/src/specialty.jl b/src/specialty.jl index 6377067..867a79d 100644 --- a/src/specialty.jl +++ b/src/specialty.jl @@ -1,34 +1,15 @@ ## Laplacian -function _imfilter_inbounds!{N}(out, A::OffsetArray, kern::Laplacian{N}, R) - ΔI = -CartesianIndex(A.offsets) - _imfilter_inbounds!(out, (parent(A), ΔI), kern, R) -end -function _imfilter_inbounds!(out, A::AbstractArray, L::Laplacian, R) +function _imfilter_inbounds!(r, out, A::AbstractArray, L::Laplacian, border::NoPad, R) TT = eltype(out) # accumtype(eltype(out), eltype(A)) n = 2*length(L.offsets) - for I in R - @inbounds tmp = convert(TT, - n * A[I]) - @inbounds for J in L.offsets + @unsafe for I in R + tmp = convert(TT, - n * A[I]) + for J in L.offsets tmp += A[I+J] tmp += A[I-J] end - @inbounds out[I] = tmp - end - out -end -function _imfilter_inbounds!(out, Ashift::Tuple{AbstractArray,CartesianIndex}, L::Laplacian, R) - A, ΔI = Ashift - TT = eltype(out) # accumtype(eltype(out), eltype(A)) - n = 2*length(L.offsets) - for I in R - Ishift = I + ΔI - @inbounds tmp = convert(TT, - n * A[Ishift]) - @inbounds for J in L.offsets - tmp += A[Ishift+J] - tmp += A[Ishift-J] - end - @inbounds out[I] = tmp + out[I] = tmp end out end diff --git a/src/utils.jl b/src/utils.jl index 9d3b317..47e2343 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -28,6 +28,8 @@ function checkextended(inds::Indices, n) end checkextended(a::AbstractArray, n) = checkextended(indices(a), n) +ranges(R::CartesianRange) = map(colon, R.start.I, R.stop.I) + _reshape{_,N}(A::OffsetArray{_,N}, ::Type{Val{N}}) = A _reshape{N}(A::OffsetArray, ::Type{Val{N}}) = OffsetArray(reshape(parent(A), Val{N}), fill_to_length(A.offsets, -1, Val{N})) _reshape{N}(A::AbstractArray, ::Type{Val{N}}) = reshape(A, Val{N}) diff --git a/test/basic.jl b/test/basic.jl index 43b5d3b..2f8bd52 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -3,7 +3,7 @@ using ImagesFiltering, OffsetArrays, Base.Test @testset "basic" begin @test eltype(KernelFactors.IIRGaussian(3)) == Float64 @test eltype(KernelFactors.IIRGaussian(Float32, 3)) == Float32 - @test isa(KernelFactors.IIRGaussian([1,2.0f0]), Tuple{KernelFactors.TriggsSdika,KernelFactors.TriggsSdika}) + @test isa(KernelFactors.IIRGaussian([1,2.0f0]), Tuple{KernelFactors.ReshapedVector,KernelFactors.ReshapedVector}) @test KernelFactors.kernelfactors(([0,3], [1,7])) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2)) @test KernelFactors.kernelfactors(([0,3], [1,7]')) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2)) From 9be7a257844c0c3ab7dfea853ab0fe3d4cef960c Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 23 Aug 2016 09:23:27 -0500 Subject: [PATCH 07/21] Fix and test LoG --- src/kernel.jl | 10 +++++----- test/specialty.jl | 9 +++++++++ 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/kernel.jl b/src/kernel.jl index 50a1028..9a77c71 100644 --- a/src/kernel.jl +++ b/src/kernel.jl @@ -162,13 +162,13 @@ function LoG{N}(σs::NTuple{N}) σ = SVector(σs) C = 1/(prod(σ)*(2π)^(N/2)) σ2 = σ.^2 - iσ4 = sum(1./σ2.^2) - function df(I::CartesianIndex, σ2, iσ4) + σ2i = sum(1./σ2) + function df(I::CartesianIndex, σ2, σ2i) x = SVector(I.I) - xσ = sum(x.^2./σ2) - (xσ - iσ4) * exp(-xσ/2) + xσ = x.^2./σ2 + (sum(xσ./σ2) - σ2i) * exp(-sum(xσ)/2) end - centered([C*df(I, σ2, iσ4) for I in R]) + centered([C*df(I, σ2, σ2i) for I in R]) end LoG(σ::Real) = LoG((σ,σ)) diff --git a/test/specialty.jl b/test/specialty.jl index 947d302..a3a098b 100644 --- a/test/specialty.jl +++ b/test/specialty.jl @@ -186,6 +186,15 @@ using Base.Test @test indices(Kernel.DoG((5,), (7,), (21,))) == (-10:10,) end + @testset "LoG" begin + img = rand(20,21) + σs = (2.5, 3.2) + kernel1 = Kernel.LoG(σs) + kernel2 = (KernelFactors.IIRGaussian(σs)..., Kernel.Laplacian()) + imgf1 = imfilter(img, kernel1) + imgf2 = imfilter(img, kernel2) + @test cor(vec(imgf1), vec(imgf2)) > 0.8 + end end nothing From 543f307398d0661b828e2bb0f2a8183bee7c8c8e Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 23 Aug 2016 09:49:00 -0500 Subject: [PATCH 08/21] Reorganize the tests a bit No actual tests changed --- test/2d.jl | 51 ------- test/border.jl | 343 ++++++++++++++++++++++++----------------------- test/runtests.jl | 1 + test/triggs.jl | 56 ++++++++ 4 files changed, 230 insertions(+), 221 deletions(-) create mode 100644 test/triggs.jl diff --git a/test/2d.jl b/test/2d.jl index 7e7adac..5765bc6 100644 --- a/test/2d.jl +++ b/test/2d.jl @@ -196,55 +196,4 @@ using Base.Test @test !(b === a) end -@testset "TriggsSdika 1d" begin - l = 1000 - x = [1, 2, 3, 5, 10, 20, l>>1, l-19, l-9, l-4, l-2, l-1, l] - σs = [3.1, 5, 10.0, 20, 50.0, 100.0] - # For plotting: - # aσ = Array{Any}(length(x)) - # for i = 1:length(x) - # aσ[i] = Array{Float64}(l, 2*length(σs)) - # end - for (n,σ) in enumerate(σs) - kernel = KernelFactors.IIRGaussian(σ) - A = [kernel.a'; eye(2) zeros(2)] - B = [kernel.b'; eye(2) zeros(2)] - I1 = zeros(3,3); I1[1,1] = 1 - @test_approx_eq kernel.M (I1+B*kernel.M*A) # Triggs & Sdika, Eq. 8 - for (i,c) in enumerate(x) - a = zeros(l) - a[c] = 1 - af = exp(-((1:l)-c).^2/(2*σ^2))/(σ*√(2π)) - @inferred(imfilter!(a, a, (kernel,), Fill(0))) - @test norm(a-af) < 0.1*norm(af) - # aσ[i][:,2n-1:2n] = [af a] - end - end -end - -@testset "TriggsSdika images" begin - imgf = zeros(5, 7); imgf[3,4] = 1 - imgg = fill(Gray{Float32}(0), 5, 7); imgg[3,4] = 1 - imgc = fill(RGB{Float64}(0,0,0), 5, 7); imgc[3,4] = RGB(1,0,0) - σ = 5 - x = -2:2 - y = (-3:3)' - kernel = KernelFactors.IIRGaussian((σ, σ)) - for img in (imgf, imgg, imgc) - imgcmp = img[3,4]*exp(-(x.^2 .+ y.^2)/(2*σ^2))/(σ^2*2*pi) - border = Fill(zero(eltype(img))) - img0 = copy(img) - imgf = @inferred(imfilter(img, kernel, border)) - @test sumabs2(imgcmp - imgf) < 0.2^2*sumabs2(imgcmp) - @test img == img0 - @test imgf != img - end - ret = imfilter(imgf, KernelFactors.IIRGaussian((0,σ)), Pad{:replicate}()) - @test ret != imgf - out = similar(ret) - imfilter!(CPU1(Algorithm.IIR()), out, imgf, KernelFactors.IIRGaussian(σ), 2, Pad{:replicate}()) - @test out == ret -end - - nothing diff --git a/test/border.jl b/test/border.jl index 06cfcb0..1720e9c 100644 --- a/test/border.jl +++ b/test/border.jl @@ -1,185 +1,188 @@ using ImagesFiltering, OffsetArrays, Colors using Base.Test -@testset "padarray" begin - A = reshape(1:25, 5, 5) - @test @inferred(padarray(A, Fill(0,(2,2),(2,2)))) == OffsetArray( - [0 0 0 0 0 0 0 0 0; - 0 0 0 0 0 0 0 0 0; - 0 0 1 6 11 16 21 0 0; - 0 0 2 7 12 17 22 0 0; - 0 0 3 8 13 18 23 0 0; - 0 0 4 9 14 19 24 0 0; - 0 0 5 10 15 20 25 0 0; - 0 0 0 0 0 0 0 0 0; - 0 0 0 0 0 0 0 0 0], (-2,-2)) - @test @inferred(padarray(A, Pad{:replicate}((2,2),(2,2)))) == OffsetArray( - [ 1 1 1 6 11 16 21 21 21; - 1 1 1 6 11 16 21 21 21; - 1 1 1 6 11 16 21 21 21; - 2 2 2 7 12 17 22 22 22; - 3 3 3 8 13 18 23 23 23; - 4 4 4 9 14 19 24 24 24; - 5 5 5 10 15 20 25 25 25; - 5 5 5 10 15 20 25 25 25; - 5 5 5 10 15 20 25 25 25], (-2,-2)) - @test @inferred(padarray(A, Pad{:circular}((2,2),(2,2)))) == OffsetArray( - [19 24 4 9 14 19 24 4 9; - 20 25 5 10 15 20 25 5 10; - 16 21 1 6 11 16 21 1 6; - 17 22 2 7 12 17 22 2 7; - 18 23 3 8 13 18 23 3 8; - 19 24 4 9 14 19 24 4 9; - 20 25 5 10 15 20 25 5 10; - 16 21 1 6 11 16 21 1 6; - 17 22 2 7 12 17 22 2 7], (-2,-2)) - @test @inferred(padarray(A, Pad{:symmetric}((2,2),(2,2)))) == OffsetArray( - [ 7 2 2 7 12 17 22 22 17; - 6 1 1 6 11 16 21 21 16; - 6 1 1 6 11 16 21 21 16; - 7 2 2 7 12 17 22 22 17; - 8 3 3 8 13 18 23 23 18; - 9 4 4 9 14 19 24 24 19; - 10 5 5 10 15 20 25 25 20; - 10 5 5 10 15 20 25 25 20; - 9 4 4 9 14 19 24 24 19], (-2,-2)) - @test @inferred(padarray(A, Pad{:reflect}((2,2),(2,2)))) == OffsetArray( - [13 8 3 8 13 18 23 18 13; - 12 7 2 7 12 17 22 17 12; - 11 6 1 6 11 16 21 16 11; - 12 7 2 7 12 17 22 17 12; - 13 8 3 8 13 18 23 18 13; - 14 9 4 9 14 19 24 19 14; - 15 10 5 10 15 20 25 20 15; - 14 9 4 9 14 19 24 19 14; - 13 8 3 8 13 18 23 18 13], (-2,-2)) - ret = @test_throws ArgumentError padarray(A, Fill(0)) - @test contains(ret.value.msg, "lacks the proper padding") - for Style in (:replicate, :circular, :symmetric, :reflect) - ret = @test_throws ArgumentError padarray(A, Pad{Style}((1,1,1),(1,1,1))) +@testset "Border" begin + @testset "padarray" begin + A = reshape(1:25, 5, 5) + @test @inferred(padarray(A, Fill(0,(2,2),(2,2)))) == OffsetArray( + [0 0 0 0 0 0 0 0 0; + 0 0 0 0 0 0 0 0 0; + 0 0 1 6 11 16 21 0 0; + 0 0 2 7 12 17 22 0 0; + 0 0 3 8 13 18 23 0 0; + 0 0 4 9 14 19 24 0 0; + 0 0 5 10 15 20 25 0 0; + 0 0 0 0 0 0 0 0 0; + 0 0 0 0 0 0 0 0 0], (-2,-2)) + @test @inferred(padarray(A, Pad{:replicate}((2,2),(2,2)))) == OffsetArray( + [ 1 1 1 6 11 16 21 21 21; + 1 1 1 6 11 16 21 21 21; + 1 1 1 6 11 16 21 21 21; + 2 2 2 7 12 17 22 22 22; + 3 3 3 8 13 18 23 23 23; + 4 4 4 9 14 19 24 24 24; + 5 5 5 10 15 20 25 25 25; + 5 5 5 10 15 20 25 25 25; + 5 5 5 10 15 20 25 25 25], (-2,-2)) + @test @inferred(padarray(A, Pad{:circular}((2,2),(2,2)))) == OffsetArray( + [19 24 4 9 14 19 24 4 9; + 20 25 5 10 15 20 25 5 10; + 16 21 1 6 11 16 21 1 6; + 17 22 2 7 12 17 22 2 7; + 18 23 3 8 13 18 23 3 8; + 19 24 4 9 14 19 24 4 9; + 20 25 5 10 15 20 25 5 10; + 16 21 1 6 11 16 21 1 6; + 17 22 2 7 12 17 22 2 7], (-2,-2)) + @test @inferred(padarray(A, Pad{:symmetric}((2,2),(2,2)))) == OffsetArray( + [ 7 2 2 7 12 17 22 22 17; + 6 1 1 6 11 16 21 21 16; + 6 1 1 6 11 16 21 21 16; + 7 2 2 7 12 17 22 22 17; + 8 3 3 8 13 18 23 23 18; + 9 4 4 9 14 19 24 24 19; + 10 5 5 10 15 20 25 25 20; + 10 5 5 10 15 20 25 25 20; + 9 4 4 9 14 19 24 24 19], (-2,-2)) + @test @inferred(padarray(A, Pad{:reflect}((2,2),(2,2)))) == OffsetArray( + [13 8 3 8 13 18 23 18 13; + 12 7 2 7 12 17 22 17 12; + 11 6 1 6 11 16 21 16 11; + 12 7 2 7 12 17 22 17 12; + 13 8 3 8 13 18 23 18 13; + 14 9 4 9 14 19 24 19 14; + 15 10 5 10 15 20 25 20 15; + 14 9 4 9 14 19 24 19 14; + 13 8 3 8 13 18 23 18 13], (-2,-2)) + ret = @test_throws ArgumentError padarray(A, Fill(0)) @test contains(ret.value.msg, "lacks the proper padding") - ret = @test_throws ArgumentError padarray(A, Pad{Style}) - @test contains(ret.value.msg, "must supply padding") - end - # arrays smaller than the padding - A = [1 2; 3 4] - @test @inferred(padarray(A, Pad{:replicate}((3,3),(3,3)))) == OffsetArray( - [ 1 1 1 1 2 2 2 2; - 1 1 1 1 2 2 2 2; - 1 1 1 1 2 2 2 2; - 1 1 1 1 2 2 2 2; - 3 3 3 3 4 4 4 4; - 3 3 3 3 4 4 4 4; - 3 3 3 3 4 4 4 4; - 3 3 3 3 4 4 4 4], (-3,-3)) - @test @inferred(padarray(A, Pad{:circular}((3,3),(3,3)))) == OffsetArray( - [ 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1; - 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1; - 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1; - 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1], (-3,-3)) - @test @inferred(padarray(A, Pad{:symmetric}((3,3),(3,3)))) == OffsetArray( - [ 4 4 3 3 4 4 3 3; - 4 4 3 3 4 4 3 3; - 2 2 1 1 2 2 1 1; - 2 2 1 1 2 2 1 1; - 4 4 3 3 4 4 3 3; - 4 4 3 3 4 4 3 3; - 2 2 1 1 2 2 1 1; - 2 2 1 1 2 2 1 1], (-3,-3)) - @test @inferred(padarray(A, Pad{:reflect}((3,3),(3,3)))) == OffsetArray( - [ 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1; - 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1; - 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1; - 4 3 4 3 4 3 4 3; - 2 1 2 1 2 1 2 1], (-3,-3)) - for Style in (:replicate, :circular, :symmetric, :reflect) - @test @inferred(padarray(A, Pad{Style}((0,0), (0,0)))) == A + for Style in (:replicate, :circular, :symmetric, :reflect) + ret = @test_throws ArgumentError padarray(A, Pad{Style}((1,1,1),(1,1,1))) + @test contains(ret.value.msg, "lacks the proper padding") + ret = @test_throws ArgumentError padarray(A, Pad{Style}) + @test contains(ret.value.msg, "must supply padding") + end + # arrays smaller than the padding + A = [1 2; 3 4] + @test @inferred(padarray(A, Pad{:replicate}((3,3),(3,3)))) == OffsetArray( + [ 1 1 1 1 2 2 2 2; + 1 1 1 1 2 2 2 2; + 1 1 1 1 2 2 2 2; + 1 1 1 1 2 2 2 2; + 3 3 3 3 4 4 4 4; + 3 3 3 3 4 4 4 4; + 3 3 3 3 4 4 4 4; + 3 3 3 3 4 4 4 4], (-3,-3)) + @test @inferred(padarray(A, Pad{:circular}((3,3),(3,3)))) == OffsetArray( + [ 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1; + 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1; + 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1; + 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1], (-3,-3)) + @test @inferred(padarray(A, Pad{:symmetric}((3,3),(3,3)))) == OffsetArray( + [ 4 4 3 3 4 4 3 3; + 4 4 3 3 4 4 3 3; + 2 2 1 1 2 2 1 1; + 2 2 1 1 2 2 1 1; + 4 4 3 3 4 4 3 3; + 4 4 3 3 4 4 3 3; + 2 2 1 1 2 2 1 1; + 2 2 1 1 2 2 1 1], (-3,-3)) + @test @inferred(padarray(A, Pad{:reflect}((3,3),(3,3)))) == OffsetArray( + [ 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1; + 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1; + 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1; + 4 3 4 3 4 3 4 3; + 2 1 2 1 2 1 2 1], (-3,-3)) + for Style in (:replicate, :circular, :symmetric, :reflect) + @test @inferred(padarray(A, Pad{Style}((0,0), (0,0)))) == A + end + @test @inferred(padarray(A, Fill(0, (0,0), (0,0)))) == A + @test @inferred(padarray(A, Pad{:replicate}((1,2), (2,0)))) == OffsetArray([1 1 1 2; 1 1 1 2; 3 3 3 4; 3 3 3 4; 3 3 3 4], (-1,-2)) + @test @inferred(padarray(A, Pad{:circular}((2,1), (0,2)))) == OffsetArray([2 1 2 1 2; 4 3 4 3 4; 2 1 2 1 2; 4 3 4 3 4], (-2,-1)) + @test @inferred(padarray(A, Pad{:symmetric}((1,2), (2,0)))) == OffsetArray([2 1 1 2; 2 1 1 2; 4 3 3 4; 4 3 3 4; 2 1 1 2], (-1,-2)) + @test @inferred(padarray(A, Fill(-1,(1,2), (2,0)))) == OffsetArray([-1 -1 -1 -1; -1 -1 1 2; -1 -1 3 4; -1 -1 -1 -1; -1 -1 -1 -1], (-1,-2)) + A = [1 2 3; 4 5 6] + @test @inferred(padarray(A, Pad{:reflect}((1,2), (2,0)))) == OffsetArray([6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6], (-1,-2)) + A = [1 2; 3 4] + @test @inferred(padarray(A, Pad(1,1))) == OffsetArray([1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4], (-1,-1)) + @test @inferred(padarray(A, Pad{:circular}((1,1), ()))) == OffsetArray([4 3 4; 2 1 2; 4 3 4], (-1,-1)) + @test @inferred(padarray(A, Pad{:symmetric}((), (1,1)))) == [1 2 2; 3 4 4; 3 4 4] + A = ["a" "b"; "c" "d"] + @test @inferred(padarray(A, Pad(1,1))) == OffsetArray(["a" "a" "b" "b"; "a" "a" "b" "b"; "c" "c" "d" "d"; "c" "c" "d" "d"], (-1,-1)) + @test_throws MethodError padarray(A, Pad{:unknown}(1,1)) + A = trues(3,3) + @test isa(parent(padarray(A, Pad((1,2), (2,1)))), BitArray{2}) + A = falses(10,10,10) + B = view(A,1:8,1:8,1:8) + @test isa(parent(padarray(A, Pad((1,1,1)))), BitArray{3}) + @test isa(parent(padarray(B, Pad((1,1,1)))), BitArray{3}) + A = reshape(1:15, 3, 5) + B = @inferred(padarray(A, Inner((1,1)))) + @test B == A + # test that it's a copy + B[1,1] = 0 + @test B != A + A = rand(RGB{U8}, 3, 5) + ret = @test_throws ErrorException padarray(A, Fill(0, (0,0), (0,0))) + @test contains(ret.value.msg, "element type ColorTypes.RGB") + A = bitrand(3, 5) + ret = @test_throws ErrorException padarray(A, Fill(7, (0,0), (0,0))) + @test contains(ret.value.msg, "element type Bool") + @test isa(parent(padarray(A, Fill(false, (1,1), (1,1)))), BitArray) end - @test @inferred(padarray(A, Fill(0, (0,0), (0,0)))) == A - @test @inferred(padarray(A, Pad{:replicate}((1,2), (2,0)))) == OffsetArray([1 1 1 2; 1 1 1 2; 3 3 3 4; 3 3 3 4; 3 3 3 4], (-1,-2)) - @test @inferred(padarray(A, Pad{:circular}((2,1), (0,2)))) == OffsetArray([2 1 2 1 2; 4 3 4 3 4; 2 1 2 1 2; 4 3 4 3 4], (-2,-1)) - @test @inferred(padarray(A, Pad{:symmetric}((1,2), (2,0)))) == OffsetArray([2 1 1 2; 2 1 1 2; 4 3 3 4; 4 3 3 4; 2 1 1 2], (-1,-2)) - @test @inferred(padarray(A, Fill(-1,(1,2), (2,0)))) == OffsetArray([-1 -1 -1 -1; -1 -1 1 2; -1 -1 3 4; -1 -1 -1 -1; -1 -1 -1 -1], (-1,-2)) - A = [1 2 3; 4 5 6] - @test @inferred(padarray(A, Pad{:reflect}((1,2), (2,0)))) == OffsetArray([6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6], (-1,-2)) - A = [1 2; 3 4] - @test @inferred(padarray(A, Pad(1,1))) == OffsetArray([1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4], (-1,-1)) - @test @inferred(padarray(A, Pad{:circular}((1,1), ()))) == OffsetArray([4 3 4; 2 1 2; 4 3 4], (-1,-1)) - @test @inferred(padarray(A, Pad{:symmetric}((), (1,1)))) == [1 2 2; 3 4 4; 3 4 4] - A = ["a" "b"; "c" "d"] - @test @inferred(padarray(A, Pad(1,1))) == OffsetArray(["a" "a" "b" "b"; "a" "a" "b" "b"; "c" "c" "d" "d"; "c" "c" "d" "d"], (-1,-1)) - @test_throws MethodError padarray(A, Pad{:unknown}(1,1)) - A = trues(3,3) - @test isa(parent(padarray(A, Pad((1,2), (2,1)))), BitArray{2}) - A = falses(10,10,10) - B = view(A,1:8,1:8,1:8) - @test isa(parent(padarray(A, Pad((1,1,1)))), BitArray{3}) - @test isa(parent(padarray(B, Pad((1,1,1)))), BitArray{3}) - A = reshape(1:15, 3, 5) - B = @inferred(padarray(A, Inner((1,1)))) - @test B == A - # test that it's a copy - B[1,1] = 0 - @test B != A - A = rand(RGB{U8}, 3, 5) - ret = @test_throws ErrorException padarray(A, Fill(0, (0,0), (0,0))) - @test contains(ret.value.msg, "element type ColorTypes.RGB") - A = bitrand(3, 5) - ret = @test_throws ErrorException padarray(A, Fill(7, (0,0), (0,0))) - @test contains(ret.value.msg, "element type Bool") - @test isa(parent(padarray(A, Fill(false, (1,1), (1,1)))), BitArray) -end -@testset "Pad" begin - @test Pad{:replicate}([1,2], [5,3]) == Pad{:replicate}((1,2), (5,3)) - @test @inferred(Pad{:replicate,2}([1,2], [5,3])) == Pad{:replicate}((1,2), (5,3)) - @test_throws TypeError Pad{:replicate,3}([1,2], [5,3]) - @test @inferred(Pad{:circular}()(rand(3,5))) == Pad{:circular,2}((0,0),(3,5)) - @test @inferred(Pad{:circular}()(centered(rand(3,5)))) == Pad{:circular,2}((1,2),(1,2)) - @test @inferred(Pad{:symmetric}()(Kernel.Laplacian())) == Pad{:symmetric,2}((1,1),(1,1)) - @test @inferred(Pad{:symmetric}()(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Pad{:symmetric,2}((1,1),(1,1)) + @testset "Pad" begin + @test Pad{:replicate}([1,2], [5,3]) == Pad{:replicate}((1,2), (5,3)) + @test @inferred(Pad{:replicate,2}([1,2], [5,3])) == Pad{:replicate}((1,2), (5,3)) + @test_throws TypeError Pad{:replicate,3}([1,2], [5,3]) + @test @inferred(Pad{:circular}()(rand(3,5))) == Pad{:circular,2}((0,0),(3,5)) + @test @inferred(Pad{:circular}()(centered(rand(3,5)))) == Pad{:circular,2}((1,2),(1,2)) + @test @inferred(Pad{:symmetric}()(Kernel.Laplacian())) == Pad{:symmetric,2}((1,1),(1,1)) + @test @inferred(Pad{:symmetric}()(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Pad{:symmetric,2}((1,1),(1,1)) - a = reshape(1:15, 3, 5) - targetinds = (OffsetArray([1,1,2,3,3], 0:4), OffsetArray([1:5;], 1:5)) - ERROLD = STDERR - tmpfile = tempname() - try - open(tmpfile, "w") do f - redirect_stderr(f) - @test @inferred(ImagesFiltering.padindices(rand(3,5), Pad{:replicate}((1,)))) == targetinds + a = reshape(1:15, 3, 5) + targetinds = (OffsetArray([1,1,2,3,3], 0:4), OffsetArray([1:5;], 1:5)) + ERROLD = STDERR + tmpfile = tempname() + try + open(tmpfile, "w") do f + redirect_stderr(f) + @test @inferred(ImagesFiltering.padindices(rand(3,5), Pad{:replicate}((1,)))) == targetinds + end + finally + redirect_stderr(ERROLD) end - finally - redirect_stderr(ERROLD) + @test contains(readstring(tmpfile), "lacks the proper padding sizes") end - @test contains(readstring(tmpfile), "lacks the proper padding sizes") -end -@testset "Inner" begin - @test @inferred(Inner((1,2))) == Inner((1,2), (1,2)) - @test @inferred(Inner((1,2), ())) == Inner((1,2), (0,0)) - @test @inferred(Inner((), (1,2))) == Inner((0,0), (1,2)) - @test @inferred(Inner{2}([1,2],[5,6])) == Inner((1,2), (5,6)) - @test Inner([1,2],[5,6]) == Inner((1,2), (5,6)) - @test @inferred(Inner()(rand(3,5))) == Inner((0,0),(3,5)) -end -@testset "Fill" begin - @test Fill(1,[1,2],[5,6]) == Fill(1, (1,2), (5,6)) - @test @inferred(Fill(-1, rand(3,5))) == Fill(-1, (0,0), (3,5)) - @test @inferred(Fill(2)(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Fill(2, (1,1),(1,1)) -end + @testset "Inner" begin + @test @inferred(Inner((1,2))) == Inner((1,2), (1,2)) + @test @inferred(Inner((1,2), ())) == Inner((1,2), (0,0)) + @test @inferred(Inner((), (1,2))) == Inner((0,0), (1,2)) + @test @inferred(Inner{2}([1,2],[5,6])) == Inner((1,2), (5,6)) + @test Inner([1,2],[5,6]) == Inner((1,2), (5,6)) + @test @inferred(Inner()(rand(3,5))) == Inner((0,0),(3,5)) + end -@testset "misc" begin - a0 = reshape([1]) # 0-dimensional - @test ImagesFiltering.accumulate_padding((), a0) == () - @test ImagesFiltering.accumulate_padding((0:1, -1:1), a0) == (0:1, -1:1) + @testset "Fill" begin + @test Fill(1,[1,2],[5,6]) == Fill(1, (1,2), (5,6)) + @test @inferred(Fill(-1, rand(3,5))) == Fill(-1, (0,0), (3,5)) + @test @inferred(Fill(2)(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Fill(2, (1,1),(1,1)) + end + + @testset "misc" begin + a0 = reshape([1]) # 0-dimensional + @test ImagesFiltering.accumulate_padding((), a0) == () + @test ImagesFiltering.accumulate_padding((0:1, -1:1), a0) == (0:1, -1:1) + end end nothing diff --git a/test/runtests.jl b/test/runtests.jl index 23981e9..ceea667 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ asa = detect_ambiguities(StaticArrays, Base) include("border.jl") include("2d.jl") +include("triggs.jl") include("cascade.jl") include("specialty.jl") include("basic.jl") diff --git a/test/triggs.jl b/test/triggs.jl new file mode 100644 index 0000000..a840b7d --- /dev/null +++ b/test/triggs.jl @@ -0,0 +1,56 @@ +using ImagesFiltering, Colors, ComputationalResources +using Base.Test + +@testset "TriggsSdika" begin + @testset "TriggsSdika 1d" begin + l = 1000 + x = [1, 2, 3, 5, 10, 20, l>>1, l-19, l-9, l-4, l-2, l-1, l] + σs = [3.1, 5, 10.0, 20, 50.0, 100.0] + # For plotting: + # aσ = Array{Any}(length(x)) + # for i = 1:length(x) + # aσ[i] = Array{Float64}(l, 2*length(σs)) + # end + for (n,σ) in enumerate(σs) + kernel = KernelFactors.IIRGaussian(σ) + A = [kernel.a'; eye(2) zeros(2)] + B = [kernel.b'; eye(2) zeros(2)] + I1 = zeros(3,3); I1[1,1] = 1 + @test_approx_eq kernel.M (I1+B*kernel.M*A) # Triggs & Sdika, Eq. 8 + for (i,c) in enumerate(x) + a = zeros(l) + a[c] = 1 + af = exp(-((1:l)-c).^2/(2*σ^2))/(σ*√(2π)) + @inferred(imfilter!(a, a, (kernel,), Fill(0))) + @test norm(a-af) < 0.1*norm(af) + # aσ[i][:,2n-1:2n] = [af a] + end + end + end + + @testset "TriggsSdika images" begin + imgf = zeros(5, 7); imgf[3,4] = 1 + imgg = fill(Gray{Float32}(0), 5, 7); imgg[3,4] = 1 + imgc = fill(RGB{Float64}(0,0,0), 5, 7); imgc[3,4] = RGB(1,0,0) + σ = 5 + x = -2:2 + y = (-3:3)' + kernel = KernelFactors.IIRGaussian((σ, σ)) + for img in (imgf, imgg, imgc) + imgcmp = img[3,4]*exp(-(x.^2 .+ y.^2)/(2*σ^2))/(σ^2*2*pi) + border = Fill(zero(eltype(img))) + img0 = copy(img) + imgf = @inferred(imfilter(img, kernel, border)) + @test sumabs2(imgcmp - imgf) < 0.2^2*sumabs2(imgcmp) + @test img == img0 + @test imgf != img + end + ret = imfilter(imgf, KernelFactors.IIRGaussian((0,σ)), Pad{:replicate}()) + @test ret != imgf + out = similar(ret) + imfilter!(CPU1(Algorithm.IIR()), out, imgf, KernelFactors.IIRGaussian(σ), 2, Pad{:replicate}()) + @test out == ret + end +end + +nothing From d41f8ac661dd6da3b7a3ca8a23c424e2298a95e8 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 24 Aug 2016 07:52:23 -0500 Subject: [PATCH 09/21] Track the interior during successive stages of filtering Previously, filter cascades would result in undefined values on the edges of intermediate states. This tracks the edge size at each stage of the cascade to ensure that no invalid values are used. --- src/border.jl | 10 ++++- src/imfilter.jl | 91 ++++++++++++++++++++++++++------------------ src/kernelfactors.jl | 5 ++- src/specialty.jl | 3 +- test/2d.jl | 7 ++-- test/specialty.jl | 5 +++ 6 files changed, 77 insertions(+), 44 deletions(-) diff --git a/src/border.jl b/src/border.jl index 962fd60..11b2bb0 100644 --- a/src/border.jl +++ b/src/border.jl @@ -282,8 +282,14 @@ interior(A::AbstractArray, kernel::Union{ArrayLike,Laplacian}) = _interior(indic interior(A, factkernel::Tuple) = _interior(indices(A), accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) function _interior{N}(indsA::NTuple{N}, indsk) indskN = fill_to_length(indsk, 0:0, Val{N}) - CartesianRange(CartesianIndex(map((ia,ik)->first(ia) + lo(ik), indsA, indskN)), - CartesianIndex(map((ia,ik)->last(ia) - hi(ik), indsA, indskN))) + map((ia,ik)->first(ia) + lo(ik) : last(ia) - hi(ik), indsA, indskN) +end + +next_interior(inds::Indices, ::Tuple{}) = inds +function next_interior(inds::Indices, kernel::Tuple) + kern = first(kernel) + iscopy(kern) && return next_interior(inds, tail(kernel)) + _interior(inds, indices(kern)) end # end diff --git a/src/imfilter.jl b/src/imfilter.jl index 5c0d6f3..cb94478 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -26,8 +26,7 @@ end # Step 4: if necessary, allocate the ouput @inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::Inner{0}, args...) - R = interior(img, kernel) - inds = to_ranges(R) + inds = interior(img, kernel) imfilter!(similar(img, T, inds), img, kernel, border, args...) end @inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny, args...) @@ -49,8 +48,7 @@ function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel: imfilter(r, T, img, kernel, Pad{:replicate}()) # supply the default border end function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::Inner{0}) - R = interior(img, kernel) - inds = to_ranges(R) + inds = interior(img, kernel) imfilter!(r, similar(img, T, inds), img, kernel, border) end function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny) @@ -218,8 +216,8 @@ function imfilter!{S,T,N}(r::AbstractResource, imfilter!(r, out, A, kernel, NoPad(border)) end -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{}, ::NoPad) - R = CartesianRange(indices(out)) +function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{}, ::NoPad, inds::Indices=indices(out)) + R = CartesianRange(inds) copy!(out, R, A, R) end @@ -229,6 +227,8 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke imfilter!(r, out, A, samedims(out, kern), border) end +const fillbuf_nan = Ref(false) + function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) kern = kernel[1] iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border) @@ -236,30 +236,39 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke # and swap them at each stage. The first of the two is the one # that holds the most recent result. A2 = similar(A) # for type-stability, let's hope it's really the *same* type... - _imfilter!(r, out, A, A2, kernel, border) + if fillbuf_nan[] + fill!(A2, NaN) # for testing purposes + end + inds = interior(r, A, kern) + _imfilter!(r, out, A, A2, kernel, border, inds) return out end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, border::NoPad) - imfilter!(r, out, A1, kernel, border) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, border::NoPad, inds::Indices) + imfilter!(r, out, A1, kernel, border, inds) end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, border::NoPad) - imfilter!(r, out, A1, kernel, border) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, border::NoPad, inds::Indices) + imfilter!(r, out, A1, kernel[1], border, inds) end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, inds::Indices) kern = kernel[1] - iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), border) + iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), border, inds) kernN = samedims(A2, kern) - imfilter!(r, A2, A1, kernN, border, interior(r, A2, kernN)) # store result in A2 - _imfilter!(r, out, A2, A1, tail(kernel), border) # swap the buffers + imfilter!(r, A2, A1, kernN, border, inds) # store result in A2 + if fillbuf_nan[] + @show A2 + end + kernelt = tail(kernel) + newinds = next_interior(inds, kernelt) + _imfilter!(r, out, A2, A1, tail(kernel), border, newinds) # swap the buffers end ## FIR filtering """ - imfilter!(::AbstractResource{FIR}, imgfilt, img, kernel, NoPad(), [R=CartesianRange(indices(imfilt))]) + imfilter!(::AbstractResource{FIR}, imgfilt, img, kernel, NoPad(), [inds=indices(imgfilt)]) Filter an array `img` with kernel `kernel` by computing their correlation, storing the result in `imgfilt`, using a finite-impulse @@ -274,11 +283,11 @@ with a specific `border`, or use for default padding. -If the CartesianRange `R` is supplied, only the elements of `imgfilt` -with indices in the domain of `R` will be calculated. This can be -particularly useful for "cascaded filters" where you pad over a larger -area and then calculate the result over just the necessary region at -each stage. +If `inds` is supplied, only the elements of `imgfilt` with indices in +the domain of `inds` will be calculated. This can be particularly +useful for "cascaded filters" where you pad over a larger area and +then calculate the result over just the necessary region at each +stage. See also: imfilter. """ @@ -287,42 +296,45 @@ function imfilter!{S,T,N}(r::CPU1{FIR}, A::AbstractArray{T,N}, kern::NDimKernel{N}, border::NoPad, - R::CartesianRange=CartesianRange(indices(out))) + inds::Indices{N}=indices(out)) (isempty(A) || isempty(kern)) && return out indso, indsA, indsk = indices(out), indices(A), indices(kern) if iscopy(kern) return copy!(out, R, A, R) end for i = 1:N - # Check that R is inbounds for out - if R.start[i] < first(indso[i]) || R.stop[i] > last(indso[i]) - throw(DimensionMismatch("output indices $indso disagrees with requested range $R")) + # Check that inds is inbounds for out + indsi, indsoi, indsAi, indski = inds[i], indso[i], indsA[i], indsk[i] + if first(indsi) < first(indsoi) || last(indsi) > last(indsoi) + throw(DimensionMismatch("output indices $indso disagrees with requested indices $inds")) end # Check that input A is big enough not to throw a BoundsError - if first(indsA[i]) > R.start[i] + first(indsk[i]) || - last(indsA[i]) < R.stop[i] + last(indsk[i]) - throw(DimensionMismatch("requested range $R and kernel indices $indsk do not agree with indices of padded input, $indsA")) + if first(indsAi) > first(indsi) + first(indski) || + last(indsA[i]) < last(indsi) + last(indski) + throw(DimensionMismatch("requested indices $inds and kernel indices $indsk do not agree with indices of padded input, $indsA")) end end - # For performance reasons, we now dispatch to a method that may strip - # off index-shift containers (in Julia 0.5 we shouldn't remove - # boundschecks from OffsetArray). - _imfilter_inbounds!(r, out, A, kern, border, R) + _imfilter_inbounds!(r, out, A, kern, border, inds) +end + +function _imfilter_inbounds!{T,N,Npre,V<:TriggsSdika}(r::AbstractResource, out, A::AbstractArray, kern::ReshapedVector{T,N,Npre,V}, border::NoPad, inds) + indspre, ind, indspost = iterdims(inds, kern) + _imfilter_dim!(r, out, A, kern.data, CartesianRange(indspre), ind, CartesianRange(indspost), border[]) end -function _imfilter_inbounds!{T,N,Npre,V<:TriggsSdika}(r::AbstractResource, out, A::AbstractArray, kern::ReshapedVector{T,N,Npre,V}, border::NoPad, R) - indspre, ind, indspost = iterdims(ranges(R), kern) +function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractVector, kern::TriggsSdika, border::NoPad, inds) + indspre, ind, indspost = iterdims(inds, kern) _imfilter_dim!(r, out, A, kern.data, CartesianRange(indspre), ind, CartesianRange(indspost), border[]) end -function _imfilter_inbounds!(r, out, A::AbstractArray, kern, border::NoPad, R) +function _imfilter_inbounds!(r, out, A::AbstractArray, kern, border::NoPad, inds) indsk = indices(kern) - Rk = CartesianRange(indsk) + R, Rk = CartesianRange(inds), CartesianRange(indsk) p = A[first(R)+first(Rk)] * first(kern) TT = typeof(p+p) for I in R tmp = zero(TT) - @unsafe for J in CartesianRange(indsk) + @unsafe for J in Rk tmp += A[I+J]*kern[J] end @unsafe out[I] = tmp @@ -463,6 +475,11 @@ function imfilter!(r::AbstractResource, out, img, kernel::TriggsSdika, dim::Inte _imfilter_dim!(r, out, img, kernel, Rbegin, inds[dim], Rend, border) end +function imfilter!(r::AbstractResource, out, A::AbstractVector, kern::TriggsSdika, border::NoPad, inds::Indices=indices(out)) + indspre, ind, indspost = iterdims(inds, kern) + _imfilter_dim!(r, out, A, kern, CartesianRange(indspre), ind, CartesianRange(indspost), border[]) +end + # Lispy and type-stable inplace (currently just Triggs-Sdika) filtering over each dimension function _imfilter_inplace_tuple!(r, out, img, kernel, Rbegin, inds, Rend, border) ind = first(inds) diff --git a/src/kernelfactors.jl b/src/kernelfactors.jl index 251f32f..8bac143 100644 --- a/src/kernelfactors.jl +++ b/src/kernelfactors.jl @@ -222,9 +222,12 @@ Base.vec(kernel::TriggsSdika) = kernel Base.ndims(kernel::TriggsSdika) = 1 Base.ndims{T<:TriggsSdika}(::Type{T}) = 1 Base.indices1(kernel::TriggsSdika) = 0:0 -Base.indices(kernel::TriggsSdika) = (indices1(kernel),) +Base.indices(kernel::TriggsSdika) = (Base.indices1(kernel),) Base.isempty(kernel::TriggsSdika) = false +iterdims(inds::Indices{1}, kern::TriggsSdika) = (), inds[1], () +_reshape(kern::TriggsSdika, ::Type{Val{1}}) = kern + # Note that there's a sign reversal between Young & Triggs. """ IIRGaussian([T], σ; emit_warning::Bool=true) diff --git a/src/specialty.jl b/src/specialty.jl index 867a79d..c0a114b 100644 --- a/src/specialty.jl +++ b/src/specialty.jl @@ -1,8 +1,9 @@ ## Laplacian -function _imfilter_inbounds!(r, out, A::AbstractArray, L::Laplacian, border::NoPad, R) +function _imfilter_inbounds!(r, out, A::AbstractArray, L::Laplacian, border::NoPad, inds) TT = eltype(out) # accumtype(eltype(out), eltype(A)) n = 2*length(L.offsets) + R = CartesianRange(inds) @unsafe for I in R tmp = convert(TT, - n * A[I]) for J in L.offsets diff --git a/test/2d.jl b/test/2d.jl index 5765bc6..6a3b631 100644 --- a/test/2d.jl +++ b/test/2d.jl @@ -12,21 +12,22 @@ using Base.Test imgi = zeros(Int, 5, 7); imgi[3,4] = 1 imgg = fill(Gray(0), 5, 7); imgg[3,4] = 1 imgc = fill(RGB(0,0,0), 5, 7); imgc[3,4] = RGB(1,0,0) - # Dense kernel + # Dense inseparable kernel kern = [0.1 0.2; 0.4 0.5] kernel = OffsetArray(kern, -1:0, 1:2) for img in (imgf, imgi, imgg, imgc) targetimg = zeros(typeof(img[1]*kern[1]), size(img)) targetimg[3:4,2:3] = rot180(kern)*img[3,4] - @test (ret = @inferred(imfilter(img, kernel))) ≈ targetimg + @test @inferred(imfilter(img, kernel)) ≈ targetimg @test @inferred(imfilter(img, (kernel,))) ≈ targetimg @test @inferred(imfilter(f32type(img), img, kernel)) ≈ float32(targetimg) + ret = imfilter(img, kernel) fill!(ret, zero(eltype(ret))) @test @inferred(imfilter!(ret, img, kernel)) ≈ targetimg fill!(ret, zero(eltype(ret))) @test @inferred(imfilter!(CPU1(Algorithm.FIR()), ret, img, kernel)) ≈ targetimg for border in (Pad{:replicate}(), Pad{:circular}(), Pad{:symmetric}(), Pad{:reflect}(), Fill(zero(eltype(img)))) - @test (ret = @inferred(imfilter(img, kernel, border))) ≈ targetimg + @test @inferred(imfilter(img, kernel, border)) ≈ targetimg @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32(targetimg) fill!(ret, zero(eltype(ret))) @test @inferred(imfilter!(ret, img, kernel, border)) ≈ targetimg diff --git a/test/specialty.jl b/test/specialty.jl index a3a098b..8f23367 100644 --- a/test/specialty.jl +++ b/test/specialty.jl @@ -194,6 +194,11 @@ using Base.Test imgf1 = imfilter(img, kernel1) imgf2 = imfilter(img, kernel2) @test cor(vec(imgf1), vec(imgf2)) > 0.8 + # Ensure that edge-trimming under successive stages of filtering works correctly + ImagesFiltering.fillbuf_nan[] = true + kernel3 = (Kernel.Laplacian(), KernelFactors.IIRGaussian(σs)...) + @test !any(isnan, imfilter(img, kernel3)) + ImagesFiltering.fillbuf_nan[] = false end end From 7b5baf8c27803e2db7383e5b624780e72a102e23 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 24 Aug 2016 09:48:51 -0500 Subject: [PATCH 10/21] A few bug fixes and more tests --- src/border.jl | 3 ++- src/imfilter.jl | 6 ++--- test/border.jl | 6 +++++ test/nd.jl | 57 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/triggs.jl | 4 ++-- 6 files changed, 70 insertions(+), 7 deletions(-) create mode 100644 test/nd.jl diff --git a/src/border.jl b/src/border.jl index 11b2bb0..6468948 100644 --- a/src/border.jl +++ b/src/border.jl @@ -165,7 +165,7 @@ end padarray(img, border::Inner) = padarray(eltype(img), img, border) padarray{T}(::Type{T}, img::AbstractArray{T}, border::Inner) = copy(img) -padarray{T}(::Type{T}, img::AbstractArray, border::Inner) = convert(Array{T}, img) +padarray{T}(::Type{T}, img::AbstractArray, border::Inner) = copy!(similar(Array{T}, indices(img)), img) """ Fill(val) @@ -192,6 +192,7 @@ Fill(value, kernel::AbstractArray) = Fill(value, indices(kernel)) Fill(value, kernel::Laplacian) = Fill(value, indices(kernel)) Fill(value, factkernel::Tuple) = Fill(value, accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) +(p::Fill)(kernel) = Fill(p.value, kernel) (p::Fill)(kernel, img, ::FIR) = Fill(p.value, kernel) function (p::Fill)(factkernel::Tuple, img, ::FFT) inds = accumulate_padding(indices(factkernel[1]), tail(factkernel)...) diff --git a/src/imfilter.jl b/src/imfilter.jl index cb94478..3ca2b4e 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -208,7 +208,7 @@ function imfilter!{S,T,N}(r::AbstractResource, img::AbstractArray{T,N}, kernel::ProcessedKernel, border::BorderSpec) - bord = border(kernel, img, Alg(r)) + bord = border(kernel, img, Alg(r)) # if it's FFT, the size of img is also relevant A = padarray(S, img, bord) # By specifying NoPad(), we ensure that dispatch will never # accidentally "go back" to an earlier routine and apply more @@ -257,9 +257,6 @@ function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{ iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), border, inds) kernN = samedims(A2, kern) imfilter!(r, A2, A1, kernN, border, inds) # store result in A2 - if fillbuf_nan[] - @show A2 - end kernelt = tail(kernel) newinds = next_interior(inds, kernelt) _imfilter!(r, out, A2, A1, tail(kernel), border, newinds) # swap the buffers @@ -300,6 +297,7 @@ function imfilter!{S,T,N}(r::CPU1{FIR}, (isempty(A) || isempty(kern)) && return out indso, indsA, indsk = indices(out), indices(A), indices(kern) if iscopy(kern) + R = CartesianRange(inds) return copy!(out, R, A, R) end for i = 1:N diff --git a/test/border.jl b/test/border.jl index 1720e9c..50fa0cd 100644 --- a/test/border.jl +++ b/test/border.jl @@ -170,11 +170,17 @@ using Base.Test @test @inferred(Inner{2}([1,2],[5,6])) == Inner((1,2), (5,6)) @test Inner([1,2],[5,6]) == Inner((1,2), (5,6)) @test @inferred(Inner()(rand(3,5))) == Inner((0,0),(3,5)) + @test @inferred(Inner()(centered(rand(3,5)))) == Inner((1,2),(1,2)) end @testset "Fill" begin @test Fill(1,[1,2],[5,6]) == Fill(1, (1,2), (5,6)) + @test @inferred(Fill(-1)(rand(3,5))) == Fill(-1, (0,0), (3,5)) + @test @inferred(Fill(-1)(centered(rand(3,5)))) == Fill(-1, (1,2), (1,2)) @test @inferred(Fill(-1, rand(3,5))) == Fill(-1, (0,0), (3,5)) + @test @inferred(Fill(-1, centered(rand(3,5)))) == Fill(-1, (1,2), (1,2)) + @test @inferred(Fill(-1)(rand(3,5))) == Fill(-1, (0,0), (3,5)) + @test @inferred(Fill(-1)(centered(rand(3,5)))) == Fill(-1, (1,2), (1,2)) @test @inferred(Fill(2)(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Fill(2, (1,1),(1,1)) end diff --git a/test/nd.jl b/test/nd.jl new file mode 100644 index 0000000..b04bcde --- /dev/null +++ b/test/nd.jl @@ -0,0 +1,57 @@ +using ImagesFiltering, OffsetArrays +using Base.Test + +@testset "1d" begin + # Cascades don't result in out-of-bounds values + img = 1:8 + k1 = centered([0.25, 0.5, 0.25]) + k2 = OffsetArray([0.5, 0.5], 1:2) + casc = imfilter(img, (k1, k2, k1)) + A0 = padarray(img, Pad{:replicate}((2,),(4,))) + A1 = imfilter(A0, k1, Inner()) + @test_approx_eq A1 OffsetArray([1.0,1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75,8.0,8.0,8.0], 0:11) + A2 = imfilter(A1, k2, Inner()) + @test_approx_eq A2 OffsetArray([1.625,2.5,3.5,4.5,5.5,6.5,7.375,7.875,8.0,8.0], 0:9) + A3 = imfilter(A2, k1, Inner()) + @test_approx_eq casc A3 + @test size(casc) == size(img) + # copy! kernels, presence/order doesn't matter + kc = centered([1]) + @test ImagesFiltering.iscopy(kc) + @test_approx_eq imfilter(img, (k1,)) [1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75] + @test_approx_eq imfilter(img, (kc, k1)) [1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75] + @test_approx_eq imfilter(img, (k1, kc)) [1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75] +end + +@testset "3d" begin + img = trues(10,10,10) + kernel = centered(trues(3,3,3)/27) + for border in (Pad{:replicate}(), Pad{:circular}(), Pad{:symmetric}(), Pad{:reflect}(), Fill(true)) + for alg in (Algorithm.FIR(), Algorithm.FFT()) + @test_approx_eq imfilter(img, kernel, border) img + end + end + # Fill(0) + target = convert(Array{Float64}, img) + for i in (1,10) + target[:,:,i] = target[:,i,:] = target[i,:,:] = 2/3 + end + for i in (1,10) + for j in (1,10) + target[:,i,j] = target[i,:,j] = target[i,j,:] = (2/3)^2 + end + end + for i in (1,10) + for j in (1,10) + for k in (1,10) + target[i,j,k] = (2/3)^3 + end + end + end + @test_approx_eq imfilter(img, kernel, Fill(0)) target + @test_approx_eq imfilter(img, kernel, Fill(0), Algorithm.FFT()) target + # Inner + @test_approx_eq imfilter(img, kernel, Inner()) OffsetArray(trues(8,8,8), 2:9, 2:9, 2:9) +end + +nothing diff --git a/test/runtests.jl b/test/runtests.jl index ceea667..2a9d902 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ asa = detect_ambiguities(StaticArrays, Base) @test isempty(setdiff(aif, asa)) include("border.jl") +include("nd.jl") include("2d.jl") include("triggs.jl") include("cascade.jl") diff --git a/test/triggs.jl b/test/triggs.jl index a840b7d..6e3b0b4 100644 --- a/test/triggs.jl +++ b/test/triggs.jl @@ -2,7 +2,7 @@ using ImagesFiltering, Colors, ComputationalResources using Base.Test @testset "TriggsSdika" begin - @testset "TriggsSdika 1d" begin + @testset "1d" begin l = 1000 x = [1, 2, 3, 5, 10, 20, l>>1, l-19, l-9, l-4, l-2, l-1, l] σs = [3.1, 5, 10.0, 20, 50.0, 100.0] @@ -28,7 +28,7 @@ using Base.Test end end - @testset "TriggsSdika images" begin + @testset "images" begin imgf = zeros(5, 7); imgf[3,4] = 1 imgg = fill(Gray{Float32}(0), 5, 7); imgg[3,4] = 1 imgc = fill(RGB{Float64}(0,0,0), 5, 7); imgc[3,4] = RGB(1,0,0) From eefdadb0cfe9a48c44b6e803db67554152e6558d Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 24 Aug 2016 13:50:37 -0500 Subject: [PATCH 11/21] Ensure commutivity when using a mixture of FIR and TriggsSdika kernels From an efficiency standpoint it's better to do the TriggsSdika kernels first. --- src/ImagesFiltering.jl | 17 ++++++++++++++++- src/border.jl | 35 ++++++++++++++++++++++++++++++----- src/imfilter.jl | 31 ++++++++++++++++++------------- test/triggs.jl | 9 +++++++++ 4 files changed, 73 insertions(+), 19 deletions(-) diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index 2b6bac0..43c0e1d 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -23,14 +23,29 @@ Alg{A<:Alg}(r::AbstractResource{A}) = r.settings include("utils.jl") include("kernelfactors.jl") using .KernelFactors: TriggsSdika, IIRFilter, ReshapedVector, iterdims + +typealias ArrayLike{T} Union{AbstractArray{T}, IIRFilter{T}, ReshapedVector{T}} +typealias ReshapedTriggsSdika{T,N,Npre,V<:TriggsSdika} ReshapedVector{T,N,Npre,V} +typealias AnyTriggs Union{TriggsSdika, ReshapedTriggsSdika} + include("kernel.jl") using .Kernel using .Kernel: Laplacian -typealias ArrayLike{T} Union{AbstractArray{T}, IIRFilter{T}, ReshapedVector{T}} +typealias NDimKernel{N,K} Union{AbstractArray{K,N},ReshapedVector{K,N},Laplacian{N}} include("border.jl") + +typealias BorderSpec{Style,T} Union{Pad{Style,0}, Fill{T,0}, Inner{0}} +typealias BorderSpecRF{T} Union{Pad{:replicate,0}, Fill{T,0}} +typealias BorderSpecNoNa{T} Union{Pad{:replicate,0}, Pad{:circular,0}, Pad{:symmetric,0}, + Pad{:reflect,0}, Fill{T,0}, Inner{0}} +typealias BorderSpecAny Union{BorderSpec,NoPad} + include("deprecated.jl") + +typealias ProcessedKernel Tuple + include("imfilter.jl") include("specialty.jl") diff --git a/src/border.jl b/src/border.jl index 6468948..44f3724 100644 --- a/src/border.jl +++ b/src/border.jl @@ -80,13 +80,13 @@ Given a filter array `kernel`, determine the amount of padding from the `indices """ (p::Pad{Style,0}){Style}(kernel::AbstractArray) = Pad{Style}(indices(kernel)) (p::Pad{Style,0}){Style}(kernel::Laplacian) = Pad{Style}(indices(kernel)) -(p::Pad{Style,0}){Style}(factkernel::Tuple) = Pad{Style}(accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) +(p::Pad{Style,0}){Style}(factkernel::Tuple) = Pad{Style}(calculate_padding(factkernel)) (p::Pad{Style,0}){Style}(factkernel::Tuple, img, ::FIR) = p(factkernel) (p::Pad{Style,0}){Style}(kernel::Laplacian, img, ::FIR) = p(kernel) # Padding for FFT: round up to next size expressible as 2^m*3^n function (p::Pad{Style,0}){Style}(factkernel::Tuple, img, ::FFT) - inds = accumulate_padding(indices(factkernel[1]), tail(factkernel)...) + inds = calculate_padding(factkernel) newinds = map(padfft, inds, map(length, indices(img))) Pad{Style}(newinds) end @@ -160,7 +160,7 @@ end (p::Inner{0})(factkernel::Tuple, img, ::FIR) = p(factkernel) (p::Inner{0})(factkernel::Tuple, img, ::FFT) = p(factkernel) -(p::Inner{0})(factkernel::Tuple) = Inner(accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) +(p::Inner{0})(factkernel::Tuple) = Inner(calculate_padding(factkernel)) (p::Inner{0})(kernel::AbstractArray) = Inner(indices(kernel)) padarray(img, border::Inner) = padarray(eltype(img), img, border) @@ -190,12 +190,12 @@ Fill(value, lo::AbstractVector, hi::AbstractVector) = Fill(value, (lo...,), (hi. Fill{T,N}(value::T, inds::Base.Indices{N}) = Fill{T,N}(value, map(lo,inds), map(hi,inds)) Fill(value, kernel::AbstractArray) = Fill(value, indices(kernel)) Fill(value, kernel::Laplacian) = Fill(value, indices(kernel)) -Fill(value, factkernel::Tuple) = Fill(value, accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) +Fill(value, factkernel::Tuple) = Fill(value, calculate_padding(factkernel)) (p::Fill)(kernel) = Fill(p.value, kernel) (p::Fill)(kernel, img, ::FIR) = Fill(p.value, kernel) function (p::Fill)(factkernel::Tuple, img, ::FFT) - inds = accumulate_padding(indices(factkernel[1]), tail(factkernel)...) + inds = calculate_padding(factkernel) newinds = map(padfft, inds, map(length, indices(img))) Fill(p.value, newinds) end @@ -260,6 +260,31 @@ end # @inline _flatten(t1, t...) = (t1, flatten(t)...) # _flatten() = () +@inline function calculate_padding(kernel::Tuple{Any, Vararg{Any}}) + inds = accumulate_padding(indices(kernel[1]), tail(kernel)...) + if hastriggs(kernel) && hasfir(kernel) + return map(doublepadding, inds) + end + inds +end + +hastriggs(kernel) = _hastriggs(false, kernel...) +_hastriggs(ret) = ret +_hastriggs(ret, kern, kernel...) = _hastriggs(ret, kernel...) +_hastriggs(ret, kern::AnyTriggs, kernel...) = true + +hasfir(kernel) = _hasfir(false, kernel...) +_hasfir(ret) = ret +_hasfir(ret, kern, kernel...) = true +_hasfir(ret, kern::AnyTriggs, kernel...) = _hasfir(ret, kernel...) + +function doublepadding(ind::AbstractUnitRange) + f, l = first(ind), last(ind) + f = f < 0 ? 2f : f + l = l > 0 ? 2l : l + f:l +end + accumulate_padding(inds::Indices, kernel1, kernels...) = accumulate_padding(_accumulate_padding(inds, indices(kernel1)), kernels...) accumulate_padding(inds::Indices, kernel1::TriggsSdika, kernels...) = diff --git a/src/imfilter.jl b/src/imfilter.jl index 3ca2b4e..228d107 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -1,12 +1,3 @@ -typealias BorderSpec{Style,T} Union{Pad{Style,0}, Fill{T,0}, Inner{0}} -typealias BorderSpecNoNa{T} Union{Pad{:replicate,0}, Pad{:circular,0}, Pad{:symmetric,0}, - Pad{:reflect,0}, Fill{T,0}, Inner{0}} -typealias BorderSpecAny Union{BorderSpec,NoPad} -typealias NDimKernel{N,K} Union{AbstractArray{K,N},ReshapedVector{K,N},Laplacian{N}} -typealias ReshapedTriggsSdika{T,N,Npre,V<:TriggsSdika} ReshapedVector{T,N,Npre,V} - -typealias ProcessedKernel Tuple - # see below for imfilter docstring # Step 1: if necessary, determine the output's element type @@ -244,12 +235,28 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke return out end +# We deliberately drop `inds` in the next methods: when `kernel` is a tuple +# that has both TriggsSdika and FIR filters, the overall padding +# gets doubled, yet we only trim off the minimum at each +# stage. Consequently, `inds` might be optimistic about the range +# available in `out`. function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, border::NoPad, inds::Indices) - imfilter!(r, out, A1, kernel, border, inds) + imfilter!(r, out, A1, kernel, border) end function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, border::NoPad, inds::Indices) - imfilter!(r, out, A1, kernel[1], border, inds) + imfilter!(r, out, A1, kernel[1], border) +end + +# However, here it's important to filter over the whole passed-in +# range, and then copy! to out +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{AnyTriggs}, border::NoPad, inds::Indices) + if inds != indices(out) + imfilter!(r, A2, A1, kernel[1], border, inds) + R = CartesianRange(indices(out)) + return copy!(out, R, A2, R) + end + imfilter!(r, out, A1, kernel[1], border) end function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, inds::Indices) @@ -433,8 +440,6 @@ end # Note this is safe for inplace use, i.e., out === img -typealias BorderSpecRF{T} Union{Pad{:replicate,0}, Fill{T,0}} - function imfilter!{S,T,N}(r::AbstractResource, out::AbstractArray{S,N}, img::AbstractArray{T,N}, diff --git a/test/triggs.jl b/test/triggs.jl index 6e3b0b4..52111c7 100644 --- a/test/triggs.jl +++ b/test/triggs.jl @@ -28,6 +28,15 @@ using Base.Test end end + @testset "commutivity" begin + img = 1:8 + k1 = KernelFactors.IIRGaussian(2) + k2 = centered(ones(3)/3) + for border in (Pad{:replicate}(), Fill(0)) + @test_approx_eq imfilter(img, (k1, k2), border) imfilter(img, (k2, k1), border) + end + end + @testset "images" begin imgf = zeros(5, 7); imgf[3,4] = 1 imgg = fill(Gray{Float32}(0), 5, 7); imgg[3,4] = 1 From a635994a95413a2aaf72574a861e773d5ead49ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Wed, 24 Aug 2016 19:55:44 -0500 Subject: [PATCH 12/21] Add imgradients Taken from Images pull request 544 --- src/specialty.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/src/specialty.jl b/src/specialty.jl index c0a114b..6f69e47 100644 --- a/src/specialty.jl +++ b/src/specialty.jl @@ -14,3 +14,78 @@ function _imfilter_inbounds!(r, out, A::AbstractArray, L::Laplacian, border::NoP end out end + +## imgradients +""" + imgradients(img, [points], [method], [border]) -> g1, g2, ... + +Performs edge detection filtering of the N-dimensional array `img`. +Gradients are computed at specified `points` (or indexes) in the +array or everywhere. Available methods: `"sobel"` and `"ando3"`. +Border options:`"replicate"`, `"circular"`, `"reflect"`, `"symmetric"`. +Returns a 2D array `G` with the gradients as rows. The number of rows +is the number of points at which the gradient was computed and the +number of columns is the dimensionality of the array. +""" +function imgradients{T,N}(img::AbstractArray{T,N}, points::AbstractVector; + method::AbstractString="ando3", border::AbstractString="replicate") + extent = size(img) + ndirs = length(extent) + npoints = length(points) + + # pad input image only on appropriate directions + imgpad = _gradientpad(img, border) + + # gradient matrix + Tret = typeof(zero(T)*zero(Float64)) + G = zeros(Tret, npoints, ndirs) + + for dir in 1:ndirs + # kernel = centered difference + perpendicular smoothing + if extent[dir] > 1 + kern = _directional_kernel(dir, extent, method) + + # compute gradient at specified points + A = zeros(kern) + shape = size(kern) + for (k, p) in enumerate(points) + icenter = CartesianIndex(ind2sub(extent, p)) + i1 = CartesianIndex(tuple(ones(Int, ndirs)...)) + for ii in CartesianRange(shape) + A[ii] = imgpad[ii + icenter - i1] + end + + G[k,dir] = sum(kern .* A) + end + end + end + + G +end + +function imgradients{T,N}(img::AbstractArray{T,N}; method::AbstractString="ando3", border::AbstractString="replicate") + extent = size(img) + ndirs = length(extent) + + # pad input image only on appropriate directions + imgpad = _gradientpad(img, border) + + # gradient tuple + Tret = typeof(zero(T)*zero(Float64)) + G = Array(Array{Tret,N}, ndirs) + + for dir in 1:ndirs + # kernel = centered difference + perpendicular smoothing + if extent[dir] > 1 + kern = _directional_kernel(dir, extent, method) + G[dir] = imfilter(imgpad, kern, "inner") + end + end + + result = (G...,) + if ndims(img) == 2 && spatialorder(img) == yx + result = (result[2], result[1]) + end + + result +end From bc236d0dd8274ea312e34d5c18f8dfd9f039f466 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Thu, 25 Aug 2016 11:02:17 -0500 Subject: [PATCH 13/21] Modernize the imfiltering version of imgradients --- src/ImagesFiltering.jl | 9 +++++- src/imfilter.jl | 2 +- src/kernel.jl | 26 ++++++++++++++-- src/kernelfactors.jl | 71 +++++++++++++++++++++++++++++++++++------- src/specialty.jl | 38 +++++++++------------- test/gradient.jl | 46 +++++++++++++++++++++++++++ 6 files changed, 154 insertions(+), 38 deletions(-) create mode 100644 test/gradient.jl diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index 43c0e1d..39c71a7 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -4,9 +4,16 @@ using Colors, FixedPointNumbers, ImagesCore, MappedArrays, FFTViews, OffsetArray using ColorVectorSpace # in case someone filters RGB arrays using Base: Indices, tail, fill_to_length, @pure -export Kernel, KernelFactors, Pad, Fill, Inner, NoPad, Algorithm, imfilter, imfilter!, padarray, centered +export Kernel, KernelFactors, Pad, Fill, Inner, NoPad, Algorithm, imfilter, imfilter!, imgradients, padarray, centered typealias FixedColorant{T<:UFixed} Colorant{T} +typealias StaticOffsetArray{T,N,A<:StaticArray} OffsetArray{T,N,A} + +# Needed for type-stability +function Base.transpose{T}(A::StaticOffsetArray{T,2}) + inds1, inds2 = indices(A) + OffsetArray(transpose(parent(A)), inds2, inds1) +end module Algorithm # deliberately don't export these, but it's expected that they diff --git a/src/imfilter.jl b/src/imfilter.jl index 228d107..ad4b35d 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -689,7 +689,7 @@ function factorstridedkernel(inds, kernel::StridedMatrix) end prod_kernel(kern::AbstractArray) = kern -prod_kernel(kern::AbstractArray, kern1, kerns...) = broadcast(.*, kern, kern1, kerns...) #prod_kernel(kern.*kern1, kerns...) +prod_kernel(kern::AbstractArray, kern1, kerns...) = broadcast(*, kern, kern1, kerns...) #prod_kernel(kern.*kern1, kerns...) prod_kernel{N}(::Type{Val{N}}, args...) = prod_kernel(Val{N}, prod_kernel(args...)) prod_kernel{_,N}(::Type{Val{N}}, kernel::AbstractArray{_,N}) = kernel function prod_kernel{N}(::Type{Val{N}}, kernel::AbstractArray) diff --git a/src/kernel.jl b/src/kernel.jl index 9a77c71..ebdb414 100644 --- a/src/kernel.jl +++ b/src/kernel.jl @@ -20,6 +20,8 @@ See also: KernelFactors.sobel, Kernel.prewitt, Kernel.ando. """ sobel() = product2d(KernelFactors.sobel()) +sobel{N}(::Type{Val{N}},d,extent=trues(N)) = broadcast(*, KernelFactors.sobel(Val{N},d,extent)...) + """ diff1, diff2 = prewitt() @@ -32,6 +34,8 @@ See also: KernelFactors.prewitt, Kernel.sobel, Kernel.ando. """ prewitt() = product2d(KernelFactors.prewitt()) +prewitt{N}(::Type{Val{N}},d,extent=trues(N)) = broadcast(*, KernelFactors.prewitt(Val{N},d,extent)...) + """ diff1, diff2 = ando3() @@ -48,6 +52,8 @@ See also: KernelFactors.ando3, Kernel.ando4, Kernel.ando5. """ ando3() = product2d(KernelFactors.ando3()) +ando3{N}(::Type{Val{N}},d,extent=trues(N)) = broadcast(*, KernelFactors.ando3(Val{N},d,extent)...) + """ diff1, diff2 = ando4() @@ -70,6 +76,14 @@ function ando4() return f', f end +function ando4{N}(::Type{Val{N}}, d, extent=trues(N)) + if N == 2 && all(extent) + return ando4()[d] + else + error("dimensions other than 2 are not yet supported") + end +end + """ diff1, diff2 = ando5() @@ -93,6 +107,14 @@ function ando5() return f', f end +function ando5{N}(::Type{Val{N}}, d, extent) + if N == 2 && all(extent) + return ando5()[d] + else + error("dimensions other than 2 are not yet supported") + end +end + """ gaussian((σ1, σ2, ...), [(l1, l2, ...]) -> g gaussian(σ) -> g @@ -107,13 +129,13 @@ constructed. See also: KernelFactors.gaussian. """ @inline gaussian{N}(σs::NTuple{N,Real}, ls::NTuple{N,Integer}) = - broadcast(.*, KernelFactors.gaussian(σs, ls)...) + broadcast(*, KernelFactors.gaussian(σs, ls)...) gaussian(σ::Tuple{Real}, l::Tuple{Integer}) = KernelFactors.gaussian(σ[1], l[1]) gaussian(σ::Tuple{}, l::Tuple{}) = reshape([1]) # 0d gaussian{T<:Real,I<:Integer}(σs::AbstractVector{T}, ls::AbstractVector{I}) = gaussian((σs...,), (ls...,)) -@inline gaussian{N}(σs::NTuple{N,Real}) = broadcast(.*, KernelFactors.gaussian(σs)...) +@inline gaussian{N}(σs::NTuple{N,Real}) = broadcast(*, KernelFactors.gaussian(σs)...) gaussian{T<:Real}(σs::AbstractVector{T}) = gaussian((σs...,)) gaussian(σ::Tuple{Real}) = KernelFactors.gaussian(σ[1]) gaussian(σ::Tuple{}) = reshape([1]) diff --git a/src/kernelfactors.jl b/src/kernelfactors.jl index 8bac143..17762d7 100644 --- a/src/kernelfactors.jl +++ b/src/kernelfactors.jl @@ -81,20 +81,43 @@ end ## gradients -"`kern1, kern2 = sobel()` returns factored Sobel filters for dimensions 1 and 2 of your image" +""" + kern1, kern2 = sobel() + +Factored Sobel filters for dimensions 1 and 2 of a two-dimensional +image. Each is a 2-tuple of one-dimensional filters. +""" function sobel() - f1 = centered(SVector( 1.0, 2.0, 1.0)) - f2 = centered(SVector(-1.0, 0.0, 1.0)) + f1 = centered(SVector( 1.0, 2.0, 1.0)/4) + f2 = centered(SVector(-1.0, 0.0, 1.0)/2) return kernelfactors((f2, f1)), kernelfactors((f1, f2)) end +""" + kern = sobel(N, d) + +Return a factored Sobel filter for computing the gradient in N dimensions along axis d. +""" +function sobel{N}(::Type{Val{N}}, d, extended=trues(N)) + gradfactors(Val{N}, d, [-1, 0, 1]/2, [1, 2, 1]/4, extended) +end + "`kern1, kern2 = prewitt()` returns factored Prewitt filters for dimensions 1 and 2 of your image" function prewitt() - f1 = centered(SVector( 1.0, 1.0, 1.0)) - f2 = centered(SVector(-1.0, 0.0, 1.0)) + f1 = centered(SVector( 1.0, 1.0, 1.0)/3) + f2 = centered(SVector(-1.0, 0.0, 1.0)/2) return kernelfactors((f2, f1)), kernelfactors((f1, f2)) end +""" + kern = prewitt(N, d) + +Return a factored Prewitt filter for computing the gradient in N dimensions along axis d. +""" +function prewitt{N}(::Type{Val{N}}, d, extended=trues(N)) + gradfactors(Val{N}, d, [-1, 0, 1]/2, [1, 1, 1]/3, extended) +end + # Consistent Gradient Operators # Ando Shigeru # IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March 2000 @@ -111,11 +134,15 @@ Ando Shigeru, IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March 2000. See also: `ando4`, `ando5`. """ function ando3() - f1 = centered(SVector(0.112737, 0.274526, 0.112737)) - f2 = centered(SVector(-1.0, 0.0, 1.0)) + f1 = centered(2*SVector(0.112737, 0.274526, 0.112737)) + f2 = centered(SVector(-1.0, 0.0, 1.0)/2) return kernelfactors((f2, f1)), kernelfactors((f1, f2)) end +function ando3{N}(::Type{Val{N}},d,extended=trues(N)) + gradfactors(Val{N}, d, [-1.0, 0.0, 1.0]/2, 2*[0.112737, 0.274526, 0.112737], extended) +end + """ `kern1, kern2 = ando4()` returns separable approximations of the optimal 4x4 filters for dimensions 1 and 2 of your image, as defined @@ -125,11 +152,19 @@ March 2000. See also: `Kernel.ando4`. """ function ando4() - f1 = centered(SVector(0.025473821998749126, 0.11299599504060115, 0.11299599504060115, 0.025473821998749126)) - f2 = centered(SVector(-0.2254400431590164, -1.0, 1.0, 0.2254400431590164)) + f1 = centered(SVector( 0.0919833, 0.408017, 0.408017, 0.0919833)) + f2 = centered(1.46205884*SVector(-0.0919833, -0.408017, 0.408017, 0.0919833)) return kernelfactors((f2, f1)), kernelfactors((f1, f2)) end +function ando4{N}(::Type{Val{N}}, d, extended=trues(N)) + if N == 2 && all(extended) + return ando4()[d] + else + error("dimensions other than 2 are not yet supported") + end +end + """ `kern1, kern2 = ando5_sep()` returns separable approximations of the optimal 5x5 gradient filters for dimensions 1 and 2 of your image, as defined @@ -139,11 +174,19 @@ March 2000. See also: `Kernel.ando5`. """ function ando5() - f1 = centered(SVector(0.03136697678461958, 0.21844976784066258, 0.37816313370270255, 0.21844976784066258, 0.03136697678461958)) - f2 = centered(SVector(-0.12288050911244743, -0.3242040200682835, 0.0, 0.3242040200682835, 0.12288050911244743)) + f1 = centered(SVector( 0.0357338, 0.248861, 0.43081, 0.248861, 0.0357338)) + f2 = centered(0.784406*SVector(-0.137424, -0.362576, 0.0, 0.362576, 0.137424)) return kernelfactors((f2, f1)), kernelfactors((f1, f2)) end +function ando5{N}(::Type{Val{N}}, d, extended=trues(N)) + if N == 2 && all(extended) + return ando5()[d] + else + error("dimensions other than 2 are not yet supported") + end +end + ## Gaussian """ @@ -313,4 +356,10 @@ end # A variant in which we just need to fill out to N dimensions kernelfactors{N}(factors::NTuple{N,AbstractArray}) = map(f->_reshape(f, Val{N}), factors) +function gradfactors{N}(::Type{Val{N}}, d::Int, k1, k2, extended=trues(N)) + kernelfactors(ntuple(i -> kdim(extended[i], i==d ? k1 : k2), Val{N})) +end + +kdim(keep::Bool, k) = keep ? centered(k) : OffsetArray([one(eltype(k))], 0:0) + end diff --git a/src/specialty.jl b/src/specialty.jl index 6f69e47..25e417a 100644 --- a/src/specialty.jl +++ b/src/specialty.jl @@ -63,29 +63,21 @@ function imgradients{T,N}(img::AbstractArray{T,N}, points::AbstractVector; G end -function imgradients{T,N}(img::AbstractArray{T,N}; method::AbstractString="ando3", border::AbstractString="replicate") - extent = size(img) - ndirs = length(extent) - - # pad input image only on appropriate directions - imgpad = _gradientpad(img, border) - - # gradient tuple - Tret = typeof(zero(T)*zero(Float64)) - G = Array(Array{Tret,N}, ndirs) - - for dir in 1:ndirs - # kernel = centered difference + perpendicular smoothing - if extent[dir] > 1 - kern = _directional_kernel(dir, extent, method) - G[dir] = imfilter(imgpad, kern, "inner") - end - end +function imgradients{T,N}(img::AbstractArray{T,N}, kernelfun::Function, border=Pad{:replicate}()) + extent = map(length, indices(img)) + kern1 = kernelfun(Val{N}, 1, map(x->x>1, extent)) + S = filter_type(img, kern1) + bord = border(kern1) + imgpad = padarray(S, img, bord) + _imgradients((), imgpad, kernelfun, extent, Inner()) +end - result = (G...,) - if ndims(img) == 2 && spatialorder(img) == yx - result = (result[2], result[1]) - end +# When all N gradients have been calculated, return the result +_imgradients{T,N}(G::NTuple{N}, imgpad::AbstractArray{T,N}, kernelfun::Function, extent, border::Inner) = G - result +# Add the next dimension to G +function _imgradients{T,M,N}(G::NTuple{M}, imgpad::AbstractArray{T,N}, kernelfun::Function, extent, border::Inner) + d = M+1 # the dimension we're working on now + kern = kernelfun(Val{N}, d, map(x->x>1, extent)) + _imgradients((G..., imfilter(imgpad, kern, border)), imgpad, kernelfun, extent, border) end diff --git a/test/gradient.jl b/test/gradient.jl new file mode 100644 index 0000000..367676c --- /dev/null +++ b/test/gradient.jl @@ -0,0 +1,46 @@ +using ImagesFiltering, Colors, ColorVectorSpace +using Base.Test + +@testset "gradient" begin + y, x = 1:5, (1:7)' + for (img, ey, ex) in ((y .* ones(x), 1, 0), + (ones(y) .* x, 0, 1), + (map(Gray, y.*ones(Float64, x)), Gray(1), Gray(0)), + (map(v->RGB(v,0,0), y.*ones(Float64, x)), RGB(1,0,0), RGB(0,0,0))) + for kernelfunc in (KernelFactors.ando3, KernelFactors.sobel, KernelFactors.prewitt, + KernelFactors.ando4, KernelFactors.ando5, + Kernel.ando3, Kernel.sobel, Kernel.prewitt, + Kernel.ando4, Kernel.ando5) + gy, gx = @inferred(imgradients(img, kernelfunc, Inner())) + for val in gy + @test abs(val - ey) < 1e-4 + end + for val in gx + @test abs(val - ex) < 1e-4 + end + gy, gx = imgradients(img, kernelfunc, Pad{:replicate}()) + @test indices(gy) == indices(gx) == indices(img) + end + end + # 3d + y, x, z = 1:5, (1:7)', reshape(1:6, 1, 1, 6) + for (img, ey, ex, ez) in ((y .* ones(x) .* ones(z), 1, 0, 0), + (ones(y) .* x .* ones(z), 0, 1, 0), + (ones(y) .* ones(x) .* z, 0, 0, 1)) + for kernelfunc in (KernelFactors.ando3, KernelFactors.sobel, KernelFactors.prewitt, + Kernel.ando3, Kernel.sobel, Kernel.prewitt) + gy, gx, gz = @inferred(imgradients(img, kernelfunc, Inner())) + for val in gy + @test abs(val - ey) < 1e-4 + end + for val in gx + @test abs(val - ex) < 1e-4 + end + for val in gz + @test abs(val - ez) < 1e-4 + end + end + end +end + +nothing From ebf19edc441ac4d2ed44c591612c2b1f60c51cb4 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 26 Aug 2016 04:48:25 -0500 Subject: [PATCH 14/21] Give ReshapedVector more AbstractArray functionality --- src/kernelfactors.jl | 36 +++++++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/src/kernelfactors.jl b/src/kernelfactors.jl index 17762d7..6b9f99d 100644 --- a/src/kernelfactors.jl +++ b/src/kernelfactors.jl @@ -3,7 +3,7 @@ module KernelFactors using StaticArrays, OffsetArrays using ..ImagesFiltering: centered, dummyind import ..ImagesFiltering: _reshape, _vec, nextendeddims -using Base: tail, Indices +using Base: tail, Indices, @pure, checkbounds_indices, throw_boundserror abstract IIRFilter{T} @@ -23,7 +23,7 @@ AbstractArray. ReshapedVectors allow one to specify a "filtering dimension" for a 1-dimensional filter. """ -immutable ReshapedVector{T,N,Npre,V} # note not <: AbstractArray{T,N} (more general) +immutable ReshapedVector{T,N,Npre,V} # not <: AbstractArray{T,N} (more general, incl. IIR) data::V function ReshapedVector(data::V) @@ -41,14 +41,44 @@ end end _reshapedvector{N,Npre}(::NTuple{N,Bool}, ::NTuple{Npre,Bool}, data) = ReshapedVector{eltype(data),N,Npre,typeof(data)}(data) +# Give ReshapedVector many of the characteristics of AbstractArray Base.eltype{T}(A::ReshapedVector{T}) = T Base.ndims{_,N}(A::ReshapedVector{_,N}) = N Base.isempty(A::ReshapedVector) = isempty(A.data) -@inline Base.indices{_,N,Npre}(A::ReshapedVector{_,N,Npre}) = Base.fill_to_length((Base.ntuple(d->0:0, Val{Npre})..., Base.indices1(A.data)), 0:0, Val{N}) +@inline Base.indices{_,N,Npre}(A::ReshapedVector{_,N,Npre}) = Base.fill_to_length((Base.ntuple(d->0:0, Val{Npre})..., UnitRange(Base.indices1(A.data))), 0:0, Val{N}) + +Base.start(A::ReshapedVector) = start(A.data) +Base.next(A::ReshapedVector, state) = next(A.data, state) +Base.done(A::ReshapedVector, state) = done(A.data, state) + +@inline function Base.getindex(A::ReshapedVector, i::Int) + @boundscheck checkbounds(A.data, i) + @inbounds ret = A.data[i] + ret +end +@inline function Base.getindex{T,N,Npre}(A::ReshapedVector{T,N,Npre}, I::Vararg{Int,N}) + @boundscheck checkbounds_indices(Bool, indices(A), I) || throw_boundserror(A, I) + @inbounds ret = A.data[I[Npre+1]] + ret +end +@inline function Base.getindex(A::ReshapedVector, I::Union{CartesianIndex,Int}...) + A[Base.IteratorsMD.flatten(I)...] +end + +Base.convert{AA<:AbstractArray}(::Type{AA}, A::ReshapedVector) = convert(AA, reshape(A.data, indices(A))) + +import Base: == +==(A::ReshapedVector, B::ReshapedVector) = convert(AbstractArray, A) == convert(AbstractArray, B) +==(A::ReshapedVector, B::AbstractArray) = convert(AbstractArray, A) == B +==(A::AbstractArray, B::ReshapedVector) = A == convert(AbstractArray, B) + +# for broadcasting +@pure Base.promote_eltype_op{S,T}(op, ::ReshapedVector{S}, ::ReshapedVector{T}) = Base.promote_op(op, S, T) _reshape{T,N}(A::ReshapedVector{T,N}, ::Type{Val{N}}) = A _vec(A::ReshapedVector) = A.data +Base.vec(A::ReshapedVector) = A.data # is this OK? (note indices won't nec. start with 1) nextendeddims(a::ReshapedVector) = 1 """ From 0a97c835dd73f06daa71542d24cff43a94386f61 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 26 Aug 2016 04:50:22 -0500 Subject: [PATCH 15/21] Substantial performance improvements for imgradients Switches to returning ReshapedVectors from the factored kernels, and adds optimized methods. It's not clear we strictly need all the optimized methods, so some of these may go away sometime. --- src/border.jl | 7 +++++ src/imfilter.jl | 63 +++++++++++++++++++++++++++++++++++++------- src/kernel.jl | 10 +++---- src/kernelfactors.jl | 13 +++++---- src/specialty.jl | 9 ++++--- test/basic.jl | 4 ++- test/specialty.jl | 4 +-- 7 files changed, 81 insertions(+), 29 deletions(-) diff --git a/src/border.jl b/src/border.jl index 44f3724..55b9723 100644 --- a/src/border.jl +++ b/src/border.jl @@ -318,4 +318,11 @@ function next_interior(inds::Indices, kernel::Tuple) _interior(inds, indices(kern)) end +allocate_output{T}(::Type{T}, img, kernel, border) = similar(img, T) +function allocate_output{T}(::Type{T}, img, kernel, ::Inner{0}) + inds = interior(img, kernel) + similar(img, T, inds) +end +allocate_output(img, kernel, border) = allocate_output(filter_type(img, kernel), img, kernel, border) + # end diff --git a/src/imfilter.jl b/src/imfilter.jl index ad4b35d..0d6ca59 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -16,12 +16,8 @@ function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, arg end # Step 4: if necessary, allocate the ouput -@inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::Inner{0}, args...) - inds = interior(img, kernel) - imfilter!(similar(img, T, inds), img, kernel, border, args...) -end @inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny, args...) - imfilter!(similar(img, T), img, kernel, border, args...) + imfilter!(allocate_output(T, img, kernel, border), img, kernel, border, args...) end # Now do the same steps for the case where the user supplies a Resource @@ -38,12 +34,9 @@ end function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ProcessedKernel) imfilter(r, T, img, kernel, Pad{:replicate}()) # supply the default border end -function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::Inner{0}) - inds = interior(img, kernel) - imfilter!(r, similar(img, T, inds), img, kernel, border) -end + function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny) - imfilter!(r, similar(img, T), img, kernel, border) + imfilter!(r, allocate_output(T, img, kernel, border), img, kernel, border) end """ @@ -347,6 +340,56 @@ function _imfilter_inbounds!(r, out, A::AbstractArray, kern, border::NoPad, inds out end +function _imfilter_inbounds!(r, out, A::AbstractArray, kern::ReshapedVector, border::NoPad, inds) + Rpre, ind, Rpost = iterdims(inds, kern) + k = kern.data + R, Rk = CartesianRange(inds), CartesianRange(indices(kern)) + p = A[first(R)+first(Rk)] * first(k) + TT = typeof(p+p) + _imfilter_inbounds(r, TT, out, A, k, Rpre, ind, Rpost) +end + +function _imfilter_inbounds{TT}(r, ::Type{TT}, out, A::AbstractArray, k::AbstractVector, Rpre::CartesianRange, ind, Rpost::CartesianRange) + indsk = indices(k, 1) + if length(indsk) <= 5 + ks, jo = to_static(k) + return @time _imfilter_inbounds(r, TT, out, A, (ks, jo), Rpre, ind, Rpost) + end + for Ipost in Rpost + for Ipre in Rpre + for i in ind + tmp = zero(TT) + @unsafe for j in indsk + tmp += A[Ipre,i+j,Ipost]*k[j] + end + @unsafe out[Ipre,i,Ipost] = tmp + end + end + end + out +end + +to_static(a::OffsetArray) = _to_static(parent(a), a.offsets[1]) +to_static(a) = _to_static(a, first(indices(a,1))-1) +_to_static(a, o) = convert(SVector{length(a)}, a), o + +function _imfilter_inbounds{TT,L}(r, ::Type{TT}, out, A::AbstractArray, ko::Tuple{SVector{L},Int}, Rpre::CartesianRange, ind, Rpost::CartesianRange) + k, jo = ko + for Ipost in Rpost + for Ipre in Rpre + for i in ind + io = i+jo + tmp = zero(TT) + @unsafe for j = 1:L + tmp += A[Ipre,io+j,Ipost]*k[j] + end + @unsafe out[Ipre,i,Ipost] = tmp + end + end + end + out +end + ### FFT filtering """ diff --git a/src/kernel.jl b/src/kernel.jl index ebdb414..a2fe55f 100644 --- a/src/kernel.jl +++ b/src/kernel.jl @@ -20,7 +20,7 @@ See also: KernelFactors.sobel, Kernel.prewitt, Kernel.ando. """ sobel() = product2d(KernelFactors.sobel()) -sobel{N}(::Type{Val{N}},d,extent=trues(N)) = broadcast(*, KernelFactors.sobel(Val{N},d,extent)...) +sobel{N}(::Type{Val{N}},d,extent=trues(N)) = (broadcast(*, KernelFactors.sobel(Val{N},d,extent)...),) """ diff1, diff2 = prewitt() @@ -34,7 +34,7 @@ See also: KernelFactors.prewitt, Kernel.sobel, Kernel.ando. """ prewitt() = product2d(KernelFactors.prewitt()) -prewitt{N}(::Type{Val{N}},d,extent=trues(N)) = broadcast(*, KernelFactors.prewitt(Val{N},d,extent)...) +prewitt{N}(::Type{Val{N}},d,extent=trues(N)) = (broadcast(*, KernelFactors.prewitt(Val{N},d,extent)...),) """ diff1, diff2 = ando3() @@ -52,7 +52,7 @@ See also: KernelFactors.ando3, Kernel.ando4, Kernel.ando5. """ ando3() = product2d(KernelFactors.ando3()) -ando3{N}(::Type{Val{N}},d,extent=trues(N)) = broadcast(*, KernelFactors.ando3(Val{N},d,extent)...) +ando3{N}(::Type{Val{N}},d,extent=trues(N)) = (broadcast(*, KernelFactors.ando3(Val{N},d,extent)...),) """ diff1, diff2 = ando4() @@ -78,7 +78,7 @@ end function ando4{N}(::Type{Val{N}}, d, extent=trues(N)) if N == 2 && all(extent) - return ando4()[d] + return (ando4()[d],) else error("dimensions other than 2 are not yet supported") end @@ -109,7 +109,7 @@ end function ando5{N}(::Type{Val{N}}, d, extent) if N == 2 && all(extent) - return ando5()[d] + return (ando5()[d],) else error("dimensions other than 2 are not yet supported") end diff --git a/src/kernelfactors.jl b/src/kernelfactors.jl index 6b9f99d..add1b37 100644 --- a/src/kernelfactors.jl +++ b/src/kernelfactors.jl @@ -371,16 +371,15 @@ If passed a tuple of general arrays, it is assumed that each is shaped appropriately along its "leading" dimensions; the dimensionality of each is "extended" to `N = length(factors)`, appending 1s to the size as needed. """ -kernelfactors{N}(factors::NTuple{N,AbstractVector}) = _kernelfactors((), factors) +kernelfactors{N}(factors::NTuple{N,AbstractVector}) = _kernelfactors((), ntuple(d->true,Val{N}), (), factors) -_kernelfactors(out::NTuple, ::Tuple{}) = out -@inline function _kernelfactors{L,M}(out::NTuple{L}, factors::NTuple{M,AbstractVector}) +_kernelfactors(pre, post, out::NTuple, ::Tuple{}) = out +@inline function _kernelfactors(pre, post, out, factors) # L+M=N f = factors[1] - ind1 = indices(f,1) - ind = dummyind(ind1) - newf = reshape(f, ntuple(d->ind, Val{L})..., ind1, tail(ntuple(d->ind, Val{M}))...) - _kernelfactors((out..., newf), tail(factors)) + newpost = tail(post) + rv = ReshapedVector(pre, f, newpost) + _kernelfactors((pre..., true), newpost, (out..., rv), tail(factors)) end # A variant in which we just need to fill out to N dimensions diff --git a/src/specialty.jl b/src/specialty.jl index 25e417a..88d6440 100644 --- a/src/specialty.jl +++ b/src/specialty.jl @@ -63,7 +63,7 @@ function imgradients{T,N}(img::AbstractArray{T,N}, points::AbstractVector; G end -function imgradients{T,N}(img::AbstractArray{T,N}, kernelfun::Function, border=Pad{:replicate}()) +function imgradients{T,N}(img::AbstractArray{T,N}, kernelfun::Function=KernelFactors.ando3, border=Pad{:replicate}()) extent = map(length, indices(img)) kern1 = kernelfun(Val{N}, 1, map(x->x>1, extent)) S = filter_type(img, kern1) @@ -73,11 +73,12 @@ function imgradients{T,N}(img::AbstractArray{T,N}, kernelfun::Function, border=P end # When all N gradients have been calculated, return the result -_imgradients{T,N}(G::NTuple{N}, imgpad::AbstractArray{T,N}, kernelfun::Function, extent, border::Inner) = G +_imgradients{T,N}(G::NTuple{N}, imgpad::AbstractArray{T,N}, kernelfun::Function, extent, border) = G # Add the next dimension to G -function _imgradients{T,M,N}(G::NTuple{M}, imgpad::AbstractArray{T,N}, kernelfun::Function, extent, border::Inner) +function _imgradients{T,M,N}(G::NTuple{M}, imgpad::AbstractArray{T,N}, kernelfun::Function, extent, border) d = M+1 # the dimension we're working on now kern = kernelfun(Val{N}, d, map(x->x>1, extent)) - _imgradients((G..., imfilter(imgpad, kern, border)), imgpad, kernelfun, extent, border) + out = allocate_output(imgpad, kern, border) + _imgradients((G..., imfilter!(out, imgpad, kern, NoPad(border))), imgpad, kernelfun, extent, border) end diff --git a/test/basic.jl b/test/basic.jl index 2f8bd52..e38def7 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -5,7 +5,7 @@ using ImagesFiltering, OffsetArrays, Base.Test @test eltype(KernelFactors.IIRGaussian(Float32, 3)) == Float32 @test isa(KernelFactors.IIRGaussian([1,2.0f0]), Tuple{KernelFactors.ReshapedVector,KernelFactors.ReshapedVector}) - @test KernelFactors.kernelfactors(([0,3], [1,7])) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2)) + @test KernelFactors.kernelfactors(([0,3], [1,7])) == (reshape([0,3], 1:2, 0:0), reshape([1,7], 0:0, 1:2)) @test KernelFactors.kernelfactors(([0,3], [1,7]')) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2)) # Warnings @@ -30,3 +30,5 @@ using ImagesFiltering, OffsetArrays, Base.Test end @test contains(readstring(fname), "too small for accuracy") end + +nothing diff --git a/test/specialty.jl b/test/specialty.jl index 8f23367..1a8f015 100644 --- a/test/specialty.jl +++ b/test/specialty.jl @@ -155,8 +155,8 @@ using Base.Test end @test_approx_eq KernelFactors.gaussian(2, 9) gaussiancmp(2, -4:4) k = KernelFactors.gaussian((2,3), (9,7)) - @test_approx_eq k[1] gaussiancmp(2, -4:4) - @test_approx_eq k[2] gaussiancmp(3, -3:3) + @test_approx_eq vec(k[1]) gaussiancmp(2, -4:4) + @test_approx_eq vec(k[2]) gaussiancmp(3, -3:3) @test_approx_eq sum(KernelFactors.gaussian(5)) 1 for k = (KernelFactors.gaussian((2,3)), KernelFactors.gaussian([2,3]), KernelFactors.gaussian([2,3], [9,7])) @test_approx_eq sum(k[1]) 1 From d8addaf4f96fbb1a4fe4e963321319d26b315bfa Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 27 Aug 2016 02:50:22 -0500 Subject: [PATCH 16/21] Get rid of size-dependent conversion to StaticArray --- src/imfilter.jl | 25 ------------------------- 1 file changed, 25 deletions(-) diff --git a/src/imfilter.jl b/src/imfilter.jl index 0d6ca59..8e63b6b 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -351,10 +351,6 @@ end function _imfilter_inbounds{TT}(r, ::Type{TT}, out, A::AbstractArray, k::AbstractVector, Rpre::CartesianRange, ind, Rpost::CartesianRange) indsk = indices(k, 1) - if length(indsk) <= 5 - ks, jo = to_static(k) - return @time _imfilter_inbounds(r, TT, out, A, (ks, jo), Rpre, ind, Rpost) - end for Ipost in Rpost for Ipre in Rpre for i in ind @@ -369,27 +365,6 @@ function _imfilter_inbounds{TT}(r, ::Type{TT}, out, A::AbstractArray, k::Abstrac out end -to_static(a::OffsetArray) = _to_static(parent(a), a.offsets[1]) -to_static(a) = _to_static(a, first(indices(a,1))-1) -_to_static(a, o) = convert(SVector{length(a)}, a), o - -function _imfilter_inbounds{TT,L}(r, ::Type{TT}, out, A::AbstractArray, ko::Tuple{SVector{L},Int}, Rpre::CartesianRange, ind, Rpost::CartesianRange) - k, jo = ko - for Ipost in Rpost - for Ipre in Rpre - for i in ind - io = i+jo - tmp = zero(TT) - @unsafe for j = 1:L - tmp += A[Ipre,io+j,Ipost]*k[j] - end - @unsafe out[Ipre,i,Ipost] = tmp - end - end - end - out -end - ### FFT filtering """ From 4713eed09eebfcce987763901dc5877a60f8472f Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 27 Aug 2016 04:00:39 -0500 Subject: [PATCH 17/21] Type and dispatch cleanups --- src/ImagesFiltering.jl | 7 +-- src/border.jl | 42 +++++++--------- src/imfilter.jl | 109 ++++++++++++++++++++--------------------- src/specialty.jl | 2 +- 4 files changed, 76 insertions(+), 84 deletions(-) diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index 39c71a7..0c111a6 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -22,8 +22,9 @@ module Algorithm immutable FFT <: Alg end immutable FIR <: Alg end immutable IIR <: Alg end + immutable Mixed <: Alg end end -using .Algorithm: Alg, FFT, FIR, IIR +using .Algorithm: Alg, FFT, FIR, IIR, Mixed Alg{A<:Alg}(r::AbstractResource{A}) = r.settings @@ -32,8 +33,8 @@ include("kernelfactors.jl") using .KernelFactors: TriggsSdika, IIRFilter, ReshapedVector, iterdims typealias ArrayLike{T} Union{AbstractArray{T}, IIRFilter{T}, ReshapedVector{T}} -typealias ReshapedTriggsSdika{T,N,Npre,V<:TriggsSdika} ReshapedVector{T,N,Npre,V} -typealias AnyTriggs Union{TriggsSdika, ReshapedTriggsSdika} +typealias ReshapedIIR{T,N,Npre,V<:IIRFilter} ReshapedVector{T,N,Npre,V} +typealias AnyIIR Union{IIRFilter, ReshapedIIR} include("kernel.jl") using .Kernel diff --git a/src/border.jl b/src/border.jl index 55b9723..6e4b87c 100644 --- a/src/border.jl +++ b/src/border.jl @@ -78,15 +78,12 @@ Pad the input image by `lo` pixels at the lower edge, and `hi` pixels at the upp Given a filter array `kernel`, determine the amount of padding from the `indices` of `kernel`. """ -(p::Pad{Style,0}){Style}(kernel::AbstractArray) = Pad{Style}(indices(kernel)) -(p::Pad{Style,0}){Style}(kernel::Laplacian) = Pad{Style}(indices(kernel)) -(p::Pad{Style,0}){Style}(factkernel::Tuple) = Pad{Style}(calculate_padding(factkernel)) -(p::Pad{Style,0}){Style}(factkernel::Tuple, img, ::FIR) = p(factkernel) -(p::Pad{Style,0}){Style}(kernel::Laplacian, img, ::FIR) = p(kernel) +(p::Pad{Style,0}){Style}(kernel) = Pad{Style}(calculate_padding(kernel)) +(p::Pad{Style,0}){Style}(kernel, img, ::Alg) = p(kernel) # Padding for FFT: round up to next size expressible as 2^m*3^n -function (p::Pad{Style,0}){Style}(factkernel::Tuple, img, ::FFT) - inds = calculate_padding(factkernel) +function (p::Pad{Style,0}){Style}(kernel, img, ::FFT) + inds = calculate_padding(kernel) newinds = map(padfft, inds, map(length, indices(img))) Pad{Style}(newinds) end @@ -158,10 +155,8 @@ end (::Type{Inner{N}}){N}(lo::AbstractVector, hi::AbstractVector) = Inner{N}((lo...,), (hi...,)) (::Type{Inner})(lo::AbstractVector, hi::AbstractVector) = Inner((lo...,), (hi...,)) # not inferrable -(p::Inner{0})(factkernel::Tuple, img, ::FIR) = p(factkernel) -(p::Inner{0})(factkernel::Tuple, img, ::FFT) = p(factkernel) -(p::Inner{0})(factkernel::Tuple) = Inner(calculate_padding(factkernel)) -(p::Inner{0})(kernel::AbstractArray) = Inner(indices(kernel)) +(p::Inner{0})(kernel, img, ::Alg) = p(kernel) +(p::Inner{0})(kernel) = Inner(calculate_padding(kernel)) padarray(img, border::Inner) = padarray(eltype(img), img, border) padarray{T}(::Type{T}, img::AbstractArray{T}, border::Inner) = copy(img) @@ -188,14 +183,12 @@ Fill{T}(value::T) = Fill{T,0}(value) Fill{T,N}(value::T, lo::Dims{N}, hi::Dims{N}) = Fill{T,N}(value, lo, hi) Fill(value, lo::AbstractVector, hi::AbstractVector) = Fill(value, (lo...,), (hi...,)) Fill{T,N}(value::T, inds::Base.Indices{N}) = Fill{T,N}(value, map(lo,inds), map(hi,inds)) -Fill(value, kernel::AbstractArray) = Fill(value, indices(kernel)) -Fill(value, kernel::Laplacian) = Fill(value, indices(kernel)) -Fill(value, factkernel::Tuple) = Fill(value, calculate_padding(factkernel)) +Fill(value, kernel) = Fill(value, calculate_padding(kernel)) (p::Fill)(kernel) = Fill(p.value, kernel) -(p::Fill)(kernel, img, ::FIR) = Fill(p.value, kernel) -function (p::Fill)(factkernel::Tuple, img, ::FFT) - inds = calculate_padding(factkernel) +(p::Fill)(kernel, img, ::Alg) = Fill(p.value, kernel) +function (p::Fill)(kernel, img, ::FFT) + inds = calculate_padding(kernel) newinds = map(padfft, inds, map(length, indices(img))) Fill(p.value, newinds) end @@ -260,23 +253,24 @@ end # @inline _flatten(t1, t...) = (t1, flatten(t)...) # _flatten() = () +calculate_padding(kernel) = indices(kernel) @inline function calculate_padding(kernel::Tuple{Any, Vararg{Any}}) inds = accumulate_padding(indices(kernel[1]), tail(kernel)...) - if hastriggs(kernel) && hasfir(kernel) + if hasiir(kernel) && hasfir(kernel) return map(doublepadding, inds) end inds end -hastriggs(kernel) = _hastriggs(false, kernel...) -_hastriggs(ret) = ret -_hastriggs(ret, kern, kernel...) = _hastriggs(ret, kernel...) -_hastriggs(ret, kern::AnyTriggs, kernel...) = true +hasiir(kernel) = _hasiir(false, kernel...) +_hasiir(ret) = ret +_hasiir(ret, kern, kernel...) = _hasiir(ret, kernel...) +_hasiir(ret, kern::AnyIIR, kernel...) = true hasfir(kernel) = _hasfir(false, kernel...) _hasfir(ret) = ret _hasfir(ret, kern, kernel...) = true -_hasfir(ret, kern::AnyTriggs, kernel...) = _hasfir(ret, kernel...) +_hasfir(ret, kern::AnyIIR, kernel...) = _hasfir(ret, kernel...) function doublepadding(ind::AbstractUnitRange) f, l = first(ind), last(ind) @@ -304,7 +298,7 @@ arraytype(A::BitArray, ::Type{Bool}) = BitArray interior(r::AbstractResource, A::AbstractArray, kernel) = interior(A, kernel) -interior(A::AbstractArray, kernel::Union{ArrayLike,Laplacian}) = _interior(indices(A), indices(kernel)) +interior(A, kernel) = _interior(indices(A), indices(kernel)) interior(A, factkernel::Tuple) = _interior(indices(A), accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) function _interior{N}(indsA::NTuple{N}, indsk) indskN = fill_to_length(indsk, 0:0, Val{N}) diff --git a/src/imfilter.jl b/src/imfilter.jl index 8e63b6b..db9f25d 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -211,7 +211,7 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke imfilter!(r, out, A, samedims(out, kern), border) end -const fillbuf_nan = Ref(false) +const fillbuf_nan = Ref(false) # used only for testing purposes function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) kern = kernel[1] @@ -243,7 +243,7 @@ end # However, here it's important to filter over the whole passed-in # range, and then copy! to out -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{AnyTriggs}, border::NoPad, inds::Indices) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{AnyIIR}, border::NoPad, inds::Indices) if inds != indices(out) imfilter!(r, A2, A1, kernel[1], border, inds) R = CartesianRange(indices(out)) @@ -265,30 +265,30 @@ end ## FIR filtering """ - imfilter!(::AbstractResource{FIR}, imgfilt, img, kernel, NoPad(), [inds=indices(imgfilt)]) + imfilter!(::AbstractResource, imgfilt, img, kernel, NoPad(), [inds=indices(imgfilt)]) Filter an array `img` with kernel `kernel` by computing their -correlation, storing the result in `imgfilt`, using a finite-impulse +correlation, storing the result in `imgfilt`, defaulting to a finite-impulse response (FIR) algorithm. Any necessary padding must have already been supplied to `img`. If you want padding applied, instead call - imfilter!(::AbstractResource{FIR}, imgfilt, img, kernel, border) + imfilter!([r::AbstractResource,] imgfilt, img, kernel, border) with a specific `border`, or use - imfilter!(imgfilt, img, kernel, Algorithm.FIR()) + imfilter!(imgfilt, img, kernel, [Algorithm.FIR()]) for default padding. If `inds` is supplied, only the elements of `imgfilt` with indices in the domain of `inds` will be calculated. This can be particularly -useful for "cascaded filters" where you pad over a larger area and -then calculate the result over just the necessary region at each -stage. +useful for "cascaded FIR filters" where you pad over a larger area and +then calculate the result over just the necessary/well-defined region +at each successive stage. See also: imfilter. """ -function imfilter!{S,T,N}(r::CPU1{FIR}, +function imfilter!{S,T,N}(r::AbstractResource, out::AbstractArray{S,N}, A::AbstractArray{T,N}, kern::NDimKernel{N}, @@ -315,17 +315,12 @@ function imfilter!{S,T,N}(r::CPU1{FIR}, _imfilter_inbounds!(r, out, A, kern, border, inds) end -function _imfilter_inbounds!{T,N,Npre,V<:TriggsSdika}(r::AbstractResource, out, A::AbstractArray, kern::ReshapedVector{T,N,Npre,V}, border::NoPad, inds) - indspre, ind, indspost = iterdims(inds, kern) - _imfilter_dim!(r, out, A, kern.data, CartesianRange(indspre), ind, CartesianRange(indspost), border[]) -end - -function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractVector, kern::TriggsSdika, border::NoPad, inds) +function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, kern::ReshapedIIR, border::NoPad, inds) indspre, ind, indspost = iterdims(inds, kern) _imfilter_dim!(r, out, A, kern.data, CartesianRange(indspre), ind, CartesianRange(indspost), border[]) end -function _imfilter_inbounds!(r, out, A::AbstractArray, kern, border::NoPad, inds) +function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, kern, border::NoPad, inds) indsk = indices(kern) R, Rk = CartesianRange(inds), CartesianRange(indsk) p = A[first(R)+first(Rk)] * first(kern) @@ -340,7 +335,7 @@ function _imfilter_inbounds!(r, out, A::AbstractArray, kern, border::NoPad, inds out end -function _imfilter_inbounds!(r, out, A::AbstractArray, kern::ReshapedVector, border::NoPad, inds) +function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, kern::ReshapedVector, border::NoPad, inds) Rpre, ind, Rpost = iterdims(inds, kern) k = kern.data R, Rk = CartesianRange(inds), CartesianRange(indices(kern)) @@ -349,7 +344,7 @@ function _imfilter_inbounds!(r, out, A::AbstractArray, kern::ReshapedVector, bor _imfilter_inbounds(r, TT, out, A, k, Rpre, ind, Rpost) end -function _imfilter_inbounds{TT}(r, ::Type{TT}, out, A::AbstractArray, k::AbstractVector, Rpre::CartesianRange, ind, Rpost::CartesianRange) +function _imfilter_inbounds{TT}(r::AbstractResource, ::Type{TT}, out, A::AbstractArray, k::AbstractVector, Rpre::CartesianRange, ind, Rpost::CartesianRange) indsk = indices(k, 1) for Ipost in Rpost for Ipre in Rpre @@ -650,6 +645,43 @@ function rightΔu{T,l}(img, uplus, Ibegin, i, Iend, kernel::TriggsSdika{T,3,l}) ret end +### Pad{:na}() boundary conditions + +function imfilter_na_inseparable!{T}(r, out::AbstractArray{T}, img, nanflag, kernel::Tuple{Vararg{TriggsSdika}}) + fc, fn = Fill(zero(T)), Fill(zero(eltype(T))) # color, numeric + copy!(out, img) + out[nanflag] = zero(T) + validpixels = copy!(similar(Array{eltype(T)}, indices(img)), mappedarray(x->!x, nanflag)) + # TriggsSdika is safe for inplace operations + imfilter!(r, out, out, kernel, fc) + imfilter!(r, validpixels, validpixels, kernel, fn) + for I in eachindex(out) + out[I] /= validpixels[I] + end + out[nanflag] = nan(T) + out +end + +function imfilter_na_inseparable!{T}(r, out::AbstractArray{T}, img, nanflag, kernel::Tuple) + fc, fn = Fill(zero(T)), Fill(zero(eltype(T))) # color, numeric + imgtmp = copy!(similar(out, indices(img)), img) + imgtmp[nanflag] = zero(T) + validpixels = copy!(similar(Array{eltype(T)}, indices(img)), mappedarray(x->!x, nanflag)) + imfilter!(r, out, imgtmp, kernel, fc) + vp = imfilter(r, validpixels, kernel, fn) + for I in eachindex(out) + out[I] /= vp[I] + end + out[nanflag] = nan(T) + out +end + +function imfilter_na_separable!{T}(r, out::AbstractArray{T}, img, kernel::Tuple) + fc, fn = Fill(zero(T)), Fill(zero(eltype(T))) # color, numeric + imfilter!(r, out, img, kernel, fc) + normalize_separable!(r, out, kernel, fn) +end + ### Utilities filter_type{S}(img::AbstractArray{S}, kernel) = filter_type(S, kernel) @@ -733,44 +765,9 @@ end filter_algorithm(out, img, kernel) = FIR() filter_algorithm(out, img, kernel::Tuple{IIRFilter,Vararg{IIRFilter}}) = IIR() -isseparable(kernels::Tuple{Vararg{TriggsSdika}}) = true +isseparable(kernels::Tuple{Vararg{AnyIIR}}) = true isseparable(kernels::Tuple) = all(x->nextendeddims(x)==1, kernels) -function imfilter_na_inseparable!{T}(r, out::AbstractArray{T}, img, nanflag, kernel::Tuple{Vararg{TriggsSdika}}) - fc, fn = Fill(zero(T)), Fill(zero(eltype(T))) # color, numeric - copy!(out, img) - out[nanflag] = zero(T) - validpixels = copy!(similar(Array{eltype(T)}, indices(img)), mappedarray(x->!x, nanflag)) - # TriggsSdika is safe for inplace operations - imfilter!(r, out, out, kernel, fc) - imfilter!(r, validpixels, validpixels, kernel, fn) - for I in eachindex(out) - out[I] /= validpixels[I] - end - out[nanflag] = nan(T) - out -end - -function imfilter_na_inseparable!{T}(r, out::AbstractArray{T}, img, nanflag, kernel::Tuple) - fc, fn = Fill(zero(T)), Fill(zero(eltype(T))) # color, numeric - imgtmp = copy!(similar(out, indices(img)), img) - imgtmp[nanflag] = zero(T) - validpixels = copy!(similar(Array{eltype(T)}, indices(img)), mappedarray(x->!x, nanflag)) - imfilter!(r, out, imgtmp, kernel, fc) - vp = imfilter(r, validpixels, kernel, fn) - for I in eachindex(out) - out[I] /= vp[I] - end - out[nanflag] = nan(T) - out -end - -function imfilter_na_separable!{T}(r, out::AbstractArray{T}, img, kernel::Tuple) - fc, fn = Fill(zero(T)), Fill(zero(eltype(T))) # color, numeric - imfilter!(r, out, img, kernel, fc) - normalize_separable!(r, out, kernel, fn) -end - normalize_separable!(r::AbstractResource, A, ::Tuple{}, border) = error("this shouldn't happen") function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,TriggsSdika}, border) inds = indices(A) @@ -780,7 +777,7 @@ function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,Trigg filtdims = ntuple(d->imfilter_inplace!(r, similar(dims->ones(dims), inds[d]), kernels[d], border), Val{N}) normalize_dims!(A, filtdims) end -function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,ReshapedTriggsSdika}, border) +function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,ReshapedIIR}, border) normalize_separable!(r, A, map(_vec, kernels), border) end diff --git a/src/specialty.jl b/src/specialty.jl index 88d6440..3a88bd2 100644 --- a/src/specialty.jl +++ b/src/specialty.jl @@ -1,6 +1,6 @@ ## Laplacian -function _imfilter_inbounds!(r, out, A::AbstractArray, L::Laplacian, border::NoPad, inds) +function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, L::Laplacian, border::NoPad, inds) TT = eltype(out) # accumtype(eltype(out), eltype(A)) n = 2*length(L.offsets) R = CartesianRange(inds) From 38970f8689ced6d1254f863b1e78a3320c6ed83a Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 27 Aug 2016 04:15:01 -0500 Subject: [PATCH 18/21] Choose FFT/FIR algorithm based on kernel size --- src/imfilter.jl | 20 ++++++++++++++++++-- test/basic.jl | 5 ++++- 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/imfilter.jl b/src/imfilter.jl index db9f25d..9c88402 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -762,8 +762,24 @@ function kernelshift(inds::Indices, A) A end -filter_algorithm(out, img, kernel) = FIR() -filter_algorithm(out, img, kernel::Tuple{IIRFilter,Vararg{IIRFilter}}) = IIR() +# Note this is not type-stable. Fortunately, all the outputs are +# allocated by the time this gets called. +function filter_algorithm(out, img, kernel::Union{AbstractArray,Tuple{Vararg{AbstractArray}}}) + L = maxlen(kernel) + if L > 30 + return FFT() + end + FIR() +end +filter_algorithm(out, img, kernel::Tuple{AnyIIR,Vararg{AnyIIR}}) = IIR() +filter_algorithm(out, img, kernel) = Mixed() + +maxlen(A::AbstractArray) = _length(A) +@inline maxlen(kernel::Tuple) = _maxlen(0, kernel...) +_maxlen(len, kernel1, kernel...) = _maxlen(max(len, _length(kernel1)), kernel...) +_maxlen(len) = len + +_length(A::AbstractArray) = length(linearindices(A)) isseparable(kernels::Tuple{Vararg{AnyIIR}}) = true isseparable(kernels::Tuple) = all(x->nextendeddims(x)==1, kernels) diff --git a/test/basic.jl b/test/basic.jl index e38def7..257d4ab 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -3,7 +3,10 @@ using ImagesFiltering, OffsetArrays, Base.Test @testset "basic" begin @test eltype(KernelFactors.IIRGaussian(3)) == Float64 @test eltype(KernelFactors.IIRGaussian(Float32, 3)) == Float32 - @test isa(KernelFactors.IIRGaussian([1,2.0f0]), Tuple{KernelFactors.ReshapedVector,KernelFactors.ReshapedVector}) + kern = KernelFactors.IIRGaussian([1,2.0f0]) + @test isa(kern, Tuple{KernelFactors.ReshapedVector,KernelFactors.ReshapedVector}) + @test isa(ImagesFiltering.filter_algorithm([1],[1],kern), Algorithm.IIR) + @test isa(ImagesFiltering.filter_algorithm([1],[1],(kern...,Kernel.Laplacian())), Algorithm.Mixed) @test KernelFactors.kernelfactors(([0,3], [1,7])) == (reshape([0,3], 1:2, 0:0), reshape([1,7], 0:0, 1:2)) @test KernelFactors.kernelfactors(([0,3], [1,7]')) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2)) From a72ed22b29fed0c9cca452dca4849d3440908d66 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 27 Aug 2016 07:28:28 -0500 Subject: [PATCH 19/21] Introduce `expand` and `shrink` for more intuitive handling of filter cascades --- src/border.jl | 69 +++++++++++++++++++++++++++++++++---------------- src/imfilter.jl | 4 +-- 2 files changed, 49 insertions(+), 24 deletions(-) diff --git a/src/border.jl b/src/border.jl index 6e4b87c..47cd296 100644 --- a/src/border.jl +++ b/src/border.jl @@ -248,16 +248,11 @@ function extend(lo::Integer, inds::AbstractUnitRange, hi::Integer) OffsetArray(newind, newind) end -# @inline flatten(t::Tuple) = _flatten(t...) -# @inline _flatten(t1::Tuple, t...) = (flatten(t1)..., flatten(t)...) -# @inline _flatten(t1, t...) = (t1, flatten(t)...) -# _flatten() = () - calculate_padding(kernel) = indices(kernel) @inline function calculate_padding(kernel::Tuple{Any, Vararg{Any}}) inds = accumulate_padding(indices(kernel[1]), tail(kernel)...) if hasiir(kernel) && hasfir(kernel) - return map(doublepadding, inds) + inds = map(doublepadding, inds) end inds end @@ -280,15 +275,8 @@ function doublepadding(ind::AbstractUnitRange) end accumulate_padding(inds::Indices, kernel1, kernels...) = - accumulate_padding(_accumulate_padding(inds, indices(kernel1)), kernels...) -accumulate_padding(inds::Indices, kernel1::TriggsSdika, kernels...) = - accumulate_padding(inds, kernels...) -accumulate_padding(inds) = inds -_accumulate_padding(inds1, inds2) = (__accumulate_padding(inds1[1], inds2[1]), _accumulate_padding(tail(inds1), tail(inds2))...) -_accumulate_padding(::Tuple{}, ::Tuple{}) = () -_accumulate_padding(::Tuple{}, inds2) = inds2 -_accumulate_padding(inds1, ::Tuple{}) = inds1 -__accumulate_padding(ind1, ind2) = first(ind1)+min(0,first(ind2)):last(ind1)+max(0,last(ind2)) + accumulate_padding(expand(inds, indices(kernel1)), kernels...) +accumulate_padding(inds::Indices) = inds modrange(x, r::AbstractUnitRange) = mod(x-first(r), length(r))+first(r) modrange(A::AbstractArray, r::AbstractUnitRange) = map(x->modrange(x, r), A) @@ -296,22 +284,42 @@ modrange(A::AbstractArray, r::AbstractUnitRange) = map(x->modrange(x, r), A) arraytype{T}(A::AbstractArray, ::Type{T}) = Array{T} # fallback arraytype(A::BitArray, ::Type{Bool}) = BitArray -interior(r::AbstractResource, A::AbstractArray, kernel) = interior(A, kernel) - interior(A, kernel) = _interior(indices(A), indices(kernel)) interior(A, factkernel::Tuple) = _interior(indices(A), accumulate_padding(indices(factkernel[1]), tail(factkernel)...)) function _interior{N}(indsA::NTuple{N}, indsk) indskN = fill_to_length(indsk, 0:0, Val{N}) - map((ia,ik)->first(ia) + lo(ik) : last(ia) - hi(ik), indsA, indskN) + map(intersect, indsA, shrink(indsA, indsk)) end -next_interior(inds::Indices, ::Tuple{}) = inds -function next_interior(inds::Indices, kernel::Tuple) +next_shrink(inds::Indices, ::Tuple{}) = inds +function next_shrink(inds::Indices, kernel::Tuple) kern = first(kernel) - iscopy(kern) && return next_interior(inds, tail(kernel)) - _interior(inds, indices(kern)) + iscopy(kern) && return next_shrink(inds, tail(kernel)) + shrink(inds, kern) end +""" + expand(inds::Indices, kernel) + expand(inds::Indices, indskernel::Indices) + +Expand an image region `inds` to account for necessary padding by `kernel`. +""" +expand(inds::Indices, A) = expand(inds, indices(A)) +expand(inds::Indices, pad::Indices) = firsttype(map_copytail(expand, inds, pad)) +expand(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)+first(pad):last(ind)+last(pad)) +expand(ind::Base.OneTo, pad::AbstractUnitRange) = expand(UnitRange(ind), pad) + +""" + shrink(inds::Indices, kernel) + shrink(inds::Indices, indskernel) + +Remove edges from an image region `inds` that correspond to padding needed for `kernel`. +""" +shrink(inds::Indices, A) = shrink(inds, indices(A)) +shrink(inds::Indices, pad::Indices) = firsttype(map_copytail(shrink, inds, pad)) +shrink(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)-first(pad):last(ind)-last(pad)) +shrink(ind::Base.OneTo, pad::AbstractUnitRange) = shrink(UnitRange(ind), pad) + allocate_output{T}(::Type{T}, img, kernel, border) = similar(img, T) function allocate_output{T}(::Type{T}, img, kernel, ::Inner{0}) inds = interior(img, kernel) @@ -319,4 +327,21 @@ function allocate_output{T}(::Type{T}, img, kernel, ::Inner{0}) end allocate_output(img, kernel, border) = allocate_output(filter_type(img, kernel), img, kernel, border) +""" + map_copytail(f, a::Tuple, b::Tuple) + +Apply `f` to paired elements of `a` and `b`, copying any tail elements +when the lengths of `a` and `b` are not equal. +""" +@inline map_copytail(f, a::Tuple, b::Tuple) = (f(a[1], b[1]), map_copytail(f, tail(a), tail(b))...) +map_copytail(f, a::Tuple{}, b::Tuple{}) = () +map_copytail(f, a::Tuple, b::Tuple{}) = a +map_copytail(f, a::Tuple{}, b::Tuple) = b + +function firsttype(t::Tuple) + T = typeof(t[1]) + map(T, t) +end +firsttype(::Tuple{}) = () + # end diff --git a/src/imfilter.jl b/src/imfilter.jl index 9c88402..0f7df98 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -223,7 +223,7 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke if fillbuf_nan[] fill!(A2, NaN) # for testing purposes end - inds = interior(r, A, kern) + inds = shrink(indices(A), kern) _imfilter!(r, out, A, A2, kernel, border, inds) return out end @@ -258,7 +258,7 @@ function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{ kernN = samedims(A2, kern) imfilter!(r, A2, A1, kernN, border, inds) # store result in A2 kernelt = tail(kernel) - newinds = next_interior(inds, kernelt) + newinds = next_shrink(inds, kernelt) _imfilter!(r, out, A2, A1, tail(kernel), border, newinds) # swap the buffers end From d11b60b311689b96e3508b338d33b11d2689ae3c Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 28 Aug 2016 06:40:23 -0500 Subject: [PATCH 20/21] Add tiled, multithreaded FIR filtering The performance is very good; at this point the main avenue for improvement is to implement lazy padding. --- src/ImagesFiltering.jl | 2 +- src/border.jl | 4 +- src/imfilter.jl | 124 ++++++++++++++++++++++++++++++++++++++++- src/kernel.jl | 2 +- src/kernelfactors.jl | 8 +++ test/2d.jl | 1 + 6 files changed, 136 insertions(+), 5 deletions(-) diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index 0c111a6..ffb5994 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -1,6 +1,6 @@ module ImagesFiltering -using Colors, FixedPointNumbers, ImagesCore, MappedArrays, FFTViews, OffsetArrays, StaticArrays, ComputationalResources +using Colors, FixedPointNumbers, ImagesCore, MappedArrays, FFTViews, OffsetArrays, StaticArrays, ComputationalResources, TiledIteration using ColorVectorSpace # in case someone filters RGB arrays using Base: Indices, tail, fill_to_length, @pure diff --git a/src/border.jl b/src/border.jl index 47cd296..938bf77 100644 --- a/src/border.jl +++ b/src/border.jl @@ -304,7 +304,7 @@ end Expand an image region `inds` to account for necessary padding by `kernel`. """ -expand(inds::Indices, A) = expand(inds, indices(A)) +expand(inds::Indices, A) = expand(inds, calculate_padding(A)) expand(inds::Indices, pad::Indices) = firsttype(map_copytail(expand, inds, pad)) expand(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)+first(pad):last(ind)+last(pad)) expand(ind::Base.OneTo, pad::AbstractUnitRange) = expand(UnitRange(ind), pad) @@ -315,7 +315,7 @@ expand(ind::Base.OneTo, pad::AbstractUnitRange) = expand(UnitRange(ind), pad) Remove edges from an image region `inds` that correspond to padding needed for `kernel`. """ -shrink(inds::Indices, A) = shrink(inds, indices(A)) +shrink(inds::Indices, A) = shrink(inds, calculate_padding(A)) shrink(inds::Indices, pad::Indices) = firsttype(map_copytail(shrink, inds, pad)) shrink(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)-first(pad):last(ind)-last(pad)) shrink(ind::Base.OneTo, pad::AbstractUnitRange) = shrink(UnitRange(ind), pad) diff --git a/src/imfilter.jl b/src/imfilter.jl index 0f7df98..1b167a0 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -262,7 +262,114 @@ function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{ _imfilter!(r, out, A2, A1, tail(kernel), border, newinds) # swap the buffers end -## FIR filtering +### When everything is FIR, we can use a tiled algorithm +function imfilter!(r::AbstractCPU{FIR}, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) + kern = kernel[1] + iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border) + tmp = tile_allocate(filter_type(A, kernel), kernel) + _imfilter_tiled!(r, out, A, kernel, border, tmp) + out +end + +# Single-threaded, pair of kernels (with only one temporary buffer required) +function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}) + k1, k2 = kernel + tile = tiles[1] + indsk1, indstile = indices(k1), indices(tile) + sz = map(length, indstile) + chunksz = map(length, shrink(indstile, indsk1)) + inds = expand(indices(out), k2) + for tileinds in TileIterator(inds, chunksz) + tileb = TileBuffer(tile, tileinds) + imfilter!(r, tileb, A, samedims(tileb, k1), border, tileinds) + imfilter!(r, out, tileb, samedims(out, k2), border, shrink(tileinds, k2)) + end + out +end + +# Multithreaded, pair of kernels +function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}) + k1, k2 = kernel + tile = tiles[1] + indsk1, indstile = indices(k1), indices(tile) + sz = map(length, indstile) + chunksz = map(length, shrink(indstile, indsk1)) + tileinds_all = collect(TileIterator(expand(indices(out), k2), chunksz)) + _imfilter_tiled_threads!(CPU1(r), out, A, samedims(out, k1), samedims(out, k2), border, tileinds_all, tiles) +end +# This must be in a separate function due to #15276 +@noinline function _imfilter_tiled_threads!{AA<:AbstractArray}(r1, out, A, k1, k2, border, tileinds_all, tile::Vector{AA}) + Threads.@threads for i = 1:length(tileinds_all) + id = Threads.threadid() + tileinds = tileinds_all[i] + tileb = TileBuffer(tile[id], tileinds) + imfilter!(r1, tileb, A, k1, border, tileinds) + imfilter!(r1, out, tileb, k2, border, shrink(tileinds, k2)) + end + out +end + +# Single-threaded, multiple kernels (requires two tile buffers, swapping on each iteration) +function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}) + k1, kt = kernel[1], tail(kernel) + tilepair = tiles[1] + indsk1, indstile = indices(k1), indices(tilepair[1]) + sz = map(length, indstile) + chunksz = map(length, shrink(indstile, indsk1)) + inds = expand(indices(out), kt) + for tileinds in TileIterator(inds, chunksz) + tileb1 = TileBuffer(tilepair[1], tileinds) + imfilter!(r, tileb1, A, samedims(tileb1, k1), border, tileinds) + _imfilter_tiled_swap!(r, out, kt, border, (tileb1, tilepair[2])) + end +end + +# Multithreaded, multiple kernels +function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}) + k1, kt = kernel[1], tail(kernel) + tilepair = tiles[1] + indsk1, indstile = indices(k1), indices(tilepair[1]) + sz = map(length, indstile) + chunksz = map(length, shrink(indstile, indsk1)) + tileinds_all = collect(TileIterator(expand(indices(out), kt), chunksz)) + _imfilter_tiled_threads!(CPU1(r), out, A, samedims(out, k1), kt, border, tileinds_all, tiles) +end +# This must be in a separate function due to #15276 +@noinline function _imfilter_tiled_threads!{AA<:AbstractArray}(r1, out, A, k1, kt, border, tileinds_all, tiles::Vector{Tuple{AA,AA}}) + Threads.@threads for i = 1:length(tileinds_all) + tileinds = tileinds_all[i] + id = Threads.threadid() + tile1, tile2 = tiles[id] + tileb1 = TileBuffer(tile1, tileinds) + imfilter!(r1, tileb1, A, k1, border, tileinds) + _imfilter_tiled_swap!(r1, out, kt, border, (tileb1, tile2)) + end + out +end + +# The first of the pair in `tmp` has the current data. We also make +# the second a plain array so there's no doubt about who's holding the +# proper indices. +function _imfilter_tiled_swap!(r, out, kernel::Tuple{Any,Any,Vararg{Any}}, border, tmp::Tuple{TileBuffer,Array}) + tileb1, tile2 = tmp + k1, kt = kernel[1], tail(kernel) + parentinds = indices(tileb1) + tileinds = shrink(parentinds, k1) + tileb2 = TileBuffer(tile2, tileinds) + imfilter!(r, tileb2, tileb1, samedims(tileb2, k1), border, tileinds) + _imfilter_tiled_swap!(r, out, kt, border, (tileb2, parent(tileb1))) +end + +# on the last call we write to `out` instead of one of the buffers +function _imfilter_tiled_swap!(r, out, kernel::Tuple{Any}, border, tmp::Tuple{TileBuffer,Array}) + tileb1 = tmp[1] + k1 = kernel[1] + parentinds = indices(tileb1) + tileinds = shrink(parentinds, k1) + imfilter!(r, out, tileb1, samedims(out, k1), border, tileinds) +end + +### FIR filtering """ imfilter!(::AbstractResource, imgfilt, img, kernel, NoPad(), [inds=indices(imgfilt)]) @@ -818,3 +925,18 @@ iscopy(kernel::AbstractArray) = all(x->x==0:0, indices(kernel)) && first(kernel) iscopy(kernel::Laplacian) = false iscopy(kernel::TriggsSdika) = all(x->x==0, kernel.a) && all(x->x==0, kernel.b) && kernel.scale == 1 iscopy(kernel::ReshapedVector) = iscopy(kernel.data) + +## Tiling utilities + +function tile_allocate{T}(::Type{T}, kernel::Tuple{Any,Any}) + sz = map(length, calculate_padding(kernel)) + tsz = padded_tilesize(T, sz) + [Array{T}(tsz) for i = 1:Threads.nthreads()] +end + +function tile_allocate{T}(::Type{T}, kernel::Tuple{Any,Any,Vararg{Any}}) + # Allocate a pair of tiles and swap buffers at each stage + sz = map(length, calculate_padding(kernel)) + tsz = padded_tilesize(T, sz) + [(Array{T}(tsz), Array{T}(tsz)) for i = 1:Threads.nthreads()] +end diff --git a/src/kernel.jl b/src/kernel.jl index a2fe55f..d711387 100644 --- a/src/kernel.jl +++ b/src/kernel.jl @@ -237,7 +237,7 @@ end Base.indices(L::Laplacian) = map(f->f ? (-1:1) : (0:0), L.flags) Base.isempty(L::Laplacian) = false function Base.convert{N}(::Type{AbstractArray}, L::Laplacian{N}) - A = zeros(Int, indices(L)...) + A = fill!(OffsetArray{Int}(indices(L)), 0) for I in L.offsets A[I] = A[-I] = 1 end diff --git a/src/kernelfactors.jl b/src/kernelfactors.jl index add1b37..fb12bf4 100644 --- a/src/kernelfactors.jl +++ b/src/kernelfactors.jl @@ -73,6 +73,14 @@ import Base: == ==(A::ReshapedVector, B::AbstractArray) = convert(AbstractArray, A) == B ==(A::AbstractArray, B::ReshapedVector) = A == convert(AbstractArray, B) +import Base: .+, .-, .*, ./ +for op in (:+, :-, :*, :/) + opdot = Symbol(:., op) + @eval begin + ($opdot)(A::ReshapedVector, B::ReshapedVector) = broadcast($op, A, B) + @inline ($opdot)(As::ReshapedVector...) = broadcast($op, As...) + end +end # for broadcasting @pure Base.promote_eltype_op{S,T}(op, ::ReshapedVector{S}, ::ReshapedVector{T}) = Base.promote_op(op, S, T) diff --git a/test/2d.jl b/test/2d.jl index 6a3b631..302218b 100644 --- a/test/2d.jl +++ b/test/2d.jl @@ -55,6 +55,7 @@ using Base.Test for img in (imgf, imgi, imgg, imgc) targetimg = zeros(typeof(img[1]*kern[1]), size(img)) targetimg[3:4,2:3] = rot180(kern)*img[3,4] + ret = similar(targetimg) @test @inferred(imfilter(img, kernel)) ≈ targetimg @test @inferred(imfilter(f32type(img), img, kernel)) ≈ float32(targetimg) for border in (Pad{:replicate}(), Pad{:circular}(), Pad{:symmetric}(), Pad{:reflect}(), Fill(zero(eltype(img)))) From 4f58febc1207eaa03c183e7925ba1f9ad10823eb Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 30 Aug 2016 11:48:17 -0500 Subject: [PATCH 21/21] Add lazy-padding (with edge-handling) for FIR filtering This means that padding is no longer necessary for pure-FIR filters (and pure-IIR filters). The most complex part of this was getting the indices-handling correct. --- src/ImagesFiltering.jl | 5 +- src/border.jl | 4 +- src/imfilter.jl | 168 +++++++++++++++++++++++++++++------------ 3 files changed, 124 insertions(+), 53 deletions(-) diff --git a/src/ImagesFiltering.jl b/src/ImagesFiltering.jl index ffb5994..498c580 100644 --- a/src/ImagesFiltering.jl +++ b/src/ImagesFiltering.jl @@ -46,8 +46,9 @@ include("border.jl") typealias BorderSpec{Style,T} Union{Pad{Style,0}, Fill{T,0}, Inner{0}} typealias BorderSpecRF{T} Union{Pad{:replicate,0}, Fill{T,0}} -typealias BorderSpecNoNa{T} Union{Pad{:replicate,0}, Pad{:circular,0}, Pad{:symmetric,0}, - Pad{:reflect,0}, Fill{T,0}, Inner{0}} +typealias PadNoNa Union{Pad{:replicate,0}, Pad{:circular,0}, Pad{:symmetric,0}, + Pad{:reflect, 0}} +typealias BorderSpecNoNa{T} Union{PadNoNa, Fill{T,0}, Inner{0}} typealias BorderSpecAny Union{BorderSpec,NoPad} include("deprecated.jl") diff --git a/src/border.jl b/src/border.jl index 938bf77..a245823 100644 --- a/src/border.jl +++ b/src/border.jl @@ -304,7 +304,7 @@ end Expand an image region `inds` to account for necessary padding by `kernel`. """ -expand(inds::Indices, A) = expand(inds, calculate_padding(A)) +expand(inds::Indices, kernel) = expand(inds, calculate_padding(kernel)) expand(inds::Indices, pad::Indices) = firsttype(map_copytail(expand, inds, pad)) expand(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)+first(pad):last(ind)+last(pad)) expand(ind::Base.OneTo, pad::AbstractUnitRange) = expand(UnitRange(ind), pad) @@ -315,7 +315,7 @@ expand(ind::Base.OneTo, pad::AbstractUnitRange) = expand(UnitRange(ind), pad) Remove edges from an image region `inds` that correspond to padding needed for `kernel`. """ -shrink(inds::Indices, A) = shrink(inds, calculate_padding(A)) +shrink(inds::Indices, kernel) = shrink(inds, calculate_padding(kernel)) shrink(inds::Indices, pad::Indices) = firsttype(map_copytail(shrink, inds, pad)) shrink(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)-first(pad):last(ind)-last(pad)) shrink(ind::Base.OneTo, pad::AbstractUnitRange) = shrink(UnitRange(ind), pad) diff --git a/src/imfilter.jl b/src/imfilter.jl index 1b167a0..d27b5b1 100644 --- a/src/imfilter.jl +++ b/src/imfilter.jl @@ -134,14 +134,23 @@ end """ imfilter!(imgfilt, img, kernel, [border=Pad{:replicate}()], [alg]) - imfilter!(r, imgfilt, img, kernel, border) + imfilter!(r, imgfilt, img, kernel, border, [inds]) + imfilter!(r, imgfilt, img, kernel, border::NoPad, [inds=indices(imgfilt)]) Filter an array `img` with kernel `kernel` by computing their correlation, storing the result in `imgfilt`. The indices of `imgfilt` determine the region over which the filtered -image is computed---you can use this fact to select just a specific region -of interest. +image is computed---you can use this fact to select just a specific +region of interest, although be aware that the input `img` might still +get padded. Alteratively, explicitly provide the indices `inds` of +`imgfilt` that you want to calculate, and use `NoPad` boundary +conditions. In such cases, you are responsible for supplying +appropriate padding: `img` must be indexable for all of the locations +needed for calculating the output. This syntax is best-supported for +FIR filtering; in particular, that that IIR filtering can lead to +results that are inconsistent with respect to filtering the entire +array. See also: imfilter. """ @@ -200,22 +209,47 @@ function imfilter!{S,T,N}(r::AbstractResource, imfilter!(r, out, A, kernel, NoPad(border)) end +# An optimized case that performs only "virtual padding" +function imfilter!{S,T,N}(r::AbstractCPU{FIR}, + out::AbstractArray{S,N}, + img::AbstractArray{T,N}, + kernel::ProcessedKernel, + border::PadNoNa) + # The fast path: handle the points that don't need padding + iinds = map(intersect, interior(img, kernel), indices(out)) + imfilter!(r, out, img, kernel, NoPad(border), iinds) + # The not-so-fast path: handle the edges + padded = view(img, padindices(img, border(kernel))...) + pkernel = kernelconv(kernel...) + _imfilter_iter!(r, out, padded, pkernel, EdgeIterator(indices(out), iinds)) +end + +### "Scheduler" methods (all with NoPad) + +# These methods handle much of what Halide calls "the schedule." +# Together they handle the order-of-operations for separable and/or +# cascaded kernels, and even implement multithreadable tiling for FIR +# filtering. + +# Trivial kernel (a copy operation) function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{}, ::NoPad, inds::Indices=indices(out)) R = CartesianRange(inds) copy!(out, R, A, R) end -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any}, border::NoPad) +# A single kernel +function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any}, border::NoPad, inds::Indices=indices(out)) kern = kernel[1] - iscopy(kern) && return imfilter!(r, out, A, (), border) - imfilter!(r, out, A, samedims(out, kern), border) + iscopy(kern) && return imfilter!(r, out, A, (), border, inds) + imfilter!(r, out, A, samedims(out, kern), border, inds) end const fillbuf_nan = Ref(false) # used only for testing purposes -function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) +# A filter cascade (2 or more filters) +function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, inds=indices(out)) kern = kernel[1] - iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border) + iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border, inds) # For multiple stages of filtering, we introduce a second buffer # and swap them at each stage. The first of the two is the one # that holds the most recent result. @@ -223,62 +257,70 @@ function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, ke if fillbuf_nan[] fill!(A2, NaN) # for testing purposes end - inds = shrink(indices(A), kern) - _imfilter!(r, out, A, A2, kernel, border, inds) + indsstep = shrink(expand(inds, calculate_padding(kernel)), kern) + _imfilter!(r, out, A, A2, kernel, border, indsstep) return out end -# We deliberately drop `inds` in the next methods: when `kernel` is a tuple -# that has both TriggsSdika and FIR filters, the overall padding -# gets doubled, yet we only trim off the minimum at each -# stage. Consequently, `inds` might be optimistic about the range -# available in `out`. -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, border::NoPad, inds::Indices) - imfilter!(r, out, A1, kernel, border) +### When everything is FIR, we can instead use a tiled algorithm for the cascaded case +function imfilter!{S,T,N}(r::AbstractCPU{FIR}, out::AbstractArray{S,N}, A::AbstractArray{T,N}, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, inds=indices(out)) + kern = kernel[1] + iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border, inds) + tmp = tile_allocate(filter_type(A, kernel), kernel) + _imfilter_tiled!(r, out, A, kernel, border, tmp, inds) + out +end + +### Scheduler support methods + +## No tiling, filter cascade + +# For these (internal) methods, `indsstep` refers to `inds` for this +# step of filtering, not the indices of `out` that we want to finally +# target. + +# When `kernel` is (originally) a tuple that has both TriggsSdika and +# FIR filters, the overall padding gets doubled, yet we only trim off +# the minimum at each stage. Consequently, `indsstep` might be +# optimistic about the range available in `out`; therefore we use +# `intersect`. +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, border::NoPad, indsstep::Indices) + imfilter!(r, out, A1, kernel, border, map(intersect, indsstep, indices(out))) end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, border::NoPad, inds::Indices) - imfilter!(r, out, A1, kernel[1], border) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, border::NoPad, indsstep::Indices) + imfilter!(r, out, A1, kernel[1], border, map(intersect, indsstep, indices(out))) end -# However, here it's important to filter over the whole passed-in -# range, and then copy! to out -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{AnyIIR}, border::NoPad, inds::Indices) - if inds != indices(out) - imfilter!(r, A2, A1, kernel[1], border, inds) - R = CartesianRange(indices(out)) +# For IIR, it's important to filter over the whole passed-in range, +# and then copy! to out +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{AnyIIR}, border::NoPad, indsstep::Indices) + if indsstep != indices(out) + imfilter!(r, A2, A1, kernel[1], border, indsstep) + R = CartesianRange(map(intersect, indsstep, indices(out))) return copy!(out, R, A2, R) end - imfilter!(r, out, A1, kernel[1], border) + imfilter!(r, out, A1, kernel[1], border, indsstep) end -function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, inds::Indices) +function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, indsstep::Indices) kern = kernel[1] - iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), border, inds) + iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), border, indsstep) kernN = samedims(A2, kern) - imfilter!(r, A2, A1, kernN, border, inds) # store result in A2 + imfilter!(r, A2, A1, kernN, border, indsstep) # store result in A2 kernelt = tail(kernel) - newinds = next_shrink(inds, kernelt) + newinds = next_shrink(indsstep, kernelt) _imfilter!(r, out, A2, A1, tail(kernel), border, newinds) # swap the buffers end -### When everything is FIR, we can use a tiled algorithm -function imfilter!(r::AbstractCPU{FIR}, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad) - kern = kernel[1] - iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border) - tmp = tile_allocate(filter_type(A, kernel), kernel) - _imfilter_tiled!(r, out, A, kernel, border, tmp) - out -end - # Single-threaded, pair of kernels (with only one temporary buffer required) -function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}) +function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}, indsout) k1, k2 = kernel tile = tiles[1] indsk1, indstile = indices(k1), indices(tile) sz = map(length, indstile) chunksz = map(length, shrink(indstile, indsk1)) - inds = expand(indices(out), k2) + inds = expand(indsout, k2) for tileinds in TileIterator(inds, chunksz) tileb = TileBuffer(tile, tileinds) imfilter!(r, tileb, A, samedims(tileb, k1), border, tileinds) @@ -288,13 +330,13 @@ function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any, end # Multithreaded, pair of kernels -function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}) +function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}, indsout) k1, k2 = kernel tile = tiles[1] indsk1, indstile = indices(k1), indices(tile) sz = map(length, indstile) chunksz = map(length, shrink(indstile, indsk1)) - tileinds_all = collect(TileIterator(expand(indices(out), k2), chunksz)) + tileinds_all = collect(TileIterator(expand(indsout, k2), chunksz)) _imfilter_tiled_threads!(CPU1(r), out, A, samedims(out, k1), samedims(out, k2), border, tileinds_all, tiles) end # This must be in a separate function due to #15276 @@ -310,13 +352,13 @@ end end # Single-threaded, multiple kernels (requires two tile buffers, swapping on each iteration) -function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}) +function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}, indsout) k1, kt = kernel[1], tail(kernel) tilepair = tiles[1] indsk1, indstile = indices(k1), indices(tilepair[1]) sz = map(length, indstile) chunksz = map(length, shrink(indstile, indsk1)) - inds = expand(indices(out), kt) + inds = expand(indsout, kt) for tileinds in TileIterator(inds, chunksz) tileb1 = TileBuffer(tilepair[1], tileinds) imfilter!(r, tileb1, A, samedims(tileb1, k1), border, tileinds) @@ -325,13 +367,13 @@ function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any, end # Multithreaded, multiple kernels -function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}) +function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}, indsout) k1, kt = kernel[1], tail(kernel) tilepair = tiles[1] indsk1, indstile = indices(k1), indices(tilepair[1]) sz = map(length, indstile) chunksz = map(length, shrink(indstile, indsk1)) - tileinds_all = collect(TileIterator(expand(indices(out), kt), chunksz)) + tileinds_all = collect(TileIterator(expand(indsout, kt), chunksz)) _imfilter_tiled_threads!(CPU1(r), out, A, samedims(out, k1), kt, border, tileinds_all, tiles) end # This must be in a separate function due to #15276 @@ -467,6 +509,20 @@ function _imfilter_inbounds{TT}(r::AbstractResource, ::Type{TT}, out, A::Abstrac out end +function _imfilter_iter!(r::AbstractResource, out, padded, kernel::AbstractArray, iter) + p = first(padded) * first(kernel) + TT = typeof(p+p) + Rk = CartesianRange(indices(kernel)) + for I in iter + tmp = zero(TT) + for J in Rk + tmp += padded[I+J]*kernel[J] + end + out[I] = tmp + end + out +end + ### FFT filtering """ @@ -515,7 +571,7 @@ function _imfilter_fft!{S,T,N}(r::AbstractCPU{FFT}, A::AbstractArray{T,N}, kernel::Tuple{AbstractArray,Vararg{AbstractArray}}, border::NoPad) - kern = prod_kernel(Val{N}, kernel...) + kern = samedims(A, kernelconv(kernel...)) krn = FFTView(zeros(eltype(kern), map(length, indices(A)))) for I in CartesianRange(indices(kern)) krn[I] = kern[I] @@ -560,7 +616,7 @@ end # Note this is safe for inplace use, i.e., out === img -function imfilter!{S,T,N}(r::AbstractResource, +function imfilter!{S,T,N}(r::AbstractResource{IIR}, out::AbstractArray{S,N}, img::AbstractArray{T,N}, kernel::Tuple{TriggsSdika, Vararg{TriggsSdika}}, @@ -926,6 +982,20 @@ iscopy(kernel::Laplacian) = false iscopy(kernel::TriggsSdika) = all(x->x==0, kernel.a) && all(x->x==0, kernel.b) && kernel.scale == 1 iscopy(kernel::ReshapedVector) = iscopy(kernel.data) +kernelconv(kernel) = kernel +function kernelconv(k1, k2, kernels...) + out = similar(Array{filter_type(eltype(k1), k2)}, calculate_padding((k1, k2))) + fill!(out, zero(eltype(out))) + k1N, k2N = samedims(out, k1), samedims(out, k2) + R1, R2 = CartesianRange(indices(k1N)), CartesianRange(indices(k2N)) + for I1 in R1 + for I2 in R2 + out[I1+I2] += k1N[I1]*k2N[I2] + end + end + kernelconv(out, kernels...) +end + ## Tiling utilities function tile_allocate{T}(::Type{T}, kernel::Tuple{Any,Any})