rcpd

Coherent point drift in Rust, with Python bindings: fit a transform between two point clouds when you do not know — and do not need to know — which point matches which.

The clouds may have different numbers of points. The fit is by EM over a Gaussian mixture (Myronenko & Song 2010), and what comes back is a closed-form transform, so it applies exactly to points that took no part in fitting it.

import rcpd

tr = rcpd.register_rigid(source, target)     # rotation, translation, scale
tr = rcpd.register_deform(source, target)    # a smooth non-rigid warp
moved = tr.apply(source)

use rcpd::{rigid_register, RigidOpts};

let fit = rigid_register(source.view(), target.view(), RigidOpts::default(), None)?;
let moved = fit.transform.apply(source.view())?;

What is here

• Rigid / similarity registration — rotation, translation and optionally one uniform scale, which can be confined to a range you consider plausible.
• Deformable registration — a smooth warp, regularised by the motion-coherence prior and solved through the low-rank formulation from the original paper.

Either as a single pair or a whole grid of pairs at once, one registration per core. Both return a function of position rather than a set of moved points, so either applies exactly to points that took no part in the fit — which is what makes fitting on a subsample free of approximation, and what saves a grid of registrations from being hundreds of gigabytes of displaced clouds.

Install

pip install rcpd

or, for the Rust crate,

cargo add rcpd

Rigid, then deformable

The two fits are not alternatives. Running both is the usual workflow: the rigid fit takes out the pose, and the deformable fit is asked only for what is left.

rigid = rcpd.register_rigid(source, target)
moved = rigid.apply(source)

deform = rcpd.register_deform(moved, target)   # note: fitted on the *moved* points
final = deform.apply(moved)

rigid.nrms, deform.nrms                        # the two are on the same scale

Points that took no part in either fit go through both, in the same order — deform.apply(rigid.apply(pts)).

Fitting the deformation on moved rather than on source is the load-bearing part. The motion-coherence prior penalises displacement, so any pose offset left for the deformation to undo is paid for out of the same budget as the warp — and past a large enough offset it is not undone at all: the fit settles into a wrong local optimum whose residual reads like two unrelated clouds. On a 400-point cloud rotated 70 degrees, shifted, and warped by something no rigid transform can express, that is exactly what happens:

Where the offset stops being recoverable is data-dependent; on that cloud it is somewhere between 45 and 70 degrees.

Judging a fit

Nearly every fit converges, and that is not the same as being right — coherent point drift will settle a cloud onto an unrelated one perfectly happily. So the transform carries the residual it left behind:

tr.rms      # residual distance, in the units that went in
tr.nrms     # the same, as a fraction of the target cloud's own radius

nrms is the one to compare between fits. It is calibrated rather than merely ordered — displace every target point by 1% of the cloud's radius and it reads 0.01 — so a cutoff can be picked by hand. Measured on five neuron skeletons:

A rigid fit's and a deformable fit's nrms are on the same scale, which is the point: run both and the pair of numbers says whether the extra freedom bought anything. On two example skeletons, rigid gives 0.087 and deformable 0.058.

Below ~1e-3 there is nothing left to measure: the E-step runs in f32, which resolves the variance to ~1e-7 relative and hence the residual — its square root — to ~3e-4.

Why not pycpd?

pycpd is the reference Python implementation and has been unmaintained since 2021. It also writes the E-step as

P = np.sum((X[None, :, :] - TY[:, None, :]) ** 2, axis=2)

which materialises an M x N x 3 array, squares it into a second, then builds three more M x N temporaries — every iteration, of which there are typically 50-100. Two ~4,500-point clouds:

The memory is the more interesting column. The M-step only ever reads four reductions of the correspondence matrix, and its normaliser is a reduction over rows — so a block of columns is self-contained, and the whole E-step fits in one exp pass over a few hundred KB. Peak memory stops depending on M x N at all, which is what makes the batch entry point possible: at a gigabyte apiece you cannot hold fourteen registrations at once.

Across a grid the gap widens, because pycpd has no parallelism to give: a 12 x 12 pairwise alignment of the same clouds takes 30 s here against ~26 min.

Deformable is where pycpd stops being usable at all, and for an algorithmic reason rather than an implementation one: its M-step solves a dense M x M system once per iteration, which is cubic in the point count in time and quadratic in memory. On the same two ~4,500-point clouds, 50 iterations:

The fix is the one in the original paper — approximate the kernel by its leading eigenpairs and apply the Woodbury identity, so each iteration is a K x K solve. Those eigenpairs come from a randomised subspace iteration that never forms the kernel, rather than the dense np.linalg.eigh(G) that cycpd uses; see core/src/lowrank.rs for why, and for why there is no fast Gauss transform here.

Held at the same iteration count — taking the deliberate convergence difference out of the comparison — the two implementations move the points to within 3e-6 of the cloud's radius at five iterations, and 3e-4 by sixty.

Deliberate differences from pycpd

Four, all fixes.

Convergence is relative. pycpd stops when its objective moves by less than tolerance in absolute terms — on a quantity that scales with the data, so the same clouds converge differently in nanometres and in microns. Here the test is on the relative change in the fitted variance, so it is dimensionless.

A deformable fit gets a second stopping rule, because it has enough freedom to drive the residual below what the f32 kernel can resolve — and there the variance stops descending and rattles by a few percent an iteration around a value it has already reached, so tolerance is never met however long it runs. Over a window of iterations the fit compares how far the variance moved with how far it travelled, and stops when it is travelling without arriving. On a well-matched pair that is 48 iterations rather than 200, for a residual 2% different; a fit that is merely slow is untouched, since a descending variance has the two distances equal.

scale is honoured. pycpd.RigidRegistration takes no scale argument in any released version, including master, so a caller that passes one has it land in **kwargs and a scale fitted regardless. scale=False here holds it at exactly 1, and scale_bounds holds it inside a range:

tr = rcpd.register_rigid(source, target, scale_bounds=(0.8, 1.25))

The limits are imposed at every EM step, so what comes back is the best alignment within them — rotation and translation re-fitted against the limited scale — rather than a free fit squashed into range afterwards. A pair that wants a scale outside the range comes back sitting exactly on the nearer limit, which is how you tell that the constraint bound. This is worth having on a grid of pairs: two clouds of the same shape at very different sizes will otherwise be talked into a flattering nrms by a scale you know to be impossible.

w is scale-free. The outlier term weighs a Gaussian density against a uniform one over some volume V; pycpd drops V, leaving a uniform density of 1, which is only right for data spanning about one unit. On a cloud spanning 1e4 units every point is declared an outlier and the fit collapses — at w=0.1 pycpd returns a scale of 0.19 against a true 1.0, and at w=0.01 it does not converge in 500 iterations. Here V is carried explicitly, as the bounding volume of the fixed cloud.

At w=0 — the default — the two agree to the precision of the f32 kernel: on real skeletons, transformed points match to 4e-8 of the cloud's extent.

beta is a fraction of the cloud, not a length. The motion-coherence kernel's width in data units silently means something different for the same object in nanometres and in microns, and the failure is quiet rather than loud: too wide a kernel makes the kernel matrix nearly all-ones, the deformation collapses to a global translation, and what comes back looks like a fit rather than like an error. Here beta is a fraction of the cloud's radius, so it means the same thing whatever the units. Note it is coupled to num_modes — the kernel's rank grows roughly as (1/beta)**3.

There is also no LinAlgError to work around. The 3x3 SVD here is a Jacobi eigensolve that cannot fail to converge, and a rank-deficient system (collinear points) completes the rotation basis rather than raising.

Development

cargo test --workspace                 # Rust core and bindings
cd py && maturin develop --release     # build the extension in place
cd py && pytest tests/                 # Python suite (checks against pycpd where installed)

Releasing

The version lives in the workspace Cargo.toml and nowhere else — the crate, the extension and the Python package all inherit it. Bump it, land the changelog entry, then tag:

git tag v0.1.0 && git push origin v0.1.0

.github/workflows/release.yml builds wheels for macOS, Linux (glibc and musl, x86-64 and aarch64) and Windows, plus an sdist that it builds back into a wheel to check it is not broken, publishes those to PyPI, then publishes the crate to crates.io. It refuses to start if the tag and the manifest disagree, since neither registry lets a version be replaced. Running it from the Actions tab instead builds everything and publishes nothing, which is how to test a change to it.

Licence

GPL-3.0-or-later.