diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status new file mode 100644 index 0000000..3b616fe --- /dev/null +++ b/.claude/freshen-package-status @@ -0,0 +1,11 @@ +DONE: design review +DONE: API review +TODO: update .gitignore +TODO: format with runic +TODO: add Aqua.jl +TODO: remove deprecations +TODO: add ExplicitImports.jl +TODO: limit struct mutability +TODO: improve test coverage +TODO: add and improve docstrings +TODO: add or improve documentation diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 71003c2..160a54c 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,7 +13,7 @@ jobs: fail-fast: false matrix: version: - - '1.6' + - 'min' - '1' # - 'nightly' os: @@ -23,7 +23,7 @@ jobs: - x64 steps: - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} @@ -47,6 +47,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@v5 with: file: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} \ No newline at end of file diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index 2c4f30b..0cf0937 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -28,8 +28,5 @@ jobs: steps: - uses: JuliaRegistries/TagBot@v1 with: - token: ${{ secrets.GITHUB_TOKEN }} - # Edit the following line to reflect the actual name of the GitHub Secret containing your private key - ssh: ${{ secrets.DOCUMENTER_KEY }} - # ssh: ${{ secrets.NAME_OF_MY_SSH_PRIVATE_KEY_SECRET }} + token: ${{ secrets.TAGBOT_TOKEN }} registry: HolyLab/HolyLabRegistry diff --git a/NEWS.md b/NEWS.md index 8101d8f..e2c9f24 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,25 @@ -# version 3.0 +# version 1.0.0 + +## Breaking changes + +- `rotation_gridsearch`: the `SD` argument is now a keyword argument. Calls that + passed `SD` positionally (`rotation_gridsearch(fixed, moving, maxshift, maxradians, + rgridsz, mySD)`) must be updated to keyword syntax + (`rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD=mySD)`). + +## Non-breaking changes + +- `qsmooth(img)` now infers the output element type from the input: `float(eltype(img))` + instead of always returning `Float32`. To keep `Float32` output, pass it explicitly: + `qsmooth(Float32, img)`. +- `default_minrot` now accepts an `AbstractArray` directly in addition to + `CartesianIndices`. +- The old `qd_rigid` and `qd_affine` deprecated stubs (which always threw an error) + have been removed; mismatched calls now produce a clearer `MethodError`. + +--- + +# version 0.3 ## Breaking changes diff --git a/Project.toml b/Project.toml index aa6da12..dbc9013 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "RegisterQD" uuid = "ac24ea0c-1830-11e9-18d4-81f172323054" -version = "0.3.3" +version = "1.0.0" [deps] CenterIndexedArrays = "46a7138f-0d70-54e1-8ada-fb8296f91f24" @@ -22,7 +22,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] AxisArrays = "0.3, 0.4" -CenterIndexedArrays = "0.2" +CenterIndexedArrays = "1" CoordinateTransformations = "0.5, 0.6" Distributions = "0.20, 0.21, 0.22, 0.23, 0.24, 0.25" ImageCore = "0.10" @@ -35,14 +35,14 @@ MappedArrays = "0.2, 0.3, 0.4" OffsetArrays = "0.11, 1" PaddedViews = "0.4, 0.5" QuadDIRECT = "0.1" -RegisterCore = "0.1, 0.2" -RegisterDeformation = "0.3, 0.4" -RegisterMismatch = "0.3, 0.4" +RegisterCore = "1" +RegisterDeformation = "1" +RegisterMismatch = "1" Rotations = "0.12, 0.13, 1" StaticArrays = "0.11, 0.12, 1" TestImages = "0.5, 0.6, 1" Unitful = "0.17, 0.18, 1" -julia = "1.6" +julia = "1.10" [extras] AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" diff --git a/src/RegisterQD.jl b/src/RegisterQD.jl index c4db561..90ee40d 100644 --- a/src/RegisterQD.jl +++ b/src/RegisterQD.jl @@ -29,21 +29,4 @@ export qd_translate, getSD, qsmooth -# Deprecations -function qd_rigid(fixed, moving, mxshift::VecLike, mxrot::Union{Number,VecLike}, minwidth_rot::VecLike, SD::AbstractMatrix=I; kwargs...) - error(""" - `qd_rigid` has a new syntax, see the help (`?qd_rigid`) and `NEWS.md`. - """) -end - -function qd_affine(fixed, moving, mxshift, linmins, linmaxs, SD; - thresh=0.5*sum(abs2.(fixed[.!(isnan.(fixed))])), - initial_tfm=IdentityTransformation(), - print_interval=100, - kwargs...) - error(""" - `qd_affine` has a new syntax, see the help (`?qd_affine`) and `NEWS.md`. - """) -end - end # module diff --git a/src/affine.jl b/src/affine.jl index d6d5224..871513f 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -69,8 +69,8 @@ end function affine_mm_slow(params, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation()) tfm = arrayscale(aff(params, moving, initial_tfm), SD) moving, fixed = warp_and_intersect(moving, fixed, tfm) - mm = mismatch0(fixed, moving; normalization=:intensity) - return ratio(mm, thresh, Inf) + mm = mismatch_zeroshift(fixed, moving; normalization=:intensity) + return ratio(mm, thresh; fillval=Inf) end """ @@ -113,6 +113,12 @@ Compute an affine transformation `tfm` that minimizes the mismatch between This only allows small translations (just a couple of pixels). As a consequence, any large translations must be supplied with reasonable accuracy in `initial_tfm` (see [`RegisterQD.qd_affine_coarse`](@ref)). + +`minwidth_mat` optionally specifies the lower limit of resolution for the *linear-map* +parameters only; the translation resolution is fixed internally. The `_mat` suffix +distinguishes this from the full-search-space `minwidth` accepted by +[`RegisterQD.qd_affine_coarse`](@ref): here it names just the linear-map subspace, which +is combined internally with the (fixed) translation resolution. """ function qd_affine_fine(fixed, moving, linmins, linmaxs; SD=I, @@ -180,10 +186,16 @@ If you have a good initial guess at the solution, pass it with the `initial_tfm` Use `SD` if your axes are not uniformly sampled, for example `SD = diagm(voxelspacing)` where `voxelspacing` is a vector encoding the spacing along all axes of the image. `thresh` enforces a certain amount of sum-of-squared-intensity overlap between the two images; with non-zero `thresh`, it is not permissible to "align" the images by shifting one entirely out of the way of the other. +The default value for `thresh` is 50% of the sum-of-squared-intensity of `fixed`. This is higher than the 10% default +used by `qd_translate` and `qd_rigid` because affine transformations have additional degrees of freedom (scaling and shear) +that make degenerate low-overlap solutions more likely. """ function qd_affine(fixed, moving, mxshift, linmins, linmaxs; presmoothed=false, SD=I, + # Affine's extra degrees of freedom (scaling/shear) make degenerate + # low-overlap solutions more likely, so the default requires more overlap + # (50%) than the 10% used by qd_translate/qd_rigid. thresh=0.5*sum(abs2.(fixed[.!(isnan.(fixed))])), initial_tfm=IdentityTransformation(), print_interval=100, diff --git a/src/gridsearch.jl b/src/gridsearch.jl index 583d40b..078be26 100644 --- a/src/gridsearch.jl +++ b/src/gridsearch.jl @@ -53,7 +53,7 @@ function grid_rotations(maxradians, rgridsz, SD) end """ -`best_tform, best_mm = rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD =Matrix{Float64}(I,ndims(fixed),ndims(fixed))))` +`best_tform, best_mm = rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD=Matrix{Float64}(I,ndims(fixed),ndims(fixed)))` Tries a grid of rotations to align `moving` to `fixed`. Also calculates the best translation (`maxshift` pixels or less) to align the images after performing the rotation. Returns an AffineMap that captures both the best rotation and shift out of the values searched, along with the mismatch value after applying that transform (`best_mm`). @@ -61,7 +61,7 @@ best rotation and shift out of the values searched, along with the mismatch valu For more on how the arguments `maxradians`, `rgridsz`, and `SD` influence the search, see the documentation for `grid_rotations`. """ -function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD = Matrix{Float64}(I,ndims(fixed),ndims(fixed))) +function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD = Matrix{Float64}(I,ndims(fixed),ndims(fixed))) rgridsz = [rgridsz...] nd = ndims(moving) @assert nd == ndims(fixed) @@ -75,7 +75,7 @@ function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD = #mm = mismatch(fixed, new_moving, maxshift; normalization=:pixels) mm = mismatch(fixed, new_moving, maxshift) thresh = 0.1*maximum(x->x.denom, mm) - best_i = indmin_mismatch(mm, thresh) + best_i = argmin_mismatch(mm, thresh) cur_best =ratio(mm[best_i], 0.0) if cur_best < best_mm best_mm = cur_best diff --git a/src/rigid.jl b/src/rigid.jl index fc6ccd8..9965761 100644 --- a/src/rigid.jl +++ b/src/rigid.jl @@ -79,8 +79,8 @@ end function rigid_mm_slow(params, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation()) tfm = arrayscale(initial_tfm ∘ tfmrigid(params, moving), SD) moving, fixed = warp_and_intersect(moving, fixed, tfm) - mm = mismatch0(fixed, moving; normalization=:intensity) - return ratio(mm, thresh, Inf) + mm = mismatch_zeroshift(fixed, moving; normalization=:intensity) + return ratio(mm, thresh; fillval=Inf) end function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot; @@ -137,8 +137,12 @@ as a vector or tuple. `mxrot` is the maximum-allowed rotation, in radians for 2d or [quaternion-units](https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation) for 3d. See [`RegisterQD.rot`](@ref) for more information. -`minwidth_rot` optionally specifies the lower limit of resolution for the rotation; -the default is a rotation that moves corner elements by 0.1 pixel. +`minwidth_rot` optionally specifies the lower limit of resolution for the *rotation* +parameters only; the translation resolution is fixed internally and is not user-adjustable. +The default is a rotation that moves corner elements by 0.1 pixel. The `_rot` suffix +distinguishes this from the full-search-space `minwidth` accepted by [`qd_translate`](@ref): +here it names just the rotation subspace, which is combined internally with the +(fixed) translation resolution. `kwargs...` can include any keyword argument that can be passed to `QuadDIRECT.analyze`. It's recommended that you pass your own stopping criteria when possible diff --git a/src/translations.jl b/src/translations.jl index 8ccf7e2..219d3f7 100644 --- a/src/translations.jl +++ b/src/translations.jl @@ -9,8 +9,8 @@ end function translate_mm_slow(params, fixed, moving, thresh; initial_tfm=IdentityTransformation()) tfm = initial_tfm ∘ tfmshift(params, moving) moving, fixed = warp_and_intersect(moving, fixed, tfm) - mm = mismatch0(fixed, moving; normalization=:intensity) - return ratio(mm, thresh, Inf) + mm = mismatch_zeroshift(fixed, moving; normalization=:intensity) + return ratio(mm, thresh; fillval=Inf) end function qd_translate_fine(fixed, moving; @@ -49,6 +49,7 @@ Do not smooth `moving`. If you have a good initial guess at the solution, pass it with the `initial_tfm` kwarg to jump-start the search. `thresh` enforces a certain amount of sum-of-squared-intensity overlap between the two images; with non-zero `thresh`, it is not permissible to "align" the images by shifting one entirely out of the way of the other. +The default value for `thresh` is 10% of the sum-of-squared-intensity of `fixed`. If the `crop` keyword arg is `true` then `fixed` is cropped by `mxshift` (after the optional `initial_tfm`) on all sides so that there will be complete overlap between `fixed` and `moving` for any evaluated shift. This avoids edge effects diff --git a/src/util.jl b/src/util.jl index 75d5a77..5833615 100644 --- a/src/util.jl +++ b/src/util.jl @@ -27,9 +27,9 @@ One can compute an overall transformation by composing `initial_tfm` with the re function best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=IdentityTransformation()) moving, fixed = warp_and_intersect(moving, fixed, initial_tfm) mms = mismatch(fixed, moving, mxshift; normalization=normalization) - best_i = indmin_mismatch(mms, thresh) + best_i = argmin_mismatch(mms, thresh) mm = mms[best_i] - return best_i.I, ratio(mm, thresh, typemax(eltype(mm))) + return best_i.I, ratio(mm, thresh; fillval=typemax(eltype(mm))) end """ @@ -116,9 +116,11 @@ end """ θ = default_minrot(ci::CartesianIndices, SD=I; Δc=0.1) + θ = default_minrot(img::AbstractArray, SD=I; Δc=0.1) Compute the rotation `θ` that results in largest change in coordinates -of size `Δc` (in pixels) for any index in `ci`. +of size `Δc` (in pixels) for any index in `ci`. When passed an array `img`, +its `CartesianIndices` are used. """ function default_minrot(ci::CartesianIndices, SD=I; Δc=0.1) L = -Inf @@ -137,6 +139,8 @@ function default_minrot(ci::CartesianIndices, SD=I; Δc=0.1) return 2*asin(ℓ/(2*L)) end +default_minrot(img::AbstractArray, SD=I; Δc=0.1) = default_minrot(CartesianIndices(img), SD; Δc) + default_minwidth_rot(ci::CartesianIndices{2}, SD=I; kwargs...) = [default_minrot(ci, SD; kwargs...)] function default_minwidth_rot(ci::CartesianIndices{3}, SD=I; kwargs...) @@ -188,13 +192,13 @@ Create a smoothed version of `img`, smoothing with the kernel of a quadratic B-s Use this on your `fixed` image in preparation for registration, and pass `presmoothed` as an option. (Do not smooth `moving`.) -`T` allows you to specify the output eltype (default `Float32`). +`T` allows you to specify the output eltype (default `float(eltype(img))`). """ function qsmooth(::Type{T}, img::AbstractArray{T2,N}) where {T,N,T2} kern1 = centered(T[1/8, 3/4, 1/8]) # quadratic B-spline kernel return imfilter(img, kernelfactors(ntuple(i->kern1, Val(N)))) end -qsmooth(img::AbstractArray) = qsmooth(Float32, img) +qsmooth(img::AbstractArray) = qsmooth(float(eltype(img)), img) function qinterp(::Type{T}, moving) where T widen1(ax) = first(ax)-1:last(ax)+1 diff --git a/test/gridsearch.jl b/test/gridsearch.jl index d157f5f..83dcb29 100644 --- a/test/gridsearch.jl +++ b/test/gridsearch.jl @@ -1,6 +1,7 @@ using RegisterQD using RegisterQD.CoordinateTransformations using RegisterQD.RegisterDeformation +using LinearAlgebra: I using Test @testset "Grid search rigid registration" begin @@ -10,7 +11,9 @@ using Test b = transform(a, tformtranslate([2.0;0.0]) ∘ tformrotate(pi/6)) tfm0 = tformtranslate([-2.0;0.0]) ∘ tformrotate(-pi/6) #note: maxshift must be GREATER than the true shift in order to find the true shift - tfm, mm = RegisterQD.rotation_gridsearch(a, b, (11,11), [pi/6], [11]) + # SD is now a keyword argument (was positional in v0.x) + SD2 = Matrix{Float64}(I, 2, 2) + tfm, mm = RegisterQD.rotation_gridsearch(a, b, (11,11), [pi/6], [11]; SD=SD2) @test tfm.translation == tfm0.translation @test tfm.linear == tfm0.linear diff --git a/test/util.jl b/test/util.jl index af0140a..0fa11f6 100644 --- a/test/util.jl +++ b/test/util.jl @@ -22,6 +22,11 @@ using Unitful: μm, mm, cm, km, s ci = CartesianIndices(img) θ = RegisterQD.default_minrot(ci) @test θ ≈ 0.1/sqrt(3^2 + 10^2 + 5^2) rtol=1e-3 + # array overload matches the CartesianIndices form + SD = [1 0 0; 0 2 0; 0 0 3] + @test RegisterQD.default_minrot(img) == RegisterQD.default_minrot(CartesianIndices(img)) + @test RegisterQD.default_minrot(img, SD) == RegisterQD.default_minrot(CartesianIndices(img), SD) + @test RegisterQD.default_minrot(img, SD; Δc=0.2) == RegisterQD.default_minrot(CartesianIndices(img), SD; Δc=0.2) end @testset "getSD" begin @@ -58,7 +63,14 @@ end @test size(getSD(timedarray)) == (2,2) end - +@testset "qsmooth eltype" begin + img64 = rand(Float64, 8, 8) + @test eltype(qsmooth(img64)) === Float64 + img32 = rand(Float32, 8, 8) + @test eltype(qsmooth(img32)) === Float32 + # explicit T still widens the compute/output eltype + @test eltype(qsmooth(Float64, img32)) === Float64 +end #TODO add a testset for other support functions #rotations