Landmark transforms¶

When no image registration exists between two template spaces, what you usually have instead is a set of matched landmarks — points a human (or another algorithm) has identified as corresponding. Two classical methods turn those into a warp:

• TpsTransform — a thin-plate spline. Interpolates the landmarks exactly and, between them, minimises the integral bending norm: the smoothest warp consistent with the data.
• MlsTransform — moving least squares (the affine flavour of Schaefer et al. 2006). Gives every point its own affine, solved on the fly from all landmarks weighted by inverse squared distance.

These mirror navis.transforms.TPStransform and navis.transforms.MovingLeastSquaresTransform, and agree with them to ~1e-14 relative.

import navis_fastcore as fastcore

tps = fastcore.TpsTransform(landmarks_source, landmarks_target)
xf = tps.xform(points)          # (N, 3) -> (N, 3)

mls = fastcore.MlsTransform(landmarks_source, landmarks_target)
xf = mls.xform(points)

Landmarks and points may be (N, 3) arrays or DataFrames with x/y/z columns. A single (3,) point is accepted and returns a (3,) point.

Which one?¶

They give similar but not identical results, and which suits a given pair of spaces is worth checking empirically. As a rule of thumb:

Both fall back to a sensible global affine far outside the landmark hull, exposed as .matrix_affine.

Direction¶

Neither method has a closed-form inverse, so "backwards" means fitting the other way round, not inverting. Both spell it -transform:

back = (-tps).xform(points)     # refits target -> source; costs another fit
back = (-mls).xform(points)     # free - MLS has no fit to redo

For MLS you can also ask at construction: MlsTransform(src, trg, direction="inverse").

No batch_size¶

The navis versions take a batch_size because both reference implementations build an intermediate sized by points × landmarks. That is what makes them expensive, and for MLS it is a hard ceiling: molesq allocates (3, M, N) arrays, so 3,400 landmarks at the default batch size wants ~23 GB — which means MovingLeastSquaresTransform cannot run at the landmark counts real registrations use.

Here that intermediate never exists. Each output row depends only on its own row of the distance matrix, so the distance and its contribution to the result are fused into one streaming pass over landmarks. Peak memory is the output array, independent of the landmark count, and there is nothing to tune.

Transforming 1M points through a 3,390-landmark spline:

Speed¶

Measured on 14 cores against the reference implementations (morphops for TPS, molesq for MLS — these are what navis calls), 1M points:

500 landmarks¶

3,390 landmarks (a real flybrains mirroring registration)¶

Worth reading honestly: TPS's win is almost entirely parallelism. Single-threaded we are only ~1.3× ahead, because cdist plus a BLAS matmul is already good code. What we add is scaling across cores and the removal of the memory ceiling.

MLS is a genuine algorithmic win — 3.7× before any parallelism — because molesq expresses the reduction as a chain of einsums over materialised per-landmark arrays.

The one place we are slower: fitting a spline¶

The TPS fit solves an (M+4) square system, so it is cubic in the landmark count. numpy sends that to hardware LAPACK (on Apple silicon, the AMX coprocessor); we use a blocked LU in portable Rust, which cannot match it:

This is a one-off cost per registration, against an xform that is 14× faster and can be called thousands of times, so it is very rarely the thing to optimise. If it is, fit with numpy and hand the coefficients over:

tps = fastcore.TpsTransform.from_coefs(source, W, A)

W and A follow the same convention as morphops.tps_coefs and navis.transforms.TPStransform.W / .A, so a navis transform converts directly:

fast = fastcore.TpsTransform.from_coefs(navis_tps.source, navis_tps.W, navis_tps.A)

MLS has no fit at all, so none of this applies to it.

Pickling¶

Both transforms pickle. TpsTransform ships its coefficients rather than just the landmarks, so unpickling in a multiprocessing worker does not repeat the fit.

API¶

            Source landmarks as x/y/z coordinates. A pandas ``DataFrame`` with
            x/y/z columns is also accepted.

            Target landmarks, one per source landmark.

>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
>>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
>>> tr = fastcore.TpsTransform(src, trg)
>>> tr.xform(np.array([[0, 0, 0], [50, 50, 50]]))
array([[ 1.        , 15.        ,  5.        ],
       [40.55555556, 54.        , 65.        ]])
>>> # Landmarks are reproduced exactly
>>> np.allclose(tr.xform(src), trg)
True
>>> # Negation refits in the opposite direction
>>> np.allclose((-tr).xform(trg), src)
True

            The source landmarks the coefficients belong to.

            Weights of the non-affine part.

            Coefficients of the affine part.

            The target landmarks. Not needed to transform points, but
            without them :meth:`__neg__` cannot refit the inverse.

>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
>>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
>>> tr = fastcore.TpsTransform(src, trg)
>>> same = fastcore.TpsTransform.from_coefs(tr.source, tr.W, tr.A)
>>> np.allclose(same.xform(src), trg)
True

    Coordinates to transform. A single ``(3,)`` point is accepted and
    returns a ``(3,)`` point; a ``DataFrame`` with x/y/z columns also works.

    Number of threads. ``None`` (default) uses all cores.

            Source landmarks as x/y/z coordinates. A pandas ``DataFrame`` with
            x/y/z columns is also accepted.

            Target landmarks, one per source landmark.

            ``"inverse"`` treats the target as the source and vice versa.
            Note this fits the warp in the opposite direction; it is not an
            exact inverse, which moving least squares does not have.

>>> import navis_fastcore as fastcore
>>> import numpy as np
>>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
>>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
>>> tr = fastcore.MlsTransform(src, trg)
>>> # A point sitting on a landmark maps to its partner
>>> np.allclose(tr.xform(src), trg)
True
>>> # Negation swaps the direction
>>> np.allclose((-tr).xform(trg), src)
True

    Coordinates to transform. A single ``(3,)`` point is accepted and
    returns a ``(3,)`` point; a ``DataFrame`` with x/y/z columns also works.

    Number of threads. ``None`` (default) uses all cores.

    Override this transform's ``direction`` for this call only:
    ``False`` maps source -> target, ``True`` target -> source.
    ``None`` (default) uses ``direction``. The underlying fit is
    direction-agnostic, so this saves rebuilding the object just
    to run it the other way.