brain_deer.warp.rbf

Variable-radius Gaussian RBF scattered-data interpolation.

Given correspondences source_i -> target_i, each with its own radius, this builds a displacement field d such that phi(x) = x + d(x) carries each target_i back onto its source_i.

The basis functions are centred on the targets, and the field is always rebuilt from scratch from the full point set – there is no incremental state.

REFERENCE BEHAVIOUR

This reimplements Slicer’s FiducialRegistrationVariableRBF CLI (netstim/ SlicerNetstim), whose source was never vendored into Lead-DBS. The conventions below were recovered empirically by probing the compiled binary; see scripts/warp_golden.py and tests/test_warp_rbf.py.

basis phi_j(x) = exp(-||x - t_j||^2 / r_j^2) centred on TARGETS rhs source - target solve (K + lam*G) w = rhs, K[i][j] = phi_j(t_i), G = diag(C/r_i) evaluate d(x) = sum_j w_j phi_j(x) C = 15 * pi**1.5 / (2*sqrt(2)) ~ 29.5305

C is not a fudge factor: it is the exact bending energy of a unit Gaussian. The regulariser is INT sum_ab (d2phi_i/dx_a dx_b)(d2phi_j/dx_a dx_b) dV, which by Parseval equals (2 pi)^-3 INT |k|^4 phihat_i conj(phihat_j) d3k; for a Gaussian of width a the diagonal evaluates to 15 (2 pi)^1.5 / (8 a) = C / a. Verified against the binary to 16 significant figures.

DELIBERATE DIVERGENCES FROM THE REFERENCE

  1. The reference builds its solve matrix as K[i][j] = exp(-||t_i - t_j||^2 / r_i^2) – indexing the radius by ROW – while evaluating with the radius indexed by COLUMN. The two disagree whenever radii differ, so the reference does not actually interpolate: a landmark dragged 4mm can land ~2mm away from where the user put it. Under uniform radii (every real session inspected, and the 15mm default) the two indexings coincide and the schemes are identical. We use the consistent column indexing, which restores exact interpolation.

  2. The reference’s regulariser picks up off-diagonal structure under mixed radii that is neither symmetric nor the true bending-energy matrix, and is not identifiable from its output. We use the analytic diagonal.

Where the problem is well-posed these agree with the reference to ~3e-07 mm. Where it is not – dense smudge trails put cond(K) around 1e6, so the regulariser dominates – they agree to <0.2mm peak and ~0.003mm mean under realistic coherent drags, well under imaging resolution.

Functions

solve_weights(→ FloatArray)

Solve for RBF weights carrying each target back onto its source.

evaluate_on_points(→ FloatArray)

Evaluate d(x) = sum_j w_j exp(-||x - c_j||^2 / r_j^2) at points.

evaluate_on_grid(→ numpy.typing.NDArray[numpy.float32])

Accumulate the field over a regular grid, one bounding box per centre.

Module Contents

brain_deer.warp.rbf.solve_weights(source: FloatArray, target: FloatArray, radii: FloatArray | float, stiffness: float = 0.1) → FloatArray

Solve for RBF weights carrying each target back onto its source.

Parameters:

  • source ((N,3) arrays of correspondences, in consistent units.)

  • target ((N,3) arrays of correspondences, in consistent units.)

  • radii (scalar, or (N,) per-landmark basis widths in the same units.)

  • stiffness (regularisation factor. 0 gives exact interpolation, which is) – ill-posed when landmarks nearly coincide – and they routinely do, since a dragged trail is sampled every ~1mm. Non-zero is strongly advised.

Returns:

brain_deer.warp.rbf.evaluate_on_points(centres: FloatArray, weights: FloatArray, radii: FloatArray | float, points: FloatArray) → FloatArray

Evaluate d(x) = sum_j w_j exp(-||x - c_j||^2 / r_j^2) at points.

centres are the targets used in solve_weights().

brain_deer.warp.rbf.evaluate_on_grid(centres: FloatArray, weights: FloatArray, radii: FloatArray | float, *, affine: FloatArray, shape: tuple[int, int, int], cutoff_radii: float = DEFAULT_CUTOFF_RADII) → numpy.typing.NDArray[numpy.float32]

Accumulate the field over a regular grid, one bounding box per centre.

affine maps voxel index -> world; shape is the grid size. Returns an (I,J,K,3) float32 displacement array.

Each Gaussian is evaluated only inside cutoff_radii * r_j of its centre. At 3 radii the kernel is ~1.2e-4 of its peak, below float32 resolution on a millimetre displacement, so truncation is free. A bounding box answers “which voxels are near centre j” with no data structure and no query overhead; a KD-tree only pays off when centres are dense enough that the boxes overlap heavily, which does not happen at these radii. Pass inf to disable.