API summary

ϕ = GridDeformation(u::Array{<:SVector}, axs) creates a deformation ϕ for an array with axes axs. u specifies the "pixel-wise" displacement at a series of nodes that are evenly-spaced over the domain specified by axs (i.e., using node-vectors range(first(axs[d]), stop=last(axs[d]), length=size(u,d))). In particular, each corner of the array is the site of one node.

ϕ = GridDeformation(u::Array{<:SVector}, nodes) specifies the node-vectors manually. u must have dimensions equal to (length(nodes[1]), length(nodes[2]), ...).

ϕ = GridDeformation(u::Array{T<:Real}, ...) constructs the deformation from a "plain" u array. For a deformation in N dimensions, u must have N+1 dimensions, where the first dimension corresponds to the displacement along each axis (and therefore size(u,1) == N).

Finally, ϕ = GridDeformation((u1, u2, ...), ...) allows you to construct the deformation using an N-tuple of shift-arrays, each with N dimensions.

Example

To represent a two-dimensional deformation over a spatial region 1:100 × 1:200 (e.g., for an image of that size),

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

ϕ = 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) constructs a vector ϕs of seqeuential deformations. The last dimension of the array u should correspond 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ϕ(ϕref, coefs)

Create a deformation with the same nodes as ϕref but using coefs for the data. This is primarily useful for testing purposes, e.g., 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.

A WarpedArray W is an AbstractArray for which W[x] = A[ϕ(x)] for some parent array A and some deformation ϕ. The object is created lazily, meaning that computation of the displaced values occurs only when you ask for them explicitly.

Create a WarpedArray like this:

W = WarpedArray(A, ϕ)

where

• The first argument A is an AbstractExtrapolation that can be evaluated anywhere. See the Interpolations package.
• ϕ is an AbstractDeformation

wimg = warp(img, ϕ) warps the array img according to the deformation ϕ.

warp!(dest, src::WarpedArray) instantiates a WarpedArray in the output dest.

warp!(dest, img, ϕ) warps img using the deformation ϕ. The result is stored in dest.

warp!(dest, img, tform, ϕ) warps img using a combination of the affine transformation tform followed by deformation with ϕ. The result is stored in dest.

warp!(T, io, img, ϕs; [nworkers=1]) writes warped images to disk. io is an IO object or HDF5/JLD2 dataset (the latter must be pre-allocated using d_create to be of the proper size). img is an image sequence, and ϕs is a vector of deformations, one per image in img. If nworkers is greater than one, it will spawn additional processes to perform the deformation.

An alternative syntax is warp!(T, io, img, uarray; [nworkers=1]), where uarray is an array of u values with size(uarray)[end] == nimages(img).

Atrans = translate(A, displacement) shifts A by an amount specified by displacement. Specifically, in simple cases Atrans[i, j, ...] = A[i+displacement[1], j+displacement[2], ...]. More generally, displacement is applied only to the spatial coordinates of A.

NaN is filled in for any missing pixels.

ϕ_c = ϕ_old(ϕ_new) computes the composition of two deformations, yielding a deformation for which ϕ_c(x) ≈ ϕ_old(ϕ_new(x)). ϕ_old must be interpolating (see interpolate(ϕ_old)).

ϕ_c, g = compose(ϕ_old, ϕ_new) also yields the gradient g of ϕ_c with respect to u_new. g[i,j,...] is the Jacobian matrix at grid position (i,j,...).

You can use _, g = compose(identity, ϕ_new) if you need the gradient for when ϕ_old is equal to the identity transformation.

ϕsi_old and ϕs_new will generate ϕs_c vector and g vector ϕsi_old is interpolated `ϕs_old: e.g) ϕsi_old = map(Interpolations.interpolate!, copy(ϕs_old))

ϕs′ = medfilt(ϕs)

Perform temporal median-filtering on a sequence of deformations. This is a form of smoothing that does not "round the corners" on sudden (but persistent) shifts.

ϕs = tinterpolate(ϕsindex, tindex, nstack) uses linear interpolation/extrapolation in time to "fill out" to times 1:nstack a deformation defined intermediate times tindex . Note that ϕs[tindex] == ϕsindex.

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(ϕ)

Create an iterator for visiting all the nodes of ϕ.

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.