JuliaImages / ImageFiltering.jl

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

kernelfactors for Gaussian kernel #77

Closed ghost closed 5 years ago

ghost commented 5 years ago

Currently, Kernel.gaussian(1)returns 2-dimensional array:

julia> kernel = Kernel.gaussian(1)
OffsetArray(::Array{Float64,2}, -2:2, -2:2) with eltype Float64 with indices -2:2×-2:2:
 0.00296902  0.0133062  0.0219382  0.0133062  0.00296902
 0.0133062   0.0596343  0.0983203  0.0596343  0.0133062
 0.0219382   0.0983203  0.162103   0.0983203  0.0219382
 0.0133062   0.0596343  0.0983203  0.0596343  0.0133062
 0.00296902  0.0133062  0.0219382  0.0133062  0.00296902

But, as is explained in Factored kernels, one can use kernelfactors for Gaussian kernel instead of 2d-array. So, why don't you use kernelfactors by default?

For speed comparison, 2d-array filtering took about 24ms on average:

julia> using Images, ImageFiltering
julia> using TestImages
julia> img = testimage("cameraman")  #512x512 Gray{N0f}

julia> kernel = Kernel.gaussian(3)
OffsetArray(::Array{Float64,2}, -6:6, -6:6) with eltype Float64 with indices -6:6×-6:6:
 0.000343881  0.000633593  0.00104462  …  0.00104462  0.000633593  0.000343881
 0.000633593  0.00116738   0.00192468     0.00192468  0.00116738   0.000633593
 ...
julia> imfilter(img, kernel); # precompile

julia> using BenchmarkTools
julia> @benchmark imfilter(img, kernel)
BenchmarkTools.Trial:
  memory estimate:  4.17 MiB
  allocs estimate:  518
  --------------
  minimum time:     23.101 ms (0.00% GC)
  median time:      23.833 ms (0.00% GC)
  mean time:        24.960 ms (3.80% GC)
  maximum time:     97.076 ms (75.65% GC)
  --------------
  samples:          201
  evals/sample:     1

The kernelfactors version took about just 5 or 6ms on average.

julia> k = KernelFactors.gaussian(3)
OffsetArray(::Array{Float64,1}, -6:6) with eltype Float64 with indices -6:6:
 0.01854402167877106
 0.03416694194322774
 ...

julia> kernelf = kernelfactors((k, k))
(ImageFiltering.KernelFactors.ReshapedOneD{Float64,2,0,OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}}([0.018544, 0.0341669, 0.0563318, 0.0831085, 0.109719, 0.129618, 0.137023, 0.129618, 0.109719, 0.0831085, 0.0563318, 0.0341669, 0.018544]), ImageFiltering.KernelFactors.ReshapedOneD{Float64,2,1,OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}}([0.018544, 0.0341669, 0.0563318, 0.0831085, 0.109719, 0.129618, 0.137023, 0.129618, 0.109719, 0.0831085, 0.0563318, 0.0341669, 0.018544]))

julia> broadcast(*, kernelf...) == kernel
true

julia> imgk = imfilter(img, kernel);
julia> imgf = imfilter(img, kernelf);
julia> imgk ≈ imgf  # verify the results
true

julia> @benchmark imfilter(img, kernelf)
BenchmarkTools.Trial:
  memory estimate:  4.16 MiB
  allocs estimate:  494
  --------------
  minimum time:     3.750 ms (0.00% GC)
  median time:      5.560 ms (0.00% GC)
  mean time:        6.353 ms (10.82% GC)
  maximum time:     294.362 ms (98.35% GC)
  --------------
  samples:          823
  evals/sample:     1

Following is my environment.

julia> versioninfo()
Julia Version 1.0.0
Commit 5d4eaca0c9 (2018-08-08 20:58 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i5-6360U CPU @ 2.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 2
timholy commented 5 years ago

The intention is that by calling Kernel.gaussian you're explicitly asking for the non-factored version---you might want to use the matrix that it returns (e.g., to plot it), rather than "just" being something you're about to pass to imfilter. If you know you want the factored version, you call KernelFactors.gaussian. In other words, you get what you ask for.

Consequently I don't see a problem, but take that with a grain of salt: since I designed this, naturally I think it makes sense, but for all I know I might be the only one. It's also possible that this is under-documented.

ghost commented 5 years ago

Thanks for your comment. Now I understand the difference between Kernel and KernelFactors modules. As you pointed out, my intension was this one:

kernelf = KernelFactors.gaussian((3, 3))
imgf = imfilter(img,  kernelf)

I agree to your design (separating normal kernel and factored ones), so I will close this PR not to violate it.

Thanks.