API summary

RegisterDeformation

A deformation (or warp) of space is represented by a function ϕ(x). For an image, the warped version of the image is specified by "looking up" the pixel value at a location ϕ(x) = x + u(x). u(x) thus expresses the displacement, in pixels, at position x. Note that a constant deformation, u(x) = x0, corresponds to a shift of the coordinates by x0, and therefore a shift of the image in the opposite direction.

In reality, deformations will be represented on a grid, and interpolation is implied at locations between grid points. For a deformation defined directly from an array, make it interpolating using ϕi = interpolate(ϕ).

The major functions/types exported by RegisterDeformation are:

- `GridDeformation`: create a deformation
- `tform2deformation`: convert an `AffineMap` to a deformation
- `ϕ_old(ϕ_new)` and `compose`: composition of two deformations
- `warp` and `warp!`: deform an image
- `WarpedArray`: create a deformed array lazily
- `warpgrid`: visualize a deformation

AbstractDeformation{T,N}

Supertype for N-dimensional deformations with displacement element type T.

A deformation maps a point x to x + u(x), where u is the displacement field. The concrete type GridDeformation represents u on a regular grid with interpolation between grid points.

Use eltype(ϕ) and ndims(ϕ) to query the element type and dimensionality.

GridDeformation(u::AbstractArray{<:SVector}, axs) -> GridDeformation
GridDeformation(u::AbstractArray{<:SVector}, nodes) -> GridDeformation
GridDeformation(u::AbstractArray{<:Real}, nodes) -> GridDeformation
GridDeformation((u1, u2, ...), nodes) -> GridDeformation

Create an N-dimensional deformation on a regular grid.

In the first form, nodes are evenly-spaced over the domain given by axs (a tuple of unit ranges, one per dimension), placing one node at each corner. In the second form, nodes is a tuple of vectors that specify node positions explicitly; u must have size (length(nodes[1]), length(nodes[2]), ...).

In the third form, u is a "plain" Real array with N+1 dimensions, where size(u, 1) == N indexes the displacement component along each axis. In the fourth form, an N-tuple of N-dimensional shift arrays is accepted instead.

Example

To represent a 2D deformation over a 1:100 × 1:200 image domain:

gridsize = (3, 5)             # a coarse grid
u = 10*randn(2, gridsize...)  # 2D displacements, ~10 pixels each
nodes = (range(1, stop=100, length=gridsize[1]), range(1, stop=200, length=gridsize[2]))
ϕ = GridDeformation(u, nodes) # "naive" deformation (not yet ready for interpolation)

NodeIterator{K,N}

Iterator returned by eachnode that visits each node of a deformation grid as an SVector. Supports length, size, and collect.

A = TransformedArray(etp, tfm)

Lazy array that applies an affine coordinate transformation on access: A[i,j] evaluates the parent array etp at the transformed coordinates [x, y] = tfm([i, j]).

etp can be an AbstractExtrapolation, AbstractInterpolation, or plain AbstractArray (automatically wrapped). tfm is an AffineMap.

See also transform, transform!.

W = WarpedArray(A, ϕ)

Lazy array for which W[x] = A[ϕ(x)], where A is the source array and ϕ is an AbstractDeformation. Indexing is evaluated on demand.

A can be any array; it is automatically wrapped in an extrapolation object so that out-of-bounds accesses return NaN rather than throwing an error.

See also warp, warp!, ImageAxes.getindex!.

ϕ = tform2deformation(tform, imgaxes, gridsize)

Construct a deformation ϕ from the affine transform tform suitable for warping arrays with axes imgaxes. The array of grid points defining ϕ has size specified by gridsize. The dimensionality of tform must match that specified by arraysize and gridsize.

Note it's more accurate to warp(img, tform) directly; the main use of this function is to initialize a GridDeformation for later optimization.

ϕs = griddeformations(u, nodes)

Construct a vector ϕs of sequential deformations. The last dimension of u corresponds to time; ϕs[i] is produced from u[:, ..., i].

ϕnew = regrid(ϕ, gridsize)

Reparametrize the deformation ϕ so that it has a grid size gridsize.

Example

ϕnew = regrid(ϕ, (13,11,7))

for a 3-dimensional deformation ϕ.

ϕ = similar_deformation(ϕref, coefs)

Create a deformation with the same nodes as ϕref but using coefs for the data. This is primarily useful for computing gradients with ForwardDiff where the elements are ForwardDiff.Dual numbers.

If ϕref is interpolating, coefs will be used for the interpolation coefficients, not the node displacements. Typically you want to create ϕref with

ϕref = interpolate!(copy(ϕ0))

rather than interpolate(ϕ0) (see interpolate!).

ϕi = interpolate(ϕ, BC=Flat(OnCell()))

Create a deformation ϕi suitable for interpolation, matching the displacements of ϕ.u at the nodes. A quadratic interpolation scheme is used, with default flat boundary conditions.

ϕi = interpolate!(ϕ, BC=InPlace(OnCell()))

Create a deformation ϕi suitable for interpolation, matching the displacements of ϕ.u at the nodes. ϕ is destroyed in the process.

A quadratic interpolation scheme is used, with the default being to use "InPlace" boundary conditions.

ϕi = extrapolate(ϕ, BC=Flat(OnCell()))

Create a deformation ϕi suitable for interpolation, matching the displacements of ϕ.u at the nodes. A quadratic interpolation scheme is used, with default flat boundary conditions and Line extrapolation.

ϕi = extrapolate!(ϕ, BC=InPlace(OnCell()))

Create a deformation ϕi suitable for interpolation, matching the displacements of ϕ.u at the nodes. ϕ is destroyed in the process.

A quadratic interpolation scheme is used, with the default being to use "InPlace" boundary conditions.

warp(img, ϕ) -> Array

Warp the array img according to the deformation ϕ. Returns an array of the same size and axes as img.

warp!(dest, src::WarpedArray) -> dest

Instantiate the WarpedArray src into the pre-allocated array dest. Returns dest.

warp!(dest, img, ϕ) -> dest

Warp img using the deformation ϕ, storing the result in the pre-allocated array dest. Returns dest.

warp!(dest, img, tform, ϕ) -> dest

Warp img by applying the affine transformation tform followed by the deformation ϕ, storing the result in the pre-allocated array dest. Returns dest.

warp!(io, img, ϕs; eltype=Float32, nworkers=1)
warp!(io, img, uarray; eltype=Float32, nworkers=1)

Write warped images to disk. io is an IO object or a pre-allocated HDF5/JLD2 dataset. img is an image sequence and ϕs is a vector of deformations, one per image in img. eltype controls the element type written to disk. If nworkers > 1, additional worker processes are spawned to parallelize the deformation.

In the second form, uarray is an array of displacement values with size(uarray)[end] == nimages(img).

translate(A, displacement) -> Array

Shift A by displacement applied to the spatial coordinates. In simple cases, result[i, j, ...] = A[i+displacement[1], j+displacement[2], ...]. Missing pixels are filled with NaN.

getindex!(dest, W::WarpedArray, coords...)

Fill dest with values from the WarpedArray W at the Cartesian product of coords. Each element of coords specifies the indices along one dimension. Returns dest.

ϕ_old ∘ ϕ_new       -> @NamedTuple{ϕ::GridDeformation, gradient}
ϕ_old(ϕ_new)         -> GridDeformation
compose(ϕ_old, ϕ_new) -> @NamedTuple{ϕ::GridDeformation, gradient}

Compute the composition of two deformations, yielding ϕ_c such that ϕ_c(x) ≈ ϕ_old(ϕ_new(x)). ϕ_old must be interpolating (see interpolate).

The functor form ϕ_old(ϕ_new) returns only the composed GridDeformation. The ∘ / compose form returns a named tuple (; ϕ, gradient) where gradient[i,j,...] is the Jacobian matrix of ϕ_c with respect to u_new at grid position (i,j,...). The result can be destructured as ϕ_c, g = compose(...).

Use compose(identity, ϕ_new) when ϕ_old is the identity transformation.

result = compose(ϕsi_old, ϕs_new)

Compose two equal-length vectors of deformations element-wise. ϕsi_old must be a vector of interpolating deformations (see interpolate!). Returns a named tuple (; ϕ, gradient) where ϕ is the vector of composed deformations and gradient is the corresponding vector of Jacobian arrays.

See also the two-argument scalar form compose(ϕ_old, ϕ_new).

tformeye(m) -> AffineMap

Return the m-dimensional identity AffineMap (identity linear part, zero translation).

tformtranslate(trans) -> AffineMap

Return an AffineMap that translates by the vector trans (identity linear part).

tformrotate(angle) -> AffineMap
tformrotate(axis, angle) -> AffineMap
tformrotate(v) -> AffineMap

Construct a rotation AffineMap (zero translation).

• tformrotate(angle): 2D rotation by angle radians.
• tformrotate(axis, angle): 3D rotation around axis (3-vector) by angle radians.
• tformrotate(v): 3D rotation where norm(v) is the angle and v/norm(v) is the axis.

See also rotation2, rotation3.

rotation2(angle) -> RotMatrix

Construct a 2D rotation matrix from angle (in radians). Returns a RotMatrix (from Rotations.jl), not an AffineMap. Use tformrotate(angle) to get an AffineMap.

rotation3(axis, angle) -> AngleAxis
rotation3(axis) -> AngleAxis

Construct a 3D rotation. Returns an AngleAxis rotation object (from Rotations.jl), not an AffineMap. Use tformrotate to get an AffineMap.

The two-argument form uses axis (a 3-vector) and angle (in radians). The one-argument form treats norm(axis) as the angle and axis/norm(axis) as the axis (angle-axis representation where the angle is encoded in the magnitude).

rotationparameters(R) -> Vector

Extract the rotation parameters from a square rotation matrix R.

• For a 2×2 matrix: returns a 1-element vector [angle] in radians.
• For a 3×3 matrix: returns a 3-element axis-angle vector angle * axis (magnitude encodes the angle, direction encodes the axis).

transform(A::TransformedArray; origin_dest=center(A), origin_src=center(A)) -> Array
transform(A, tfm::AffineMap; origin_dest=center(A), origin_src=center(A)) -> Array

Materialize the transformed array A over its entire domain. By default the transformation is assumed to operate around the center of the input array, and output coordinates are referenced relative to the center of the output.

The two-argument form wraps A in a TransformedArray first. The getindex behavior (A[i,j]) assumes the origin at zero; to match it, pass origin_src = zeros(N) and origin_dest = zeros(N), or equivalently offset the transform by origin_src - tfm.linear * origin_dest.

See also transform!.

transform!(dest, src::TransformedArray; origin_dest=center(dest), origin_src=center(src)) -> dest
transform!(dest, src, tfm::AffineMap; origin_dest=center(dest), origin_src=center(src)) -> dest

In-place version of transform. Writes the result into the pre-allocated array dest and returns it.

ϕs′ = tmedfilt(ϕs, n)

Perform temporal median-filtering on a sequence of deformations with window size n (which must be odd). This is a form of smoothing that does not "round the corners" on sudden (but persistent) shifts.

See also tmedfilt!.

tmedfilt!(out, ϕs, n)

In-place version of tmedfilt. Writes the filtered deformations into the pre-allocated vector out, which must have the same length as ϕs.

ϕs = tinterpolate(ϕsindex, tindex, nstack)

Linearly interpolate/extrapolate a sparse set of deformations in time to fill all frames 1:nstack. ϕsindex is the vector of known deformations defined at the integer time indices tindex. Satisfies ϕs[tindex[i]] == ϕsindex[i].

img = nodegrid(T, nodes)
img = nodegrid(nodes)
img = nodegrid([T], ϕ)

Returns an image img which has value 1 at points that are "on the nodes." If nodes/ϕ is two-dimensional, the image will consist of grid lines; for three-dimensional inputs, the image will have grid planes.

The meaning of "on the nodes" is that x[d] == nodes[d][i] for some dimension d and node index i. An exception is made for the edge values, which are brought inward by one pixel to prevent complete loss under subpixel deformations.

Optionally specify the element type T of img.

See also warpgrid.

img = warpgrid(ϕ; [scale=1, showidentity=false])

Returns an image img that permits visualization of the deformation ϕ. The output is a warped rectangular grid with nodes centered on the control points as specified by the nodes of ϕ. If ϕ is a two-dimensional deformation, the image will consist of grid lines; for a three-dimensional deformation, the image will have grid planes.

scale multiplies ϕ.u, effectively making the deformation stronger (for scale > 1). This can be useful if you are trying to visualize subtle changes. If showidentity is true, an RGB image is returned that has the warped grid in magenta and the original grid in green.

See also nodegrid.

iter = eachnode(ϕ)
iter = eachnode(nodes)

Create an iterator that visits each node of ϕ (or nodes) as an SVector. The iteration order follows column-major (first index varies fastest).

Example

ϕ = GridDeformation(zeros(SVector{2,Float64}, 3, 4), (1:3, 1:4))
for x in eachnode(ϕ)
    # x is an SVector{2,Float64}
end

caxs = centeraxes(axs)

Return a set of axes centered on zero. Specifically, if axs[i] is a range, then caxs[i] is a UnitRange of the same length that is approximately symmetric around 0.

arraysize(nodes) -> NTuple

Return the array size implied by nodes, where each element of the tuple is Int(maximum(n) - minimum(n) + 1) for the corresponding node range n.

vecindex(A, x::SVector) -> element

Index into A at the position given by the SVector x, equivalent to A[x[1], x[2], ...]. For AbstractInterpolation objects the callable form is used instead of getindex.

vecgradient!(g, itp, x::SVector)

Compute the gradient of interpolation object itp at position x and store the result in g. Thin wrapper around Interpolations.gradient! that accepts an SVector position.