Metadata-Version: 2.4
Name: rcpd
Version: 0.1.0
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: GNU General Public License v3 (GPLv3)
Classifier: Natural Language :: English
Classifier: Programming Language :: Rust
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: Free Threading :: 2 - Beta
Requires-Dist: numpy>=1.21.0
License-File: LICENSE
Summary: Coherent point drift: registering point clouds with no correspondences given.
Author-email: Philipp Schlegel <pms70@cam.ac.uk>
Requires-Python: >=3.9
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Issues, https://github.com/schlegelp/rcpd/issues
Project-URL: Repository, https://github.com/schlegelp/rcpd

# 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](https://arxiv.org/pdf/0905.2635.pdf)), and what comes back is a
closed-form transform, so it applies exactly to points that took no part in fitting it.

```python
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)
```

```rust
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.

```python
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:

| | `nrms` |
|---|---|
| rigid alone | 0.053 |
| deformable alone | 0.25 — no better than two unrelated clouds |
| rigid, then deformable | 1.5e-4 — at the numerical floor |

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:

```python
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:

| `nrms` | |
|---|---|
| 2e-5 | one cloud against a rotated copy of itself |
| 0.008 | one neuron against its own other half |
| 0.06 – 0.11 | two different neurons of the same type |
| ~0.2 | two unrelated point clouds |

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

```python
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:

| | time | peak RSS |
|---|---|---|
| `pycpd` | 10.8 s | +1,759 MB |
| `rcpd`, 1 core | 2.48 s | +0.4 MB |
| `rcpd`, 14 cores | 0.29 s | +25 MB |

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:

| | time | peak RSS |
|---|---|---|
| `pycpd` | 36.6 s | +3,191 MB |
| `rcpd` | 0.80 s | +95 MB |

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.

**`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:

```python
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.

