API Reference

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

Examples

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

tform, mm = qd_rigid(fixed, moving, mxshift, mxrot;
                     presmoothed=false,
                     SD=I, minwidth_rot=default_minwidth_rot(fixed, SD),
                     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, using a fourier method to speed up the search for the best whole-pixel shift. The second step performs the rotation + shift in combination to obtain sub-pixel accuracy. Any NaN-valued pixels are not included in the mismatch; you can use this to mask out any regions of fixed that you don't want to align against.

mxshift is the maximum-allowed translation, in units of array indices. It can be passed as a vector or tuple. mxrot is the maximum-allowed rotation, in radians for 2d or quaternion-units for 3d. 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 distinguishes this from the full-search-space minwidth accepted by qd_translate: 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 (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; there is currently no way to adjust these criteria for the coarse step.

The rotation returned 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.

Use SD ("spatial displacements") 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. See arrayscale for more information about SD.

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 you have a good initial guess at the solution, pass it with the initial_tfm kwarg to jump-start the search.

Use presmoothed=true if you have called qsmooth on fixed before calling qd_rigid. 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.

Examples

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

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 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 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 and qd_rigid because affine transformations have extra degrees of freedom (scaling, shear) that make degenerate low-overlap solutions more likely.

Examples

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

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

Examples

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

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), 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 for how maxradians, rgridsz, and SD shape the search.

Examples

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);

julia> best_mm < 0.1
true

itfm = arrayscale(ptfm, SD)

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

julia> using CoordinateTransformations, LinearAlgebra

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

getSD(A::AbstractArray)

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, qd_rigid, and qd_affine when the image axes are not uniformly sampled.

Attach physical-spacing metadata to your array via ImageAxes.jl or a compatible package.

Examples

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

imgq = qsmooth(img)
imgq = qsmooth(T, img)

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.

Examples

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