diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 4fa61f5..49f7f25 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -68,3 +68,29 @@ let grp = SUITE["ROF"]
end
end
end
+
+SUITE["Gabor"] = BenchmarkGroup()
+let grp = SUITE["Gabor"]
+ for sz in ((19, 19), (256, 256), (512, 512))
+ out = Matrix{ComplexF64}(undef, sz)
+ kern = Kernel.Gabor(sz, 2, 0)
+ grp["Gabor_ComplexF64"*"_"*sz2str(sz)] = @benchmarkable copyto!($out, $kern)
+
+ out = Matrix{ComplexF32}(undef, sz)
+ kern = Kernel.Gabor(sz, 2.0f0, 0.0f0)
+ grp["Gabor_ComplexF32"*"_"*sz2str(sz)] = @benchmarkable copyto!($out, $kern)
+ end
+end
+
+SUITE["LogGabor"] = BenchmarkGroup()
+let grp = SUITE["LogGabor"]
+ for sz in ((19, 19), (256, 256), (512, 512))
+ out = Matrix{ComplexF64}(undef, sz)
+ kern = Kernel.LogGabor(sz, 1/6, 0)
+ grp["LogGabor_ComplexF64"*"_"*sz2str(sz)] = @benchmarkable copyto!($out, $kern)
+
+ out = Matrix{ComplexF32}(undef, sz)
+ kern = Kernel.LogGabor(sz, 1.0f0/6, 0.0f0)
+ grp["LogGabor_ComplexF32"*"_"*sz2str(sz)] = @benchmarkable copyto!($out, $kern)
+ end
+end
diff --git a/docs/Project.toml b/docs/Project.toml
index f0b4b9c..68d33e0 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/docs/demos/filters/gabor_filter.jl b/docs/demos/filters/gabor_filter.jl
new file mode 100644
index 0000000..dc15c8a
--- /dev/null
+++ b/docs/demos/filters/gabor_filter.jl
@@ -0,0 +1,178 @@
+# ---
+# 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. A similar
+# example is [Log Gabor filter](@ref demo_log_gabor_filter).
+
+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)
diff --git a/docs/demos/filters/loggabor_filter.jl b/docs/demos/filters/loggabor_filter.jl
new file mode 100644
index 0000000..8d5e30a
--- /dev/null
+++ b/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
diff --git a/docs/src/function_reference.md b/docs/src/function_reference.md
index 217c564..9ace6f5 100644
--- a/docs/src/function_reference.md
+++ b/docs/src/function_reference.md
@@ -23,7 +23,9 @@ Kernel.gaussian
Kernel.DoG
Kernel.LoG
Kernel.Laplacian
-Kernel.gabor
+Kernel.Gabor
+Kernel.LogGabor
+Kernel.LogGaborComplex
Kernel.moffat
```
diff --git a/src/gaborkernels.jl b/src/gaborkernels.jl
new file mode 100644
index 0000000..b882629
--- /dev/null
+++ b/src/gaborkernels.jl
@@ -0,0 +1,323 @@
+module GaborKernels
+
+export Gabor, LogGabor, LogGaborComplex
+
+"""
+ 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
+
+
+"""
+ 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)
+ 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
diff --git a/src/kernel.jl b/src/kernel.jl
index 5d0993b..ff76a30 100644
--- a/src/kernel.jl
+++ b/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
@@ -424,66 +427,6 @@ function laplacian2d(alpha::Number=0)
return centered([lc lb lc; lb lm lb; lc lb lc])
end
-"""
- gabor(size_x,size_y,σ,θ,λ,γ,ψ) -> (k_real,k_complex)
-
-Returns a 2 Dimensional Complex Gabor kernel contained in a tuple where
-
- - `size_x`, `size_y` denote the size of the kernel
- - `σ` denotes the standard deviation of the Gaussian envelope
- - `θ` represents the orientation of the normal to the parallel stripes of a Gabor function
- - `λ` represents the wavelength of the sinusoidal factor
- - `γ` is the spatial aspect ratio, and specifies the ellipticity of the support of the Gabor function
- - `ψ` is the phase offset
-
-#Citation
-N. Petkov and P. Kruizinga, “Computational models of visual neurons specialised in the detection of periodic and aperiodic oriented visual stimuli: bar and grating cells,” Biological Cybernetics, vol. 76, no. 2, pp. 83–96, Feb. 1997. doi.org/10.1007/s004220050323
-"""
-function gabor(size_x::Integer, size_y::Integer, σ::Real, θ::Real, λ::Real, γ::Real, ψ::Real)
-
- σx = σ
- σy = σ/γ
- nstds = 3
- c = cos(θ)
- s = sin(θ)
-
- validate_gabor(σ,λ,γ)
-
- if(size_x > 0)
- xmax = floor(Int64,size_x/2)
- else
- @warn "The input parameter size_x should be positive. Using size_x = 6 * σx + 1 (Default value)"
- xmax = round(Int64,max(abs(nstds*σx*c),abs(nstds*σy*s),1))
- end
-
- if(size_y > 0)
- ymax = floor(Int64,size_y/2)
- else
- @warn "The input parameter size_y should be positive. Using size_y = 6 * σy + 1 (Default value)"
- ymax = round(Int64,max(abs(nstds*σx*s),abs(nstds*σy*c),1))
- end
-
- xmin = -xmax
- ymin = -ymax
-
- x = [j for i in xmin:xmax,j in ymin:ymax]
- y = [i for i in xmin:xmax,j in ymin:ymax]
- xr = x*c + y*s
- yr = -x*s + y*c
-
- kernel_real = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*cos.(2*(π/λ)*xr .+ ψ))
- kernel_imag = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*sin.(2*(π/λ)*xr .+ ψ))
-
- kernel = (kernel_real,kernel_imag)
- return kernel
-end
-
-function validate_gabor(σ::Real,λ::Real,γ::Real)
- if !(σ>0 && λ>0 && γ>0)
- throw(ArgumentError("The parameters σ, λ and γ must be positive numbers."))
- end
-end
-
"""
moffat(α, β, ls) -> k
@@ -537,4 +480,6 @@ if Base.VERSION >= v"1.4.2" && ccall(:jl_generating_output, Cint, ()) == 1
end
end
+include("kernel_deprecated.jl")
+
end
diff --git a/src/kernel_deprecated.jl b/src/kernel_deprecated.jl
new file mode 100644
index 0000000..1264aca
--- /dev/null
+++ b/src/kernel_deprecated.jl
@@ -0,0 +1,62 @@
+# Deprecated
+
+"""
+ gabor(size_x,size_y,σ,θ,λ,γ,ψ) -> (k_real,k_complex)
+
+Returns a 2 Dimensional Complex Gabor kernel contained in a tuple where
+
+ - `size_x`, `size_y` denote the size of the kernel
+ - `σ` denotes the standard deviation of the Gaussian envelope
+ - `θ` represents the orientation of the normal to the parallel stripes of a Gabor function
+ - `λ` represents the wavelength of the sinusoidal factor
+ - `γ` is the spatial aspect ratio, and specifies the ellipticity of the support of the Gabor function
+ - `ψ` is the phase offset
+
+# Citation
+
+N. Petkov and P. Kruizinga, “Computational models of visual neurons specialised in the detection of periodic and aperiodic oriented visual stimuli: bar and grating cells,” Biological Cybernetics, vol. 76, no. 2, pp. 83–96, Feb. 1997. doi.org/10.1007/s004220050323
+"""
+function gabor(size_x::Integer, size_y::Integer, σ::Real, θ::Real, λ::Real, γ::Real, ψ::Real)
+ Base.depwarn("use `Kernel.Gabor` instead.", :gabor)
+ σx = σ
+ σy = σ/γ
+ nstds = 3
+ c = cos(θ)
+ s = sin(θ)
+
+ validate_gabor(σ,λ,γ)
+
+ if(size_x > 0)
+ xmax = floor(Int64,size_x/2)
+ else
+ @warn "The input parameter size_x should be positive. Using size_x = 6 * σx + 1 (Default value)"
+ xmax = round(Int64,max(abs(nstds*σx*c),abs(nstds*σy*s),1))
+ end
+
+ if(size_y > 0)
+ ymax = floor(Int64,size_y/2)
+ else
+ @warn "The input parameter size_y should be positive. Using size_y = 6 * σy + 1 (Default value)"
+ ymax = round(Int64,max(abs(nstds*σx*s),abs(nstds*σy*c),1))
+ end
+
+ xmin = -xmax
+ ymin = -ymax
+
+ x = [j for i in xmin:xmax,j in ymin:ymax]
+ y = [i for i in xmin:xmax,j in ymin:ymax]
+ xr = x*c + y*s
+ yr = -x*s + y*c
+
+ kernel_real = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*cos.(2*(π/λ)*xr .+ ψ))
+ kernel_imag = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*sin.(2*(π/λ)*xr .+ ψ))
+
+ kernel = (kernel_real,kernel_imag)
+ return kernel
+end
+
+function validate_gabor(σ::Real,λ::Real,γ::Real)
+ if !(σ>0 && λ>0 && γ>0)
+ throw(ArgumentError("The parameters σ, λ and γ must be positive numbers."))
+ end
+end
diff --git a/test/cuda/Project.toml b/test/cuda/Project.toml
index da371cd..6673f9a 100644
--- a/test/cuda/Project.toml
+++ b/test/cuda/Project.toml
@@ -1,10 +1,12 @@
[deps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e"
ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
ImageQualityIndexes = "2996bd0c-7a13-11e9-2da2-2f5ce47296a9"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
diff --git a/test/cuda/gaborkernels.jl b/test/cuda/gaborkernels.jl
new file mode 100644
index 0000000..4b217f7
--- /dev/null
+++ b/test/cuda/gaborkernels.jl
@@ -0,0 +1,66 @@
+@testset "GaborKernels" begin
+
+@testset "Gabor" begin
+ # Gray
+ img = float32.(TestImages.shepp_logan(127))
+ kern = Kernel.Gabor(size(img), 3.0f0, 0f0)
+ img_freq = fft(channelview(img))
+ kern_freq = ifftshift(fft(kern))
+ out = abs.(ifft(img_freq .* kern_freq))
+
+ # TODO(johnnychen94): currently Gabor can't be converted to CuArray directly and thus
+ # FFTW is applied in the CPU side.
+ img_cu = CuArray(img)
+ img_freq = fft(channelview(img_cu))
+ kern_freq = CuArray(ifftshift(fft(kern)))
+ out_cu = abs.(ifft(img_freq .* kern_freq))
+
+ @test out ≈ Array(out_cu)
+
+ # RGB
+ img = float32.(testimage("monarch"))
+ kern = Kernel.Gabor(size(img), 3.0f0, 0f0)
+ kern_freq = reshape(ifftshift(fft(kern)), 1, size(kern)...)
+ img_freq = fft(channelview(img), 2:3)
+ out = colorview(RGB, abs.(ifft(img_freq .* kern_freq)))
+
+ img_cu = CuArray(img)
+ kern_freq = CuArray(reshape(ifftshift(fft(kern)), 1, size(kern)...))
+ img_freq = fft(channelview(img_cu), 2:3)
+ out_cu = colorview(RGB, abs.(ifft(img_freq .* kern_freq)))
+
+ @test out ≈ Array(out_cu)
+end
+
+@testset "LogGabor" begin
+ # Gray
+ img = float32.(TestImages.shepp_logan(127))
+ kern = Kernel.LogGabor(size(img), 1.0f0/6, 0f0)
+ kern_freq = OffsetArrays.no_offset_view(ifftshift(kern))
+ img_freq = fft(channelview(img))
+ out = abs.(ifft(kern_freq .* img_freq))
+
+ # TODO(johnnychen94): remove this no_offset_view wrapper
+ img_cu = CuArray(img)
+ kern_freq = CuArray(OffsetArrays.no_offset_view(ifftshift(kern)))
+ img_freq = fft(channelview(img_cu))
+ out_cu = abs.(ifft(kern_freq .* img_freq))
+
+ @test out ≈ Array(out_cu)
+
+ # RGB
+ img = float32.(testimage("monarch"))
+ kern = Kernel.LogGabor(size(img), 1.0f0/6, 0f0)
+ kern_freq = reshape(ifftshift(kern), 1, size(kern)...)
+ img_freq = fft(channelview(img), 2:3)
+ out = colorview(RGB, abs.(ifft(img_freq .* kern_freq)))
+
+ img_cu = CuArray(img)
+ kern_freq = CuArray(reshape(ifftshift(kern), 1, size(kern)...))
+ img_freq = fft(channelview(img_cu), 2:3)
+ out_cu = colorview(RGB, abs.(ifft(img_freq .* kern_freq)))
+
+ @test out ≈ Array(out_cu)
+end
+
+end
diff --git a/test/cuda/runtests.jl b/test/cuda/runtests.jl
index 5789d00..f86b6dc 100644
--- a/test/cuda/runtests.jl
+++ b/test/cuda/runtests.jl
@@ -7,11 +7,15 @@ using ImageBase
using ImageQualityIndexes
using Test
using Random
+using FFTW
+using OffsetArrays
+using CUDA.Adapt
CUDA.allowscalar(false)
@testset "ImageFiltering" begin
if CUDA.functional()
include("models.jl")
+ include("gaborkernels.jl")
end
end
diff --git a/test/gabor.jl b/test/deprecated.jl
similarity index 95%
rename from test/gabor.jl
rename to test/deprecated.jl
index 5c77508..63eee00 100644
--- a/test/gabor.jl
+++ b/test/deprecated.jl
@@ -1,5 +1,3 @@
-using ImageFiltering, Test, Statistics
-
@testset "gabor" begin
σx = 8
σy = 12
@@ -7,7 +5,6 @@ using ImageFiltering, Test, Statistics
size_y = 6*σy+1
γ = σx/σy
# zero size forces default kernel width, with warnings
- @info "Four warnings are expected"
kernel = Kernel.gabor(0,0,σx,0,5,γ,0)
@test isequal(size(kernel[1]),(size_x,size_y))
kernel = Kernel.gabor(0,0,σx,π,5,γ,0)
diff --git a/test/gaborkernels.jl b/test/gaborkernels.jl
new file mode 100644
index 0000000..dddda84
--- /dev/null
+++ b/test/gaborkernels.jl
@@ -0,0 +1,207 @@
+@testset "GaborKernels" begin
+
+function is_symmetric(X)
+ @assert all(isodd, size(X))
+ rst = map(CartesianIndices(size(X).÷2)) do i
+ X[i] ≈ X[-i]
+ end
+ return all(rst)
+end
+
+@testset "Gabor" begin
+ @testset "API" begin
+ # Normally it just return a ComplexF64 matrix if users are careless
+ kern = @inferred Kernel.Gabor((11, 11), 2, 0)
+ @test size(kern) == (11, 11)
+ @test axes(kern) == (-5:5, -5:5)
+ @test eltype(kern) == ComplexF64
+ # ensure that this is an efficient lazy array
+ @test isbitstype(typeof(kern))
+
+ # but still allow construction of ComplexF32 matrix
+ kern = @inferred Kernel.Gabor((11, 11), 2.0f0, 0.0f0)
+ @test eltype(kern) == ComplexF32
+
+ # passing axes explicitly allows building a subregion of it
+ kern1 = @inferred Kernel.Gabor((1:10, 1:10), 2.0, 0.0)
+ kern2 = @inferred Kernel.Gabor((21, 21), 2.0, 0.0)
+ @test kern2[1:end, 1:end] == kern1
+
+ # keyword version makes it more explicit
+ kern1 = @inferred Kernel.Gabor((11, 11), 2, 0)
+ kern2 = @inferred Kernel.Gabor((11, 11); wavelength=2, orientation=0)
+ @test kern1 === kern2
+ # and currently we can't use both positional and keyword
+ @test_throws UndefKeywordError Kernel.Gabor((5, 5))
+ @test_throws MethodError Kernel.Gabor((11, 11), 2)
+ @test_throws MethodError Kernel.Gabor((11, 11), 2; orientation=0)
+ @test_throws MethodError Kernel.Gabor((11, 11), 2; wavelength=3)
+ @test_throws MethodError Kernel.Gabor((11, 11), 2, 0; wavelength=3)
+
+ # test default keyword values
+ kern1 = @inferred Kernel.Gabor((11, 11), 2, 0)
+ kern2 = @inferred Kernel.Gabor((11, 11), 2, 0; bandwidth=1, phase_offset=0, aspect_ratio=0.5)
+ @test kern1 === kern2
+ end
+
+ @testset "getindex" begin
+ kern = @inferred Kernel.Gabor((11, 11), 2, 0)
+ @test kern[0, 0] ≈ 1
+ @test kern[1] == kern[-5, -5]
+ @test kern[end] == kern[5, 5]
+ end
+
+ @testset "symmetricity" begin
+ # for some special orientation we can easily check the symmetric
+ for θ in [0, π/2, -π/2, π, -π]
+ kern = Kernel.Gabor((11, 11), 2, θ)
+ @test is_symmetric(kern)
+ end
+ # for other orientations the symmetricity still holds, just that it doesn't
+ # align along the x-y axis.
+ kern = Kernel.Gabor((11, 11), 2, π/4)
+ @test !is_symmetric(kern)
+ end
+
+ @testset "invalid inputs" begin
+ if VERSION >= v"1.7-rc1"
+ msg = "the expected kernel size is expected to be larger than"
+ @test_warn msg Kernel.Gabor((11, 11), 5, 0)
+ msg = "`wavelength` should be equal to or greater than 2"
+ @test_warn msg Kernel.Gabor((11, 11), 1, 0)
+ end
+ @test_throws ArgumentError Kernel.Gabor((-1, -1), 2, 0)
+ @test_throws ArgumentError Kernel.Gabor((11, 1), 2, 0; bandwidth=-1)
+ @test_throws ArgumentError Kernel.Gabor((11, 1), 2, 0; aspect_ratio=-1)
+ end
+
+ @testset "numeric" begin
+ function gabor(size_x, size_y, σ, θ, λ, γ, ψ)
+ # plain implementation https://en.wikipedia.org/wiki/Gabor_filter
+ # See also: https://github.com/JuliaImages/ImageFiltering.jl/pull/57
+ σx, σy = σ, σ/γ
+ s, c = sincos(θ)
+
+ xmax = floor(Int, size_x/2)
+ ymax = floor(Int, size_y/2)
+ xmin = -xmax
+ ymin = -ymax
+
+ # The original implementation transposes x-y axis
+ # x = [j for i in xmin:xmax,j in ymin:ymax]
+ # y = [i for i in xmin:xmax,j in ymin:ymax]
+ x = [i for i in xmin:xmax,j in ymin:ymax]
+ y = [j for i in xmin:xmax,j in ymin:ymax]
+ xr = x*c + y*s
+ yr = -x*s + y*c
+
+ kernel_real = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*cos.(2*(π/λ)*xr .+ ψ))
+ kernel_imag = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*sin.(2*(π/λ)*xr .+ ψ))
+
+ kernel = (kernel_real, kernel_imag)
+ return kernel
+ end
+
+ kern1 = Kernel.Gabor((11, 11), 2, 0)
+ σ, θ, λ, γ, ψ = kern1.σ, kern1.θ, kern1.λ, kern1.γ, kern1.ψ
+ kern2_real, kern2_imag = gabor(map(length, kern1.ax)..., σ, θ, λ, γ, ψ)
+ @test collect(real.(kern1)) ≈ kern2_real
+ @test collect(imag.(kern1)) ≈ kern2_imag
+ end
+
+ @testset "applications" begin
+ img = TestImages.shepp_logan(127)
+ kern = Kernel.Gabor((19, 19), 2, 0)
+ img_out1 = imfilter(img, real.(kern))
+
+ kern_freq = freqkernel(real.(kern), size(img))
+ img_out2 = Gray.(real.(ifft(fft(channelview(img)) .* kern_freq)))
+
+ # Except for the boundary, these two methods produce the same result
+ @test img_out1[10:end-10, 10:end-10] ≈ img_out2[10:end-10, 10:end-10]
+ 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
diff --git a/test/runtests.jl b/test/runtests.jl
index cfb5ae5..a36ad8c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,8 @@ using TestImages
using ImageQualityIndexes
import StaticArrays
using Random
+using FFTW
+using Statistics
@testset "Project meta quality checks" begin
# Ambiguity test
@@ -43,7 +45,7 @@ include("gradient.jl")
include("mapwindow.jl")
include("extrema.jl")
include("basic.jl")
-include("gabor.jl")
+include("gaborkernels.jl")
include("models.jl")
@@ -70,4 +72,7 @@ if CUDA_FUNCTIONAL
else
@warn "CUDA test: disabled"
end
+
+@info "Beginning deprecation tests, warnings are expected"
+include("deprecated.jl")
nothing