diff --git a/src/RegisterMismatchCommon.jl b/src/RegisterMismatchCommon.jl index 1d76ae0..3ea2770 100644 --- a/src/RegisterMismatchCommon.jl +++ b/src/RegisterMismatchCommon.jl @@ -8,21 +8,21 @@ export padranges, shiftrange, checksize_maxshift, register_translate const DimsLike = Union{AbstractVector{Int}, Dims} # try to avoid these and just use Dims tuples for sizes -const WidthLike = Union{AbstractVector,Tuple} -FFTPROD = [2,3] +const WidthLike = Union{AbstractVector, Tuple} +FFTPROD = [2, 3] set_FFTPROD(v) = global FFTPROD = v -mismatch(fixed::AbstractArray{T}, moving::AbstractArray{T}, maxshift::DimsLike; normalization = :intensity) where {T<:AbstractFloat} = mismatch(T, fixed, moving, maxshift; normalization=normalization) -mismatch(fixed::AbstractArray, moving::AbstractArray, maxshift::DimsLike; normalization = :intensity) = mismatch(Float32, fixed, moving, maxshift; normalization=normalization) +mismatch(fixed::AbstractArray{T}, moving::AbstractArray{T}, maxshift::DimsLike; normalization = :intensity) where {T <: AbstractFloat} = mismatch(T, fixed, moving, maxshift; normalization = normalization) +mismatch(fixed::AbstractArray, moving::AbstractArray, maxshift::DimsLike; normalization = :intensity) = mismatch(Float32, fixed, moving, maxshift; normalization = normalization) -mismatch_apertures(fixed::AbstractArray{T}, moving::AbstractArray{T}, args...; kwargs...) where {T<:AbstractFloat} = mismatch_apertures(T, fixed, moving, args...; kwargs...) +mismatch_apertures(fixed::AbstractArray{T}, moving::AbstractArray{T}, args...; kwargs...) where {T <: AbstractFloat} = mismatch_apertures(T, fixed, moving, args...; kwargs...) mismatch_apertures(fixed::AbstractArray, moving::AbstractArray, args...; kwargs...) = mismatch_apertures(Float32, fixed, moving, args...; kwargs...) -function mismatch_apertures(::Type{T}, fixed::AbstractArray, moving::AbstractArray, gridsize::DimsLike, maxshift::DimsLike; kwargs...) where T +function mismatch_apertures(::Type{T}, fixed::AbstractArray, moving::AbstractArray, gridsize::DimsLike, maxshift::DimsLike; kwargs...) where {T} cs = coords_spatial(fixed) - aperture_centers = aperture_grid(map(d->size(fixed, d), cs), gridsize) + aperture_centers = aperture_grid(map(d -> size(fixed, d), cs), gridsize) aperture_width = default_aperture_width(fixed, gridsize) - mismatch_apertures(T, fixed, moving, aperture_centers, aperture_width, maxshift; kwargs...) + return mismatch_apertures(T, fixed, moving, aperture_centers, aperture_width, maxshift; kwargs...) end """ @@ -36,52 +36,52 @@ whenever `i` or `j` is zero. Data are imputed by averaging the adjacent non-suspect values. This function works in-place, overwriting the original `mm`. """ -function correctbias!(mm::MismatchArray{ND,N}, w = correctbias_weight(mm)) where {ND,N} +function correctbias!(mm::MismatchArray{ND, N}, w = correctbias_weight(mm)) where {ND, N} T = eltype(ND) mxshift = maxshift(mm) Imax = CartesianIndex(mxshift) - Imin = CartesianIndex(map(x->-x,mxshift)::NTuple{N,Int}) - I1 = CartesianIndex(ntuple(d->d>2 ? 0 : 1, N)::NTuple{N,Int}) # only first 2 dims + Imin = CartesianIndex(map(x -> -x, mxshift)::NTuple{N, Int}) + I1 = CartesianIndex(ntuple(d -> d > 2 ? 0 : 1, N)::NTuple{N, Int}) # only first 2 dims for I in eachindex(mm) if w[I] == 0 - mms = NumDenom{T}(0,0) + mms = NumDenom{T}(0, 0) ws = zero(T) - strt = max(Imin, I-I1) - stop = min(Imax, I+I1) - for J in CartesianIndices(ntuple(d->strt[d]:stop[d],N)) + strt = max(Imin, I - I1) + stop = min(Imax, I + I1) + for J in CartesianIndices(ntuple(d -> strt[d]:stop[d], N)) wJ = w[J] if wJ != 0 - mms += wJ*mm[J] + mms += wJ * mm[J] ws += wJ end end - mm[I] = mms/ws + mm[I] = mms / ws end end - mm + return mm end "`correctbias!(mms)` runs `correctbias!` on each element of an array-of-MismatchArrays." -function correctbias!(mms::AbstractArray{M}) where M<:MismatchArray +function correctbias!(mms::AbstractArray{M}) where {M <: MismatchArray} for mm in mms correctbias!(mm) end - mms + return mms end -function correctbias_weight(mm::MismatchArray{ND,N}) where {ND,N} +function correctbias_weight(mm::MismatchArray{ND, N}) where {ND, N} T = eltype(ND) w = CenterIndexedArray(ones(T, size(mm))) for I in eachindex(mm) anyzero = false - for d = 1:min(N,2) # only first 2 dims + for d in 1:min(N, 2) # only first 2 dims anyzero |= I[d] == 0 end if anyzero w[I] = 0 end end - w + return w end """ @@ -94,12 +94,12 @@ function nanpad(fixed, moving) if size(fixed) == size(moving) return fixed, moving end - rng = map(d->1:max(size(fixed,d), size(moving,d)), 1:ndims(fixed)) + rng = map(d -> 1:max(size(fixed, d), size(moving, d)), 1:ndims(fixed)) T = promote_type(eltype(fixed), eltype(moving)) - get(fixed, rng, nanval(T)), get(moving, rng, nanval(T)) + return get(fixed, rng, nanval(T)), get(moving, rng, nanval(T)) end -nanval(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN) +nanval(::Type{T}) where {T <: AbstractFloat} = convert(T, NaN) nanval(::Type{T}) where {T} = convert(Float32, NaN) """ @@ -107,18 +107,18 @@ nanval(::Type{T}) where {T} = convert(Float32, NaN) "as-is" mismatch between `fixed` and `moving`, without any shift. `normalization` may be either `:intensity` (the default) or `:pixels`. """ -function mismatch0(fixed::AbstractArray{Tf,N}, moving::AbstractArray{Tm,N}; normalization = :intensity) where {Tf,Tm,N} +function mismatch0(fixed::AbstractArray{Tf, N}, moving::AbstractArray{Tm, N}; normalization = :intensity) where {Tf, Tm, N} size(fixed) == size(moving) || throw(DimensionMismatch("Size $(size(fixed)) of fixed is not equal to size $(size(moving)) of moving")) - _mismatch0(zero(Float64), zero(Float64), fixed, moving; normalization=normalization) + return _mismatch0(zero(Float64), zero(Float64), fixed, moving; normalization = normalization) end -function _mismatch0(num::T, denom::T, fixed::AbstractArray{Tf,N}, moving::AbstractArray{Tm,N}; normalization = :intensity) where {T,Tf,Tm,N} +function _mismatch0(num::T, denom::T, fixed::AbstractArray{Tf, N}, moving::AbstractArray{Tm, N}; normalization = :intensity) where {T, Tf, Tm, N} if normalization == :intensity for i in eachindex(fixed, moving) vf = T(fixed[i]) vm = T(moving[i]) if isfinite(vf) && isfinite(vm) - num += (vf-vm)^2 + num += (vf - vm)^2 denom += vf^2 + vm^2 end end @@ -127,14 +127,14 @@ function _mismatch0(num::T, denom::T, fixed::AbstractArray{Tf,N}, moving::Abstra vf = T(fixed[i]) vm = T(moving[i]) if isfinite(vf) && isfinite(vm) - num += (vf-vm)^2 + num += (vf - vm)^2 denom += 1 end end else error("Normalization $normalization not recognized") end - NumDenom(num, denom) + return NumDenom(num, denom) end """ @@ -143,14 +143,14 @@ mismatch between `fixed` and `moving`, without any shift. The mismatch is represented in `mms` as an aperture-wise Arrays-of-MismatchArrays. """ -function mismatch0(mms::AbstractArray{M}) where M<:MismatchArray +function mismatch0(mms::AbstractArray{M}) where {M <: MismatchArray} mm0 = eltype(M)(0, 0) cr = eachindex(first(mms)) - z = first(cr)+last(cr) # all-zeros CartesianIndex + z = first(cr) + last(cr) # all-zeros CartesianIndex for mm in mms mm0 += mm[z] end - mm0 + return mm0 end """ @@ -159,19 +159,19 @@ grid of aperture centers. The grid has size `gridsize`, and is constructed for an image of spatial size `ssize`. Along each dimension the first and last elements are at the image corners. """ -function aperture_grid(ssize::Dims{N}, gridsize) where N +function aperture_grid(ssize::Dims{N}, gridsize) where {N} if length(gridsize) != N - if length(gridsize) == N-1 + if length(gridsize) == N - 1 @info("ssize and gridsize disagree; possible fix is to use a :time axis (AxisArrays) for the image") end error("ssize and gridsize must have the same length, got $ssize and $gridsize") end - grid = Array{NTuple{N,Float64},N}(undef, (gridsize...,)) - centers = map(i-> gridsize[i] > 1 ? collect(range(1, stop=ssize[i], length=gridsize[i])) : [(ssize[i]+1)/2], 1:N) + grid = Array{NTuple{N, Float64}, N}(undef, (gridsize...,)) + centers = map(i -> gridsize[i] > 1 ? collect(range(1, stop = ssize[i], length = gridsize[i])) : [(ssize[i] + 1) / 2], 1:N) for I in CartesianIndices(size(grid)) - grid[I] = ntuple(i->centers[i][I[i]], N) + grid[I] = ntuple(i -> centers[i][I[i]], N) end - grid + return grid end """ @@ -190,10 +190,10 @@ array-of-vectors) or an `N+1`-dimensional array with the center positions specified along the first dimension. (But you may find the `gridsize` syntax to be simpler.) """ -function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{C}, maxshift) where {T,C<:Union{AbstractVector,Tuple}} +function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{C}, maxshift) where {T, C <: Union{AbstractVector, Tuple}} isempty(aperture_centers) && error("aperture_centers is empty") N = length(first(aperture_centers)) - sz = map(x->2*x+1, maxshift) + sz = map(x -> 2 * x + 1, maxshift) mm = MismatchArray(T, sz...) mms = Array{typeof(mm)}(undef, size(aperture_centers)) f = true @@ -205,26 +205,26 @@ function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{C}, maxshi mms[i] = MismatchArray(T, sz...) end end - mms + return mms end -function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{R}, maxshift) where {T,R<:Real} - N = ndims(aperture_centers)-1 - mms = Array{MismatchArray{T,N}}(undef, size(aperture_centers)[2:end]) - sz = map(x->2*x+1, maxshift) +function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{R}, maxshift) where {T, R <: Real} + N = ndims(aperture_centers) - 1 + mms = Array{MismatchArray{T, N}}(undef, size(aperture_centers)[2:end]) + sz = map(x -> 2 * x + 1, maxshift) for i in eachindex(mms) mms[i] = MismatchArray(T, sz...) end - mms + return mms end -function allocate_mmarrays(::Type{T}, gridsize::NTuple{N,Int}, maxshift) where {T<:Real,N} - mms = Array{MismatchArray{NumDenom{T},N}}(undef, gridsize) - sz = map(x->2*x+1, maxshift) +function allocate_mmarrays(::Type{T}, gridsize::NTuple{N, Int}, maxshift) where {T <: Real, N} + mms = Array{MismatchArray{NumDenom{T}, N}}(undef, gridsize) + sz = map(x -> 2 * x + 1, maxshift) for i in eachindex(mms) mms[i] = MismatchArray(T, sz...) end - mms + return mms end struct ContainerIterator{C} @@ -234,22 +234,22 @@ end Base.iterate(iter::ContainerIterator) = iterate(iter.data) Base.iterate(iter::ContainerIterator, state) = iterate(iter.data, state) -struct FirstDimIterator{A<:AbstractArray,R<:CartesianIndices} +struct FirstDimIterator{A <: AbstractArray, R <: CartesianIndices} data::A rng::R - FirstDimIterator{A,R}(data::A) where {A,R} = new{A,R}(data, CartesianIndices(Base.tail(size(data)))) + FirstDimIterator{A, R}(data::A) where {A, R} = new{A, R}(data, CartesianIndices(Base.tail(size(data)))) end -FirstDimIterator(A::AbstractArray) = FirstDimIterator{typeof(A),typeof(CartesianIndices(Base.tail(size(A))))}(A) +FirstDimIterator(A::AbstractArray) = FirstDimIterator{typeof(A), typeof(CartesianIndices(Base.tail(size(A))))}(A) function Base.iterate(iter::FirstDimIterator) isempty(iter.rng) && return nothing index, state = iterate(iter.rng) - iter.data[:, index], state + return iter.data[:, index], state end function Base.iterate(iter::FirstDimIterator, state) state == last(iter.rng) && return nothing index, state = iterate(iter.rng, state) - iter.data[:, index], state + return iter.data[:, index], state end """ @@ -259,9 +259,9 @@ AbstractArray-of-tuples or -AbstractVectors, or may be an `AbstractArray` where each point is represented along the first dimension (e.g., columns of a matrix). """ -each_point(aperture_centers::AbstractArray{C}) where {C<:Union{AbstractVector,Tuple}} = ContainerIterator(aperture_centers) +each_point(aperture_centers::AbstractArray{C}) where {C <: Union{AbstractVector, Tuple}} = ContainerIterator(aperture_centers) -each_point(aperture_centers::AbstractArray{R}) where {R<:Real} = FirstDimIterator(aperture_centers) +each_point(aperture_centers::AbstractArray{R}) where {R <: Real} = FirstDimIterator(aperture_centers) """ `rng = aperture_range(center, width)` returns a tuple of @@ -270,7 +270,7 @@ and has width `width[d]`. """ function aperture_range(center, width) length(center) == length(width) || error("center and width must have the same length") - ntuple(d->leftedge(center[d], width[d]):rightedge(center[d], width[d]), length(center)) + return ntuple(d -> leftedge(center[d], width[d]):rightedge(center[d], width[d]), length(center)) end """ @@ -284,34 +284,34 @@ collection of apertures will yield full coverage of the image. function default_aperture_width(img, gridsize::DimsLike, overlap::DimsLike = zeros(Int, sdims(img))) sc = coords_spatial(img) length(sc) == length(gridsize) == length(overlap) || error("gridsize and overlap must have length equal to the number of spatial dimensions in img") - for i = 1:length(sc) + for i in 1:length(sc) if gridsize[i] > size(img, sc[i]) error("gridsize $gridsize is too large, given the size $(size(img)[sc]) of the image") end end - gsz1 = max.(1, [gridsize...].-1) - gflag = [gridsize...].>1 - tuple((([map(d->size(img,d),sc)...]-gflag)./gsz1+2*[overlap...].*gflag)...) + gsz1 = max.(1, [gridsize...] .- 1) + gflag = [gridsize...] .> 1 + return tuple((([map(d -> size(img, d), sc)...] - gflag) ./ gsz1 + 2 * [overlap...] .* gflag)...) end """ `truncatenoise!(mm, thresh)` zeros out any entries of the MismatchArray `mm` whose `denom` values are less than `thresh`. """ -function truncatenoise!(mm::AbstractArray{NumDenom{T}}, thresh::Real) where T<:Real +function truncatenoise!(mm::AbstractArray{NumDenom{T}}, thresh::Real) where {T <: Real} for I in eachindex(mm) if mm[I].denom <= thresh - mm[I] = NumDenom{T}(0,0) + mm[I] = NumDenom{T}(0, 0) end end - mm + return mm end -function truncatenoise!(mms::AbstractArray{A}, thresh::Real) where A<:MismatchArray - for i = 1:length(denoms) +function truncatenoise!(mms::AbstractArray{A}, thresh::Real) where {A <: MismatchArray} + for i in 1:length(denoms) truncatenoise!(mms[i], thresh) end - nothing + return nothing end """ @@ -321,36 +321,36 @@ computes the integer-valued translation which best aligns images Optionally specify `thresh`, the fraction (0<=thresh<=1) of overlap required between `fixed` and `moving` (default 0.25). """ -function register_translate(fixed, moving, maxshift, thresh=nothing) +function register_translate(fixed, moving, maxshift, thresh = nothing) mm = mismatch(fixed, moving, maxshift) _, denom = separate(mm) - if thresh==nothing + if thresh == nothing thresh = 0.25maximum(denom) end - indmin_mismatch(mm, thresh) + return indmin_mismatch(mm, thresh) end function checksize_maxshift(A::AbstractArray, maxshift) ndims(A) == length(maxshift) || error("Array is $(ndims(A))-dimensional, but maxshift has length $(length(maxshift))") - for i = 1:ndims(A) - size(A,i) == 2*maxshift[i]+1 || error("Along dimension $i, the output size $(size(A,i)) does not agree with maxshift[$i] = $(maxshift[i])") + for i in 1:ndims(A) + size(A, i) == 2 * maxshift[i] + 1 || error("Along dimension $i, the output size $(size(A, i)) does not agree with maxshift[$i] = $(maxshift[i])") end - nothing + return nothing end function padranges(blocksize, maxshift) padright = [maxshift...] - transformdims = findall(padright.>0) - paddedsz = [blocksize...] + 2*padright + transformdims = findall(padright .> 0) + paddedsz = [blocksize...] + 2 * padright for i in transformdims # Pick a size for which one can efficiently calculate ffts padright[i] += padsize(blocksize, maxshift, i) - paddedsz[i] end - rng = UnitRange{Int}[ 1-maxshift[i]:blocksize[i]+padright[i] for i = 1:length(blocksize) ] + return rng = UnitRange{Int}[ (1 - maxshift[i]):(blocksize[i] + padright[i]) for i in 1:length(blocksize) ] end -padsize(blocksize::Dims{N}, maxshift::Dims{N}) where N = map(padsize, blocksize, maxshift, ntuple(identity, Val(N))) +padsize(blocksize::Dims{N}, maxshift::Dims{N}) where {N} = map(padsize, blocksize, maxshift, ntuple(identity, Val(N))) function padsize(blocksize::Int, maxshift::Int, dim::Int) p = blocksize + 2maxshift @@ -358,7 +358,7 @@ function padsize(blocksize::Int, maxshift::Int, dim::Int) end function assertsamesize(A, B) - if !issamesize(A,B) + return if !issamesize(A, B) error("Arrays are not the same size") end end @@ -366,31 +366,31 @@ end function issamesize(A::AbstractArray, B::AbstractArray) n = ndims(A) ndims(B) == n || return false - for i = 1:n + for i in 1:n axes(A, i) == axes(B, i) || return false end - true + return true end function issamesize(A::AbstractArray, indices) n = ndims(A) length(indices) == n || return false - for i = 1:n + for i in 1:n size(A, i) == length(indices[i]) || return false end - true + return true end # This yields the _effective_ overlap, i.e., sets to zero if gridsize==1 along a coordinate # imgssz = image spatial size function computeoverlap(imgssz, blocksize, gridsize) - gsz1 = max(1,[gridsize...].-1) - tmp = [imgssz...]./gsz1 - blocksize - [ceil(Int, x) for x in tmp] + gsz1 = max(1, [gridsize...] .- 1) + tmp = [imgssz...] ./ gsz1 + return blocksize - [ceil(Int, x) for x in tmp] end -leftedge(center, width) = ceil(Int, center-width/2) -rightedge(center, width) = leftedge(center+width, width) - 1 +leftedge(center, width) = ceil(Int, center - width / 2) +rightedge(center, width) = leftedge(center + width, width) - 1 # These avoid making a copy if it's not necessary tovec(v::AbstractVector) = v @@ -402,20 +402,22 @@ shiftrange(r, s) = r .+ s # TODO: redesign this whole thing to be safer? using Base: ViewIndex, to_indices, unsafe_length, index_shape, tail -@inline function extraunsafe_view(V::SubArray{T,N}, I::Vararg{ViewIndex,N}) where {T,N} +@inline function extraunsafe_view(V::SubArray{T, N}, I::Vararg{ViewIndex, N}) where {T, N} idxs = unsafe_reindex(V, V.indices, to_indices(V, I)) - SubArray(V.parent, idxs) + return SubArray(V.parent, idxs) end function get_index_wo_boundcheck(r::AbstractUnitRange, s::AbstractUnitRange{<:Integer}) # BlockRegistration issue #36 f = first(r) - st = oftype(f, f + first(s)-1) - range(st, length=length(s)) + st = oftype(f, f + first(s) - 1) + return range(st, length = length(s)) end unsafe_reindex(V, idxs::Tuple{UnitRange, Vararg{Any}}, subidxs::Tuple{UnitRange, Vararg{Any}}) = - (Base.@_propagate_inbounds_meta; @inbounds new1 = get_index_wo_boundcheck(idxs[1],subidxs[1]); - (new1, unsafe_reindex(V, tail(idxs), tail(subidxs))...)) + ( + Base.@_propagate_inbounds_meta; @inbounds new1 = get_index_wo_boundcheck(idxs[1], subidxs[1]); + (new1, unsafe_reindex(V, tail(idxs), tail(subidxs))...) +) unsafe_reindex(V, idxs, subidxs) = Base.reindex(V, idxs, subidxs) @@ -424,15 +426,15 @@ unsafe_reindex(V, idxs, subidxs) = Base.reindex(V, idxs, subidxs) function padsize(blocksize, maxshift) Base.depwarn("padsize(::$(typeof(blocksize)), ::$(typeof(maxshift)) is deprecated, use Dims-tuples instead", :padsize) sz = Vector{Int}(undef, length(blocksize)) - padsize!(sz, blocksize, maxshift) + return padsize!(sz, blocksize, maxshift) end function padsize!(sz::Vector, blocksize, maxshift) n = length(blocksize) - for i = 1:n + for i in 1:n sz[i] = padsize(blocksize, maxshift, i) end - sz + return sz end function padsize(blocksize, maxshift, dim)