From f25ddf9c3ff4198c73bff2880321a7b9981982bf Mon Sep 17 00:00:00 2001 From: JohnnyChen Date: Tue, 2 Nov 2021 05:51:21 +0800 Subject: [PATCH 1/5] new gabor filter: lazy `Gabor` array This commit introduces an lazy array interface for the Gabor filter, benchmark shows that it's 2-3x faster than the previous version. The documentation is heavly rewrote to explain this filter and all its parameters. --- docs/Project.toml | 2 + docs/demos/filters/gabor_filter.jl | 177 ++++++++++++++++++++++++++++ docs/src/function_reference.md | 2 +- src/gaborkernels.jl | 182 +++++++++++++++++++++++++++++ src/kernel.jl | 5 +- test/gaborkernels.jl | 125 ++++++++++++++++++++ test/runtests.jl | 4 +- 7 files changed, 494 insertions(+), 3 deletions(-) create mode 100644 docs/demos/filters/gabor_filter.jl create mode 100644 src/gaborkernels.jl create mode 100644 test/gaborkernels.jl 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..a90f365 --- /dev/null +++ b/docs/demos/filters/gabor_filter.jl @@ -0,0 +1,177 @@ +# --- +# title: Gabor filter +# id: demo_gabor_filter +# cover: assets/gabor.png +# author: Johnny Chen +# date: 2021-11-01 +# --- + +# This example shows how one can apply spatial space kernesl [`Gabor`](@ref Kernel.Gabor) +# using fourier transformation and convolution theorem to extract image features. + +using ImageCore, ImageShow, ImageFiltering # or you could just `using Images` +using FFTW +using TestImages + +# ## Definition +# +# Mathematically, Gabor kernel is defined in spatial space: +# +# ```math +# g(x, y) = \exp(-\frac{x'^2 + \gamma^2y'^2}{2\sigma^2})\exp(i(2\pi\frac{x'}{\lambda} + \psi)) +# ``` +# where ``i`` is imaginary unit `Complex(0, 1)`, and +# ```math +# x' = x\cos\theta + x\sin\theta \\ +# y' = -x\sin\theta + y\cos\theta +# ``` +# + +# First of all, Gabor kernel is a complex-valued matrix: + +kern = Kernel.Gabor((10, 10), 2, 0.1) + +# !!! tip "Lazy array" +# The `Gabor` type is a lazy array, which means when you build the Gabor kernel, you +# actually don't need to allocate any memories. +# +# ```julia +# using BenchmarkTools +# kern = @btime Kernel.Gabor((64, 64), 5, 0); # 36.481 ns (0 allocations: 0 bytes) +# @btime collect($kern); # 75.278 μs (2 allocations: 64.05 KiB) +# ``` + +# To explain the parameters of Gabor filter, let's introduce some small helpers function to +# display complex-valued kernels. +show_phase(kern) = @. Gray(log(abs(imag(kern)) + 1)) +show_mag(kern) = @. Gray(log(abs(real(kern)) + 1)) +show_abs(kern) = @. Gray(log(abs(kern) + 1)) +nothing #hide + +# ## Keywords +# +# ### `wavelength` (λ) +# λ specifies the wavelength of the sinusoidal factor. + +bandwidth, orientation, phase_offset, aspect_ratio = 1, 0, 0, 0.5 +f(wavelength) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset)) +mosaic(f.((5, 10, 15)), nrow=1) + +# ### `orientation` (θ) +# θ specifies the orientation of the normal to the parallel stripes of a Gabor function. + +wavelength, bandwidth, phase_offset, aspect_ratio = 10, 1, 0, 0.5 +f(orientation) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset)) +mosaic(f.((0, π/4, π/2)), nrow=1) + +# ### `phase_offset` (ψ) + +wavelength, bandwidth, orientation, aspect_ratio = 10, 1, 0, 0.5 +f(phase_offset) = show_phase(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset)) +mosaic(f.((-π/2, 0, π/2, π)), nrow=1) + +# ### `aspect_ratio` (γ) +# γ specifies the ellipticity of the support of the Gabor function. + +wavelength, bandwidth, orientation, phase_offset = 10, 1, 0, 0 +f(aspect_ratio) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset)) +mosaic(f.((0.5, 1, 2)), nrow=1) + +# ### `bandwidth` (b) +# The half-response spatial frequency bandwidth (b) of a Gabor filter is related to the +# ratio σ / λ, where σ and λ are the standard deviation of the Gaussian factor of the Gabor +# function and the preferred wavelength, respectively, as follows: +# +# ```math +# b = \log_2\frac{\frac{\sigma}{\lambda}\pi + \sqrt{\frac{\ln 2}{2}}}{\frac{\sigma}{\lambda}\pi - \sqrt{\frac{\ln 2}{2}}} +# ``` + +wavelength, orientation, phase_offset, aspect_ratio = 10, 0, 0, 0.5 +f(bandwidth) = show_abs(Kernel.Gabor((100, 100); wavelength, bandwidth, orientation, aspect_ratio, phase_offset)) +mosaic(f.((0.5, 1, 2)), nrow=1) + +# ## Examples +# +# There are two options to apply a spatial space kernel: 1) via `imfilter`, and 2) via the +# convolution theorem. + +# ### `imfilter` +# +# [`imfilter`](@ref) does not require the kernel size to be the same as the image size. +# Usually kernel size is at least 5 times larger than the wavelength. + +img = TestImages.shepp_logan(127) +kern = Kernel.Gabor((19, 19), 3, 0) +out = imfilter(img, real.(kern)) +mosaic(img, show_abs(kern), show_mag(out); nrow=1) + +# ### convolution theorem + +# The [convolution theorem](https://en.wikipedia.org/wiki/Convolution_theorem) tells us +# that `fft(conv(X, K))` is equivalent to `fft(X) .* fft(K)`. Because Gabor kernel is +# defined with regard to its center point (0, 0), we need to do `ifftshift` first so that +# the frequency centers of both `fft(X)` and `fft(kern)` align well. + +kern = Kernel.Gabor(size(img), 3, 0) +out = ifft(fft(channelview(img)) .* ifftshift(fft(kern))) +mosaic(img, show_abs(kern), show_mag(out); nrow=1) + +# As you may have notice, using convolution theorem generates different results. This is +# simply because the kernel size are different. If we create a smaller kernel, we then need +# to apply [`freqkernel`](@ref) first so that we can do element-wise multiplication. + +## freqkernel = zero padding + fftshift + fft +kern = Kernel.Gabor((19, 19), 3, 0) +kern_freq = freqkernel(real.(kern), size(img)) +out = ifft(fft(channelview(img)) .* kern_freq) +mosaic(img, show_abs(kern), show_mag(out); nrow=1) + +# !!! note "Performance on different kernel size" +# When the kernel size is small, `imfilter` works more efficient than fft-based +# convolution. This benchmark isn't backed by CI so the result might be outdated, but +# you get the idea. +# +# ```julia +# using BenchmarkTools +# img = TestImages.shepp_logan(127); +# +# kern = Kernel.Gabor((19, 19), 3, 0); +# fft_conv(img, kern) = ifft(fft(channelview(img)) .* freqkernel(real.(kern), size(img))) +# @btime imfilter($img, real.($kern)); # 236.813 μs (118 allocations: 418.91 KiB) +# @btime fft_conv($img, $kern) # 1.777 ms (127 allocations: 1.61 MiB) +# +# kern = Kernel.Gabor(size(img), 3, 0) +# fft_conv(img, kern) = ifft(fft(channelview(img)) .* ifftshift(fft(kern))) +# @btime imfilter($img, real.($kern)); # 5.318 ms (163 allocations: 5.28 MiB) +# @btime fft_conv($img, $kern); # 2.218 ms (120 allocations: 1.73 MiB) +# ``` +# + +## Filter bank + +# A filter bank is just a list of filter kernels, applying the filter bank generates +# multiple outputs: + +filters = [Kernel.Gabor(size(img), 3, θ) for θ in -π/2:π/4:π/2]; +X_freq = fft(channelview(img)) +out = map(filters) do kern + ifft(X_freq .* ifftshift(fft(kern))) +end +mosaic( + map(show_abs, filters)..., + map(show_abs, out)...; + nrow=2, rowmajor=true +) + +## save covers #src +using FileIO #src +mkpath("assets") #src +filters = [Kernel.Gabor((32, 32), 5, θ) for θ in range(-π/2, stop=π/2, length=9)] #src +save("assets/gabor.png", mosaic(map(show_abs, filters); nrow=3, npad=2, fillvalue=Gray(1)); fps=2) #src + +# # References +# +# - [1] [Wikipedia: Gabor filter](https://en.wikipedia.org/wiki/Gabor_filter) +# - [2] [Wikipedia: Gabor transformation](https://en.wikipedia.org/wiki/Gabor_transform) +# - [3] [Wikipedia: Gabor atom](https://en.wikipedia.org/wiki/Gabor_atom) +# - [4] [Gabor filter for image processing and computer vision](http://matlabserver.cs.rug.nl/edgedetectionweb/web/edgedetection_params.html) diff --git a/docs/src/function_reference.md b/docs/src/function_reference.md index 217c564..3edea3c 100644 --- a/docs/src/function_reference.md +++ b/docs/src/function_reference.md @@ -23,7 +23,7 @@ Kernel.gaussian Kernel.DoG Kernel.LoG Kernel.Laplacian -Kernel.gabor +Kernel.Gabor Kernel.moffat ``` diff --git a/src/gaborkernels.jl b/src/gaborkernels.jl new file mode 100644 index 0000000..46d26e8 --- /dev/null +++ b/src/gaborkernels.jl @@ -0,0 +1,182 @@ +module GaborKernels + +export Gabor + +""" + Gabor(size_or_axes, wavelength, orientation; kwargs...) + Gabor(size_or_axes; wavelength, orientation, kwargs...) + +Generate the 2-D Gabor kernel in the spatial space. + +# Arguments + +## `kernel_size::Dims{2}` or `kernel_axes::NTuple{2,<:AbstractUnitRange}` + +Specifies either the size or axes of the output kernel. The axes at each dimension will be +`-r:r` if the size is odd. + +## `wavelength::Real`(λ>=2) + +This is the wavelength of the sinusoidal factor of the Gabor filter kernel and herewith the +preferred wavelength of this filter. Its value is specified in pixels. + +The value `λ=2` should not be used in combination with phase offset `ψ=-π/2` or `ψ=π/2` +because in these cases the Gabor function is sampled in its zero crossings. + +In order to prevent the occurence of undesired effects at the image borders, the wavelength +value should be smaller than one fifth of the input image size. + +## `orientation::Real`(θ∈[0, 2pi]): + +This parameter specifies the orientation of the normal to the parallel stripes of a Gabor +function [3]. + +# Keywords + +## `bandwidth=1` (b>0) + +The half-response spatial frequency bandwidth b (in octaves) of a Gabor filter is related to +the ratio σ / λ, where σ and λ are the standard deviation of the Gaussian factor of the +Gabor function and the preferred wavelength, respectively, as follows: + +```math +b = \\log_2\\frac{\\frac{\\sigma}{\\lambda}\\pi + \\sqrt{\\frac{\\ln 2}{2}}}{\\frac{\\sigma}{\\lambda}\\pi - \\sqrt{\\frac{\\ln 2}{2}}} +``` + +## `aspect_ratio=0.5`(γ>0) + +This parameter, called more precisely the spatial aspect ratio, specifies the ellipticity of +the support of the Gabor function [3]. For `γ = 1`, the support is circular. For γ < 1 the +support is elongated in orientation of the parallel stripes of the function. + +## `phase_offset=0` (ψ∈[-π, π]) + +The values `0` and `π` correspond to center-symmetric 'center-on' and 'center-off' +functions, respectively, while -π/2 and π/2 correspond to anti-symmetric functions. All +other cases correspond to asymmetric functions. + +# Examples + +There are two ways to use gabor filters: 1) via [`imfilter`](@ref), 2) via `fft` and +convolution theorem. Usually `imfilter` requires a small kernel to work efficiently, while +using `fft` can be more efficient when the kernel has the same size of the image. + +## imfilter + +```jldoctest gabor_example +julia> using ImageFiltering, FFTW, TestImages, ImageCore + +julia> img = TestImages.shepp_logan(256); + +julia> kern = Kernel.Gabor((19, 19), 3, 0); # usually a small kernel size + +julia> g_img = imfilter(img, real.(kern)); + +``` + +## convolution theorem + +The convolution theorem is stated as `fft(conv(A, K)) == fft(A) .* fft(K)`: + +```jldoctest gabor_example +julia> kern = Kernel.Gabor(size(img), 3, 0); + +julia> g_img = ifft(fft(channelview(img)) .* ifftshift(fft(kern))); # apply convolution theorem + +julia> @. Gray(abs(g_img)); + +``` + +See the [gabor filter demo](@ref demo_gabor_filter) for more examples with images. + +# Extended help + +Mathematically, gabor filter is defined in spatial space: + +```math +g(x, y) = \\exp(-\\frac{x'^2 + \\gamma^2y'^2}{2\\sigma^2})\\exp(i(2\\pi\\frac{x'}{\\lambda} + \\psi)) +``` + +where ``x' = x\\cos\\theta + y\\sin\\theta`` and ``y' = -x\\sin\\theta + y\\cos\\theta``. + +# References + +- [1] [Wikipedia: Gabor filter](https://en.wikipedia.org/wiki/Gabor_filter) +- [2] [Wikipedia: Gabor transformation](https://en.wikipedia.org/wiki/Gabor_transform) +- [3] [Wikipedia: Gabor atom](https://en.wikipedia.org/wiki/Gabor_atom) +- [4] [Gabor filter for image processing and computer vision](http://matlabserver.cs.rug.nl/edgedetectionweb/web/edgedetection_params.html) +""" +struct Gabor{T<:Complex, TP<:Real, R<:AbstractUnitRange} <: AbstractMatrix{T} + ax::Tuple{R,R} + λ::TP + θ::TP + ψ::TP + σ::TP + γ::TP + + # cache values + sc::Tuple{TP, TP} # sincos(θ) + λ_scaled::TP # 2π/λ + function Gabor{T,TP,R}(ax::Tuple{R,R}, λ::TP, θ::TP, ψ::TP, σ::TP, γ::TP) where {T,TP,R} + λ > 0 || throw(ArgumentError("`λ` should be positive: $λ")) + new{T,TP,R}(ax, λ, θ, ψ, σ, γ, sincos(θ), 2π/λ) + end +end +function Gabor(size_or_axes::Tuple; wavelength, orientation, kwargs...) + Gabor(size_or_axes, wavelength, orientation; kwargs...) +end +function Gabor( + size_or_axes::Tuple, + wavelength::Real, + orientation::Real; + bandwidth::Real=1.0f0, + phase_offset::Real=0.0f0, + aspect_ratio::Real=0.5f0, + ) + bandwidth > 0 || throw(ArgumentError("`bandwidth` should be positive: $bandwidth")) + aspect_ratio > 0 || throw(ArgumentError("`aspect_ratio` should be positive: $aspect_ratio")) + wavelength >= 2 || @warn "`wavelength` should be equal to or greater than 2" wavelength + # we still follow the math symbols in the implementation + λ, γ, ψ = wavelength, aspect_ratio, phase_offset + # for orientation: Julia follow column-major order, thus here we compensate the orientation by π/2 + θ = convert(float(typeof(orientation)), mod2pi(orientation+π/2)) + # follows reference [4] + σ = convert(float(typeof(λ)), (λ/π)*sqrt(0.5log(2)) * (2^bandwidth+1)/(2^bandwidth-1)) + + params = float.(promote(λ, θ, ψ, σ, γ)) + T = typeof(params[1]) + + ax = _to_axes(size_or_axes) + all(map(length, ax) .> 0) || throw(ArgumentError("Kernel size should be positive: $size_or_axes")) + if any(r->5λ > length(r), ax) + expected_size = @. 5λ * length(ax) + @warn "for wavelength `λ=$λ` the expected kernel size is expected to be larger than $expected_size." + end + + Gabor{Complex{T}, T, typeof(first(ax))}(ax, params...) +end + +@inline Base.size(kern::Gabor) = map(length, kern.ax) +@inline Base.axes(kern::Gabor) = kern.ax +@inline function Base.getindex(kern::Gabor, x::Int, y::Int) + ψ, σ, γ = kern.ψ, kern.σ, kern.γ + s, c = kern.sc # sincos(θ) + λ_scaled = kern.λ_scaled # 2π/λ + + xr = x*c + y*s + yr = -x*s + y*c + yr = γ * yr + return exp(-(xr^2 + yr^2)/(2σ^2)) * cis(xr*λ_scaled + ψ) +end + +# Utils + +function _to_axes(sz::Dims) + map(sz) do n + r=n÷2 + isodd(n) ? UnitRange(-r, n-r-1) : UnitRange(-r+1, n-r) + end +end +_to_axes(ax::NTuple{N, T}) where {N, T<:AbstractUnitRange} = ax + +end diff --git a/src/kernel.jl b/src/kernel.jl index 5d0993b..aadd21e 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 diff --git a/test/gaborkernels.jl b/test/gaborkernels.jl new file mode 100644 index 0000000..ac15b1e --- /dev/null +++ b/test/gaborkernels.jl @@ -0,0 +1,125 @@ +@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 + +end diff --git a/test/runtests.jl b/test/runtests.jl index cfb5ae5..e4cf368 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ using TestImages using ImageQualityIndexes import StaticArrays using Random +using FFTW @testset "Project meta quality checks" begin # Ambiguity test @@ -43,7 +44,8 @@ include("gradient.jl") include("mapwindow.jl") include("extrema.jl") include("basic.jl") -include("gabor.jl") +# include("gabor.jl") +include("gaborkernels.jl") include("models.jl") From 1e3e2f3c6dff88175330b7afabf78977c07864ce Mon Sep 17 00:00:00 2001 From: JohnnyChen Date: Tue, 2 Nov 2021 05:53:44 +0800 Subject: [PATCH 2/5] introduce log gabor filter: `LogGabor` and `LogGaborComplex` --- docs/demos/filters/gabor_filter.jl | 3 +- docs/demos/filters/loggabor_filter.jl | 107 +++++++++++++++++++ docs/src/function_reference.md | 2 + src/gaborkernels.jl | 143 +++++++++++++++++++++++++- test/gaborkernels.jl | 82 +++++++++++++++ 5 files changed, 335 insertions(+), 2 deletions(-) create mode 100644 docs/demos/filters/loggabor_filter.jl diff --git a/docs/demos/filters/gabor_filter.jl b/docs/demos/filters/gabor_filter.jl index a90f365..dc15c8a 100644 --- a/docs/demos/filters/gabor_filter.jl +++ b/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 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 3edea3c..9ace6f5 100644 --- a/docs/src/function_reference.md +++ b/docs/src/function_reference.md @@ -24,6 +24,8 @@ Kernel.DoG Kernel.LoG Kernel.Laplacian Kernel.Gabor +Kernel.LogGabor +Kernel.LogGaborComplex Kernel.moffat ``` diff --git a/src/gaborkernels.jl b/src/gaborkernels.jl index 46d26e8..b882629 100644 --- a/src/gaborkernels.jl +++ b/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) diff --git a/test/gaborkernels.jl b/test/gaborkernels.jl index ac15b1e..dddda84 100644 --- a/test/gaborkernels.jl +++ b/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 From 2dbafa603bdaa9d8edab6a207ff395cb6d4d72c2 Mon Sep 17 00:00:00 2001 From: JohnnyChen Date: Tue, 2 Nov 2021 06:38:21 +0800 Subject: [PATCH 3/5] deprecate `Kernel.gabor` --- src/kernel.jl | 62 ++------------------------------ src/kernel_deprecated.jl | 62 ++++++++++++++++++++++++++++++++ test/{gabor.jl => deprecated.jl} | 3 -- test/runtests.jl | 5 ++- 4 files changed, 68 insertions(+), 64 deletions(-) create mode 100644 src/kernel_deprecated.jl rename test/{gabor.jl => deprecated.jl} (95%) diff --git a/src/kernel.jl b/src/kernel.jl index aadd21e..ff76a30 100644 --- a/src/kernel.jl +++ b/src/kernel.jl @@ -427,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 @@ -540,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/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/runtests.jl b/test/runtests.jl index e4cf368..a36ad8c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using ImageQualityIndexes import StaticArrays using Random using FFTW +using Statistics @testset "Project meta quality checks" begin # Ambiguity test @@ -44,7 +45,6 @@ include("gradient.jl") include("mapwindow.jl") include("extrema.jl") include("basic.jl") -# include("gabor.jl") include("gaborkernels.jl") include("models.jl") @@ -72,4 +72,7 @@ if CUDA_FUNCTIONAL else @warn "CUDA test: disabled" end + +@info "Beginning deprecation tests, warnings are expected" +include("deprecated.jl") nothing From e0f182abfdf3e6c9c621badcb04d3a75c10b6f05 Mon Sep 17 00:00:00 2001 From: JohnnyChen Date: Wed, 3 Nov 2021 18:48:11 +0800 Subject: [PATCH 4/5] add benchmark group for `Gabor` and `LogGabor` --- benchmark/benchmarks.jl | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) 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 From b1981ea0cec964ad9e4bd39146525fa0f6800910 Mon Sep 17 00:00:00 2001 From: Johnny Chen Date: Wed, 3 Nov 2021 04:40:22 +0800 Subject: [PATCH 5/5] add CUDA testset for `Gabor` and `LogGabor` Currently we don't have very native CUDA support for these two kernels, so the test set mainly serves as "okay, we can convert the lazy kernel array to CuArray by following the steps in tests and then use GPU to do multiplication and fft" --- test/cuda/Project.toml | 2 ++ test/cuda/gaborkernels.jl | 66 +++++++++++++++++++++++++++++++++++++++ test/cuda/runtests.jl | 4 +++ 3 files changed, 72 insertions(+) create mode 100644 test/cuda/gaborkernels.jl 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