Module registration 

Rigid-body (6-DOF) registration between two volumes: recovering the transform that geometry can already apply.

Multi-orientation QSM needs two things from it — the orientations on one common grid, and the B0 direction each one contributes to that grid. The second comes out of the recovered rotation, because the object rotates and B0 does not; see registration::RigidTransform::b0_direction_in_fixed. Rigid-body (6-DOF) volume registration: estimating the transform between two volumes.

crate::geometry can already apply an arbitrary transform — crate::geometry::GridMap maps one voxel grid onto another given two affines, and crate::geometry::resample_onto, crate::geometry::resample_mask_onto and crate::geometry::resample_complex_onto move continuous data, labels and wrapped phase across it. What was missing is finding the transform in the first place, which is what this module does.

Why QSM needs it

Multi-orientation QSM (COSMOS, [crate::inversion::cosmos]) and susceptibility tensor imaging ([crate::inversion::sti]) both take N local field maps on one common grid plus one B0 direction per orientation in that grid’s voxel frame. Registration produces both, and the second is not a by-product:

The object rotates and B0 does not. So the rigid rotation that aligns orientation t to the reference is where B0 pointed during orientation t, relative to the common grid.

That matters most for the dataset shape where the affine is useless: identically prescribed slabs, where every orientation has the same affine and the anatomy moved in voxel space instead. crate::geometry::b0_direction_from_affine returns the same direction for all of them, and N identical directions quietly collapse COSMOS to an unregularized single-orientation inversion. Registration is the only thing that can recover the real answer there — see [RigidTransform::b0_direction_in_fixed], which is where that convention lives so no caller has to re-derive it.

Register on magnitude, resample phase through the complex domain

Wrapped phase cannot be interpolated directly (halfway between +3.0 and -3.0 rad a linear interpolator returns 0.0; see the crate::geometry module docs), so it must not be used as the registration input either. Register magnitude against magnitude, then move phase with [RigidTransform::resample_complex], which goes through mag·e^{iφ}.

Metric: normalised cross-correlation, not mutual information

Every use case here is mono-modal — GRE magnitude against GRE magnitude of the same subject, same sequence. NCC is invariant to an affine intensity change (gain and offset), which is what actually differs between two orientations of one acquisition, and it needs no joint histogram, no binning choice and no parzen window. Mutual information buys nothing on same-contrast data and is noisier at the small sample counts a coarse pyramid level uses. Sum-of-squared-differences would be cheaper still but is not gain-invariant, so a receive-gain difference between runs would bias it.

The one thing NCC does not absorb is spatially varying shading: the receive field is fixed to the coil, so when the head rotates the bias field rotates with the coil rather than with the anatomy. That is a multiplicative field, not a scalar gain, and no global metric (MI included) is invariant to it. Two mitigations, both available here: pass a brain mask so the metric ignores everything outside the head, and/or hand in magnitude that has already been bias-corrected (crate::homogeneity).

Optimiser, pyramid and cost

Powell’s direction-set method with a bracketed Brent line search, over a coarse-to-fine Gaussian pyramid. Powell needs no derivatives of the metric — convenient, because the metric’s derivative runs through a trilinear interpolator — and it is what the established intensity registration tools use for 6-DOF problems. No linear-algebra crate is involved; the whole module is f64 arrays, consistent with the rest of qsm-core.

The pyramid is deep on purpose. It is what sets the capture range — how far apart two volumes can start and still converge — because smoothing and decimating flattens the cost surface’s ripples, and the ripples are what a search walks into instead of the true minimum. See [RigidParams::levels] for the measured numbers; the short version is that three levels converge to 40° and five to beyond 70°, and 45° of repositioning is ordinary for STI.

Work per metric evaluation is bounded independently of volume size by [RigidParams::max_samples]: each level picks the smallest integer voxel stride whose sample count fits the cap. Coarse levels therefore see smoothed, decimated data (a smooth cost landscape, which is the point of the pyramid), while the finest level strides the unsmoothed full-resolution volume (full spatial detail at the same cost, which is what gives the final precision). Peak extra memory is the pyramid above the finest level — about 1/7 of one input volume — plus the sample list, so a 32-bit WASM host with a 4 GB heap is not at risk.

Structs

Acc 🔒

Running sums for one chunk of samples, so the reduction is order-independent.

CostFn 🔒

The cost being minimised: 1 - NCC over the overlapping, masked samples.

RigidParams

Registration parameters.

RigidTransform

A recovered rigid transform, together with everything needed to use it.

Sample 🔒

One fixed-grid sample the metric is evaluated on: its voxel coordinate and its intensity.

Volume 🔒

One pyramid level: either a borrowed view of the caller’s data (level 0) or a decimated copy.

Constants

COST_DEGENERATE 🔒

Costs assigned where NCC is not a meaningful number. Both are above the worst real cost (1 - (-1) = 2), so a degenerate candidate can never win, and the overlap case keeps a gradient that pushes the search back into the valid region.

COST_NO_OVERLAP 🔒

IDENTITY4 🔒

SAMPLE_CHUNK 🔒

Fixed samples per parallel chunk. Fixed size, so the reduction combines in index order and the metric is bit-identical regardless of thread count.

Functions

binomial_blur3 🔒

Separable [1 4 6 4 1]/16 blur with edge replication.

decimate_mask 🔒

Decimate a mask the same way Volume::halved decimates data: keep the even samples, with no blur, so the coarse mask stays binary and sits on the coarse sample positions exactly.

direction_cosines 🔒

The direction-cosine matrix of an affine: its 3×3 with the voxel sizes divided out, so what is left is the pure rotation taking voxel axes to world axes.

euler_rotation 🔒

Rz(rz)·Ry(ry)·Rx(rx), angles in radians.

grid_center_world 🔒

World coordinate of a grid’s centre voxel.

line_minimise 🔒

Minimise cost along p + x·dir: bracket the minimum, then Brent.

mat3_transpose_vec 🔒

mat3_vec 🔒

mat4_mul 🔒

matrix_from_params 🔒

The rigid 4×4 for a parameter vector: x_moving = center + R·(x_fixed - center) + t.

params_from_matrix 🔒

Recover a parameter vector from a rigid 4×4, about center. The inverse of matrix_from_params, used by RigidTransform::inverse.

powell 🔒

Powell’s direction-set method.

register_rigid

Estimate the rigid transform that aligns moving to fixed.

rigid_inverse 🔒

Inverse of a rigid 4×4: transpose the rotation, and carry the translation back through it.

rotation_of 🔒

The 3×3 block of a 4×4.

sample_set 🔒

Pick the fixed voxels the metric looks at: the smallest integer stride whose masked sample count fits max_samples.

to_mat4 🔒