diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index a5c6a36..70a4a34 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -41,6 +41,7 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v4 with: file: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} diff --git a/Project.toml b/Project.toml index b9785a4..021aca7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RegisterMismatchCommon" uuid = "abb2e897-52bf-5d28-a379-6ca321e3b878" -authors = ["Tim Holy "] version = "1.0.0" +authors = ["Tim Holy "] [deps] CenterIndexedArrays = "46a7138f-0d70-54e1-8ada-fb8296f91f24" @@ -9,13 +9,20 @@ ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" RegisterCore = "67712758-55e7-5c3c-8e85-dda1d7758434" [compat] -CenterIndexedArrays = "0.2" +Aqua = "0.8" +CenterIndexedArrays = "0.2, 1" +Documenter = "1" +ExplicitImports = "1" ImageCore = "0.8.1, 0.9, 0.10" -RegisterCore = "0.2, 1" +RegisterCore = "1" +Test = "1" julia = "1.10" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Aqua", "Documenter", "ExplicitImports", "Test"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..2dce4b9 --- /dev/null +++ b/README.md @@ -0,0 +1,187 @@ +# RegisterMismatchCommon.jl + +[![CI](https://github.com/HolyLab/RegisterMismatchCommon.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/RegisterMismatchCommon.jl/actions/workflows/CI.yml) +[![codecov](https://codecov.io/gh/HolyLab/RegisterMismatchCommon.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HolyLab/RegisterMismatchCommon.jl) + +RegisterMismatchCommon provides the shared types, utilities, and aperture-based +workflow helpers used for image-registration mismatch computation in the HolyLab +ecosystem. Concrete `mismatch` implementations live in downstream packages: + +- [`RegisterMismatch`](https://github.com/HolyLab/RegisterMismatch.jl) — CPU/FFTW +- [`RegisterMismatchCuda`](https://github.com/HolyLab/RegisterMismatchCuda.jl) — GPU/CUFFT + +In practice you will `using RegisterMismatch` (or its CUDA variant) rather than +this package directly. RegisterMismatchCommon is an explicit dependency for code +that needs to type-annotate or dispatch on `MismatchArray` / `NumDenom` without +caring which backend is loaded. + +## Installation + +This package lives in the [HolyLab registry](https://github.com/HolyLab/HolyLabRegistry). +Add the registry once, then install normally: + +```julia +using Pkg +pkg"registry add General https://github.com/HolyLab/HolyLabRegistry.git" +Pkg.add("RegisterMismatchCommon") +``` + +## Key concepts + +### `NumDenom` and `MismatchArray` + +Mismatch is stored as a ratio: a *numerator* (sum of squared pixel differences) +over a *denominator* (normalization factor). Keeping them separate lets +downstream code combine apertures, apply thresholds, and interpolate without +losing information. `NumDenom{T}` is a two-field struct from +[`RegisterCore`](https://github.com/HolyLab/RegisterCore.jl), and a +`MismatchArray` is a `CenterIndexedArray` of `NumDenom` values indexed from +`-maxshift` to `+maxshift` along each dimension. + +### Aperture workflow + +Rather than computing a single mismatch over the whole image, the aperture +workflow divides the image into overlapping sub-regions (*apertures*) and +computes one `MismatchArray` per aperture. This yields a spatially-resolved +shift field, which is the starting point for non-rigid registration. + +The typical sequence is: + +``` +aperture_grid → center coordinates for each aperture +allocate_mmarrays → pre-allocated output array +mismatch_apertures → fill the output (needs a backend loaded) +correctbias! → remove pixel-bias artifacts +truncatenoise! → zero out low-signal entries +``` + +### `mismatch` requires a backend + +`mismatch` and functions that call it (`mismatch_apertures`, `register_translate`) +are stubs defined here but not implemented. Loading `RegisterMismatch` or +`RegisterMismatchCuda` extends them with a concrete method. Calling them without +a backend will throw a `MethodError`. + +## Usage examples + +### Aperture grid + +```julia +using RegisterMismatchCommon + +# 64×64 image, 4×4 grid of apertures +ag = aperture_grid((64, 64), (4, 4)) +size(ag) # (4, 4) +ag[1, 1] # (1.0, 1.0) — top-left corner +ag[4, 4] # (64.0, 64.0) — bottom-right corner +``` + +### Aperture width + +```julia +img = zeros(Float32, 64, 64) +default_aperture_width(img, (4, 4)) # (21.0, 21.0) +``` + +### Allocating output arrays + +```julia +# Pre-allocate a 4×4 grid of MismatchArrays, each covering ±5 pixels +mms = allocate_mmarrays(Float32, (4, 4), (5, 5)) +size(mms) # (4, 4) +size(mms[1, 1]) # (11, 11) — 2*5+1 per dimension +``` + +### Zero-shift mismatch (no backend needed) + +`mismatch0` measures the mismatch at zero shift directly, without FFTs: + +```julia +fixed = [1.0 2.0; 3.0 4.0] +moving = [1.0 2.0; 3.0 4.0] +mismatch0(fixed, moving) # NumDenom(0.0, 60.0) — perfect match + +moving2 = [2.0 3.0; 4.0 5.0] +mm0 = mismatch0(fixed, moving2) # NumDenom(4.0, 84.0) +mm0.num / mm0.denom # ≈ 0.048 — normalized mismatch +``` + +### Padding mismatched arrays + +```julia +a = [1.0 2.0; 3.0 4.0] # 2×2 +b = [5.0 6.0 7.0; 8.0 9.0 10.0] # 2×3 +ap, bp = nanpad(a, b) +# ap is 2×3, padded with NaN in the third column +# bp is returned unchanged +``` + +### Post-processing + +```julia +# Zero out entries whose denominator is too small to be reliable +truncatenoise!(mms, 0.5f0) + +# Impute zero-shift entries that are corrupted by pixel-bias +correctbias!(mms) +``` + +### Full-image translation (with a backend) + +```julia +using RegisterMismatch # or RegisterMismatchCuda + +fixed = rand(Float32, 64, 64) +moving = rand(Float32, 64, 64) + +shift = register_translate(fixed, moving, (10, 10)) +# CartesianIndex of best integer shift +``` + +## API reference + +### Core computation + +| Function | Description | +|---|---| +| `mismatch` | Full-array mismatch (requires backend) | +| `mismatch_apertures` | Aperture-wise mismatch on a grid (requires backend) | +| `mismatch0` | Zero-shift mismatch without FFTs | +| `register_translate` | Best integer translation (requires backend) | + +### Aperture workflow + +| Function | Description | +|---|---| +| `aperture_grid` | Uniformly-spaced grid of aperture center coordinates | +| `allocate_mmarrays` | Pre-allocate an array of `MismatchArray`s | +| `default_aperture_width` | Compute aperture width for a given grid | +| `aperture_range` | `UnitRange` indices for one aperture | +| `each_point` | Iterate over aperture centers in any layout | + +### Post-processing + +| Function | Description | +|---|---| +| `correctbias!` | Impute pixel-bias-corrupted zero-shift entries | +| `truncatenoise!` | Zero out low-denominator (unreliable) entries | + +### Utilities + +| Function | Description | +|---|---| +| `nanpad` | Pad two arrays to the same size with `NaN` | +| `padsize` | FFT-friendly padded size | +| `padranges` | Padded index ranges for FFT cross-correlation | +| `checksize_maxshift` | Validate a mismatch array's size | +| `assertsamesize` | Throw if two arrays differ in size | +| `tovec` | Convert a tuple to a `Vector` | +| `shiftrange` | Shift a range by a scalar | +| `set_FFTPROD` | Set the allowed FFT prime factors | + +### Types + +| Name | Description | +|---|---| +| `DimsLike` | `Union{AbstractVector{Int}, Dims}` — dimension-size argument | +| `WidthLike` | `Union{AbstractVector, Tuple}` — aperture-width argument | diff --git a/src/RegisterMismatchCommon.jl b/src/RegisterMismatchCommon.jl index 3ea2770..7042cad 100644 --- a/src/RegisterMismatchCommon.jl +++ b/src/RegisterMismatchCommon.jl @@ -1,17 +1,83 @@ +""" +RegisterMismatchCommon provides shared types, utilities, and aperture-based +workflow helpers for image-registration mismatch computation. Concrete +`mismatch` implementations are supplied by downstream packages: +`RegisterMismatch` (CPU/FFTW) and `RegisterMismatchCuda` (GPU/CUFFT). + +Main entry points: +- [`mismatch`](@ref): full-array mismatch between two images +- [`mismatch_apertures`](@ref): aperture-wise mismatch on a grid +- [`register_translate`](@ref): find the best integer translation + +Aperture workflow helpers: [`aperture_grid`](@ref), [`allocate_mmarrays`](@ref), +[`default_aperture_width`](@ref), [`aperture_range`](@ref). + +Post-processing: [`correctbias!`](@ref), [`truncatenoise!`](@ref), [`mismatch0`](@ref). +""" module RegisterMismatchCommon -using RegisterCore, CenterIndexedArrays, ImageCore +using CenterIndexedArrays: CenterIndexedArrays, CenterIndexedArray +using ImageCore: ImageCore, coords_spatial, sdims +using RegisterCore: RegisterCore, MismatchArray, NumDenom, argmin_mismatch, maxshift, separate export correctbias!, nanpad, mismatch0, aperture_grid, allocate_mmarrays, default_aperture_width, truncatenoise! export DimsLike, WidthLike, each_point, aperture_range, assertsamesize, tovec, mismatch, padsize, set_FFTPROD 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 DimsLike = Union{AbstractVector{Int}, Dims} + +Type alias accepted wherever dimension sizes can be passed. + +!!! note + Prefer `Dims` tuples (e.g. `(128, 128)`) over `AbstractVector{Int}` where + possible; tuple inputs carry dimensionality in the type, which enables + more precise dispatch and avoids runtime length checks. +""" +const DimsLike = Union{AbstractVector{Int}, Dims} + +""" + const WidthLike = Union{AbstractVector, Tuple} + +Type alias for values that specify aperture widths — accepted as either a +`Tuple` or an `AbstractVector`. +""" const WidthLike = Union{AbstractVector, Tuple} + FFTPROD = [2, 3] + +""" + set_FFTPROD(v) + +Set the global list of prime factors used when selecting FFT-friendly array +sizes. The default is `[2, 3]`, meaning padded sizes are products of powers +of 2 and 3. + +This setting affects `padsize` and `padranges`. + +# Example +```julia +set_FFTPROD([2, 3, 5]) # allow sizes of the form 2^a * 3^b * 5^c +``` +""" set_FFTPROD(v) = global FFTPROD = v +""" + mismatch(fixed, moving, maxshift; normalization=:intensity) -> MismatchArray + +Compute the mismatch between `fixed` and `moving` for all integer shifts up to +`maxshift`. + +Returns a `MismatchArray` of `NumDenom` values indexed from `-maxshift` to +`+maxshift` along each dimension. `normalization` may be `:intensity` +(default) or `:pixels`. + +!!! note + This function dispatches to a concrete implementation that must be provided + by a downstream package such as `RegisterMismatch` (CPU) or + `RegisterMismatchCuda` (GPU). +""" 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) @@ -26,15 +92,19 @@ function mismatch_apertures(::Type{T}, fixed::AbstractArray, moving::AbstractArr end """ -`correctbias!(mm::MismatchArray)` replaces "suspect" mismatch -data with imputed data. If each pixel in your camera has a different -bias, then matching that bias becomes an incentive to avoid -shifts. Likewise, CMOS cameras tend to have correlated row/column -noise. These two factors combine to imply that `mm[i,j,...]` is unreliable -whenever `i` or `j` is zero. + correctbias!(mm::MismatchArray[, w]) -> mm + correctbias!(mms::AbstractArray{<:MismatchArray}) -> mms -Data are imputed by averaging the adjacent non-suspect values. This -function works in-place, overwriting the original `mm`. +Replace "suspect" mismatch entries with values imputed from their neighbors. + +Camera pixel-to-pixel bias and CMOS row/column noise create a spurious +incentive to avoid shifts, because `mm[i,j,...]` is unreliable whenever +`i` or `j` is zero (along the first two spatial dimensions). Suspect +entries are identified via weight array `w` (zero = suspect; default +computed by `correctbias_weight`), and replaced by a weighted average of +adjacent non-suspect values. + +Both forms modify their first argument in place and return it. """ function correctbias!(mm::MismatchArray{ND, N}, w = correctbias_weight(mm)) where {ND, N} T = eltype(ND) @@ -61,7 +131,6 @@ function correctbias!(mm::MismatchArray{ND, N}, w = correctbias_weight(mm)) wher return mm end -"`correctbias!(mms)` runs `correctbias!` on each element of an array-of-MismatchArrays." function correctbias!(mms::AbstractArray{M}) where {M <: MismatchArray} for mm in mms correctbias!(mm) @@ -85,9 +154,19 @@ function correctbias_weight(mm::MismatchArray{ND, N}) where {ND, N} end """ -`fixedpad, movingpad = nanpad(fixed, moving)` will pad `fixed` and/or -`moving` with NaN as needed to ensure that `fixedpad` and `movingpad` -have the same size. + fixedpad, movingpad = nanpad(fixed, moving) + +Pad `fixed` and/or `moving` with `NaN` so that both outputs have the same +size. + +If `fixed` and `moving` already have the same size they are returned unchanged. +Otherwise each is expanded along every dimension to +`max(size(fixed, d), size(moving, d))`, padding out-of-bounds positions with +`NaN`. + +The element type of the outputs is `promote_type(eltype(fixed), eltype(moving))`. +For non-floating-point inputs this promotes to at least `Float32` (the smallest +float type that can represent `NaN`). """ function nanpad(fixed, moving) ndims(fixed) == ndims(moving) || error("fixed and moving must have the same dimensionality") @@ -103,9 +182,16 @@ nanval(::Type{T}) where {T <: AbstractFloat} = convert(T, NaN) nanval(::Type{T}) where {T} = convert(Float32, NaN) """ -`mm0 = mismatch0(fixed, moving, [normalization])` computes the -"as-is" mismatch between `fixed` and `moving`, without any shift. -`normalization` may be either `:intensity` (the default) or `:pixels`. + mm0 = mismatch0(fixed, moving; normalization=:intensity) -> NumDenom{Float64} + +Compute the "as-is" mismatch between `fixed` and `moving` at zero shift. + +`normalization` may be `:intensity` (default, normalizes by `vf² + vm²`) or +`:pixels` (normalizes by the count of finite pixel pairs). Returns a +`NumDenom{Float64}`; the ratio `mm0.num / mm0.denom` gives the normalized +mismatch. + +See also: [`mismatch0(mms)`](@ref mismatch0(::AbstractArray)). """ 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")) @@ -138,10 +224,13 @@ function _mismatch0(num::T, denom::T, fixed::AbstractArray{Tf, N}, moving::Abstr end """ -`mm0 = mismatch0(mms)` computes the "as-is" -mismatch between `fixed` and `moving`, without any shift. The -mismatch is represented in `mms` as an aperture-wise -Arrays-of-MismatchArrays. + mm0 = mismatch0(mms::AbstractArray{<:MismatchArray}) -> NumDenom + +Extract and sum the zero-shift `NumDenom` from each element of an +aperture-wise array-of-`MismatchArray`s, returning a single `NumDenom` +representing the overall as-is mismatch. + +See also: [`mismatch0(fixed, moving)`](@ref mismatch0(::AbstractArray, ::AbstractArray)). """ function mismatch0(mms::AbstractArray{M}) where {M <: MismatchArray} mm0 = eltype(M)(0, 0) @@ -154,10 +243,29 @@ function mismatch0(mms::AbstractArray{M}) where {M <: MismatchArray} end """ -`ag = aperture_grid(ssize, gridsize)` constructs a uniformly-spaced -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. + ag = aperture_grid(ssize, gridsize) -> Array{NTuple{N,Float64},N} + +Construct a uniformly-spaced grid of aperture centers for an image of spatial +size `ssize`. + +The returned array has size `gridsize` and element type `NTuple{N,Float64}`. +Along each dimension, centers are linearly spaced between 1 and `ssize[d]` +(inclusive). When `gridsize[d] == 1`, the single center is placed at the +midpoint `(ssize[d] + 1) / 2`. + +# Examples +```jldoctest +julia> ag = aperture_grid((256, 256), (4, 4)); + +julia> size(ag) +(4, 4) + +julia> ag[1, 1] +(1.0, 1.0) + +julia> ag[4, 4] +(256.0, 256.0) +``` """ function aperture_grid(ssize::Dims{N}, gridsize) where {N} if length(gridsize) != N @@ -175,20 +283,19 @@ function aperture_grid(ssize::Dims{N}, gridsize) where {N} end """ -`mms = allocate_mmarrays(T, gridsize, maxshift)` allocates storage for -aperture-wise mismatch computation. `mms` will be an -Array-of-MismatchArrays with element type `NumDenom{T}` and half-size -`maxshift`. `mms` will be an array of size `gridsize`. This syntax is -recommended when your apertures are centered at points of a grid. + mms = allocate_mmarrays(T, gridsize::NTuple, maxshift) -> Array{MismatchArray{NumDenom{T},N},N} + mms = allocate_mmarrays(T, aperture_centers, maxshift) -> Array{MismatchArray{NumDenom{T},N}} + +Allocate an array of `MismatchArray{NumDenom{T}}` objects, each with half-size `maxshift`. + +**`gridsize` form** (recommended for regular grids): `mms` is an `N`-dimensional +array of size `gridsize`. -`mms = allocate_mmarrays(T, aperture_centers, maxshift)` returns `mms` -with a shape that matches that of `aperture_centers`. The centers can -in general be provided as an vector-of-tuples, vector-of-vectors, or a -matrix with each point in a column. If your centers are arranged in a -rectangular grid, you can use an `N`-dimensional array-of-tuples (or -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.) +**`aperture_centers` form**: `mms` matches the shape of `aperture_centers`. Centers +may be provided as: +- An array of tuples or `AbstractVector`s: `mms` has the same shape. +- An `AbstractMatrix{<:Real}` where each column encodes one point: `mms` has shape + `size(aperture_centers)[2:end]`, i.e., one entry per column. """ function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{C}, maxshift) where {T, C <: Union{AbstractVector, Tuple}} isempty(aperture_centers) && error("aperture_centers is empty") @@ -208,9 +315,9 @@ function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{C}, maxshi return mms end -function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{R}, maxshift) where {T, R <: Real} +function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{R}, maxshift) where {T <: Real, R <: Real} N = ndims(aperture_centers) - 1 - mms = Array{MismatchArray{T, N}}(undef, size(aperture_centers)[2:end]) + mms = Array{MismatchArray{NumDenom{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...) @@ -233,14 +340,18 @@ end Base.iterate(iter::ContainerIterator) = iterate(iter.data) Base.iterate(iter::ContainerIterator, state) = iterate(iter.data, state) +Base.length(iter::ContainerIterator) = length(iter.data) +Base.eltype(::Type{ContainerIterator{C}}) where {C} = eltype(C) 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(size(data)[2:end])) end -FirstDimIterator(A::AbstractArray) = FirstDimIterator{typeof(A), typeof(CartesianIndices(Base.tail(size(A))))}(A) +FirstDimIterator(A::AbstractArray) = FirstDimIterator{typeof(A), typeof(CartesianIndices(size(A)[2:end]))}(A) +Base.length(iter::FirstDimIterator) = length(iter.rng) +Base.eltype(::Type{FirstDimIterator{A, R}}) where {A, R} = Vector{eltype(A)} function Base.iterate(iter::FirstDimIterator) isempty(iter.rng) && return nothing index, state = iterate(iter.rng) @@ -253,11 +364,14 @@ function Base.iterate(iter::FirstDimIterator, state) end """ -`iter = each_point(points)` yields an iterator `iter` over all the -points in `points`. `points` may be represented as an -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). + iter = each_point(points) + +Return an iterator over the points in `points`. + +`points` may be: +- An `AbstractArray` of tuples or `AbstractVector`s: each element is yielded as-is. +- An `AbstractArray{<:Real}` where points are laid out along the first dimension + (e.g., columns of a matrix): each point is yielded as a `Vector`. """ each_point(aperture_centers::AbstractArray{C}) where {C <: Union{AbstractVector, Tuple}} = ContainerIterator(aperture_centers) @@ -274,19 +388,24 @@ function aperture_range(center, width) end """ -`aperturesize = default_aperture_width(img, gridsize, [overlap])` -calculates the aperture width for a regularly-spaced grid of aperture -centers with size `gridsize`. Apertures that are adjacent along -dimension `d` may overlap by a number pixels specified by -`overlap[d]`; the default value is 0. For non-negative `overlap`, the -collection of apertures will yield full coverage of the image. + aperturesize = default_aperture_width(img, gridsize[, overlap]) -> NTuple{N,Float64} + +Compute the aperture width for a regularly-spaced grid of aperture centers +with size `gridsize` over image `img`. + +`overlap` is a `DimsLike` giving the number of pixels by which adjacent +apertures overlap along each spatial dimension; it defaults to zero in every +dimension. For non-negative `overlap`, the collection of apertures provides +full coverage of the image. + +Returns a `Tuple` of `Float64` widths, one per spatial dimension. """ 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 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") + error("gridsize $gridsize is too large, given the size $(map(d -> size(img, d), sc)) of the image") end end gsz1 = max.(1, [gridsize...] .- 1) @@ -295,8 +414,13 @@ function default_aperture_width(img, gridsize::DimsLike, overlap::DimsLike = zer end """ -`truncatenoise!(mm, thresh)` zeros out any entries of the -MismatchArray `mm` whose `denom` values are less than `thresh`. + truncatenoise!(mm::AbstractArray{NumDenom{T}}, thresh) -> mm + truncatenoise!(mms::AbstractArray{<:MismatchArray}, thresh) -> mms + +Zero out entries whose `denom` is ≤ `thresh`, replacing them with +`NumDenom(0, 0)`. + +Both forms modify their first argument in place and return it. """ function truncatenoise!(mm::AbstractArray{NumDenom{T}}, thresh::Real) where {T <: Real} for I in eachindex(mm) @@ -308,18 +432,28 @@ function truncatenoise!(mm::AbstractArray{NumDenom{T}}, thresh::Real) where {T < end function truncatenoise!(mms::AbstractArray{A}, thresh::Real) where {A <: MismatchArray} - for i in 1:length(denoms) + for i in 1:length(mms) truncatenoise!(mms[i], thresh) end - return nothing + return mms end """ -`shift = register_translate(fixed, moving, maxshift, [thresh])` -computes the integer-valued translation which best aligns images -`fixed` and `moving`. All shifts up to size `maxshift` are considered. -Optionally specify `thresh`, the fraction (0<=thresh<=1) of overlap -required between `fixed` and `moving` (default 0.25). + shift = register_translate(fixed, moving, maxshift[, thresh]) -> CartesianIndex + +Compute the integer-valued translation that best aligns `fixed` and `moving`. +All shifts up to `maxshift` (in each dimension) are considered. + +`thresh` sets the minimum `denom` value for a mismatch entry to be considered +reliable. Entries with `denom ≤ thresh` are excluded from the search. The +default is `0.25 * maximum(denom)`, i.e., entries whose denominator falls +below 25% of the peak denominator are excluded. + +Returns a `CartesianIndex` of the best integer shift. + +!!! note + Requires a concrete `mismatch` implementation to be loaded, e.g. from + `RegisterMismatch` (CPU) or `RegisterMismatchCuda` (GPU). """ function register_translate(fixed, moving, maxshift, thresh = nothing) mm = mismatch(fixed, moving, maxshift) @@ -327,10 +461,19 @@ function register_translate(fixed, moving, maxshift, thresh = nothing) if thresh == nothing thresh = 0.25maximum(denom) end - return indmin_mismatch(mm, thresh) + return argmin_mismatch(mm, thresh) end +""" + checksize_maxshift(A, maxshift) + +Validate that array `A` has the size expected for a mismatch array with the +given `maxshift`. Checks that `size(A, d) == 2maxshift[d] + 1` for every +dimension `d`. + +Throws an `ErrorException` on failure; returns `nothing` on success. +""" function checksize_maxshift(A::AbstractArray, maxshift) ndims(A) == length(maxshift) || error("Array is $(ndims(A))-dimensional, but maxshift has length $(length(maxshift))") for i in 1:ndims(A) @@ -339,6 +482,16 @@ function checksize_maxshift(A::AbstractArray, maxshift) return nothing end +""" + padranges(blocksize, maxshift) -> Vector{UnitRange{Int}} + +Compute padded index ranges for an FFT-based cross-correlation. + +For each dimension `d`, returns a `UnitRange{Int}` that covers the block +`1:blocksize[d]` and is extended symmetrically by `maxshift[d]`, then rounded +up to an FFT-friendly length via `padsize`. Ranges start at `1 - maxshift[d]` +to accommodate negative shifts. +""" function padranges(blocksize, maxshift) padright = [maxshift...] transformdims = findall(padright .> 0) @@ -350,6 +503,20 @@ function padranges(blocksize, maxshift) return rng = UnitRange{Int}[ (1 - maxshift[i]):(blocksize[i] + padright[i]) for i in 1:length(blocksize) ] end +""" + padsize(blocksize, maxshift, dim) -> Int + padsize(blocksize::Dims, maxshift::Dims) -> Dims + +Compute the padded FFT size. + +The single-dimension form returns the smallest integer ≥ `blocksize + 2maxshift` +that is efficient for FFT computation: a power of 2 along dimension 1, and a +product of `FFTPROD` primes along other dimensions. When `maxshift == 0`, +the input size is returned unchanged (that dimension will not be transformed). + +The multi-dimension form applies `padsize` independently to each dimension and +returns a tuple of padded sizes. +""" 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) @@ -357,6 +524,12 @@ function padsize(blocksize::Int, maxshift::Int, dim::Int) return maxshift > 0 ? (dim == 1 ? nextpow(2, p) : nextprod(FFTPROD, p)) : p # we won't FFT along dimensions with maxshift 0 end +""" + assertsamesize(A, B) + +Throw an `ErrorException` if `A` and `B` do not have the same axes along every +dimension. Returns `nothing` on success. +""" function assertsamesize(A, B) return if !issamesize(A, B) error("Arrays are not the same size") @@ -384,25 +557,36 @@ 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) + gsz1 = max.(1, [gridsize...] .- 1) tmp = [imgssz...] ./ gsz1 - return blocksize - [ceil(Int, x) for x in tmp] + 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 -# These avoid making a copy if it's not necessary +""" + tovec(v) -> AbstractVector + +Convert `v` to an `AbstractVector`, returning `v` unchanged if it is already +an `AbstractVector` (no copy), or converting a `Tuple` to a `Vector`. +""" tovec(v::AbstractVector) = v tovec(v::Tuple) = [v...] +""" + shiftrange(r, s) + +Shift the range `r` by scalar `s`, returning `r .+ s`. +""" shiftrange(r, s) = r .+ s ### Utilities for unsafe indexing of views # TODO: redesign this whole thing to be safer? -using Base: ViewIndex, to_indices, unsafe_length, index_shape, tail +using Base: to_indices +_tail(t::Tuple) = t[2:end] -@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{Union{Real, AbstractArray}, N}) where {T, N} idxs = unsafe_reindex(V, V.indices, to_indices(V, I)) return SubArray(V.parent, idxs) end @@ -415,27 +599,17 @@ 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))...) + @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) - -### Deprecations - -function padsize(blocksize, maxshift) - Base.depwarn("padsize(::$(typeof(blocksize)), ::$(typeof(maxshift)) is deprecated, use Dims-tuples instead", :padsize) - sz = Vector{Int}(undef, length(blocksize)) - return padsize!(sz, blocksize, maxshift) -end - -function padsize!(sz::Vector, blocksize, maxshift) - n = length(blocksize) - for i in 1:n - sz[i] = padsize(blocksize, maxshift, i) - end - return sz -end +# Scalar indices in idxs are dropped dimensions — pass through without consuming a subindex +unsafe_reindex(V, idxs::Tuple{Real, Vararg{Any}}, subidxs::Tuple) = + (idxs[1], unsafe_reindex(V, _tail(idxs), subidxs)...) +# AbstractArray indices: map subindex through the stored index +unsafe_reindex(V, idxs::Tuple{AbstractArray, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) = + (@inbounds new1 = idxs[1][subidxs[1]]; (new1, unsafe_reindex(V, _tail(idxs), _tail(subidxs))...)) +unsafe_reindex(V, ::Tuple{}, ::Tuple{}) = () function padsize(blocksize, maxshift, dim) m = maxshift[dim] diff --git a/test/runtests.jl b/test/runtests.jl index fd9c58b..5eb128e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,298 @@ using Test +using Aqua +using Documenter +using ExplicitImports import RegisterMismatchCommon +using RegisterMismatchCommon +using RegisterCore: MismatchArray, NumDenom, maxshift, separate # Tests are in RegisterMismatch and RegisterMismatchCuda + +DocMeta.setdocmeta!(RegisterMismatchCommon, :DocTestSetup, :(using RegisterMismatchCommon); recursive=true) + +@testset "RegisterMismatchCommon" begin + @testset "Doctests" begin + doctest(RegisterMismatchCommon; manual=false) + end + Aqua.test_all(RegisterMismatchCommon) + @testset "ExplicitImports" begin + test_explicit_imports(RegisterMismatchCommon) + end + + @testset "set_FFTPROD" begin + old = copy(RegisterMismatchCommon.FFTPROD) + set_FFTPROD([2, 3, 5]) + @test RegisterMismatchCommon.FFTPROD == [2, 3, 5] + set_FFTPROD(old) + @test RegisterMismatchCommon.FFTPROD == old + end + + @testset "tovec" begin + v = [1, 2, 3] + @test tovec(v) === v + @test tovec((1, 2, 3)) == [1, 2, 3] + end + + @testset "shiftrange" begin + @test shiftrange(1:5, 3) == 4:8 + @test shiftrange(-2:2, -1) == -3:1 + end + + @testset "nanpad" begin + a = [1.0 2.0; 3.0 4.0] + b = [5.0 6.0; 7.0 8.0] + # Same size: returns originals unchanged + ap, bp = nanpad(a, b) + @test ap === a + @test bp === b + + # Different sizes: pads smaller to larger + c = [1.0 2.0 3.0; 4.0 5.0 6.0] # 2×3 + d = [7.0 8.0; 9.0 10.0] # 2×2 + cp, dp = nanpad(c, d) + @test size(cp) == size(dp) == (2, 3) + @test cp[1, 1] == 1.0 && isnan(dp[1, 3]) + + # Dimension mismatch + @test_throws ErrorException nanpad([1.0], [1.0 2.0]) + end + + @testset "mismatch0 (fixed/moving)" begin + a = [1.0 2.0; 3.0 4.0] + # Identical arrays → num==0 + nd = mismatch0(a, a) + @test nd.num == 0.0 + @test nd.denom > 0.0 + + # Different arrays → positive num + b = [2.0 3.0; 4.0 5.0] + nd2 = mismatch0(a, b) + @test nd2.num > 0.0 + + # :pixels normalization + nd3 = mismatch0(a, b; normalization = :pixels) + @test nd3.denom == 4.0 # 4 finite pixels + + # NaN elements are skipped + c = [1.0 NaN; 3.0 4.0] + nd4 = mismatch0(a, c) + @test nd4.denom ≈ a[1,1]^2 + c[1,1]^2 + a[2,1]^2 + c[2,1]^2 + a[2,2]^2 + c[2,2]^2 + + # Size mismatch error + @test_throws DimensionMismatch mismatch0(a, [1.0 2.0]) + + # Unrecognized normalization + @test_throws ErrorException mismatch0(a, a; normalization = :bogus) + end + + @testset "mismatch0 (array-of-MismatchArrays)" begin + mm = MismatchArray(Float64, 3, 3) + mm[0, 0] = NumDenom{Float64}(1.0, 2.0) + mms = [mm] + nd = mismatch0(mms) + @test nd == NumDenom{Float64}(1.0, 2.0) + end + + @testset "aperture_grid" begin + # 1D, single cell: center at midpoint + g = aperture_grid((10,), (1,)) + @test size(g) == (1,) + @test g[1] == (5.5,) + + # 1D, two cells: at corners + g2 = aperture_grid((10,), (2,)) + @test g2[1] == (1.0,) + @test g2[2] == (10.0,) + + # 2D grid + g3 = aperture_grid((10, 20), (2, 3)) + @test size(g3) == (2, 3) + @test g3[1, 1] == (1.0, 1.0) + @test g3[2, 3] == (10.0, 20.0) + + # Length mismatch error + @test_throws ErrorException aperture_grid((10, 20), (2,)) + end + + @testset "aperture_range" begin + rng = aperture_range((5.0, 10.0), (4.0, 6.0)) + @test length(rng) == 2 + @test rng[1] == 3:6 + @test rng[2] == 7:12 + + @test_throws ErrorException aperture_range((1.0,), (1.0, 2.0)) + end + + @testset "each_point" begin + # Array-of-tuples + pts = [(1.0, 2.0), (3.0, 4.0)] + collected = collect(each_point(pts)) + @test collected == pts + + # Real matrix (columns are points) + M = [1.0 3.0; 2.0 4.0] # 2×2: two 2D points + collected2 = collect(each_point(M)) + @test collected2[1] == [1.0, 2.0] + @test collected2[2] == [3.0, 4.0] + end + + @testset "checksize_maxshift" begin + mm = MismatchArray(Float64, 3, 5) + @test checksize_maxshift(mm, (1, 2)) === nothing + + # Wrong ndims + @test_throws ErrorException checksize_maxshift(mm, (1,)) + + # Wrong size along dim 1 + @test_throws ErrorException checksize_maxshift(mm, (2, 2)) + end + + @testset "padsize" begin + # dim==1 uses nextpow(2,...) + @test padsize(10, 3, 1) == nextpow(2, 16) + # dim==2 uses nextprod + @test padsize(10, 3, 2) == nextprod(RegisterMismatchCommon.FFTPROD, 16) + # maxshift==0: no padding, no FFT + @test padsize(10, 0, 1) == 10 + @test padsize(10, 0, 2) == 10 + + # Dims version + ps = padsize((8, 12), (2, 3)) + @test ps[1] == padsize(8, 2, 1) + @test ps[2] == padsize(12, 3, 2) + end + + @testset "padranges" begin + rngs = padranges((10, 12), (2, 3)) + @test rngs[1].start == 1 - 2 + @test rngs[2].start == 1 - 3 + for (i, r) in enumerate(rngs) + @test length(r) >= 10 + 2 * [2, 3][i] + end + end + + @testset "issamesize / assertsamesize" begin + a = zeros(3, 4) + b = zeros(3, 4) + c = zeros(3, 5) + + @test RegisterMismatchCommon.issamesize(a, b) + @test !RegisterMismatchCommon.issamesize(a, c) + @test !RegisterMismatchCommon.issamesize(a, zeros(2, 3, 4)) + + # Array vs. indices form + @test RegisterMismatchCommon.issamesize(a, (1:3, 1:4)) + @test !RegisterMismatchCommon.issamesize(a, (1:3, 1:5)) + @test !RegisterMismatchCommon.issamesize(a, (1:3,)) + + @test assertsamesize(a, b) === nothing + @test_throws ErrorException assertsamesize(a, c) + end + + @testset "allocate_mmarrays" begin + # Array-of-tuples centers + centers = [(1.0, 2.0), (3.0, 4.0), (5.0, 6.0)] + mms = allocate_mmarrays(Float32, centers, (2, 3)) + @test size(mms) == (3,) + @test maxshift(mms[1]) == (2, 3) + + # Real matrix centers: rows are coordinates, columns are points + # A 2×3 matrix means 3 points in 2D space → N=1 (ndims-1=1), 3 mmarrays + M = [1.0 3.0 5.0; 2.0 4.0 6.0] + mms2 = allocate_mmarrays(Float32, M, (2,)) + @test size(mms2) == (3,) + + # Gridsize tuple + mms3 = allocate_mmarrays(Float32, (2, 3), (1, 2)) + @test size(mms3) == (2, 3) + @test maxshift(mms3[1, 1]) == (1, 2) + end + + @testset "truncatenoise!" begin + mm = MismatchArray(Float64, 3, 3) + for I in eachindex(mm) + mm[I] = NumDenom{Float64}(1.0, Float64(abs(I[1]) + abs(I[2]) + 1)) + end + thresh = 2.5 + truncatenoise!(mm, thresh) + for I in eachindex(mm) + d = abs(I[1]) + abs(I[2]) + 1 + if d <= thresh + @test mm[I] == NumDenom{Float64}(0.0, 0.0) + end + end + + # Array-of-MismatchArrays form + mms = [MismatchArray(Float64, 3, 3) for _ in 1:2] + for mm in mms, I in eachindex(mm) + mm[I] = NumDenom{Float64}(1.0, 0.5) + end + result = truncatenoise!(mms, 1.0) + @test result === mms + for mm in mms, I in eachindex(mm) + @test mm[I] == NumDenom{Float64}(0.0, 0.0) + end + end + + @testset "correctbias!" begin + # 2D MismatchArray: row 0 and column 0 are suspect + mm = MismatchArray(Float64, 5, 5) # shifts -2:2 in each dim + good = NumDenom{Float64}(1.0, 1.0) + bad = NumDenom{Float64}(99.0, 99.0) + for I in eachindex(mm) + mm[I] = (I[1] == 0 || I[2] == 0) ? bad : good + end + correctbias!(mm) + # All suspect entries should now be replaced + for I in eachindex(mm) + if I[1] == 0 || I[2] == 0 + @test mm[I] != bad + end + end + + # Array-of-MismatchArrays form + mms = [MismatchArray(Float64, 5, 5) for _ in 1:2] + for mm2 in mms, I in eachindex(mm2) + mm2[I] = (I[1] == 0 || I[2] == 0) ? bad : good + end + result = correctbias!(mms) + @test result === mms + for mm2 in mms, I in eachindex(mm2) + if I[1] == 0 || I[2] == 0 + @test mm2[I] != bad + end + end + end + + @testset "default_aperture_width" begin + img = zeros(Float32, 20, 30) + # Single cell: full image width + w = default_aperture_width(img, (1, 1)) + @test w == (20.0, 30.0) + + # 2×3 grid + w2 = default_aperture_width(img, (2, 3)) + @test length(w2) == 2 + + # gridsize too large + @test_throws ErrorException default_aperture_width(img, (25, 3)) + end + + @testset "nanval" begin + @test isnan(RegisterMismatchCommon.nanval(Float64)) + @test isnan(RegisterMismatchCommon.nanval(Float32)) + # Non-float falls back to Float32 NaN + v = RegisterMismatchCommon.nanval(Int32) + @test isa(v, Float32) && isnan(v) + end + + @testset "computeoverlap" begin + # 2×2 grid on a 10×10 image with 6×6 blocksize → spacing = ceil(10/1) = 10, overlap = 6-10 = -4 + ov = RegisterMismatchCommon.computeoverlap((10, 10), (6, 6), (2, 2)) + @test length(ov) == 2 + @test ov == [6 - 10, 6 - 10] + # gridsize 1 along first dim: gsz1 = max.(1, [0, 1]) = [1, 1] + ov2 = RegisterMismatchCommon.computeoverlap((10, 10), (6, 6), (1, 2)) + @test length(ov2) == 2 + end +end