new gabor filter: lazy Gabor array

This commit introduces an lazy array interface for the Gabor filter,
benchmark shows that it's 2-3x faster than the previous version.

The documentation is heavly rewrote to explain this filter and all
its parameters.

Author: johnnychen94

docs/Project.toml

@@ -1,6 +1,8 @@
 [deps]
 DemoCards = "311a05b2-6137-4a5a-b473-18580a3d38b5"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageDistances = "51556ac3-7006-55f5-8cb3-34580c88182d"

docs/demos/filters/gabor_filter.jl

@@ -0,0 +1,177 @@
+# ---
+# title: Gabor filter
+# id: demo_gabor_filter
+# cover: assets/gabor.png
+# author: Johnny Chen
+# date: 2021-11-01
+# ---
+
+# This example shows how one can apply spatial space kernesl [`Gabor`](@ref Kernel.Gabor)
+# using fourier transformation and convolution theorem to extract image features.
+
+using ImageCore, ImageShow, ImageFiltering # or you could just `using Images`
+using FFTW
+using TestImages
+
+# ## Definition
+#
+# Mathematically, Gabor kernel is defined in spatial space:
+#
+# ```math
+# g(x, y) = \exp(-\frac{x'^2 + \gamma^2y'^2}{2\sigma^2})\exp(i(2\pi\frac{x'}{\lambda} + \psi))
+# ```
+# where ``i`` is imaginary unit `Complex(0, 1)`, and
+# ```math
+# x' = x\cos\theta + x\sin\theta \\
+# y' = -x\sin\theta + y\cos\theta
+# ```
+#
+
+# First of all, Gabor kernel is a complex-valued matrix:
+
+kern = Kernel.Gabor((10, 10), 2, 0.1)
+
+# !!! tip "Lazy array"
+#     The `Gabor` type is a lazy array, which means when you build the Gabor kernel, you
+#     actually don't need to allocate any memories.
+#
+#     ```julia
+#     using BenchmarkTools
+#     kern = @btime Kernel.Gabor((64, 64), 5, 0); # 36.481 ns (0 allocations: 0 bytes)
+#     @btime collect($kern); # 75.278 μs (2 allocations: 64.05 KiB)
+#     ```
+
+# To explain the parameters of Gabor filter, let's introduce some small helpers function to
+# display complex-valued kernels.
+show_phase(kern) = @. Gray(log(abs(imag(kern)) + 1))
+show_mag(kern) = @. Gray(log(abs(real(kern)) + 1))
+show_abs(kern) = @. Gray(log(abs(kern) + 1))
+nothing #hide
+
+# ## Keywords
+#
+# ### `wavelength` (λ)
+# λ specifies the wavelength of the sinusoidal factor.
+
+bandwidth, orientation, phase_offset, aspect_ratio = 1, 0, 0, 0.5
+f(wavelength) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset))
+mosaic(f.((5, 10, 15)), nrow=1)
+
+# ### `orientation` (θ)
+# θ specifies the orientation of the normal to the parallel stripes of a Gabor function.
+
+wavelength, bandwidth, phase_offset, aspect_ratio = 10, 1, 0, 0.5
+f(orientation) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset))
+mosaic(f.((0, π/4, π/2)), nrow=1)
+
+# ### `phase_offset` (ψ)
+
+wavelength, bandwidth, orientation, aspect_ratio = 10, 1, 0, 0.5
+f(phase_offset) = show_phase(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset))
+mosaic(f.((-π/2, 0, π/2, π)), nrow=1)
+
+# ### `aspect_ratio` (γ)
+# γ specifies the ellipticity of the support of the Gabor function.
+
+wavelength, bandwidth, orientation, phase_offset = 10, 1, 0, 0
+f(aspect_ratio) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset))
+mosaic(f.((0.5, 1, 2)), nrow=1)
+
+# ### `bandwidth` (b)
+# The half-response spatial frequency bandwidth (b) of a Gabor filter is related to the
+# ratio σ / λ, where σ and λ are the standard deviation of the Gaussian factor of the Gabor
+# function and the preferred wavelength, respectively, as follows:
+#
+# ```math
+# b = \log_2\frac{\frac{\sigma}{\lambda}\pi + \sqrt{\frac{\ln 2}{2}}}{\frac{\sigma}{\lambda}\pi - \sqrt{\frac{\ln 2}{2}}}
+# ```
+
+wavelength, orientation, phase_offset, aspect_ratio = 10, 0, 0, 0.5
+f(bandwidth) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset))
+mosaic(f.((0.5, 1, 2)), nrow=1)
+
+# ## Examples
+#
+# There are two options to apply a spatial space kernel: 1) via `imfilter`, and 2) via the
+# convolution theorem.
+
+# ### `imfilter`
+#
+# [`imfilter`](@ref) does not require the kernel size to be the same as the image size.
+# Usually kernel size is at least 5 times larger than the wavelength.
+
+img = TestImages.shepp_logan(127)
+kern = Kernel.Gabor((19, 19), 3, 0)
+out = imfilter(img, real.(kern))
+mosaic(img, show_abs(kern), show_mag(out); nrow=1)
+
+# ### convolution theorem
+
+# The [convolution theorem](https://en.wikipedia.org/wiki/Convolution_theorem) tells us
+# that `fft(conv(X, K))` is equivalent to `fft(X) .* fft(K)`. Because Gabor kernel is
+# defined with regard to its center point (0, 0), we need to do `ifftshift` first so that
+# the frequency centers of both `fft(X)` and `fft(kern)` align well.
+
+kern = Kernel.Gabor(size(img), 3, 0)
+out = ifft(fft(channelview(img)) .* ifftshift(fft(kern)))
+mosaic(img, show_abs(kern), show_mag(out); nrow=1)
+
+# As you may have notice, using convolution theorem generates different results. This is
+# simply because the kernel size are different. If we create a smaller kernel, we then need
+# to apply [`freqkernel`](@ref) first so that we can do element-wise multiplication.
+
+## freqkernel = zero padding + fftshift + fft
+kern = Kernel.Gabor((19, 19), 3, 0)
+kern_freq = freqkernel(real.(kern), size(img))
+out = ifft(fft(channelview(img)) .* kern_freq)
+mosaic(img, show_abs(kern), show_mag(out); nrow=1)
+
+# !!! note "Performance on different kernel size"
+#     When the kernel size is small, `imfilter` works more efficient than fft-based
+#     convolution. This benchmark isn't backed by CI so the result might be outdated, but
+#     you get the idea.
+#
+#     ```julia
+#     using BenchmarkTools
+#     img = TestImages.shepp_logan(127);
+#
+#     kern = Kernel.Gabor((19, 19), 3, 0);
+#     fft_conv(img, kern) = ifft(fft(channelview(img)) .* freqkernel(real.(kern), size(img)))
+#     @btime imfilter($img, real.($kern)); # 236.813 μs (118 allocations: 418.91 KiB)
+#     @btime fft_conv($img, $kern) # 1.777 ms (127 allocations: 1.61 MiB)
+#
+#     kern = Kernel.Gabor(size(img), 3, 0)
+#     fft_conv(img, kern) =  ifft(fft(channelview(img)) .* ifftshift(fft(kern)))
+#     @btime imfilter($img, real.($kern)); # 5.318 ms (163 allocations: 5.28 MiB)
+#     @btime fft_conv($img, $kern); # 2.218 ms (120 allocations: 1.73 MiB)
+#     ```
+#
+
+## Filter bank
+
+# A filter bank is just a list of filter kernels, applying the filter bank generates
+# multiple outputs:
+
+filters = [Kernel.Gabor(size(img), 3, θ) for θ in -π/2:π/4:π/2];
+X_freq = fft(channelview(img))
+out = map(filters) do kern
+    ifft(X_freq .* ifftshift(fft(kern)))
+end
+mosaic(
+    map(show_abs, filters)...,
+    map(show_abs, out)...;
+    nrow=2, rowmajor=true
+)
+
+## save covers #src
+using FileIO #src
+mkpath("assets")  #src
+filters = [Kernel.Gabor((32, 32), 5, θ) for θ in range(-π/2, stop=π/2, length=9)] #src
+save("assets/gabor.png", mosaic(map(show_abs, filters); nrow=3, npad=2, fillvalue=Gray(1)); fps=2) #src
+
+# # References
+#
+# - [1] [Wikipedia: Gabor filter](https://en.wikipedia.org/wiki/Gabor_filter)
+# - [2] [Wikipedia: Gabor transformation](https://en.wikipedia.org/wiki/Gabor_transform)
+# - [3] [Wikipedia: Gabor atom](https://en.wikipedia.org/wiki/Gabor_atom)
+# - [4] [Gabor filter for image processing and computer vision](http://matlabserver.cs.rug.nl/edgedetectionweb/web/edgedetection_params.html)

docs/src/function_reference.md

@@ -23,7 +23,7 @@ Kernel.gaussian
 Kernel.DoG
 Kernel.LoG
 Kernel.Laplacian
-Kernel.gabor
+Kernel.Gabor
 Kernel.moffat
 ```
 

src/gaborkernels.jl

@@ -0,0 +1,182 @@
+module GaborKernels
+
+export Gabor
+
+"""
+    Gabor(size_or_axes, wavelength, orientation; kwargs...)
+    Gabor(size_or_axes; wavelength, orientation, kwargs...)
+
+Generate the 2-D Gabor kernel in the spatial space.
+
+# Arguments
+
+## `kernel_size::Dims{2}` or `kernel_axes::NTuple{2,<:AbstractUnitRange}`
+
+Specifies either the size or axes of the output kernel. The axes at each dimension will be
+`-r:r` if the size is odd.
+
+## `wavelength::Real`(λ>=2)
+
+This is the wavelength of the sinusoidal factor of the Gabor filter kernel and herewith the
+preferred wavelength of this filter. Its value is specified in pixels.
+
+The value `λ=2` should not be used in combination with phase offset `ψ=-π/2` or `ψ=π/2`
+because in these cases the Gabor function is sampled in its zero crossings.
+
+In order to prevent the occurence of undesired effects at the image borders, the wavelength
+value should be smaller than one fifth of the input image size.
+
+## `orientation::Real`(θ∈[0, 2pi]):
+
+This parameter specifies the orientation of the normal to the parallel stripes of a Gabor
+function [3].
+
+# Keywords
+
+## `bandwidth=1` (b>0)
+
+The half-response spatial frequency bandwidth b (in octaves) of a Gabor filter is related to
+the ratio σ / λ, where σ and λ are the standard deviation of the Gaussian factor of the
+Gabor function and the preferred wavelength, respectively, as follows:
+
+```math
+b = \\log_2\\frac{\\frac{\\sigma}{\\lambda}\\pi + \\sqrt{\\frac{\\ln 2}{2}}}{\\frac{\\sigma}{\\lambda}\\pi - \\sqrt{\\frac{\\ln 2}{2}}}
+```
+
+## `aspect_ratio=0.5`(γ>0)
+
+This parameter, called more precisely the spatial aspect ratio, specifies the ellipticity of
+the support of the Gabor function [3]. For `γ = 1`, the support is circular. For γ < 1 the
+support is elongated in orientation of the parallel stripes of the function.
+
+## `phase_offset=0` (ψ∈[-π, π])
+
+The values `0` and `π` correspond to center-symmetric 'center-on' and 'center-off'
+functions, respectively, while -π/2 and π/2 correspond to anti-symmetric functions. All
+other cases correspond to asymmetric functions.
+
+# Examples
+
+There are two ways to use gabor filters: 1) via [`imfilter`](@ref), 2) via `fft` and
+convolution theorem. Usually `imfilter` requires a small kernel to work efficiently, while
+using `fft` can be more efficient when the kernel has the same size of the image.
+
+## imfilter
+
+```jldoctest gabor_example
+julia> using ImageFiltering, FFTW, TestImages, ImageCore
+
+julia> img = TestImages.shepp_logan(256);
+
+julia> kern = Kernel.Gabor((19, 19), 3, 0); # usually a small kernel size
+
+julia> g_img = imfilter(img, real.(kern));
+
+```
+
+## convolution theorem
+
+The convolution theorem is stated as `fft(conv(A, K)) == fft(A) .* fft(K)`:
+
+```jldoctest gabor_example
+julia> kern = Kernel.Gabor(size(img), 3, 0);
+
+julia> g_img = ifft(fft(channelview(img)) .* ifftshift(fft(kern))); # apply convolution theorem
+
+julia> @. Gray(abs(g_img));
+
+```
+
+See the [gabor filter demo](@ref demo_gabor_filter) for more examples with images.
+
+# Extended help
+
+Mathematically, gabor filter is defined in spatial space:
+
+```math
+g(x, y) = \\exp(-\\frac{x'^2 + \\gamma^2y'^2}{2\\sigma^2})\\exp(i(2\\pi\\frac{x'}{\\lambda} + \\psi))
+```
+
+where ``x' = x\\cos\\theta + y\\sin\\theta`` and ``y' = -x\\sin\\theta + y\\cos\\theta``.
+
+# References
+
+- [1] [Wikipedia: Gabor filter](https://en.wikipedia.org/wiki/Gabor_filter)
+- [2] [Wikipedia: Gabor transformation](https://en.wikipedia.org/wiki/Gabor_transform)
+- [3] [Wikipedia: Gabor atom](https://en.wikipedia.org/wiki/Gabor_atom)
+- [4] [Gabor filter for image processing and computer vision](http://matlabserver.cs.rug.nl/edgedetectionweb/web/edgedetection_params.html)
+"""
+struct Gabor{T<:Complex, TP<:Real, R<:AbstractUnitRange} <: AbstractMatrix{T}
+    ax::Tuple{R,R}
+    λ::TP
+    θ::TP
+    ψ::TP
+    σ::TP
+    γ::TP
+
+    # cache values
+    sc::Tuple{TP, TP} # sincos(θ)
+    λ_scaled::TP # 2π/λ
+    function Gabor{T,TP,R}(ax::Tuple{R,R}, λ::TP, θ::TP, ψ::TP, σ::TP, γ::TP) where {T,TP,R}
+        λ > 0 || throw(ArgumentError("`λ` should be positive: $λ"))
+        new{T,TP,R}(ax, λ, θ, ψ, σ, γ, sincos(θ), 2π/λ)
+    end
+end
+function Gabor(size_or_axes::Tuple; wavelength, orientation, kwargs...)
+    Gabor(size_or_axes, wavelength, orientation; kwargs...)
+end
+function Gabor(
+        size_or_axes::Tuple,
+        wavelength::Real,
+        orientation::Real;
+        bandwidth::Real=1.0f0,
+        phase_offset::Real=0.0f0,
+        aspect_ratio::Real=0.5f0,
+    )
+    bandwidth > 0 || throw(ArgumentError("`bandwidth` should be positive: $bandwidth"))
+    aspect_ratio > 0 || throw(ArgumentError("`aspect_ratio` should be positive: $aspect_ratio"))
+    wavelength >= 2 || @warn "`wavelength` should be equal to or greater than 2" wavelength
+    # we still follow the math symbols in the implementation
+    λ, γ, ψ = wavelength, aspect_ratio, phase_offset
+    # for orientation: Julia follow column-major order, thus here we compensate the orientation by π/2
+    θ = convert(float(typeof(orientation)), mod2pi(orientation+π/2))
+    # follows reference [4]
+    σ = convert(float(typeof(λ)), (λ/π)*sqrt(0.5log(2)) * (2^bandwidth+1)/(2^bandwidth-1))
+
+    params = float.(promote(λ, θ, ψ, σ, γ))
+    T = typeof(params[1])
+
+    ax = _to_axes(size_or_axes)
+    all(map(length, ax) .> 0) || throw(ArgumentError("Kernel size should be positive: $size_or_axes"))
+    if any(r->5λ > length(r), ax)
+        expected_size = @. 5λ * length(ax)
+        @warn "for wavelength `λ=$λ` the expected kernel size is expected to be larger than $expected_size."
+    end
+
+    Gabor{Complex{T}, T, typeof(first(ax))}(ax, params...)
+end
+
+@inline Base.size(kern::Gabor) = map(length, kern.ax)
+@inline Base.axes(kern::Gabor) = kern.ax
+@inline function Base.getindex(kern::Gabor, x::Int, y::Int)
+    ψ, σ, γ = kern.ψ, kern.σ, kern.γ
+    s, c = kern.sc # sincos(θ)
+    λ_scaled = kern.λ_scaled # 2π/λ
+
+    xr = x*c + y*s
+    yr = -x*s + y*c
+    yr = γ * yr
+    return exp(-(xr^2 + yr^2)/(2σ^2)) * cis(xr*λ_scaled + ψ)
+end
+
+# Utils
+
+function _to_axes(sz::Dims)
+    map(sz) do n
+        r=n÷2
+        isodd(n) ? UnitRange(-r, n-r-1) : UnitRange(-r+1, n-r)
+    end
+end
+_to_axes(ax::NTuple{N, T}) where {N, T<:AbstractUnitRange} = ax
+
+end

src/kernel.jl

@@ -11,7 +11,7 @@ dimensionality. The following kernels are supported:
   - `DoG` (Difference-of-Gaussian)
   - `LoG` (Laplacian-of-Gaussian)
   - `Laplacian`
-  - `gabor`
+  - `Gabor`
   - `moffat`
 
 See also: [`KernelFactors`](@ref).
@@ -23,6 +23,9 @@ using ..ImageFiltering
 using ..ImageFiltering.KernelFactors
 import ..ImageFiltering: _reshape, IdentityUnitRange
 
+include("gaborkernels.jl")
+using .GaborKernels
+
 # We would like to do `using ..ImageFiltering.imgradients` so that that
 # Documenter.jl (the documentation system) can parse a reference such as `See
 # also: [`ImageFiltering.imgradients`](@ref)`. However, imgradients is not yet