WIP: introduce log gabor filter

Author: johnnychen94

docs/demos/filters/gabor_filter.jl

@@ -7,7 +7,8 @@
 # ---
 
 # 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 fourier transformation and convolution theorem to extract image features. A similar
+# example is [Log Gabor filter](@ref demo_log_gabor_filter).
 
 using ImageCore, ImageShow, ImageFiltering # or you could just `using Images`
 using FFTW

docs/demos/filters/loggabor_filter.jl

@@ -0,0 +1,106 @@
+# ---
+# title: Log Gabor filter
+# id: demo_log_gabor_filter
+# cover: assets/log_gabor.png
+# author: Johnny Chen
+# date: 2021-11-01
+# ---
+
+# This example shows how one can apply frequency space kernesl [`LogGabor`](@ref
+# Kernel.LogGabor) and [`LogGaborComplex`](@ref Kernel.LogGaborComplex) using fourier
+# transformation and convolution theorem to extract image features. A similar example
+# is the sptaial space kernel [Gabor filter](@ref demo_gabor_filter).
+
+using ImageCore, ImageShow, ImageFiltering # or you could just `using Images`
+using FFTW
+using TestImages
+
+# ## Definition
+#
+# Mathematically, log gabor filter is defined in spatial space as the composition of
+# its frequency component `r` and angular component `a`:
+#
+# ```math
+# r(\omega, \theta) = \exp(-\frac{(\log(\omega/\omega_0))^2}{2\sigma_\omega^2}) \\
+# a(\omega, \theta) = \exp(-\frac{(\theta - \theta_0)^2}{2\sigma_\theta^2})
+# ```
+#
+# `LogGaborComplex` provides a complex-valued matrix with value `Complex(r, a)`, while
+# `LogGabor` provides real-valued matrix with value `r * a`.
+
+kern_c = Kernel.LogGaborComplex((10, 10), 1/6, 0)
+kern_r = Kernel.LogGabor((10, 10), 1/6, 0)
+kern_r == @. real(kern_c) * imag(kern_c)
+
+# !!! note "Lazy array"
+#     The `LogGabor` and `LogGaborComplex` types are lazy arrays, which means when you build
+#     the Log Gabor kernel, you actually don't need to allocate any memories. The computation
+#     does not happen until you request the value.
+#
+#     ```julia
+#     using BenchmarkTools
+#     kern = @btime Kernel.LogGabor((64, 64), 1/6, 0); # 1.711 ns (0 allocations: 0 bytes)
+#     @btime collect($kern); # 146.260 μs (2 allocations: 64.05 KiB)
+#     ```
+
+
+# To explain the parameters of Log Gabor filter, let's introduce one small helper 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
+
+# From left to right are visualization of the kernel in frequency space: frequency `r`,
+# algular `a`, `sqrt(r^2 + a^2)`, `r * a`, and its spatial space kernel.
+kern = Kernel.LogGaborComplex((32, 32), 100, 0)
+mosaic(
+    show_mag(kern),
+    show_phase(kern),
+    show_abs(kern),
+    Gray.(Kernel.LogGabor(kern)),
+    show_abs(centered(ifftshift(ifft(kern)))),
+    nrow=1
+)
+
+# ## Examples
+#
+# To apply the filter, we need to use [the convolution
+# theorem](https://en.wikipedia.org/wiki/Convolution_theorem):
+# ```math
+# \mathcal{F}(x \circledast k) = \mathcal{F}(x) \odot \mathcal{F}(k)
+# ```
+# where ``\circledast`` is convolution, ``\odot`` is pointwise-multiplication, and
+# ``\mathcal{F}`` is the fourier transformation.
+#
+# Because Log Gabor kernel is defined around center point (0, 0), we have to apply
+# `ifftshift` first before we do pointwise-multiplication. Because Log Gabor kernel
+# is defined on frequency space, we don't need to call `fft(kern)`.
+
+img = TestImages.shepp_logan(127)
+kern = Kernel.LogGaborComplex(size(img), 50, π/4)
+out = ifft(centered(fft(channelview(img))) .* ifftshift(kern))
+mosaic(img, show_abs(kern), show_mag(out); nrow=1)
+
+# A filter bank is just a list of filter kernels, applying the filter bank generates
+# multiple outputs:
+
+X_freq = centered(fft(channelview(img)))
+filters = vcat(
+    [Kernel.LogGaborComplex(size(img), 50, θ) for θ in -π/2:π/4:π/2],
+    [Kernel.LogGabor(size(img), 50, θ) for θ in -π/2:π/4:π/2]
+)
+out = map(filters) do kern
+    ifft(X_freq .* ifftshift(kern))
+end
+mosaic(
+    map(show_abs, filters)...,
+    map(show_abs, out)...;
+    nrow=4, rowmajor=true
+)
+
+## save covers #src
+using FileIO #src
+mkpath("assets")  #src
+filters = [Kernel.LogGaborComplex((32, 32), 5, θ) for θ in range(-π/2, stop=π/2, length=9)] #src
+save("assets/log_gabor.png", mosaic(map(show_abs, filters); nrow=3, npad=2, fillvalue=Gray(1)); fps=2) #src

docs/src/function_reference.md

@@ -24,6 +24,8 @@ Kernel.DoG
 Kernel.LoG
 Kernel.Laplacian
 Kernel.Gabor
+Kernel.LogGabor
+Kernel.LogGaborComplex
 Kernel.moffat
 ```
 

src/gaborkernels.jl

@@ -1,6 +1,6 @@
 module GaborKernels
 
-export Gabor
+export Gabor, LogGabor, LogGaborComplex
 
 """
     Gabor(size_or_axes, wavelength, orientation; kwargs...)
@@ -147,6 +147,130 @@ end
     return exp(-(xr^2 + yr^2)/(2σ^2)) * cis(xr*λ_scaled + ψ)
 end
 
+
+"""
+    LogGaborComplex(size_or_axes, ω, θ; σω=1, σθ=1)
+    LogGaborComplex(size_or_axes; ω, θ, σω=1, σθ=1)
+
+Generates the 2-D Log Gabor kernel in spatial space by `Complex(r, a)`, where `r` and `a`
+are the frequency and angular components, respectively.
+
+More detailed documentation can be found in the `r * a` version [`LogGabor`](@ref).
+"""
+struct LogGaborComplex{T, TP,R<:AbstractUnitRange} <: AbstractMatrix{T}
+    ax::Tuple{R,R}
+    ω::TP
+    θ::TP
+    σω::TP
+    σθ::TP
+
+    # cache values
+    ω_denom::TP # 1/(2(log(σω/ω))^2)
+    θ_denom::TP # 1/(2σθ^2)
+    function LogGaborComplex{T,TP,R}(ax::Tuple{R,R}, ω::TP, θ::TP, σω::TP, σθ::TP) where {T,TP,R}
+        σω > 0 || throw(ArgumentError("`σω` should be positive: $σω"))
+        σθ > 0 || throw(ArgumentError("`σθ` should be positive: $σθ"))
+        new{T,TP,R}(ax, ω, θ, σω, σθ, 1/(2(log(σω/ω))^2), 1/(2σθ^2))
+    end
+end
+LogGaborComplex(size_or_axes::Tuple, ω, θ; kwargs...) = LogGaborComplex(size_or_axes; ω, θ, kwargs...)
+function LogGaborComplex(size_or_axes::Tuple; ω::Real, θ::Real, σω::Real=1, σθ::Real=1)
+    params = float.(promote(ω, θ, σω, σθ))
+    T = typeof(params[1])
+    ax = _to_axes(size_or_axes)
+    LogGaborComplex{Complex{T}, T, typeof(first(ax))}(ax, params...)
+end
+
+@inline Base.size(kern::LogGaborComplex) = map(length, kern.ax)
+@inline Base.axes(kern::LogGaborComplex) = kern.ax
+
+@inline function Base.getindex(kern::LogGaborComplex, x::Int, y::Int)
+    ω_denom, θ_denom = kern.ω_denom, kern.θ_denom
+    # normalize: from reference [1]
+    x, y = (x, y)./size(kern)
+    ω = sqrt(x^2 + y^2) # this is faster than hypot(x, y)
+    θ = atan(y, x)
+    r = exp((-(log(ω/kern.ω))^2)*ω_denom) # radial component
+    a = exp((-(θ-kern.θ).^2)*θ_denom) # angular component
+    return r + a*im
+end
+
+
+"""
+    LogGabor(size_or_axes, ω, θ; σω=1, σθ=1)
+    LogGabor(size_or_axes; ω, θ, σω=1, σθ=1)
+
+Generates the 2-D Log Gabor kernel in spatial space by `r * a`, where `r` and `a` are the
+frequency and angular components, respectively.
+
+See also [`LogGaborComplex`](@ref) for the `Complex(r, a)` version.
+
+# Arguments
+
+- `kernel_size::Dims{2}`: the Log Gabor kernel size. The axes at each dimension will be
+  `-r:r` if the size is odd.
+- `kernel_axes::NTuple{2, <:AbstractUnitRange}`: the axes of the Log Gabor kernel.
+- `ω`: the center frequency.
+- `θ`: the center orientation.
+
+# Keywords
+
+- `σω=1`: scale component for `ω`. Larger `σω` makes the filter more sensitive to center region.
+- `σθ=1`: scale component for `θ`. Larger `σθ` makes the filter less sensitive to orientation.
+
+# Examples
+
+To apply log gabor filter `g` on image `X`, one need to use convolution theorem, i.e.,
+`conv(A, K) == ifft(fft(A) .* fft(K))`. Because this `LogGabor` function generates Log Gabor
+kernel directly in spatial space, we don't need to apply `fft(K)` here:
+
+```jldoctest
+julia> using ImageFiltering, FFTW, TestImages, ImageCore
+
+julia> img = TestImages.shepp_logan(256);
+
+julia> kern = Kernel.LogGabor(size(img), 1/6, 0);
+
+julia> g_img = ifft(fft(channelview(img)) .* ifftshift(fft(kern))); # apply convolution theorem
+
+julia> @. Gray(abs(g_img));
+
+```
+
+# Extended help
+
+Mathematically, log gabor filter is defined in spatial space as the product of
+its frequency component `r` and angular component `a`:
+
+```math
+r(\\omega, \\theta) = \\exp(-\\frac{(\\log(\\omega/\\omega_0))^2}{2\\sigma_\\omega^2}) \\
+a(\\omega, \\theta) = \\exp(-\\frac{(\\theta - \\theta_0)^2}{2\\sigma_\\theta^2})
+```
+
+# References
+
+- [1] [What Are Log-Gabor Filters and Why Are They
+  Good?](https://www.peterkovesi.com/matlabfns/PhaseCongruency/Docs/convexpl.html)
+- [2] Kovesi, Peter. "Image features from phase congruency." _Videre: Journal of computer
+  vision research_ 1.3 (1999): 1-26.
+- [3] Field, David J. "Relations between the statistics of natural images and the response
+  properties of cortical cells." _Josa a_ 4.12 (1987): 2379-2394.
+"""
+struct LogGabor{T, AT<:LogGaborComplex} <: AbstractMatrix{T}
+    complex_data::AT
+end
+LogGabor(complex_data::AT) where AT<:LogGaborComplex = LogGabor{eltype(AT), AT}(complex_data)
+LogGabor(size_or_axes::Tuple, ω, θ; kwargs...) = LogGabor(size_or_axes; ω, θ, kwargs...)
+LogGabor(size_or_axes::Tuple; kwargs...) = LogGabor(LogGaborComplex(size_or_axes; kwargs...))
+
+@inline Base.size(kern::LogGabor) = size(kern.complex_data)
+@inline Base.axes(kern::LogGabor) = axes(kern.complex_data)
+Base.@propagate_inbounds function Base.getindex(kern::LogGabor, inds::Int...)
+    v = kern.complex_data[inds...]
+    return real(v) * imag(v)
+end
+
+
 # Utils
 
 function _to_axes(sz::Dims)