introduce log gabor filter: LogGabor and LogGaborComplex

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,107 @@
+# ---
+# 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 small helper functions 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
+#
+# Because the filter is defined on frequency space, we can 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.
+#
+# Also, since Log Gabor kernel is defined around center point (0, 0), we have to apply
+# `ifftshift` first before we do pointwise-multiplication.
+
+img = TestImages.shepp_logan(127)
+kern = Kernel.LogGaborComplex(size(img), 50, π/4)
+## we don't need to call `fft(kern)` here because it's already on frequency space
+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...)
@@ -169,6 +169,147 @@ end
     return exp(-(xr^2 + yr^2)/(2σ^2)) * cis(xr*λ_scaled + ψ)
 end
 
+
+"""
+    LogGaborComplex(size_or_axes, ω, θ; σω=1, σθ=1, normalize=true)
+
+Generate 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 and example 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
+    normalize::Bool
+
+    # cache values
+    freq_scale::Tuple{TP, TP} # only used when normalize is true
+    ω_denom::TP # 1/(2(log(σω/ω))^2)
+    θ_denom::TP # 1/(2σθ^2)
+    function LogGaborComplex{T,TP,R}(ax::Tuple{R,R}, ω::TP, θ::TP, σω::TP, σθ::TP, normalize::Bool) where {T,TP,R}
+        σω > 0 || throw(ArgumentError("`σω` should be positive: $σω"))
+        σθ > 0 || throw(ArgumentError("`σθ` should be positive: $σθ"))
+        ω_denom = 1/(2(log(σω/ω))^2)
+        θ_denom = 1/(2σθ^2)
+        freq_scale = map(r->1/length(r), ax)
+        new{T,TP,R}(ax, ω, θ, σω, σθ, normalize, freq_scale, ω_denom, θ_denom)
+    end
+end
+function LogGaborComplex(
+        size_or_axes::Tuple, ω::Real, θ::Real;
+        σω::Real=1, σθ::Real=1, normalize::Bool=true,
+    )
+    params = float.(promote(ω, θ, σω, σθ))
+    T = typeof(params[1])
+    ax = _to_axes(size_or_axes)
+    LogGaborComplex{Complex{T}, T, typeof(first(ax))}(ax, params..., normalize)
+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
+    # Although in `getindex`, the computation is heavy enough that this runtime if-branch is
+    # harmless to the overall performance at all
+    if kern.normalize
+        # normalize: from reference [1] of LogGabor
+        # By changing division to multiplication gives about 5-10% performance boost
+        x, y = (x, y) .* kern.freq_scale
+    end
+    ω = 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 Complex(r, a)
+end
+
+
+"""
+    LogGabor(size_or_axes, ω, θ; σω=1, σθ=1, normalize=true)
+
+Generate 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.
+- `normalize=true`: whether to normalize the frequency domain into [-0.5, 0.5]x[-0.5, 0.5]
+  domain via `inds = inds./size(kern)`. For image-related tasks where the [Weber–Fechner
+  law](https://en.wikipedia.org/wiki/Weber%E2%80%93Fechner_law) applies, this usually
+  provides more stable pipeline.
+
+# 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(centered(fft(channelview(img))) .* ifftshift(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{real(eltype(AT)), AT}(complex_data)
+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...)
+    # cache the result to avoid repeated computation
+    v = kern.complex_data[inds...]
+    return real(v) * imag(v)
+end
+
+
 # Utils
 
 function _to_axes(sz::Dims)

test/gaborkernels.jl

@@ -122,4 +122,86 @@ end
     end
 end
 
+
+@testset "LogGabor" begin
+    @testset "API" begin
+        # LogGabor: r * a
+        # LogGaborComplex: Complex(r, a)
+        kern = @inferred Kernel.LogGabor((11, 11), 1/6, 0)
+        kern_c = Kernel.LogGaborComplex((11, 11), 1/6, 0)
+        @test kern == @. real(kern_c) * imag(kern_c)
+
+        # Normally it just return a ComplexF64 matrix if users are careless
+        kern = @inferred Kernel.LogGabor((11, 11), 2, 0)
+        @test size(kern) == (11, 11)
+        @test axes(kern) == (-5:5, -5:5)
+        @test eltype(kern) == Float64
+        # ensure that this is an efficient lazy array
+        @test isbitstype(typeof(kern))
+
+        # but still allow construction of ComplexF32 matrix
+        kern = @inferred Kernel.LogGabor((11, 11), 1.0f0/6, 0.0f0)
+        @test eltype(kern) == Float32
+
+        # passing axes explicitly allows building a subregion of it
+        kern1 = @inferred Kernel.LogGabor((1:10, 1:10), 1/6, 0)
+        @test axes(kern1) == (1:10, 1:10)
+        kern2 = @inferred Kernel.LogGabor((21, 21), 1/6, 0)
+        # but they are not the same in normalize=true mode
+        @test kern2[1:end, 1:end] != kern1
+
+        # when normalize=false, they're the same
+        kern1 = @inferred Kernel.LogGabor((1:10, 1:10), 1/6, 0; normalize=false)
+        kern2 = @inferred Kernel.LogGabor((21, 21), 1/6, 0; normalize=false)
+        @test kern2[1:end, 1:end] == kern1
+
+        # test default keyword values
+        kern1 = @inferred Kernel.LogGabor((11, 11), 2, 0)
+        kern2 = @inferred Kernel.LogGabor((11, 11), 2, 0; σω=1, σθ=1, normalize=true)
+        @test kern1 === kern2
+
+        @test_throws ArgumentError Kernel.LogGabor((11, 11), 2, 0; σω=-1)
+        @test_throws ArgumentError Kernel.LogGabor((11, 11), 2, 0; σθ=-1)
+        @test_throws ArgumentError Kernel.LogGaborComplex((11, 11), 2, 0; σω=-1)
+        @test_throws ArgumentError Kernel.LogGaborComplex((11, 11), 2, 0; σθ=-1)
+    end
+
+    @testset "symmetricity" begin
+        kern = Kernel.LogGaborComplex((11, 11), 1/12, 0)
+        # the real part is an even-symmetric function
+        @test is_symmetric(real.(kern))
+        # the imaginary part is a mirror along y-axis
+        rows =  [imag.(kern[i, :]) for i in axes(kern, 1)]
+        @test all(map(is_symmetric,  rows))
+
+        # NOTE(johnnychen94): I'm not sure the current implementation is the standard or
+        # correct. This numerical references are generated only to make sure future changes
+        # don't accidentally break it.
+        # Use N0f8 to ignore "insignificant" numerical changes to keep unit test happy
+        ref_real = N0f8[
+            0.741  0.796  0.82   0.796  0.741
+            0.796  0.886  0.941  0.886  0.796
+            0.82   0.941  0.0    0.941  0.82
+            0.796  0.886  0.941  0.886  0.796
+            0.741  0.796  0.82   0.796  0.741
+        ]
+        ref_imag = N0f8[
+            0.063  0.027  0.008  0.027  0.063
+            0.125  0.063  0.008  0.063  0.125
+            0.29   0.29   1.0    0.29   0.29
+            0.541  0.733  1.0    0.733  0.541
+            0.733  0.898  1.0    0.898  0.733
+        ]
+        kern = Kernel.LogGaborComplex((5, 5), 1/12, 0)
+        @test collect(N0f8.(real(kern))) == ref_real
+        @test collect(N0f8.(imag(kern))) == ref_imag
+    end
+
+    @testset "applications" begin
+        img = TestImages.shepp_logan(127)
+        kern = Kernel.LogGabor(size(img), 1/100, π/4)
+        out = @test_nowarn ifft(centered(fft(channelview(img))) .* ifftshift(kern))
+    end
+end
+
 end