RegisterQD v1.0.0 API review

.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

.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 }}

.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

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
 

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"

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

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,

src/gridsearch.jl

@@ -53,15 +53,15 @@ 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`).
 
 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

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

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

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

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
 

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