diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status deleted file mode 100644 index 5647ec9..0000000 --- a/.claude/freshen-package-status +++ /dev/null @@ -1,11 +0,0 @@ -DONE: design review -DONE: API review -DONE: update .gitignore -DONE: 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/Documenter.yml b/.github/workflows/Documenter.yml new file mode 100644 index 0000000..6a8068e --- /dev/null +++ b/.github/workflows/Documenter.yml @@ -0,0 +1,26 @@ +name: Documenter +on: + push: + branches: + - master + tags: '*' + pull_request: + +jobs: + build: + permissions: + contents: write + statuses: write + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: '1' + - name: Add HolyLab registry + run: julia -e 'using Pkg; pkg"registry add General https://github.com/HolyLab/HolyLabRegistry.git"' + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.gitignore b/.gitignore index 6991c1c..56b88fc 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,5 @@ Manifest.toml Manifest-v*.toml test/Manifest.toml .vscode/ +docs/build/ +docs/Manifest.toml diff --git a/Project.toml b/Project.toml index dbc9013..abf51da 100644 --- a/Project.toml +++ b/Project.toml @@ -21,32 +21,42 @@ Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] +Aqua = "0.8" AxisArrays = "0.3, 0.4" CenterIndexedArrays = "1" CoordinateTransformations = "0.5, 0.6" Distributions = "0.20, 0.21, 0.22, 0.23, 0.24, 0.25" +Documenter = "1" +ExplicitImports = "1" ImageCore = "0.10" ImageFiltering = "0.7" ImageMagick = "0.7, 1" ImageMetadata = "0.9" ImageTransformations = "0.10" Interpolations = "0.12, 0.13, 0.14, 0.15" +LinearAlgebra = "1" MappedArrays = "0.2, 0.3, 0.4" OffsetArrays = "0.11, 1" PaddedViews = "0.4, 0.5" QuadDIRECT = "0.1" +Random = "1" RegisterCore = "1" RegisterDeformation = "1" RegisterMismatch = "1" +RegisterMismatchCommon = "1" Rotations = "0.12, 0.13, 1" StaticArrays = "0.11, 0.12, 1" +Test = "1" TestImages = "0.5, 0.6, 1" Unitful = "0.17, 0.18, 1" julia = "1.10" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -57,4 +67,4 @@ TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["Test", "ImageMagick", "TestImages", "Random", "AxisArrays", "ImageMetadata", "Unitful", "Distributions", "LinearAlgebra", "RegisterMismatch"] +test = ["Aqua", "Documenter", "ExplicitImports", "Test", "ImageMagick", "TestImages", "Random", "AxisArrays", "ImageMetadata", "Unitful", "Distributions", "LinearAlgebra", "RegisterMismatch"] diff --git a/README.md b/README.md index 9ed7b5c..4de5e0f 100644 --- a/README.md +++ b/README.md @@ -1,23 +1,60 @@ # RegisterQD -[![CI](https://github.com/HolyLab/RegisterQD.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/RegisterQD.jl/actions/workflows/CI.yml)[![codecov](https://codecov.io/gh/HolyLab/RegisterQD.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HolyLab/RegisterQD.jl) +[![CI](https://github.com/HolyLab/RegisterQD.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/RegisterQD.jl/actions/workflows/CI.yml) +[![codecov](https://codecov.io/gh/HolyLab/RegisterQD.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HolyLab/RegisterQD.jl) +[![Stable docs](https://img.shields.io/badge/docs-stable-blue.svg)](https://HolyLab.github.io/RegisterQD.jl/stable/) +[![Dev docs](https://img.shields.io/badge/docs-dev-blue.svg)](https://HolyLab.github.io/RegisterQD.jl/dev/) +[![Aqua QA](https://juliatesting.github.io/Aqua.jl/dev/assets/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) RegisterQD performs image registration using the global optimization routine [QuadDIRECT](https://github.com/timholy/QuadDIRECT.jl). -Unlike many other registration packages, this is not "greedy" descent based on an initial guess---it attempts to find the globally-optimal alignment of your images. +Unlike many other registration packages, this is not "greedy" descent based on an initial guess — it attempts to find the globally-optimal alignment of your images. -To use this package, users must choose between using either the CPU or the GPU. For CPU processing, you must manually load the [RegisterMismatch package](https://github.com/HolyLab/RegisterMismatch.jl): `using RegisterMismatch, RegisterQD`. For GPU processing, you should instead load the [RegisterMismatchCuda package](https://github.com/HolyLab/RegisterMismatchCuda.jl): `using RegisterMismatchCuda, RegisterQD`. *Note that loading both mismatch packages in the same session will cause method conflicts.* Both mismatch packages are registered in the publicly-available [HolyLabRegistry](https://github.com/HolyLab/HolyLabRegistry), and users are advised to add that registry. -In the current absense of Github resources for GPU code, "gpu_test.jl" should be run on your personal machine as required. +## Installation + +RegisterQD and its dependencies live in the [HolyLab registry](https://github.com/HolyLab/HolyLabRegistry). +Add the registry once, then install: + +```julia +using Pkg +pkg"registry add https://github.com/HolyLab/HolyLabRegistry.git" +Pkg.add("RegisterQD") +``` + +You also need a mismatch backend. For CPU processing, load [RegisterMismatch](https://github.com/HolyLab/RegisterMismatch.jl): + +```julia +Pkg.add("RegisterMismatch") +``` + +For GPU processing, use [RegisterMismatchCuda](https://github.com/HolyLab/RegisterMismatchCuda.jl) instead. +*Do not load both in the same session — they conflict.* + +## Quick start + +```julia +using RegisterMismatch, RegisterQD + +fixed = Float64.(reshape(1:25, 5, 5)) +moving = circshift(fixed, (2, 1)) # known shift: 2 rows, 1 column + +tform, mm = qd_translate(fixed, moving, (3, 3)) +# tform.translation == [2.0, 1.0] +# mm == 0.0 +``` + +## Registration functions -This package exports the following registration functions: - `qd_translate`: register images by shifting one with respect to another (translations only) - `qd_rigid`: register images using rotations and translations - `qd_affine`: register images using arbitrary affine transformations -In general, using more degrees of freedom allows you to solve harder optimization problems, but also makes it harder to find the global optimum. Your best strategy is to permit no more degrees of freedom than needed to solve the problem. +In general, using more degrees of freedom allows you to solve harder optimization problems, but also makes it harder to find the global optimum. +Use no more degrees of freedom than your problem requires. -See the help on these functions for details about how to call them. +## Anisotropic sampling -Another important feature of this package is that it supports images that were sampled anisotropically. This is particularly common for three-dimensional biomedical imaging, where MRI and optical microscopy typically have one axis sampled at lower resolution. -A rotation (from a rigid transformation) in physical space needs to be modified before applying it to an anisotropically-sampled image; see `arrayscale` and `getSD` for more information. +This package supports images sampled anisotropically, which is common in 3-D biomedical imaging (e.g. MRI, optical sections where the axial resolution differs from the in-plane resolution). +Pass `SD = diagm(voxelspacing)` to the registration functions to account for non-uniform spacing. +See [`arrayscale`](https://HolyLab.github.io/RegisterQD.jl/stable/api/#RegisterQD.arrayscale) and [`getSD`](https://HolyLab.github.io/RegisterQD.jl/stable/api/#RegisterQD.getSD) for details, and the [User Guide](https://HolyLab.github.io/RegisterQD.jl/stable/guide/) for a full explanation. -**NOTE**: see NEWS.md for information about a recent breaking change. +**NOTE**: see NEWS.md for information about recent breaking changes. diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..bd35d4a --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,7 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +RegisterQD = "ac24ea0c-1830-11e9-18d4-81f172323054" + +[compat] +Documenter = "1" +julia = "1.10" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..168b4ee --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,25 @@ +using RegisterQD +using Documenter + +makedocs(; + modules = [RegisterQD], + authors = "HolyLab", + sitename = "RegisterQD.jl", + format = Documenter.HTML(; + prettyurls = get(ENV, "CI", "false") == "true", + canonical = "https://HolyLab.github.io/RegisterQD.jl", + ), + pages = [ + "Home" => "index.md", + "User Guide" => "guide.md", + "API Reference" => "api.md", + ], + checkdocs = :exports, + doctest = :none, # doctests are run in the test suite via doctest(RegisterQD; manual=false) +) + +deploydocs(; + repo = "github.com/HolyLab/RegisterQD.jl", + devbranch = "master", + push_preview = true, +) diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 0000000..ab8de2a --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,24 @@ +# API Reference + +## Registration functions + +```@docs +qd_translate +qd_rigid +qd_affine +``` + +## Rotation grid search + +```@docs +grid_rotations +rotation_gridsearch +``` + +## Utilities + +```@docs +arrayscale +getSD +qsmooth +``` diff --git a/docs/src/guide.md b/docs/src/guide.md new file mode 100644 index 0000000..d1d3808 --- /dev/null +++ b/docs/src/guide.md @@ -0,0 +1,140 @@ +# User Guide + +## How registration works + +Image registration is the problem of finding a geometric transformation that +aligns one image (`moving`) to another (`fixed`). RegisterQD measures how well +two images are aligned using a *mismatch* value — roughly, the sum-of-squared +differences after normalising for intensity — and minimises it with the +[QuadDIRECT](https://github.com/timholy/QuadDIRECT.jl) global optimiser. + +Because QuadDIRECT searches the entire parameter space rather than following a +local gradient, it can find the correct alignment even without a good starting +guess. The trade-off is that it is slower than greedy methods; providing a good +`initial_tfm` when one is available can speed things up significantly. + +## Mismatch backend + +The mismatch computation is implemented in a separate *backend* package that must +be loaded before calling any registration function: + +| Use case | Package to load | +|----------|----------------| +| CPU | [RegisterMismatch](https://github.comlyHolyLab/RegisterMismatch.jl) | +| GPU (CUDA) | [RegisterMismatchCuda](https://github.com/HolyLab/RegisterMismatchCuda.jl) | + +Load exactly one backend — loading both in the same session causes method +conflicts. + +## Choosing a registration function + +| Function | Degrees of freedom | When to use | +|---|---|---| +| [`qd_translate`](@ref) | N translations | images are related by a pure shift | +| [`qd_rigid`](@ref) | N translations + rotation(s) | images share the same shape but may be rotated | +| [`qd_affine`](@ref) | N translations + N² linear | images may have scaling or shear in addition to rotation | + +Start with the fewest degrees of freedom that your problem requires. +More degrees of freedom make the global search harder and increase the risk of +finding a degenerate (low-overlap) solution. + +## Anisotropic sampling and the `SD` parameter + +Three-dimensional biomedical images (MRI, optical sections) are often sampled at +different resolutions along different axes — for example, 0.5 mm in-plane but +2 mm axially. A physical-space rotation of such a volume does not look like a +rotation in array-index space. + +The `SD` ("spatial displacements") parameter encodes the physical spacing. Its +columns are the physical displacements corresponding to one array-element step +along each axis: + +```julia +# 2-D example: dim 2 sampled 4× coarser than dim 1 +SD = [1.0 0.0; + 0.0 4.0] + +tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD) +``` + +If your image carries axis metadata (e.g. from +[AxisArrays.jl](https://github.com/JuliaArrays/AxisArrays.jl)), use +[`getSD`](@ref) to extract `SD` automatically: + +```julia +SD = getSD(fixed) +tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD) +``` + +The returned transformation is always in *physical* coordinates. Before warping +the array, convert it with [`arrayscale`](@ref): + +```julia +itfm = arrayscale(tform, SD) +warped = warp(moving, itfm) +``` + +## Centre of rotation + +Rotations are defined around the *origin of coordinates*, which for a plain +Julia array is index `(1,1,…)`. For most images, rotation around the image +centre is more natural. Call `centered` (from ImageTransformations) on both +images to shift the origin to the centre before registering: + +```julia +using ImageTransformations: centered + +cfixed = centered(fixed) +cmoving = centered(moving) + +tform, mm = qd_rigid(collect(cfixed), collect(float(cmoving)), mxshift, mxrot) +``` + +Remember to call `centered` again when applying the result to a new image, or +re-encode the transformation with a different origin using +`recenter(tform, newctr)`. + +## Initial guess + +If you already have an approximate transformation, pass it as `initial_tfm` to +warm-start the search: + +```julia +guess = Translation(1.0, 0.0) +tform, mm = qd_translate(fixed, moving, mxshift; initial_tfm=guess) +``` + +## Pre-smoothing with `qsmooth` + +[`qsmooth`](@ref) smooths an image with a quadratic B-spline kernel and returns +an interpolant that can be warped cheaply. Pre-smoothing `fixed` and passing +`presmoothed=true` reduces the mismatch evaluation cost: + +```julia +fixed_s = qsmooth(fixed) +tform, mm = qd_translate(fixed_s, moving, mxshift; presmoothed=true) +``` + +Do **not** smooth `moving`. + +## Overlap threshold (`thresh`) + +`thresh` prevents the optimiser from "aligning" images by sliding one entirely +out of view. The default is 10 % of the sum-of-squared intensity of `fixed` +for translations and rigid transforms, and 50 % for affine transforms (because +the extra degrees of freedom make degenerate solutions more tempting). Increase +`thresh` if you see unexpected results. + +## Coarse rotation search + +For large rotations where QuadDIRECT may not converge, use +[`rotation_gridsearch`](@ref) to identify a good starting rotation, then refine +with [`qd_rigid`](@ref): + +```julia +best_rot, _ = rotation_gridsearch(fixed, moving, mxshift, maxradians, gridsz) +tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; initial_tfm=best_rot) +``` + +[`grid_rotations`](@ref) generates the grid of candidate rotations used +internally by `rotation_gridsearch` and can also be used standalone. diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..1a01694 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,62 @@ +# RegisterQD.jl + +```@docs +RegisterQD.RegisterQD +``` + +RegisterQD performs image registration using the global optimization routine +[QuadDIRECT](https://github.com/timholy/QuadDIRECT.jl). +Unlike greedy descent methods, it searches for the globally-optimal alignment, +making it robust to poor initial guesses. + +## Installation + +RegisterQD and its dependencies live in the +[HolyLab registry](https://github.com/HolyLab/HolyLabRegistry). +Add the registry once, then install the package: + +```julia +using Pkg +pkg"registry add https://github.com/HolyLab/HolyLabRegistry.git" +Pkg.add("RegisterQD") +``` + +RegisterQD also requires a *mismatch backend* to be loaded before calling any +registration function. For CPU processing, install and load +[RegisterMismatch](https://github.com/HolyLab/RegisterMismatch.jl): + +```julia +Pkg.add("RegisterMismatch") +``` + +For GPU processing, use +[RegisterMismatchCuda](https://github.com/HolyLab/RegisterMismatchCuda.jl) instead. +Do **not** load both in the same session — they conflict. + +## Quick start + +```julia +using RegisterMismatch, RegisterQD + +fixed = Float64.(reshape(1:25, 5, 5)) +moving = circshift(fixed, (2, 1)) # known shift: 2 rows, 1 column + +tform, mm = qd_translate(fixed, moving, (3, 3)) +# tform.translation ≈ [2.0, 1.0] +# mm ≈ 0.0 +``` + +For a rotation search: + +```julia +using RegisterMismatch, RegisterQD, CoordinateTransformations, Rotations, ImageTransformations + +fixed = Float64.(reshape(1:100, 10, 10)) +moving = warp(centered(fixed), LinearMap(RotMatrix(0.1))) + +tform, mm = qd_rigid(collect(centered(fixed)), collect(float(moving)), (2,2), (0.3,)) +# mm ≈ 0.0 (rotation recovered) +``` + +See the [User Guide](@ref) for a full explanation of concepts, and the +[API Reference](@ref) for all exported functions. diff --git a/src/RegisterQD.jl b/src/RegisterQD.jl index 3900329..98b91ff 100644 --- a/src/RegisterQD.jl +++ b/src/RegisterQD.jl @@ -1,16 +1,48 @@ +""" + RegisterQD + +Image registration using the [QuadDIRECT](https://github.com/timholy/QuadDIRECT.jl) +global optimization algorithm. + +The three main entry points are: +- [`qd_translate`](@ref): optimize a pure translation +- [`qd_rigid`](@ref): optimize a rigid transformation (rotation + translation) +- [`qd_affine`](@ref): optimize a full affine transformation + +All three return `(tform, mm)` where `tform` is a +[CoordinateTransformations.jl](https://github.com/JuliaGeometry/CoordinateTransformations.jl) +transform object and `mm` is the residual mismatch value (lower is better). + +!!! note + A mismatch backend such as + [RegisterMismatch.jl](https://github.com/HolyLab/RegisterMismatch.jl) must be + loaded (`using RegisterMismatch`) before calling any registration function. + +# Utilities +- [`arrayscale`](@ref): convert a physical-space transform to array-index space +- [`getSD`](@ref): extract the spatial-directions matrix from an annotated image +- [`qsmooth`](@ref): pre-smooth an image for registration +- [`grid_rotations`](@ref) / [`rotation_gridsearch`](@ref): coarse rotation grid search +""" module RegisterQD -using ImageCore, ImageTransformations, ImageFiltering -using CoordinateTransformations -using QuadDIRECT -using RegisterMismatchCommon -using RegisterCore -using RegisterDeformation, PaddedViews, MappedArrays -using Rotations -using Interpolations, CenterIndexedArrays, StaticArrays, OffsetArrays -using LinearAlgebra - -using ImageTransformations: CornerIterator +using CenterIndexedArrays: CenterIndexedArrays +using CoordinateTransformations: CoordinateTransformations, AbstractAffineMap, AffineMap, + IdentityTransformation, LinearMap, Translation +using ImageCore: ImageCore, spacedirections +using ImageFiltering: ImageFiltering, centered, imfilter, kernelfactors +using ImageTransformations: ImageTransformations, warp +using Interpolations: Interpolations, BSpline, Free, OnCell, Quadratic, extrapolate +using LinearAlgebra: LinearAlgebra, I, UniformScaling, eigen, norm +using MappedArrays: MappedArrays, of_eltype +using OffsetArrays: OffsetArrays, OffsetArray +using PaddedViews: PaddedViews, PaddedView +using QuadDIRECT: QuadDIRECT, value +using RegisterCore: RegisterCore, argmin_mismatch, ratio +using RegisterDeformation: RegisterDeformation, tformeye, tformrotate, tformtranslate, transform +using RegisterMismatchCommon: RegisterMismatchCommon, mismatch, mismatch_zeroshift +using Rotations: Rotations, RotMatrix, isrotation +using StaticArrays: StaticArrays, @SMatrix, SMatrix, SVector const VecLike = Union{AbstractVector{<:Number}, Tuple{Number, Vararg{Number}}} diff --git a/src/affine.jl b/src/affine.jl index 1e9a954..f30b32b 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -160,43 +160,76 @@ end # You supply one or the other, so I don't see a problem. #TODO I think that this is a tad loquatious, and not enough examples. Permission to rework it? """ -`tform, mm = qd_affine(fixed, moving, mxshift, linmins, linmaxs, SD=I; presmoothed=false, thresh, initial_tfm, kwargs...)` -`tform, mm = qd_affine(fixed, moving, mxshift, SD=I; presmoothed=false, thresh, initial_tfm, kwargs...)` -optimizes an affine transformation (linear map + translation) to minimize the mismatch between `fixed` and -`moving` using the QuadDIRECT algorithm. The algorithm is run twice: the first step samples the search space -at a coarser resolution than the second. `kwargs...` may contain any keyword argument that can be passed to -`QuadDIRECT.analyze`. It's recommended that you pass your own stopping criteria when possible -(i.e. `rtol`, `atol`, and/or `fvalue`). If you provide `rtol` and/or `atol` they will apply only to the -second (fine) step of the registration; the user may not adjust these criteria for the coarse step. - -`tform` will be centered on the origin-of-coordinates, i.e. (0,0) for a 2D image. Usually it is more natural to consider rotations -around the center of the image. If you would like `mxrot` and the returned rotation to act relative to the center of the image, then you must -move the origin to the center of the image by calling `centered(img)` from the `ImageTransformations` package. Call `centered` on both the -fixed and moving image to generate the `fixed` and `moving` that you provide as arguments. If you later want to apply the returned transform -to an image you must remember to call `centered` on that image as well. Alternatively you can re-encode the transformation in terms of a -different origin by calling `recenter(tform, newctr)` where `newctr` is the displacement of the new center from the old center. - -The `linmins` and `linmaxs` arguments set the minimum and maximum allowable values in the linear map matrix. -They can be supplied as NxN matrices or flattened vectors. If omitted then a modest default search space is chosen. -`mxshift` sets the magnitude of the largest allowable translation in each dimension (It's a vector of length N). -This default search-space allows for very little rotation. -Alternatively, you can submit `dmax` or `ndmax` values as keyword functions, which will use diagonal or non-diagonal variation from the identity matrix -to generate less modest `linmins` and `linmaxs` arguments for you. - -Use `presmoothed=true` if you have called [`qsmooth`](@ref) on `fixed` before calling `qd_affine`. -Do not smooth `moving`. - -`kwargs...` can also include any other keyword argument that can be passed to `QuadDIRECT.analyze`. -It's recommended that you pass your own stopping criteria when possible (i.e. `rtol`, `atol`, and/or `fvalue`). - -If you have a good initial guess at the solution, pass it with the `initial_tfm` kwarg to jump-start the search. - -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. + tform, mm = qd_affine(fixed, moving, mxshift, linmins, linmaxs; + presmoothed=false, SD=I, + thresh=0.5*sum(abs2,fixed), initial_tfm=IdentityTransformation(), + kwargs...) + tform, mm = qd_affine(fixed, moving, mxshift; + dmax=0.05, ndmax=0.05, + presmoothed=false, SD=I, + thresh=0.5*sum(abs2,fixed), initial_tfm=IdentityTransformation(), + kwargs...) + +Optimize an affine transformation (linear map + translation) to minimize the mismatch +between `fixed` and `moving` using the QuadDIRECT algorithm. + +Returns `(tform, mm)` where `tform` is an `AffineMap` and `mm` is the residual mismatch +value (lower is better). + +The algorithm runs in two steps: the first step samples the search space at a coarser +resolution than the second. `kwargs...` may contain any keyword argument accepted by +`QuadDIRECT.analyze`. Supplying your own stopping criteria (`rtol`, `atol`, and/or +`fvalue`) is recommended. Any `rtol`/`atol` you supply will apply only to the second +(fine) step; the coarse step uses fixed internal criteria. + +`tform` is centered on the origin of coordinates, i.e. `(0, 0)` for 2D images. To +rotate around the image center instead, call `centered(img)` (from `ImageFiltering` or +`ImageTransformations`) on both `fixed` and `moving` before calling `qd_affine`. +To re-encode the result relative to a different center, use +`recenter(tform, newctr)`. + +`linmins` and `linmaxs` bound the allowable values of the linear-map matrix entries. +They can be `N×N` matrices or flattened vectors. If omitted, a modest default search +space is used, controllable via the `dmax` (diagonal) and `ndmax` (off-diagonal) keyword +arguments; e.g. `dmax=0.05` permits diagonal entries in `[0.95, 1.05]`. + +`mxshift` sets the maximum allowable translation in each dimension. + +Use `presmoothed=true` if you have called [`qsmooth`](@ref) on `fixed` before calling +`qd_affine`. Do not smooth `moving`. + +If you have a good initial guess, pass it with `initial_tfm` to jump-start the search. + +Use `SD` if your axes are not uniformly sampled, for example +`SD = diagm(voxelspacing)` where `voxelspacing` encodes the physical spacing along each +axis. See [`arrayscale`](@ref) for details. + +`thresh` enforces a minimum sum-of-squared-intensity overlap; with non-zero `thresh`, +it is not permissible to "align" the images by shifting one entirely out of view. +The default is 50% of the sum-of-squared-intensity of `fixed` — higher than the 10% +used by [`qd_translate`](@ref) and [`qd_rigid`](@ref) because affine transformations +have extra degrees of freedom (scaling, shear) that make degenerate low-overlap +solutions more likely. + +!!! note + A mismatch backend such as + [RegisterMismatch.jl](https://github.com/HolyLab/RegisterMismatch.jl) must be + loaded before calling this function. + +# Examples + +```jldoctest +julia> using RegisterMismatch + +julia> fixed = Float64.(reshape(1:25, 5, 5)); + +julia> moving = circshift(fixed, (1, 0)); + +julia> tform, mm = qd_affine(fixed, moving, (3, 3); print_interval=typemax(Int)); + +julia> mm < 1e-6 +true +``` """ function qd_affine( fixed, moving, mxshift, linmins, linmaxs; diff --git a/src/gridsearch.jl b/src/gridsearch.jl index ff75d09..74a4051 100644 --- a/src/gridsearch.jl +++ b/src/gridsearch.jl @@ -1,19 +1,30 @@ """ -`rotations = grid_rotations(maxradians, rgridsz, SD)` generates -a set of rotations (AffineMap) useful for a gridsearch of -possible rotations to align a pair of images. - -`maxradians` is either a single maximum angle (in 2D) or a set of -Euler angles (in 3D and higher). `rgridsz` is one or more integers -specifying the number of gridpoints to search in each of the rotation -axes, corresponding with entries in `maxradians`. `SD` is a matrix -specifying the sample spacing. - -e.g. `grid_rotations((pi/8,pi/8,pi/8), (3,3,3), Matrix{Float64}(I,3,3))` would -return an array of 27 rotations with 3 possible angles for each -Euler axis: -pi/8, 0, pi/8. Passing `Matrix{Float64}(I,3,3)` for SD indicates -that the resulting transforms are meant to be applied to an image with isotropic -pixel spacing. + rotations = grid_rotations(maxradians, rgridsz, SD) + +Generate a `Vector{AffineMap}` of candidate rotation transforms for use in a +coarse rotation grid search. + +`maxradians` is either a single maximum angle (in 2D) or a tuple of maximum +Euler angles (in 3D and higher). `rgridsz` is one or more integers specifying +the number of grid points along each rotation axis; values should be odd so that +zero rotation is included. `SD` is a matrix specifying the sample spacing +(see [`arrayscale`](@ref)). + +# Examples + +```jldoctest +julia> using LinearAlgebra + +julia> rots = grid_rotations(π/8, 3, Matrix{Float64}(I, 2, 2)); + +julia> length(rots) +3 + +julia> rots3d = grid_rotations((π/8, π/8, π/8), (3,3,3), Matrix{Float64}(I, 3, 3)); + +julia> length(rots3d) +27 +``` """ function grid_rotations(maxradians, rgridsz, SD) rgridsz = [rgridsz...] @@ -53,13 +64,41 @@ 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)))` -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`). + best_tform, best_mm = rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; + SD=Matrix{Float64}(I,ndims(fixed),ndims(fixed))) + +Search a grid of rotations to coarsely align `moving` to `fixed`. + +Returns `(best_tform, best_mm)` where `best_tform` is an `AffineMap` encoding the +best-found rotation and translation, and `best_mm` is the corresponding mismatch +value (lower is better). + +For each candidate rotation (generated by [`grid_rotations`](@ref)), the best +integer-pixel translation of up to `maxshift` pixels is found via a Fourier method. +The candidate with the lowest combined mismatch is returned. + +See [`grid_rotations`](@ref) for how `maxradians`, `rgridsz`, and `SD` shape the search. + +!!! note + A mismatch backend such as + [RegisterMismatch.jl](https://github.com/HolyLab/RegisterMismatch.jl) must be + loaded before calling this function. + +# Examples + +```jldoctest +julia> using RegisterMismatch, CoordinateTransformations, Rotations, ImageFiltering, ImageTransformations + +julia> fixed = Float64.(reshape(1:100, 10, 10)); + +julia> moving = warp(centered(fixed), LinearMap(RotMatrix(π/16))); + +julia> best_tform, best_mm = rotation_gridsearch( + collect(centered(fixed)), collect(float(moving)), (3,3), π/8, 3); -For more on how the arguments `maxradians`, `rgridsz`, and `SD` influence the search, see the documentation for -`grid_rotations`. +julia> best_mm < 0.1 +true +``` """ function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD = Matrix{Float64}(I, ndims(fixed), ndims(fixed))) rgridsz = [rgridsz...] diff --git a/src/rigid.jl b/src/rigid.jl index 32309dc..c234662 100644 --- a/src/rigid.jl +++ b/src/rigid.jl @@ -133,7 +133,8 @@ end tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; presmoothed=false, SD=I, minwidth_rot=default_minwidth_rot(fixed, SD), - thresh=thresh, initial_tfm=IdentityTransformation(), kwargs...) + thresh=0.1*sum(abs2,fixed), initial_tfm=IdentityTransformation(), + kwargs...) Optimize a rigid transformation (rotation + shift) to minimize the mismatch between `fixed` and `moving` using the QuadDIRECT algorithm. The algorithm is run twice: the first step finds the optimal rotation, @@ -146,7 +147,7 @@ any regions of `fixed` that you don't want to align against. 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. +See `RegisterQD.rot` for more information. `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 @@ -188,6 +189,27 @@ Do not smooth `moving`. Both the output `tfm` and any `initial_tfm` are represented in *physical* coordinates; as long as `initial_tfm` is a rigid transformation, `tfm` will be a pure rotation+translation. If `SD` is not the identity, use `arrayscale` before applying the result to `moving`. + +!!! note + A mismatch backend such as + [RegisterMismatch.jl](https://github.com/HolyLab/RegisterMismatch.jl) must be + loaded before calling this function. + +# Examples + +```jldoctest +julia> using RegisterMismatch, CoordinateTransformations, Rotations, ImageFiltering, ImageTransformations + +julia> fixed = Float64.(reshape(1:100, 10, 10)); + +julia> moving = warp(centered(fixed), LinearMap(RotMatrix(0.1))); + +julia> tform, mm = qd_rigid(collect(centered(fixed)), collect(float(moving)), (2,2), (0.3,); + print_interval=typemax(Int)); + +julia> mm < 1e-4 +true +``` """ function qd_rigid( fixed, moving, mxshift::VecLike, mxrot::Union{Number, VecLike}; diff --git a/src/translations.jl b/src/translations.jl index cd5616c..20009be 100644 --- a/src/translations.jl +++ b/src/translations.jl @@ -30,32 +30,61 @@ function qd_translate_fine( end """ -`tform, mm = qd_translate(fixed, moving, mxshift; presmoothed=false, thresh=thresh, kwargs...)` -optimizes a simple shift (translation) to minimize the mismatch between `fixed` and -`moving` using the QuadDIRECT algorithm with the constraint that no shifts larger than -``mxshift` (after an optional `initial_tfm`) will be considered. - -Both `mxshift` and the returned translation are specified in terms of pixel units, so the -algorithm need not be aware of anisotropic sampling. - -The algorithm involves two steps: the first step uses a fourier method to speed up -the search for the best whole-pixel shift. The second step refines the search for sub-pixel accuracy. -The default precision of this step is 1% of one pixel (0.01) for each dimension of the image. -You can override the default with the `minwidth` argument. `kwargs...` can also include -any other keyword argument that can be passed to `QuadDIRECT.analyze`. -It's recommended that you pass your own stopping criteria when possible (i.e. `rtol`, `atol`, and/or `fvalue`). - -Use `presmoothed=true` if you have called [`qsmooth`](@ref) on `fixed` before calling `qd_affine`. -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 -that can occur due to normalization when the transformed `moving` doesn't fully overlap with `fixed`. + tform, mm = qd_translate(fixed, moving, mxshift; + presmoothed=false, thresh=0.1*sum(abs2,fixed), kwargs...) + +Optimize a translation to minimize the mismatch between `fixed` and `moving` using the +QuadDIRECT algorithm. No shift larger than `mxshift` (after an optional `initial_tfm`) +will be considered. + +Returns `(tform, mm)` where `tform` is a `Translation` and `mm` is the residual +mismatch value (lower is better). + +Both `mxshift` and the returned translation are in pixel units, so the algorithm does +not need to know the physical sampling. + +The algorithm runs in two steps: the first uses a Fourier method to find the best +whole-pixel shift; the second refines for sub-pixel accuracy with default precision of +1% of one pixel (`minwidth=fill(0.01, ndims(fixed))`). Override with the `minwidth` +keyword argument. `kwargs...` can include any keyword argument accepted by +`QuadDIRECT.analyze`. Supplying your own stopping criteria (`rtol`, `atol`, and/or +`fvalue`) is recommended. + +Use `presmoothed=true` if you have called [`qsmooth`](@ref) on `fixed` before calling +`qd_translate`. Do not smooth `moving`. + +If you have a good initial guess, pass it with `initial_tfm` to jump-start the search. + +`thresh` enforces a minimum sum-of-squared-intensity overlap between the two images; +with non-zero `thresh`, shifting one image entirely out of view is not a valid solution. +The default is 10% of the sum-of-squared-intensity of `fixed`. + +If `crop=true`, `fixed` is cropped by `mxshift` on all sides so that there is complete +overlap between `fixed` and `moving` for every evaluated shift. This avoids edge-effect +normalization artifacts when the transformed `moving` does not fully overlap `fixed`. + +!!! note + A mismatch backend such as + [RegisterMismatch.jl](https://github.com/HolyLab/RegisterMismatch.jl) must be + loaded before calling this function. + +# Examples + +```jldoctest +julia> using RegisterMismatch + +julia> fixed = Float64.(reshape(1:25, 5, 5)); + +julia> moving = circshift(fixed, (2, 1)); # known shift: 2 rows, 1 column + +julia> tform, mm = qd_translate(fixed, moving, (3, 3); print_interval=typemax(Int)); + +julia> println(tform.translation) +[2.0, 1.0] + +julia> mm +0.0 +``` """ function qd_translate( fixed, moving, mxshift; diff --git a/src/util.jl b/src/util.jl index 58d4d30..58b9692 100644 --- a/src/util.jl +++ b/src/util.jl @@ -1,3 +1,6 @@ +_corners(ci::CartesianIndices) = + (CartesianIndex(c) for c in Iterators.product(((first(r), last(r)) for r in ci.indices)...)) + function warp_and_intersect(moving, fixed, tfm::IdentityTransformation) if axes(moving) == axes(fixed) return moving, fixed @@ -35,19 +38,38 @@ end """ itfm = arrayscale(ptfm, SD) -Convert the physical-space transformation `ptfm` into one, `itfm`, that operates on index-space -for arrays. -For example, suppose `ptfm` is a pure 3d rotation, but you want to apply it to an array -in which the sampling along the three axes is (0.5mm, 0.5mm, 2mm). Then by setting -`SD = Diagonal(SVector(1, 1, 4))` (the ratio of scales along each axis), -one obtains an `itfm` suitable for warping the array. +Convert the physical-space transformation `ptfm` into an equivalent transformation +`itfm` that operates on array index-space. Returns an `AffineMap` suitable for +warping the array (e.g. via `ImageTransformations.warp`). + +For example, suppose `ptfm` is a pure 3D rotation, but you want to apply it to an array +sampled at (0.5 mm, 0.5 mm, 2 mm) along the three axes. By setting +`SD = Diagonal(SVector(1, 1, 4))` (the spacing ratios), you obtain an `itfm` that +correctly accounts for the anisotropy when warping. + +Any translational component of `ptfm` is interpreted in physical-space units, not +index-units. + +`SD` does not need to be diagonal; the columns of `SD` encode the physical-coordinate +displacement achieved by advancing one array element along each axis: `SD[:,j]` is +the physical displacement per voxel along dimension `j`. + +# Examples -Any translational component of `ptfm` is interpreted in physical-space units, not index-units. +```jldoctest +julia> using CoordinateTransformations, LinearAlgebra -`SD` does not even have to be diagonal, if the array sampling is skewed. -The columns of `SD` should correspond to the physical-coordinate displacement -achieved by shifting by one array element along each axis of the array. -Specifically, `SD[:,j]` should be the physical displacement of one array voxel along dimension `j`. +julia> SD = [1.0 0.0; 0.0 2.0]; # dim 2 sampled 2× coarser than dim 1 + +julia> ptfm = LinearMap([0.0 -1.0; 1.0 0.0]); # 90° CCW rotation in physical space + +julia> itfm = arrayscale(ptfm, SD); + +julia> itfm.linear +2×2 Matrix{Float64}: + 0.0 -2.0 + 0.5 0.0 +``` """ arrayscale(ptfm::AbstractAffineMap, SD::AbstractMatrix) = arrayscale(ptfm, LinearMap(SD)) @@ -124,7 +146,7 @@ its `CartesianIndices` are used. """ function default_minrot(ci::CartesianIndices, SD = I; Δc = 0.1) L = -Inf - for x in CornerIterator(ci) + for x in _corners(ci) x′ = SD * SVector(Tuple(x)) # position of corner point in physical space L = max(L, norm(x′)) end @@ -151,25 +173,32 @@ end """ getSD(A::AbstractArray) -If your image is not uniformily sampled, use this to get the `SD` matrix, which represents spacing along all axes of an image. +Return the spatial-directions matrix `SD` for array `A`. + +`SD` is an `N×N` `SMatrix` (returned as `SMatrix{N,N,Float64}`) whose columns encode +the physical-coordinate displacement corresponding to a one-element step along each +array axis: `SD[:,j]` is the physical displacement per voxel along dimension `j`. +This matrix is used by [`arrayscale`](@ref), [`qd_rigid`](@ref), and +[`qd_affine`](@ref) when the image axes are not uniformly sampled. + +Attach physical-spacing metadata to your array via +[ImageAxes.jl](https://github.com/JuliaImages/ImageAxes.jl) or a compatible package. # Examples -```julia-repl -julia> myimage -Normed ImageMeta with: - data: 3-dimensional AxisArray{N2f14,3,...} with axes: - :x, 0.0 μm:0.71 μm:6.39 μm - :l, 0.0 μm:0.71 μm:6.39 μm - :z, 0.0 μm:6.2 μm:55.8 μm -And data, a 10×10×10 Array{N2f14,3} with eltype Normed{UInt16,14} - properties: - imagineheader: - -julia> getSD(myimage) -3×3 SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3): - 0.71 0.0 0.0 - 0.0 0.71 0.0 - 0.0 0.0 6.2 + +```jldoctest +julia> using AxisArrays + +julia> img = AxisArray(rand(4, 5, 6), + Axis{:x}(0.0:0.5:1.5), + Axis{:y}(0.0:0.5:2.0), + Axis{:z}(0.0:2.0:10.0)); + +julia> getSD(img) +3×3 StaticArraysCore.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3): + 0.5 0.0 0.0 + 0.0 0.5 0.0 + 0.0 0.0 2.0 ``` """ getSD(A::AbstractArray) = getSD(spacedirections(A)) #takes advantage of how spacedirections automatically strips out the time component @@ -188,11 +217,26 @@ end imgq = qsmooth(img) imgq = qsmooth(T, img) -Create a smoothed version of `img`, smoothing with the kernel of a quadratic B-spline. -Use this on your `fixed` image in preparation for registration, and pass `presmoothed` -as an option. (Do not smooth `moving`.) +Return a smoothed copy of `img` with element type `T` (default: `float(eltype(img))`), +smoothed with the kernel of a quadratic B-spline. + +Apply to `fixed` before registration and pass `presmoothed=true` to the registration +function. Do not smooth `moving`. -`T` allows you to specify the output eltype (default `float(eltype(img))`). +# Examples + +```jldoctest +julia> img = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]; + +julia> qsmooth(img) +3×3 Matrix{Float64}: + 1.5 2.375 3.25 + 4.125 5.0 5.875 + 6.75 7.625 8.5 + +julia> eltype(qsmooth(Float32, Float32.(img))) +Float32 +``` """ 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 diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 17e0f22..441b4fc 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -23,7 +23,7 @@ using Test function fixedmov(img, tfm) img = float(img) img2 = warp(img, tfm) - inds = OffsetArrays.IdentityUnitRange.(intersect.(axes(img), axes(img2))) + inds = OffsetArrays.IdentityUnitRange.(intersect.(Base.axes(img), axes(img2))) fixed = img[inds...] moving = img2[inds...] return fixed, moving diff --git a/test/gridsearch.jl b/test/gridsearch.jl index 8ac01f0..12c63d2 100644 --- a/test/gridsearch.jl +++ b/test/gridsearch.jl @@ -3,6 +3,7 @@ using RegisterQD.CoordinateTransformations using RegisterQD.RegisterDeformation using LinearAlgebra: I using Test +using OffsetArrays @testset "Grid search rigid registration" begin ## 2D @@ -27,3 +28,14 @@ using Test @test tfm.translation == tfm0.translation @test tfm.linear == tfm0.linear end + +@testset "grid_rotations rounds even rgridsz to odd" begin + SD = Matrix{Float64}(I, 2, 2) + rots = @test_logs (:warn, r"rgridsz should be odd") RegisterQD.grid_rotations([pi / 6], [4], SD) + @test length(rots) == 5 # 4 rounded up to 5 +end + +@testset "grid_rotations unsupported dimensionality" begin + SD = Matrix{Float64}(I, 1, 1) + @test_throws ErrorException RegisterQD.grid_rotations([0.1], [3], SD) +end diff --git a/test/initial_tfm.jl b/test/initial_tfm.jl index 54e9601..78954bf 100644 --- a/test/initial_tfm.jl +++ b/test/initial_tfm.jl @@ -27,7 +27,7 @@ EYE = Matrix(1.0 * I, 3, 3) tform = RotXYZ(0.1, 0.1, 0.1) mytform = AffineMap(tform, [0, 0, 0]) - testimage2 = warp(testimage1, mytform, axes(testimage1)) + testimage2 = warp(testimage1, mytform, Base.axes(testimage1)) # insufficient parameters + no initiat_tfm causes large misalignment @@ -85,7 +85,7 @@ end tform = RotXYZ(0.1, 0.1, 0.1) mytform = AffineMap(tform, [0, 0, 0]) - testimage2 = warp(testimage1, mytform, axes(testimage1)) + testimage2 = warp(testimage1, mytform, Base.axes(testimage1)) mxrot2 = (0.2, 0.2, 0.2) minwidth_rot = RegisterQD.default_minwidth_rot(CartesianIndices(testimage2), EYE) @@ -193,7 +193,7 @@ end tform = RotXYZ(0.1, 0.1, 0.1) mytform = AffineMap(tform, [0.0, 0.0, 0.0]) - testimage2 = warp(testimage1, mytform, axes(testimage1)) + testimage2 = warp(testimage1, mytform, Base.axes(testimage1)) tformtest1, mm1 = qd_affine(testimage2, testimage1, mxshift; print_interval = typemax(Int)) # diff --git a/test/qd_random.jl b/test/qd_random.jl index afb951c..6f67b18 100644 --- a/test/qd_random.jl +++ b/test/qd_random.jl @@ -41,7 +41,7 @@ using Test, TestImages newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) - fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation + fixed = etp(Base.axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])) mxshift = (5, 5, 5) @@ -56,7 +56,7 @@ using Test, TestImages newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) - fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation + fixed = etp(Base.axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])) mxshift = (10, 10) mxrot = pi / 90 @@ -73,7 +73,7 @@ using Test, TestImages newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) - fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation + fixed = etp(Base.axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])) mxshift = (5, 5, 5) mxrot = [pi / 90; pi / 90; pi / 90] @@ -95,7 +95,7 @@ using Test, TestImages newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) - fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation + fixed = etp(Base.axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))])) mxshift = (5, 5) SD = SDiagonal(@SVector(ones(ndims(fixed)))) @@ -116,7 +116,7 @@ using Test, TestImages #mat = SArray{Tuple{3,3}}(eye(3) + rand(3,3)./30 + -rand(3,3)./30); #tfm0 = AffineMap(mat, shft); #ground truth #newfixed = warp(moving, tfm0); - #inds = intersect.(axes(moving), axes(newfixed)) + #inds = intersect.(Base.axes(moving), axes(newfixed)) #fixed = newfixed[inds...] #moving = moving[inds...] #thresh = 0.1 * (sum(abs2.(fixed[.!(isnan.(fixed))]))+sum(abs2.(moving[.!(isnan.(moving))]))); @@ -139,7 +139,7 @@ using Test, TestImages ##tfm0 = recenter(tfm00, center(moving)); #ground truth #tfm0 = tfm00 #ground truth #newfixed = warp(moving, tfm0); - #inds = intersect.(axes(moving), axes(newfixed)) + #inds = intersect.(Base.axes(moving), axes(newfixed)) #fixed = newfixed[inds...] #moving = moving[inds...] #thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))])); @@ -154,6 +154,12 @@ using Test, TestImages #TODO this test passes if run individually, but breaks when included in the testset. end #tests with random images +@testset "qd_translate crop error for undersized moving image" begin + fixed = rand(20, 20) + moving = rand(5, 5) # too small for mxshift = (5, 5): needs size ≥ 2*(5+1) = 12 + @test_throws ErrorException RegisterQD.qd_translate(fixed, moving, (5, 5); crop = true) +end + @testset "Suppression of printing" begin a, b = rand(5, 5), rand(5, 5) ca, cb = centered(a), centered(b) diff --git a/test/qd_standard.jl b/test/qd_standard.jl index 1e2800b..9d46264 100644 --- a/test/qd_standard.jl +++ b/test/qd_standard.jl @@ -18,7 +18,7 @@ using Random function fixedmov(img, tfm) img = float(img) img2 = warp(img, tfm) - inds = OffsetArrays.IdentityUnitRange.(intersect.(axes(img), axes(img2))) + inds = OffsetArrays.IdentityUnitRange.(intersect.(Base.axes(img), Base.axes(img2))) fixed = img[inds...] moving = img2[inds...] return fixed, moving diff --git a/test/runtests.jl b/test/runtests.jl index 5807884..5dc1b04 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,24 @@ using ImageMagick using RegisterQD, RegisterMismatch -using Test +using AxisArrays +using Aqua, Documenter, ExplicitImports, Test + +Aqua.test_all(RegisterQD) + +# BSplineInterpolation is intentionally accessed as internal API to bypass prefiltering +test_explicit_imports(RegisterQD; ignore=(:BSplineInterpolation,)) include("util.jl") include("qd_random.jl") include("qd_standard.jl") include("gridsearch.jl") include("initial_tfm.jl") + +DocMeta.setdocmeta!( + RegisterQD, :DocTestSetup, + :(using RegisterQD, RegisterMismatch, AxisArrays, CoordinateTransformations, LinearAlgebra); + recursive = true, +) +@testset "Doctests" begin + doctest(RegisterQD; manual=false) +end diff --git a/test/util.jl b/test/util.jl index 503d742..3c8c8fb 100644 --- a/test/util.jl +++ b/test/util.jl @@ -1,11 +1,13 @@ using ImageMagick using TestImages using RegisterQD +using RegisterQD.CoordinateTransformations: IdentityTransformation using LinearAlgebra using ImageMetadata import AxisArrays using AxisArrays: AxisArray, Axis using Unitful: μm, mm, cm, km, s +using OffsetArrays @testset "default_minwidth_rot" begin img = rand(3, 10) @@ -72,8 +74,10 @@ end @test eltype(qsmooth(Float64, img32)) === Float64 end -#TODO add a testset for other support functions -#rotations -#arrayscale -#pscale -#restrict +@testset "warp_and_intersect IdentityTransformation with mismatched axes" begin + moving = OffsetArray(rand(5, 5), 1:5, 1:5) + fixed = OffsetArray(rand(5, 5), 3:7, 3:7) + vm, vf = RegisterQD.warp_and_intersect(moving, fixed, IdentityTransformation()) + @test Base.axes(vm) == Base.axes(vf) + @test length.(Base.axes(vm)) == (3, 3) # intersection of 1:5 and 3:7 has 3 elements per dim +end