|
| 1 | +""" |
| 2 | +BlobLoG stores information about the location of peaks as discovered by `blob_LoG`. |
| 3 | +It has fields: |
| 4 | +
|
| 5 | +- location: the location of a peak in the filtered image (a CartesianIndex) |
| 6 | +- σ: the value of σ which lead to the largest `-LoG`-filtered amplitude at this location |
| 7 | +- amplitude: the value of the `-LoG(σ)`-filtered image at the peak |
| 8 | +
|
| 9 | +Note that the radius is equal to σ√2. |
| 10 | +
|
| 11 | +See also: [`blob_LoG`](@ref). |
| 12 | +""" |
| 13 | +struct BlobLoG{T,S,N} |
| 14 | + location::CartesianIndex{N} |
| 15 | + σ::S |
| 16 | + amplitude::T |
| 17 | +end |
| 18 | +BlobLoG(; location, σ, amplitude) = BlobLoG(location, σ, amplitude) |
| 19 | + |
| 20 | +function Base.show(io::IO, bl::BlobLoG) |
| 21 | + print(io, "BlobLoG(location=", bl.location, ", σ=", bl.σ, ", amplitude=", bl.amplitude, ")") |
| 22 | +end |
| 23 | + |
| 24 | + |
| 25 | +""" |
| 26 | + blob_LoG(img, σscales; edges::Tuple=(true, false, ...), σshape::Tuple=(1, ...), rthresh=0.001) -> Vector{BlobLoG} |
| 27 | +
|
| 28 | +Find "blobs" in an N-D image using the negative Lapacian of Gaussians |
| 29 | +with the specifed vector or tuple of σ values. The algorithm searches for places |
| 30 | +where the filtered image (for a particular σ) is at a peak compared to all |
| 31 | +spatially- and σ-adjacent voxels, where σ is `σscales[i] * σshape` for some i. |
| 32 | +By default, `σshape` is an ntuple of 1s. |
| 33 | +
|
| 34 | +The optional `edges` argument controls whether peaks on the edges are |
| 35 | +included. `edges` can be `true` or `false`, or a N+1-tuple in which |
| 36 | +the first entry controls whether edge-σ values are eligible to serve |
| 37 | +as peaks, and the remaining N entries control each of the N dimensions |
| 38 | +of `img`. |
| 39 | +
|
| 40 | +`rthresh` controls the minimum amplitude of peaks in the -LoG-filtered image, |
| 41 | +as a fraction of `maximum(abs, img)` and the volume of the Gaussian. |
| 42 | +
|
| 43 | +# Examples |
| 44 | +
|
| 45 | +While most images are 2- or 3-dimensional, it will be easier to illustrate this with |
| 46 | +a one-dimensional "image" containing two Gaussian blobs of different sizes: |
| 47 | +
|
| 48 | +```jldoctest; setup=:(using ImageFiltering), filter=r"amplitude=.*"] |
| 49 | +julia> σs = 2.0.^(1:6); |
| 50 | +
|
| 51 | +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]; |
| 52 | +
|
| 53 | +julia> blob_LoG(img, σs; edges=false) |
| 54 | +2-element Vector{BlobLoG{Float64, Tuple{Float64}, 1}}: |
| 55 | + location=CartesianIndex(25,), σ=(4.0,), amplitude=0.10453155018303673 |
| 56 | + location=CartesianIndex(65,), σ=(8.0,), amplitude=0.046175719034527364 |
| 57 | +``` |
| 58 | +
|
| 59 | +The other two are centered in their corresponding "features," and the width `σ` |
| 60 | +reflects the width of the feature itself. |
| 61 | +
|
| 62 | +`blob_LoG` tends to work best for shapes that are "Gaussian-like" but does |
| 63 | +generalize somewhat. |
| 64 | +
|
| 65 | +# Citation: |
| 66 | +
|
| 67 | +Lindeberg T (1998), "Feature Detection with Automatic Scale Selection", |
| 68 | +International Journal of Computer Vision, 30(2), 79–116. |
| 69 | +
|
| 70 | +See also: [`BlobLoG`](@ref). |
| 71 | +""" |
| 72 | +function blob_LoG(img::AbstractArray{T,N}, σscales; |
| 73 | + edges::Union{Bool,Tuple{Bool,Vararg{Bool,N}}}=(true, ntuple(d->false, Val(N))...), |
| 74 | + σshape::NTuple{N,Real}=ntuple(d->1, Val(N)), |
| 75 | + rthresh::Real=1//1000) where {T<:Union{AbstractGray,Real},N} |
| 76 | + if edges isa Bool |
| 77 | + edges = (edges, ntuple(d->edges,Val(N))...) |
| 78 | + end |
| 79 | + sigmas = sort(σscales) |
| 80 | + img_LoG = multiLoG(img, sigmas, σshape) |
| 81 | + maxima = findlocalmaxima(img_LoG; edges=edges) |
| 82 | + # The "density" should not be much smaller than 1/volume of the Gaussian |
| 83 | + if !iszero(rthresh) |
| 84 | + athresh = rthresh./(sigmas.^N .* prod(σshape)) |
| 85 | + imgmax = maximum(abs, img) |
| 86 | + [BlobLoG(CartesianIndex(tail(x.I)), sigmas[x[1]].*σshape, img_LoG[x]) for x in maxima if img_LoG[x] > athresh[x[1]]*imgmax] |
| 87 | + else |
| 88 | + [BlobLoG(CartesianIndex(tail(x.I)), sigmas[x[1]].*σshape, img_LoG[x]) for x in maxima] |
| 89 | + end |
| 90 | +end |
| 91 | + |
| 92 | +function multiLoG(img::AbstractArray{T,N}, sigmas, σshape) where {T,N} |
| 93 | + issorted(sigmas) || error("sigmas must be sorted") |
| 94 | + img_LoG = similar(img, float(eltype(T)), (Base.OneTo(length(sigmas)), axes(img)...)) |
| 95 | + colons = ntuple(d->Colon(), Val(N)) |
| 96 | + @inbounds for (isigma, σ) in enumerate(sigmas) |
| 97 | + LoG_slice = @view img_LoG[isigma, colons...] |
| 98 | + imfilter!(LoG_slice, img, Kernel.LoG(ntuple(i->σ*σshape[i], Val(N))), "reflect") |
| 99 | + LoG_slice .*= -σ |
| 100 | + end |
| 101 | + return img_LoG |
| 102 | +end |
| 103 | + |
| 104 | +default_window(img) = (cs = coords_spatial(img); ntuple(d -> d ∈ cs ? 3 : 1, ndims(img))) |
| 105 | + |
| 106 | +""" |
| 107 | + findlocalmaxima(img; window=default_window(img), edges=true) -> Vector{CartesianIndex} |
| 108 | +
|
| 109 | +Returns the coordinates of elements whose value is larger than all of |
| 110 | +their immediate neighbors. `edges` is a boolean specifying whether to include the |
| 111 | +first and last elements of each dimension, or a tuple-of-Bool |
| 112 | +specifying edge behavior for each dimension separately. |
| 113 | +
|
| 114 | +The `default_window` is 3 for each spatial dimension of `img`, and 1 otherwise, implying |
| 115 | +that maxima are detected over nearest-neighbors in each spatial "slice" by default. |
| 116 | +""" |
| 117 | +findlocalmaxima(img::AbstractArray; window=default_window(img), edges=true) = |
| 118 | + findlocalextrema(>, img, window, edges) |
| 119 | + |
| 120 | +""" |
| 121 | + findlocalminima(img; window=default_window(img), edges=true) -> Vector{CartesianIndex} |
| 122 | +
|
| 123 | +Like [`findlocalmaxima`](@ref), but returns the coordinates of the smallest elements. |
| 124 | +""" |
| 125 | +findlocalminima(img::AbstractArray; window=default_window(img), edges=true) = |
| 126 | + findlocalextrema(<, img, window, edges) |
| 127 | + |
| 128 | + |
| 129 | +findlocalextrema(f, img::AbstractArray{T,N}, window, edges::Bool) where {T,N} = findlocalextrema(f, img, window, ntuple(d->edges,Val(N))) |
| 130 | + |
| 131 | +function findlocalextrema(f::F, img::AbstractArray{T,N}, window::Dims{N}, edges::NTuple{N,Bool}) where {F,T<:Union{Gray,Number},N} |
| 132 | + extrema = Vector{CartesianIndex{N}}(undef, 0) |
| 133 | + Iedge = CartesianIndex(map(!, edges)) |
| 134 | + R0 = CartesianIndices(img) |
| 135 | + R = clippedinds(R0, Iedge) |
| 136 | + halfwindow = CartesianIndex(map(x -> x >> 1, window)) |
| 137 | + Rinterior = clippedinds(R0, halfwindow) |
| 138 | + Rwindow = _colon(-halfwindow, halfwindow) |
| 139 | + z = zero(halfwindow) |
| 140 | + for i in R |
| 141 | + isextrema = true |
| 142 | + img_i = img[i] |
| 143 | + if i ∈ Rinterior |
| 144 | + # If i is in the interior, we don't have to worry about i+j being out-of-bounds |
| 145 | + for j in Rwindow |
| 146 | + j == z && continue |
| 147 | + if !f(img_i, img[i+j]) |
| 148 | + isextrema = false |
| 149 | + break |
| 150 | + end |
| 151 | + end |
| 152 | + else |
| 153 | + for j in Rwindow |
| 154 | + (j == z || i+j ∉ R0) && continue |
| 155 | + if !f(img_i, img[i+j]) |
| 156 | + isextrema = false |
| 157 | + break |
| 158 | + end |
| 159 | + end |
| 160 | + end |
| 161 | + isextrema && push!(extrema, i) |
| 162 | + end |
| 163 | + extrema |
| 164 | +end |
| 165 | + |
| 166 | +clippedinds(Router, Iclip) = _colon(first(Router)+Iclip, last(Router)-Iclip) |
0 commit comments