Move extremization, blob_LoG from Images

Project.toml

@@ -1,7 +1,7 @@
 name = "ImageFiltering"
 uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 author = ["Tim Holy <tim.holy@gmail.com>", "Jan Weidner <jw3126@gmail.com>"]
-version = "0.6.22"
+version = "0.7.0"
 
 [deps]
 CatIndices = "aafaddc9-749c-510e-ac4f-586e18779b91"

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

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)))

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)

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

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")