diff --git a/Project.toml b/Project.toml index e23edad8..0efab21e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ImageFiltering" uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5" author = ["Tim Holy ", "Jan Weidner "] -version = "0.6.22" +version = "0.7.0" [deps] CatIndices = "aafaddc9-749c-510e-ac4f-586e18779b91" diff --git a/src/ImageFiltering.jl b/src/ImageFiltering.jl index f9703a57..27a68ca5 100644 --- a/src/ImageFiltering.jl +++ b/src/ImageFiltering.jl @@ -18,7 +18,9 @@ export Kernel, KernelFactors, imfilter, imfilter!, mapwindow, mapwindow!, imgradients, padarray, centered, kernelfactors, reflect, - freqkernel, spacekernel + freqkernel, spacekernel, + findlocalminima, findlocalmaxima, + blob_LoG, BlobLoG FixedColorant{T<:Normed} = Colorant{T} StaticOffsetArray{T,N,A<:StaticArray} = OffsetArray{T,N,A} @@ -56,6 +58,7 @@ using .Algorithm: Alg, FFT, FIR, FIRTiled, IIR, Mixed Alg(r::AbstractResource{A}) where {A<:Alg} = r.settings include("utils.jl") +include("compat.jl") include("kernelfactors.jl") using .KernelFactors: TriggsSdika, IIRFilter, ReshapedOneD, iterdims, kernelfactors @@ -85,6 +88,7 @@ include("specialty.jl") include("mapwindow.jl") using .MapWindow +include("extrema.jl") function __init__() # See ComputationalResources README for explanation diff --git a/src/compat.jl b/src/compat.jl new file mode 100644 index 00000000..d22a4480 --- /dev/null +++ b/src/compat.jl @@ -0,0 +1,11 @@ +if VERSION <= v"1.0.5" + # https://github.com/JuliaLang/julia/pull/29442 + _oneunit(::CartesianIndex{N}) where {N} = _oneunit(CartesianIndex{N}) + _oneunit(::Type{CartesianIndex{N}}) where {N} = CartesianIndex(ntuple(x -> 1, Val(N))) +else + const _oneunit = Base.oneunit +end + +# Equivalent to I:J on later Julia versions +_colon(I::CartesianIndex{N}, J::CartesianIndex{N}) where N = + CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J))) diff --git a/src/extrema.jl b/src/extrema.jl new file mode 100644 index 00000000..0eb02b5e --- /dev/null +++ b/src/extrema.jl @@ -0,0 +1,166 @@ +""" +BlobLoG stores information about the location of peaks as discovered by `blob_LoG`. +It has fields: + +- location: the location of a peak in the filtered image (a CartesianIndex) +- σ: the value of σ which lead to the largest `-LoG`-filtered amplitude at this location +- amplitude: the value of the `-LoG(σ)`-filtered image at the peak + +Note that the radius is equal to σ√2. + +See also: [`blob_LoG`](@ref). +""" +struct BlobLoG{T,S,N} + location::CartesianIndex{N} + σ::S + amplitude::T +end +BlobLoG(; location, σ, amplitude) = BlobLoG(location, σ, amplitude) + +function Base.show(io::IO, bl::BlobLoG) + print(io, "BlobLoG(location=", bl.location, ", σ=", bl.σ, ", amplitude=", bl.amplitude, ")") +end + + +""" + blob_LoG(img, σscales; edges::Tuple=(true, false, ...), σshape::Tuple=(1, ...), rthresh=0.001) -> Vector{BlobLoG} + +Find "blobs" in an N-D image using the negative Lapacian of Gaussians +with the specifed vector or tuple of σ values. The algorithm searches for places +where the filtered image (for a particular σ) is at a peak compared to all +spatially- and σ-adjacent voxels, where σ is `σscales[i] * σshape` for some i. +By default, `σshape` is an ntuple of 1s. + +The optional `edges` argument controls whether peaks on the edges are +included. `edges` can be `true` or `false`, or a N+1-tuple in which +the first entry controls whether edge-σ values are eligible to serve +as peaks, and the remaining N entries control each of the N dimensions +of `img`. + +`rthresh` controls the minimum amplitude of peaks in the -LoG-filtered image, +as a fraction of `maximum(abs, img)` and the volume of the Gaussian. + +# Examples + +While most images are 2- or 3-dimensional, it will be easier to illustrate this with +a one-dimensional "image" containing two Gaussian blobs of different sizes: + +```jldoctest; setup=:(using ImageFiltering), filter=r"amplitude=.*"] +julia> σs = 2.0.^(1:6); + +julia> img = zeros(100); img[20:30] = [exp(-x^2/(2*4^2)) for x=-5:5]; img[50:80] = [exp(-x^2/(2*8^2)) for x=-15:15]; + +julia> blob_LoG(img, σs; edges=false) +2-element Vector{BlobLoG{Float64, Tuple{Float64}, 1}}: + location=CartesianIndex(25,), σ=(4.0,), amplitude=0.10453155018303673 + location=CartesianIndex(65,), σ=(8.0,), amplitude=0.046175719034527364 +``` + +The other two are centered in their corresponding "features," and the width `σ` +reflects the width of the feature itself. + +`blob_LoG` tends to work best for shapes that are "Gaussian-like" but does +generalize somewhat. + +# Citation: + +Lindeberg T (1998), "Feature Detection with Automatic Scale Selection", +International Journal of Computer Vision, 30(2), 79–116. + +See also: [`BlobLoG`](@ref). +""" +function blob_LoG(img::AbstractArray{T,N}, σscales; + edges::Union{Bool,Tuple{Bool,Vararg{Bool,N}}}=(true, ntuple(d->false, Val(N))...), + σshape::NTuple{N,Real}=ntuple(d->1, Val(N)), + rthresh::Real=1//1000) where {T<:Union{AbstractGray,Real},N} + if edges isa Bool + edges = (edges, ntuple(d->edges,Val(N))...) + end + sigmas = sort(σscales) + img_LoG = multiLoG(img, sigmas, σshape) + maxima = findlocalmaxima(img_LoG; edges=edges) + # The "density" should not be much smaller than 1/volume of the Gaussian + if !iszero(rthresh) + athresh = rthresh./(sigmas.^N .* prod(σshape)) + imgmax = maximum(abs, img) + [BlobLoG(CartesianIndex(tail(x.I)), sigmas[x[1]].*σshape, img_LoG[x]) for x in maxima if img_LoG[x] > athresh[x[1]]*imgmax] + else + [BlobLoG(CartesianIndex(tail(x.I)), sigmas[x[1]].*σshape, img_LoG[x]) for x in maxima] + end +end + +function multiLoG(img::AbstractArray{T,N}, sigmas, σshape) where {T,N} + issorted(sigmas) || error("sigmas must be sorted") + img_LoG = similar(img, float(eltype(T)), (Base.OneTo(length(sigmas)), axes(img)...)) + colons = ntuple(d->Colon(), Val(N)) + @inbounds for (isigma, σ) in enumerate(sigmas) + LoG_slice = @view img_LoG[isigma, colons...] + imfilter!(LoG_slice, img, Kernel.LoG(ntuple(i->σ*σshape[i], Val(N))), "reflect") + LoG_slice .*= -σ + end + return img_LoG +end + +default_window(img) = (cs = coords_spatial(img); ntuple(d -> d ∈ cs ? 3 : 1, ndims(img))) + +""" + findlocalmaxima(img; window=default_window(img), edges=true) -> Vector{CartesianIndex} + +Returns the coordinates of elements whose value is larger than all of +their immediate neighbors. `edges` is a boolean specifying whether to include the +first and last elements of each dimension, or a tuple-of-Bool +specifying edge behavior for each dimension separately. + +The `default_window` is 3 for each spatial dimension of `img`, and 1 otherwise, implying +that maxima are detected over nearest-neighbors in each spatial "slice" by default. +""" +findlocalmaxima(img::AbstractArray; window=default_window(img), edges=true) = + findlocalextrema(>, img, window, edges) + +""" + findlocalminima(img; window=default_window(img), edges=true) -> Vector{CartesianIndex} + +Like [`findlocalmaxima`](@ref), but returns the coordinates of the smallest elements. +""" +findlocalminima(img::AbstractArray; window=default_window(img), edges=true) = + findlocalextrema(<, img, window, edges) + + +findlocalextrema(f, img::AbstractArray{T,N}, window, edges::Bool) where {T,N} = findlocalextrema(f, img, window, ntuple(d->edges,Val(N))) + +function findlocalextrema(f::F, img::AbstractArray{T,N}, window::Dims{N}, edges::NTuple{N,Bool}) where {F,T<:Union{Gray,Number},N} + extrema = Vector{CartesianIndex{N}}(undef, 0) + Iedge = CartesianIndex(map(!, edges)) + R0 = CartesianIndices(img) + R = clippedinds(R0, Iedge) + halfwindow = CartesianIndex(map(x -> x >> 1, window)) + Rinterior = clippedinds(R0, halfwindow) + Rwindow = _colon(-halfwindow, halfwindow) + z = zero(halfwindow) + for i in R + isextrema = true + img_i = img[i] + if i ∈ Rinterior + # If i is in the interior, we don't have to worry about i+j being out-of-bounds + for j in Rwindow + j == z && continue + if !f(img_i, img[i+j]) + isextrema = false + break + end + end + else + for j in Rwindow + (j == z || i+j ∉ R0) && continue + if !f(img_i, img[i+j]) + isextrema = false + break + end + end + end + isextrema && push!(extrema, i) + end + extrema +end + +clippedinds(Router, Iclip) = _colon(first(Router)+Iclip, last(Router)-Iclip) diff --git a/test/extrema.jl b/test/extrema.jl new file mode 100644 index 00000000..92bb20ed --- /dev/null +++ b/test/extrema.jl @@ -0,0 +1,59 @@ +@testset "extrema" begin + @testset "local extrema" begin + A = zeros(Int, 9, 9); A[[1:2;5],5].=1 + @test findlocalmaxima(A) == [CartesianIndex((5,5))] + @test findlocalmaxima(A; window=(1,3)) == [CartesianIndex((1,5)),CartesianIndex((2,5)),CartesianIndex((5,5))] + @test findlocalmaxima(A; window=(1,3), edges=false) == [CartesianIndex((2,5)),CartesianIndex((5,5))] + A = zeros(Int, 9, 9, 9); A[[1:2;5],5,5].=1 + @test findlocalmaxima(A) == [CartesianIndex((5,5,5))] + @test findlocalmaxima(A; window=(1,3,1)) == [CartesianIndex((1,5,5)),CartesianIndex((2,5,5)),CartesianIndex((5,5,5))] + @test findlocalmaxima(A, window=(1,3,1), edges=false) == [CartesianIndex((2,5,5)),CartesianIndex((5,5,5))] + A = zeros(Int, 9, 9); A[[1:2;5],5].=-1 + @test findlocalminima(A) == [CartesianIndex((5,5))] + end + + @testset "blob_LoG" begin + A = zeros(Int, 9, 9); A[5, 5] = 1 + blobs = blob_LoG(A, 2.0.^[0.5,0,1]) + @test length(blobs) == 1 + blob = blobs[1] + @test blob.amplitude ≈ 0.3183098861837907 + @test blob.σ === (1.0, 1.0) + @test blob.location == CartesianIndex((5,5)) + str = sprint(print, blob) + @test occursin("σ=$((1.0, 1.0))", str) + @test eval(Meta.parse(str)) == blob + @test blob_LoG(A, [1.0]) == blobs + @test blob_LoG(A, [1.0]; edges=(true, false, false)) == blobs + @test isempty(blob_LoG(A, [1.0]; edges=false)) + A = zeros(Int, 9, 9); A[1, 5] = 1 + blobs = blob_LoG(A, 2.0.^[0,0.5,1]) + A = zeros(Int, 9, 9); A[1,5] = 1 + blobs = blob_LoG(A, 2.0.^[0.5,0,1]) + @test all(b.amplitude < 1e-16 for b in blobs) + blobs = filter(b->b.amplitude > 0.1, blob_LoG(A, 2.0.^[0.5,0,1]; edges=true)) + @test length(blobs) == 1 + @test blobs[1].location == CartesianIndex((1,5)) + @test filter(b->b.amplitude > 0.1, blob_LoG(A, 2.0.^[0.5,0,1], edges=(true, true, false))) == blobs + @test isempty(blob_LoG(A, 2.0.^[0,1], edges=(false, true, false))) + blobs = blob_LoG(A, 2.0.^[0,0.5,1], edges=(true, false, true)) + @test all(b.amplitude < 1e-16 for b in blobs) + # stub test for N-dimensional blob_LoG: + A = zeros(Int, 9, 9, 9); A[5, 5, 5] = 1 + blobs = blob_LoG(A, 2.0.^[0.5, 0, 1]) + @test length(blobs) == 1 + @test blobs[1].location == CartesianIndex((5,5,5)) + # kinda anisotropic image + A = zeros(Int,9,9,9); A[5,4:6,5] .= 1; + blobs = blob_LoG(A,2 .^ [1.,0,0.5], σshape=(1.,3.,1.)) + @test length(blobs) == 1 + @test blobs[1].location == CartesianIndex((5,5,5)) + A = zeros(Int,9,9,9); A[1,1,4:6] .= 1; + blobs = filter(b->b.amplitude > 0.1, blob_LoG(A, 2.0.^[0.5,0,1], edges=true, σshape=(1.,1.,3.))) + @test length(blobs) == 1 + @test blobs[1].location == CartesianIndex((1,1,5)) + @test filter(b->b.amplitude > 0.1, blob_LoG(A, 2.0.^[0.5,0,1], edges=(true, true, true, false), σshape=(1.,1.,3.))) == blobs + @test isempty(blob_LoG(A, 2.0.^[0,1], edges=(false, true, false, false), σshape=(1.,1.,3.))) + @test length(blob_LoG([zeros(10); 1.0; 0.0], [4]; edges=true, rthresh=0)) > length(blob_LoG([zeros(10); 1.0; 0.0], [4]; edges=true)) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 29d6f7c9..fd407e9c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,6 +28,7 @@ include("cascade.jl") include("specialty.jl") include("gradient.jl") include("mapwindow.jl") +include("extrema.jl") include("basic.jl") include("gabor.jl")