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