API Reference

This module computes the total registration penalty, combining both "data" (the mismatch between fixed and moving images) and "regularization" (a penalty applied to deformations that do not fit some pre-conceived notion of "goodness").

The main exported types/functions are:

• AffinePenalty: regularization that penalizes deviations from an affine transformation
• DeformationPenalty: abstract supertype; subtype to implement custom regularization
• penalty!: compute the total, data, or regularization penalty
• interpolate_mm!: prepare mismatch arrays for interpolation

AffinePenalty(nodes::NTuple{N, <:AbstractVector}, λ)
AffinePenalty(nodes::AbstractVector{<:AbstractVector}, λ)
AffinePenalty(nodes::AbstractMatrix, λ)

Construct an affine-residual regularization penalty for use with penalty!. Returns an AffinePenalty{T,N} where T is the element type inferred from nodes and N is the spatial dimensionality.

The penalty measures how much a deformation's displacement field u deviates from any affine transformation:

penalty = (λ/n) * ∑_i ‖u_i - a_i‖²

where {a_i} is the least-squares affine fit to {u_i} and n = length(u).

For regular grids, pass nodes as an N-tuple of coordinate vectors (one per spatial dimension) or as a vector of such vectors. For irregular point clouds, pass an N×npoints matrix whose columns are the node coordinates.

Examples

julia> nodes = (range(-1.0, 1.0, length=5), range(-1.0, 1.0, length=5));

julia> ap = AffinePenalty(nodes, 0.1);

julia> ndims(ap), eltype(ap)
(2, Float64)

DeformationPenalty{T, N}

Abstract supertype for regularization penalties on N-dimensional deformation fields with element type T.

Subtypes implement regularization by defining a penalty! method:

penalty!(g, dp::MyPenalty, ϕ_c::AbstractDeformation) -> scalar

where g is pre-allocated gradient storage (written in-place) and ϕ_c is the (possibly composed) deformation being penalized. The built-in subtype is AffinePenalty.

Base.eltype and Base.ndims are defined for all subtypes and return T and N respectively.

mmi  = interpolate_mm!(mm)
mmi  = interpolate_mm!(mm, Quadratic)
mmi  = interpolate_mm!(mm, Linear)
mmi  = interpolate_mm!(mm, itype)
mmis = interpolate_mm!(mms)
mmis = interpolate_mm!(mms, Quadratic)
mmis = interpolate_mm!(mms, itype)

Prepare a MismatchArray (as returned by RegisterMismatch) for sub-pixel interpolation. Returns a CenterIndexedArray{NumDenom{T},N} wrapping a B-spline or gridded interpolant, ready for use with penalty!.

The order argument is a Degree subtype from Interpolations.jl: Quadratic (default) or Linear. With Quadratic, the array data are overwritten in-place with B-spline prefilter coefficients. For fine-grained control over the boundary condition or interpolation scheme, pass an InterpolationType directly:

mmi = interpolate_mm!(mm, BSpline(Quadratic(InPlace(OnCell()))))

The array-of-MismatchArrays form processes each element and returns a Vector{CenterIndexedArray{NumDenom{T},N}}.

Examples

julia> using RegisterCore

julia> mm = MismatchArray(Float64, (11, 11));

julia> mmi = interpolate_mm!(mm);

julia> eltype(mmi), ndims(mmi)
(NumDenom{Float64}, 2)

julia> mms = [MismatchArray(Float64, (5, 5)) for _ in 1:4];

julia> mmis = interpolate_mm!(mms);

julia> length(mmis), eltype(first(mmis))
(4, NumDenom{Float64})

p = penalty!(g, ϕ, ϕ_old, dp::DeformationPenalty, mmis)
p = penalty!(g, ϕ, ϕ_old, dp::DeformationPenalty, mmis, keep)

Compute the total penalty (regularization + data) for deformation ϕ, mismatch data mmis, and an optional "base" deformation ϕ_old. Returns a scalar of the same element type as ϕ.

When ϕ_old is not identity, the regularization penalty is evaluated on the composed deformation ϕ_c = ϕ_old(ϕ), while the data penalty uses ϕ and mmis (which should express the residual mismatch after applying ϕ_old). This supports an incremental registration workflow:

• Compute initial deformation ϕ_0 that partially aligns fixed and moving
• Warp moving with ϕ_0
• Compute residual mismatch between fixed and the warped image
• Optimize ϕ_1 to correct the residual; final deformation is ϕ_0(ϕ_1)

This method resolves to:

val =  penalty!(g, dp, ϕ_c, g_c)   # regularization; see [`penalty!(g, dp, ϕ_c)`](@ref)
val += penalty!(g, ϕ, mmis, keep)  # data;           see [`penalty!(g, ϕ, mmis)`](@ref)

g_c is the Jacobian of ϕ_c with respect to ϕ.u. When ϕ_old == identity, no composition is needed and ϕ is used directly.

ϕ_old must be interpolating if not identity; ϕ must not be interpolating.

g may be a flat Vector{T} or an array of SVector{N,T} with the same shape as ϕ.u. It is written by the regularization term and incremented by the data term; initialize it with zeros before calling.

p = penalty!(gs, ϕs, ϕs_old, dp, λt, mmis)
p = penalty!(gs, ϕs, ϕs_old, dp, λt, mmis, keeps)

Evaluate the total penalty for a temporal sequence of deformations ϕs and the corresponding mismatch data mmis. λt is the temporal roughness penalty coefficient. Returns a scalar of the same element type as the deformations.

This is the temporal-sequence analogue of the single-frame penalty!(g, ϕ, ϕ_old, dp, mmis): it calls the single-frame penalty for each time point and adds the temporal-roughness penalty penalty!(gs, λt, ϕs).

p = penalty!(g, λt, ϕs)

Compute the temporal-roughness penalty

(λt/2) * ∑_i ‖ϕ_{i+1} - ϕ_i‖²

for a vector ϕs of GridDeformations. Returns a non-negative scalar of the same element type as the deformations.

g, if non-nothing and non-empty, must be a flat Vector whose length equals length(ϕs) * length(first(ϕs).u). On return, g holds the gradient of the penalty with respect to all displacement vectors.

p = penalty!(g, ϕ, mmis)
p = penalty!(g, ϕ, mmis, keep)
p = penalty!(g, u, mmis)
p = penalty!(g, u, mmis, keep)

Compute the data penalty — the total mismatch between fixed and moving given deformation ϕ (or displacement array u). Returns a scalar of the same element type as the mismatch data.

mmis is an array-of-CenterIndexedArrays prepared by interpolate_mm!. keep is an optional boolean array the same size as mmis; blocks with keep[i] == false are skipped (treated as zero mismatch).

The penalty is a globally-normalized ratio:

    pnum_1 + pnum_2 + ... + pnum_n

p = –––––––––––––––– pden1 + pden2 + ... + pden_n

where each pnum_i / pden_i is the mismatch at aperture i evaluated at ϕ.u[i] (or u[i]).

g is pre-allocated gradient storage (same shape as ϕ.u or u), and may be nothing or empty to skip gradient computation. This function adds to g; initialize it with zeros or call the regularization penalty first.

p = penalty!(g, dp::DeformationPenalty, ϕ_c)
p = penalty!(g, dp::DeformationPenalty, ϕ_c, g_c)

Compute the regularization penalty for a (possibly composed) deformation ϕ_c. Returns a non-negative scalar of the same element type as ϕ_c.

The _c suffix indicates "composed": use the two-argument form when ϕ_c = ϕ_old(ϕ) and g_c is the Jacobian of that composition with respect to ϕ.u. When ϕ is not derived by composition, pass it directly as ϕ_c and omit g_c.

dp determines the type of penalty; the built-in implementation is AffinePenalty. Custom regularization can be added by subtyping DeformationPenalty and defining a new penalty! method.

If g is non-nothing and non-empty, the gradient of the penalty is written into g (unlike the data penalty, which adds to g). When g_c is supplied, g is adjusted by the chain rule so that it is the gradient with respect to ϕ.u rather than ϕ_c.u.