search

JuliaImages / ImageFiltering.jl

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

Feature request: Specify a tuple of padding styles #94

Closed aramirezreyes closed 5 years ago

aramirezreyes commented 5 years ago

Hi:

It would be nice to have the possibility to specify a tuple of padding styles in the call to imfilter: I am thinking, for example, where the array to be filtered has dimensions with different physical meanings. The particular use case I have in mind is one where the first two dimensions are periodic boundaries, so "circular" works fine, but for the 3rd dimension (temporal) is more suited for "replicate", or another style.

I'm thinking something like:

using ImageFiltering
borders = ("circular","circular","replicate")
imfilter(f,w,borders)

This can be worked around by calling twice imfilter with different kernels, each with it's own border style, so it is by no means an urgent matter.

Additionally, with some orientation, I think I would be able to dig into this.

jw3126 commented 5 years ago

Are you still interested in doing this? I think it is quite a bit of effort. Also it is nor clear, what to do e.g. at the corners of an image with a ("circular", "replicate") padding? Do you use "circular" or "replicate"? If you want to implement this, the way to go is probably define struct ProductPad{...} <: AbstractBorder.

timholy commented 5 years ago

For separable kernels you'd just filter one dimension at a time. Do you have non-separable kernels where this is desired?

aramirezreyes commented 5 years ago

The reason I ask for this is because currently imfilter! allocates an amount of memory equivalent to the filtered array for every call (related to #76) , so for example:

julia> using ImageFiltering, BenchmarkTools

julia> U = buf = rand(100,100,100,100); #Array of size 762.94 MiB each
julia> smooth_x = smooth_time=10;

julia> function kernel4d(window_h)
           kernel_h = ones(window_h)/(window_h)
           return kernelfactors(( centered(kernel_h), centered(kernel_h), centered([1.0]), centered([1.0]) ))
 end

julia> @btime imfilter!($buf,$U,kernel4d(smooth_x),"circular");
  #2.752 s (2390047 allocations: 1006.57 MiB)

So this call allocated a lot, even when calling the in-place function with a preallocated buffer.

If I call it once again using a different kernel for the other direction (with different padding option), It will do the allocation once more.

When the array size increases to GiB, this calls become a huge bottleneck, altough I don't know if this because of the allocations exclusively, it is to be expected because of the nature of the operation, or it is an example of deficient usage of the functions.

timholy commented 5 years ago

Small comment, your example doesn't run; kernel4d only takes a single argument, but you call it with 2. But a single-argument call works fine.

I built Julia with https://github.com/JuliaLang/julia/pull/31915 and ran the following:

julia> Profile.init(alloc_rate=1.0, delay=1.0)

julia> @profile imfilter!(buf,U,kern,"circular");

julia> Profile.clear()

julia> @profile imfilter!(buf,U,kern,"circular");

julia> Profile.print()
30634 ./task.jl:268; (::getfield(Revise, Symbol("##74#76")){REPL.REPLBackend})()
 30634 /home/tim/.julia/dev/Revise/src/Revise.jl:943; run_backend(::REPL.REPLBackend)
  30634 /home/tim/src/julia-master/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:86; eval_user_input(::Any, ::REPL.REPLBackend)
   30634 ./boot.jl:330; eval(::Module, ::Any)
    30634 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:597; imfilter!(::Array{Float64,4}, ::Array{Float64,4}, ::Tuple{ImageFiltering.KernelFactors.ReshapedOneD{Flo...
     1     /home/tim/.julia/dev/ImageFiltering/src/border.jl:139; borderinstance(::String)
      1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:167; Type
       1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:202; Type
        1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:126; Type
     30633 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:606; imfilter!
      5     /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:1537; filter_algorithm(::Array{Float64,4}, ::Array{Float64,4}, ::Tuple{ImageFiltering.KernelFactors.Reshape...
       5 /home/tim/.julia/dev/ImageFiltering/src/border.jl:945; calculate_padding
        4 /home/tim/.julia/dev/ImageFiltering/src/border.jl:969; accumulate_padding(::NTuple{4,UnitRange{Int64}}, ::ImageFiltering.KernelFactors.ReshapedOneD{Float64,4,...
      30628 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:612; imfilter!(::Array{Float64,4}, ::Array{Float64,4}, ::Tuple{ImageFiltering.KernelFactors.ReshapedOneD{F...
       6     /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:701; imfilter!
        6 /home/tim/.julia/dev/ImageFiltering/src/border.jl:226; Pad
         1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:225; (::Pad{0})(::Tuple{ImageFiltering.KernelFactors.ReshapedOneD{Float64,4,0,OffsetArrays.OffsetArray{Floa...
          1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:207; Type
           1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:126; Type
            1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:126; Type
         5 /home/tim/.julia/dev/ImageFiltering/src/border.jl:945; (::Pad{0})(::Tuple{ImageFiltering.KernelFactors.ReshapedOneD{Float64,4,0,OffsetArrays.OffsetArray{Floa...
          4 /home/tim/.julia/dev/ImageFiltering/src/border.jl:969; accumulate_padding(::NTuple{4,UnitRange{Int64}}, ::ImageFiltering.KernelFactors.ReshapedOneD{Float64,4...
       30622 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:702; imfilter!
        14    /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:711; imfilter!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array{Float64,4}, :...
         1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:661; padarray
          1 /home/tim/.julia/dev/ImageFiltering/src/borderarray.jl:67; Type
         4 /home/tim/.julia/dev/ImageFiltering/src/border.jl:662; padarray
          4 ./abstractarray.jl:626; similar
           4 /home/tim/.julia/dev/ImageFiltering/src/borderarray.jl:77; similar
            3 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:100; similar(::Array{Float64,4}, ::Type{Float64}, ::NTuple{4,UnitRange{Int64}})
             3 ./array.jl:318; similar
              3 ./boot.jl:416; Type
            1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:101; similar(::Array{Float64,4}, ::Type{Float64}, ::NTuple{4,UnitRange{Int64}})
             1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
              1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
         9 /home/tim/.julia/dev/ImageFiltering/src/border.jl:663; padarray
          9 /home/tim/.julia/dev/ImageFiltering/src/borderarray.jl:161; copy!(::OffsetArrays.OffsetArray{Float64,4,Array{Float64,4}}, ::BorderArray{Float64,4,Array{Float64,4}...
           8 /home/tim/.julia/dev/ImageFiltering/src/borderarray.jl:174; _copy!(::OffsetArrays.OffsetArray{Float64,4,Array{Float64,4}}, ::BorderArray{Float64,4,Array{Float64,...
            6 /home/tim/.julia/dev/ImageFiltering/src/border.jl:247; padindices
             6 /home/tim/.julia/dev/ImageFiltering/src/border.jl:254; _padindices
              6 /home/tim/.julia/dev/ImageFiltering/src/border.jl:254; _padindices
               4 /home/tim/.julia/dev/ImageFiltering/src/border.jl:254; _padindices
                2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:254; _padindices
                 2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:909; padindex(::Pad{4}, ::Int64, ::Base.OneTo{Int64}, ::Int64)
                  2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:974; modrange
                   2 ./array.jl:548; map
                    2 ./array.jl:624; _collect(::OffsetArrays.OffsetArray{Int64,1,UnitRange{Int64}}, ::Base.Generator{OffsetArrays.Off...
                     2 ./array.jl:519; _similar_for
                      1 ./boot.jl:404; similar
                      1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:101; similar
                       1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
                        1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
                2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:909; padindex(::Pad{4}, ::Int64, ::Base.OneTo{Int64}, ::Int64)
                 2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:974; modrange
                  2 ./array.jl:548; map
                   2 ./array.jl:624; _collect(::OffsetArrays.OffsetArray{Int64,1,UnitRange{Int64}}, ::Base.Generator{OffsetArrays.Off...
                    2 ./array.jl:519; _similar_for
                     1 ./boot.jl:404; similar
                     1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:101; similar
                      1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
                       1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
               2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:909; padindex(::Pad{4}, ::Int64, ::Base.OneTo{Int64}, ::Int64)
                2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:974; modrange
                 2 ./array.jl:548; map
                  2 ./array.jl:624; _collect(::OffsetArrays.OffsetArray{Int64,1,UnitRange{Int64}}, ::Base.Generator{OffsetArrays.Offs...
                   2 ./array.jl:519; _similar_for
                    1 ./boot.jl:404; similar
                    1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:101; similar
                     1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
                      1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
            2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:254; padindices
             2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:909; padindex(::Pad{4}, ::Int64, ::Base.OneTo{Int64}, ::Int64)
              2 /home/tim/.julia/dev/ImageFiltering/src/border.jl:974; modrange
               2 ./array.jl:548; map
                2 ./array.jl:624; _collect(::OffsetArrays.OffsetArray{Int64,1,UnitRange{Int64}}, ::Base.Generator{OffsetArrays.Offse...
                 2 ./array.jl:519; _similar_for
                  1 ./boot.jl:404; similar
                  1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:101; similar
                   1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
                    1 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
           1 /home/tim/.julia/dev/ImageFiltering/src/borderarray.jl:175; _copy!(::OffsetArrays.OffsetArray{Float64,4,Array{Float64,4}}, ::BorderArray{Float64,4,Array{Float64,...
            1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:685; copydata!
             1 ./tuple.jl:142; map
        30608 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:715; imfilter!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array{Float64,4}, :...
         1     /home/tim/.julia/dev/ImageFiltering/src/border.jl:8; Type
          1 /home/tim/.julia/dev/ImageFiltering/src/border.jl:8; Type
         30607 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:770; imfilter!
          6     /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:772; imfilter!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array{Float64,4}, ...
           6 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:1643; tile_allocate
            5 ./array.jl:606; collect(::Base.Generator{UnitRange{Int64},getfield(ImageFiltering, Symbol("##74#75")){Float64,NTuple...
             5 ./generator.jl:47; iterate
              5 ./none:0; #74
               4 ./boot.jl:421; Type
                4 ./boot.jl:416; Type
            1 ./array.jl:611; collect(::Base.Generator{UnitRange{Int64},getfield(ImageFiltering, Symbol("##74#75")){Float64,NTuple...
             1 ./array.jl:598; _array_for
              1 ./abstractarray.jl:671; similar
               1 ./abstractarray.jl:672; similar
                1 ./boot.jl:413; Type
                 1 ./boot.jl:404; Type
          30601 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:773; imfilter!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array{Float64,4}, ...
           1     /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:870; _imfilter_tiled!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array{Float...
            1 ./essentials.jl:201; tail
           7200  /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:877; _imfilter_tiled!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array{Float...
            5400 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:170; Type
             3600 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:193; _tileview
              3600 ./pointer.jl:84; unsafe_wrap
               3600 ./pointer.jl:84; #unsafe_wrap#62
             1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:56; OffsetArrays.OffsetArray(::Array{Float64,4}, ::NTuple{4,UnitRange{Int64}})
              1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
               1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
            1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:171; Type
             1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:157; Type
              1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:157; Type
           23400 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:879; _imfilter_tiled!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array{Float...
            1800  /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:911; _imfilter_tiled_swap!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array...
             1800 ./essentials.jl:201; tail
            7200  /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:914; _imfilter_tiled_swap!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array...
             5400 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:170; Type
              3600 ./pointer.jl:84; _tileview
              1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:56; OffsetArrays.OffsetArray(::Array{Float64,4}, ::NTuple{4,Base.IdentityUnitRange{UnitRange{Int64}}})
               1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
                1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
             1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:171; Type
              1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:157; Type
               1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:157; Type
            12600 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:916; _imfilter_tiled_swap!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array...
             1800 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:911; _imfilter_tiled_swap!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array...
              1800 ./essentials.jl:201; tail
             7200 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:914; _imfilter_tiled_swap!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array...
              5400 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:170; Type
               3600 ./pointer.jl:84; _tileview
               1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:56; OffsetArrays.OffsetArray(::Array{Float64,4}, ::NTuple{4,Base.IdentityUnitRange{UnitRange{Int64}}})
                1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:20; Type
                 1800 /home/tim/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:15; Type
              1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:171; Type
               1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:157; Type
                1800 /home/tim/.julia/dev/TiledIteration/src/TiledIteration.jl:157; Type
             1800 /home/tim/.julia/dev/ImageFiltering/src/imfilter.jl:916; _imfilter_tiled_swap!(::ComputationalResources.CPU1{ImageFiltering.Algorithm.FIRTiled{4}}, ::Array...

julia> size(buf)
(30, 30, 30, 30)

Notice that the big allocations all come in integer multiples of 30*30=900. This, and the line numbers of the allocations, basically shows you that the allocations come from creating 2d slice wrappers (e.g., the TileBuffer creation in https://github.com/JuliaImages/ImageFiltering.jl/blob/0ae7df1f31310eb967bd140f290e664da8cde6c7/src/imfilter.jl#L830-L844). See https://stackoverflow.com/questions/47590839/unexpected-memory-allocation-when-using-array-views-julia/47607539#47607539.

If you instead profile for performance, these do not show up as hotspots, and the amount of time in the gc is not a large fraction of the total time. So I don't think there's an actual performance problem here.

aramirezreyes commented 5 years ago

Hi Tim: Thanks for the comment, cleaning to get a MWE sometimes goes wrong.

For some reason I failt to get the performance of the code being used here before (it was a Matlab code). While all the other parts run faster, the filtering par overhead in my version is slower enough to make the overall program slower (this is the most expensive part). I was mistakenly thinking this was due to GC from this case. Thanks for the thorough explanation.

timholy commented 5 years ago

Well, there may be other bottlenecks, I was just addressing your direct question. If you're only filtering along 2 dimensions, why not just have the kernel comprised of those two? 

kern = kernel4d(smooth_x)
imfilter!(buf,U,kern[1:2],"circular")

because only the first two do anything.

To make it faster:

aramirezreyes commented 5 years ago

Thanks a lot. Closing this and sorry for the noise.