# Image registration Intensity- and landmark-based registration between the two loaded workspaces: three independent engines, per-run analytics, a deformation vector field you can see and export, and the option to restrict any of it to a single structure. Beside them, a registration by structures aligns the two images on the surfaces of structures contoured on both, without looking at a voxel value (*Align by structures*, below). elastix and plastimatch are C++/ITK toolboxes; nothing of either is linked - the algorithms are **re-implemented natively in Rust**. ![registration](screenshot_registration.png) *Deformable registration of two breathing phases: the module reports the recovered transform and metric improvement; the fusion overlay shows aligned anatomy in gray, residual respiratory mismatch as magenta/green fringes.* ## The three engines | | elastix rigid | elastix B-spline | plastimatch B-spline | plastimatch landmarks | by structures | |---|---|---|---|---|---| | transform | Euler 6-DOF | rigid + cubic FFD | centre-of-gravity + cubic FFD | RBF warp | Euler 6-DOF or translation, then optionally a local FFD per structure | | samples | ~3000 random, redrawn every iteration | same | every eligible voxel | none (geometric) | up to 4000 surface points per structure and side | | gradient | stochastic estimate | stochastic estimate | exact analytic | closed form | exact analytic | | optimizer | ASGD | ASGD | L-BFGS + line search | direct solve | Gauss-Newton with Levenberg-Marquardt damping | | metric | mean squares | mean squares | mean squares **or** Mattes MI | landmark residual | RMS surface distance (mm) | | regularizer | - | - | bending energy | stiffness | - | | deterministic | seeded | seeded | yes | yes | yes | | multi-modal | no | no | yes (MI) | yes | yes (no intensities) | ### elastix - stochastic sampling and ASGD A native re-implementation of elastix's own defaults - `Optimizer AdaptiveStochasticGradientDescent`, `ImageSampler RandomCoordinate`, `NewSamplesEveryIteration true` and `Metric AdvancedMeanSquares` - and the engine to reach for first: iterations are cheap, so thousands are affordable, and the estimate's noise carries the search past small local minima. * **Multi-resolution Gaussian pyramids** (`NumberOfResolutions`, default 3): [1 2 1]/4 smoothing + factor-2 decimation, voxel-centre origin bookkeeping. * **Random coordinate sampling**, fresh every iteration (`NumberOfSpatialSamples`, default 3000), within a body mask from a configurable HU threshold (default −500) and, for a local run, the region; drawn from a pre-built eligible-voxel list, so every draw is a hit. * **Metric:** mean squared difference with analytic gradients. * **Optimizer:** Adaptive Stochastic Gradient Descent (Klein et al., IJCV 2009 - elastix's default) with automatic gain estimation, the sigmoid time-adaptation rule and a trust-region step cap. * **Rigid:** 6-DOF Euler transform about the fixed-image centre, with automatic rotation/translation parameter scaling. * **Deformable:** the rigid result composed with a cubic B-spline free-form deformation (`FinalGridSpacingInPhysicalUnits`, default 32 mm), optimized coarse-to-fine across the pyramid. ### plastimatch - a dense exact gradient and L-BFGS Following plastimatch's `bspline` (Shackleford et al., *High performance deformable image registration algorithms for manycore processors*) - the opposite trade: far more work per iteration, far fewer of them. 1. `xform=align_center` - a translation matching the centres of gravity of the two thresholded images. Skipped for a local run or a refinement, which already start aligned. 2. `xform=bspline` per resolution level, coarse to fine: the cost and its **exact analytic gradient** over every eligible fixed voxel, each voxel's contribution scattered onto the 64 control points that support it. * **Metric.** `mse` is the mean squared difference divided by the fixed image's variance, so the cost is dimensionless. `mi` is **Mattes mutual information** over a 32 × 32 joint histogram, zero-order Parzen window on the fixed image and cubic B-spline window on the moving one (Mattes et al., IEEE TMI 2003) - the multi-modal (CT-MR, CT-CBCT) option. * **Regularizer.** `young_modulus` weights the discrete bending energy of the control lattice (second differences, mixed terms counted twice), made dimensionless by the lattice spacing; its gradient is exact. * **Optimizer.** L-BFGS (two-loop recursion, history 6) with an Armijo backtracking line search; plastimatch's default L-BFGS-B differs only in box constraints, which B-spline coefficients lack. "Dense" is capped at 400 000 samples per level, thinned deterministically so the set is the same on every iteration - on a 512³ study every eligible voxel is tens of millions. ### plastimatch - landmark warp `landmark_warp` interpolates paired points and reads no image intensity - for CT against MR, a post-operative cavity, anatomy that genuinely changed, or when specific anatomical points must be honoured and nothing else. | kernel | φ(r) | support | affine term | |---|---|---|---| | Thin-plate spline | `r` | global | yes | | Gaussian | `exp(−r² / 2R²)` | global, decaying | no | | Wendland ψ₃,₁ | `(1 − r/R)⁴ (4r/R + 1)` | compact, zero beyond `R` | no | The thin-plate spline minimizes bending energy over the whole domain; its affine term reproduces exactly any global shift or rotation implied by the landmarks, and it needs at least four non-coplanar pairs. The radial kernels have no affine term, so displacement decays to zero away from the landmarks; the compactly supported Wendland kernel *provably* leaves distant anatomy untouched. `stiffness` (plastimatch's regularization) is added to the diagonal of the interpolation matrix: zero passes exactly through every landmark, larger values smooth the field and tolerate inconsistent pairs. Put the crosshair on the same anatomy in both workspaces and press **➕ Add pair** in the *Landmarks* section (with *View ▶ Sync crosshairs* off, or both crosshairs move together). Each pair shows its displacement and, after a run, its residual. ## Running a registration *Modules ▶ Image registration* puts the section - the two images, method, region, parameters, landmarks, result and vector field - in the right panel. **Fixed image** and **Moving image** each name one image series of either workspace: the displayed ones, a phase of a 4DCT, or two series of one workspace (a cardiac CT and a phase of the 4DCT it arrived with). A series that is not on display is loaded for the run. The fixed image may also be *every phase of a 4D group*, which is one registration per phase (below). The run goes on a background thread with progress and a **Cancel** button. The result names both images, and **Fusion overlay on** chooses which one carries the overlay: on the fixed image the moving one is warped onto it, on the moving image the fixed one comes back through the inverse. The overlay, the vector field and the crosshair link appear in whichever workspace displays that image, so two series of one workspace show the fusion once the same folder is loaded as a second workspace too; propagation works either way. The transform maps **fixed → moving** patient coordinates, as in elastix, ITK and plastimatch; the inverse (for the crosshair link and propagation) is exact for the rigid part and a fixed-point iteration for the deformable one. **Start from** says where the search begins. The engines take steps of a few millimetres, so two images that do not overlap at the identity never find each other: a cardiac CT and a 4DCT of one patient are two acquisitions in two frames of reference, hundreds of millimetres apart in patient coordinates, and a run started from the identity has no gradient to follow (it now says so rather than returning the identity as a result). *Automatic* keeps the identity when the images overlap - every same-frame pair, so nothing changes for those - and matches the centres of gravity when they do not (elastix's `AutomaticTransformInitialization`, what the plastimatch engine always did as `align_center`). *Centroids of a structure* matches one structure contoured on both workspaces, which is the surest start for an organ: the heart on a cardiac CT and on a planning CT. A local run always starts from the identity. On the bundled data (two breathing phases of the 4DFBCT, 512 × 512 × 133 CT each): elastix rigid pre-alignment plus three B-spline resolution levels, 1800 iterations total, ≈ 20 s on a desktop CPU, driving the mean-squared HU difference from ≈ 9700 to ≈ 1800. ## The matrix, typed in by hand Under the *Register* row is a foldable **Transform matrix**: the same 4 × 4 a planning system shows. Three rows of direction cosines with the shift in the right-hand column, in patient millimetres, mapping fixed → moving like everything else here; the bottom row is `0 0 0 1` for any spatial transform, so it is shown and not editable. Sometimes the number is already known - a couch shift from the record, a transform from a planning system or a colleague, or a pure 5 mm translation as a sanity check - and typing it is faster and more honest than tuning a registration until it produces it. **Use this matrix** is the switch, and nothing reads the matrix unless it is on. While it is off the grid *follows the tool*: it shows the transform of the active registration, so ticking the switch takes those numbers over and edits them rather than starting from an identity nobody asked for. Once it is on the numbers are yours, and a later run does not reach in and change them; untick it and the grid follows again. A deformable result is not a matrix - what the grid shows of it is the rigid part it starts from, and the section says so, because using it drops the deformation. A run against a 4D group makes one transform per phase, not one. The row of buttons above the grid - **Of 0% · 10% · …** - is which of them the grid is showing; ten grids one under the other would be a wall of numbers, one grid and a picker is the same information a phase at a time. The active registration is the first entry, because that is the one *Apply as the registration* acts on. With nothing registered yet the section says so rather than showing an identity that means nothing. *Identity* puts back a matrix that moves nothing; *Invert* replaces it with the mapping the other way, and is disabled for a matrix that flattens space, because that one has no inverse; *From the result* pulls the run's transform back in after a hand edit has wandered. **📋 Copy** and **📥 Paste** move the sixteen numbers through the clipboard whitespace-separated, the common plain text form of a 4 × 4 matrix, so a matrix from elsewhere goes in without retyping (commas separate them just as well). **▶ Apply as the registration** installs it as the active registration of its two workspaces without running anything. Everything downstream reads the transform rather than the engine that made it, so the fusion overlay, the crosshair link, the vector field, propagation, the analytics and the REG export all follow it at once. The result is filed as *Given - a transform matrix, not a recovered one*: nothing was optimized, so the metric numbers a run reports would mean nothing and are not shown. The same editor, in the same shape, sits in *Structure propagation* ([propagation.md](propagation.md)) and in *Transfer by relationship* ([motion-4d.md](motion-4d.md#transfer-by-relationship-srcapptransfer_winrs)), where it stands in for that tool's own transform for the one run. ## Against a 4D group The **Fixed image** list ends with **every phase of a 4D group** of either workspace. That runs one registration per phase against the moving image: the phases of one acquisition differ by breathing, so a single transform for the group would be answering a question nobody asked. The moving image can be any series - of the group's own workspace (a planning CT or a cardiac CT beside its 4DCT) or of the other one. Each phase reports its own metric line, and the transforms are kept so that propagating structures onto the same group afterwards costs no registration ([propagation.md](propagation.md)). A propagation onto a group - anchored or not - leaves its per-phase transforms here in the same way. **The deformation field of a group run.** When a group run (a propagation onto a 4D group, anchored or not, or a registration against the group) finishes, the transform of the phase on display becomes the active registration, so the result section below and its **Vector field** settings are there as after any other run; nothing is switched on by itself. To draw the field: - in the phase table, the **Vector field** column: **👁 show** draws that phase's field on the views and in 3D (the phase is put on display first if it is not); **👁 on the views** hides it again; - or, in *Structure propagation*, the **Deformation field** row under *Last run*: one **👁** button per phase, with **deformation only** beside it; - or **Show the deformation field** under the result's *Vector field*. For a deformable transform the field starts as *deformation only* (see below), because an anchored run's rigid part is the jump from one frame of reference to the other. **Clear registration** then puts the phase away and keeps the group; **Clear group registration** drops the per-phase transforms, as does clearing a registration that was run here. ## Local registration Any method can be restricted to a **region** - an RTSTRUCT ROI or a painted segmentation of the fixed workspace, dilated by a margin. Three things change: * samples come from inside the region only; * the B-spline control lattice covers the region's bounding box, so a small structure can be aligned at a grid spacing unaffordable globally; * the centre of rotation and parameter scaling are the region's, not the patient's. A **local deformable** run skips the rigid stage on purpose: a rigid body fitted to one structure would move the whole volume. Confined to its lattice, the correction is exactly zero outside the region. A **local rigid** run instead reports how that structure moved *as a rigid body*; the transform is global. Without a margin nothing outside the structure constrains the boundary being aligned. ### Refining **▶ Refine** recovers a correction *on top of* the active registration: the moving image is sampled through the existing transform plus the new deformation, and the result is the two composed - typically a global registration, then a local refinement on the structure that matters, leaving the rest of the patient on the global result. ## Registration by structures *Align by structures*, a sub-section of the Image registration section, aligns the same two picked images on structures contoured on both of them and on nothing else: the voxel values take no part. It is the registration to reach for when the contours are trusted more than the grey values - a contrast CT against a plain one, a CT against an MR or a CBCT, two scans whose anatomy agrees and whose intensities never will - or when the alignment must follow one set of organs (the heart and the great vessels for a cardiac target, the spine for a setup check) and nothing else. `src/registration/shape.rs` is the engine; the propagation module's anchored run aligns its anchor by the same distance maps ([star-target-propagation.md](star-target-propagation.md)). **Which structures.** The table lists every structure drawn on the fixed image whose name (case-insensitive) is also drawn on the moving image: a structure set or segmentation series that references the series, or one that names no series of its workspace. Tick the ones to align on; *Pair by hand* adds a pair named differently on the two sides (`Heart` against `heart_total`). Each structure is normalised by its own number of surface points, so a large organ does not outvote a small one, and its **weight** (1 by default) then says how much it counts against the others: 0.1 for a structure that should only break a tie. **The rigid stage.** Each structure becomes a signed distance map on its own image's lattice - millimetres to the surface, negative inside, clamped at 40 mm so a point far away has no gradient to follow - computed in a box around the structure. The surface points are the midpoints of the faces between a voxel inside and a neighbour outside (a face on the edge of the image is not surface: that is a structure cut by the field of view). The cost is ``` E(θ) = Σ_k w_k [ 1/n_k Σ_i ρ(D_k^moving(T x_i)) + 1/m_k Σ_j ρ(D_k^fixed(T⁻¹ y_j)) ] ``` the fixed surface points `x_i` laid onto the moving structure's map and, with *Both ways* on (the default), the moving surface points `y_j` laid back onto the fixed structure's map through the exact inverse. The second term is what keeps a structure contoured over a shorter length on one image (a spinal cord, an oesophagus) from sliding along its partner. `ρ` is the square, or with *Robust* Huber's function of the width given, so a slice contoured differently on one side counts linearly rather than squared. The six parameters (three with *Translation only*, the rotations kept at zero) are found by Gauss-Newton with Levenberg-Marquardt damping on the exact derivative of the trilinear maps, rotations about the centre of the fixed surfaces: deterministic, a few dozen iterations, well under a second after the maps. Every sum over points is taken in fixed pieces and in order, so a run gives the same numbers on any number of threads. **Where it starts.** The section's *Start from* row applies: *Automatic* and *Centres of gravity* match the centroids of the ticked structures (the contour analogue of the centres of gravity, and what lets two images in two frames of reference find each other), *Identity* keeps the identity, and a structure there matches that structure's own centroids. **Then refine.** Optionally, a local B-spline refinement (elastix or plastimatch) runs on each structure's distance maps, one structure at a time from the largest to the smallest: its region is the structure grown by the *Margin*, its start the transform so far, so each is a local correction composed onto the rigid result and the rest of the patient keeps that result exactly. The grid spacing is the sub-section's own *Grid*; the resolutions, iterations and samples are the ones under *Parameters*. **▶ Refine active** skips the rigid stage and refines the active registration of the same two images instead - an intensity registration, say, corrected on the structures that matter. **What comes back** is an ordinary registration, installed as the active one: fusion, crosshair link, vector field, analysis, propagation and REG export work on it unchanged. Its metric is the RMS surface distance in millimetres before and after (`RMS mm 6.31 ▶ 0.42`), and *By structures* lists for every structure the mean surface distance and the Dice (the moving structure carried onto the fixed image, as *Score structures* carries it, against the fixed one) where the search started and where it ended, with each refinement's line, its 95th percentile displacement and folded fraction in the tooltip. A structure whose distance barely falls while the others meet is one whose two contours disagree: untick it, or lower its weight. The method aligns the contours it is given; a contouring difference becomes registration error by construction, which is what the per-structure table is there to show. One round structure alone does not fix the rotations (a sphere turns freely): use *Translation only*, or add a second structure. The same engine is the MCP tool `register_structures` ([mcp.md](mcp.md)) and the workflow step *Register by structures* ([workflows.md](workflows.md)). ## What the result says The result block reports method, region if any, metric before and after, and deformation model. The **Analysis** section is measured on the transform itself, on a lattice over the fixed image (or the region), so it means the same for every method: * **Image overlap** - the Dice coefficient of the two images' tissue, after the registration and before it. Every voxel at or above a tissue threshold (-300 HU for CT, a quarter of the way up its own value range for anything that carries no air) counts as tissue, and the score is the overlap of the fixed image's tissue with the moving image's, sampled on the same lattice as the rest of the analysis. It is the headline number because it is the one that says whether the result is usable at all: 0.80 and above reads as a good match (green), 0.60 to 0.80 wants a look (amber), below 0.60 is a failure to explain (red). It is an *image* score - it says the two workspaces now cover the same space, not that any one organ lines up. * **Best-fitting rigid body** - the orthogonal Procrustes fit: translation, three Euler angles in the same `Rz Ry Rx` convention as the rigid transform, and the RMS residual those six numbers do *not* explain. * **Displacements** - min / mean / p95 / max / RMS of `|T(p) − p|` in millimetres, plus the mean *vector* (systematic shift vs. scattered local motion). * **Jacobian determinant** - `det(I + ∂d/∂x)` by central differences: above 1 the tissue expanded, below 1 it compressed, at or below zero it folded. The folded fraction is reported; a regularized B-spline should show none. * **Per structure** - mean and maximum displacement over each contoured structure's own points: "the tumour moved 9 mm and the cord 0.4 mm" rather than "4 mm on average". *Score structures (Dice)* adds the anatomical half of the question: every structure of the fixed workspace is paired with the structure of the same name on the moving workspace (a contour or a segmentation, matched case-insensitively), the moving one is carried through this registration, and the overlap is scored against the fixed one - after the registration and before it, coloured by the same three bands. A good image score with a poor structure score is the case worth catching: the patient lines up, the organ does not. ## The fusion overlay and the vector field **Fusion** blends the transformed moving image into the fixed image's green channel (aligned anatomy gray, mismatch magenta/green) with a blend slider; the cross-study crosshair link maps through the recovered transform, inverse included. The **vector field** is the transform sampled onto a regular lattice - once, not per pixel on every repaint: a B-spline evaluation is 64 weighted lookups, a landmark warp a sum over every landmark. It is drawn in all three MPR views of the fixed workspace and, optionally, in the 3D window: * **Arrows** from where anatomy is to where it goes, exaggerated by an adjustable factor (millimetre motion is invisible at 1×) and coloured by magnitude; out-of-plane displacement becomes a disc sized by that component. * **Deformed grid** - the sampling lattice pushed through the deformation: warped graph paper, showing compression and expansion. * Lattice spacing, arrow scale and colouring are adjustable; changing the spacing re-samples on a worker thread. * **Deformation only** draws what the B-spline adds on top of the rigid alignment instead of the whole displacement. Between two frames of reference - a cardiac CT onto a 4DCT phase - the rigid part moves every point by the same hundreds of millimetres, and a field drawn whole shows that jump and hides the deformation. In the **3D window**, *Workspace B through the registration* meshes the other workspace's structures and maps every vertex through the recovered transform, so both anatomies stand in one frame of reference with independent opacities. The field can be overlaid as 3-D arrows in the same scene. ## DICOM interchange A rigid matrix from a DICOM **REG** object or a **Deformable Spatial Registration** object's displacement grid can be applied instead of running the optimizer; it becomes the active registration and everything downstream (fusion, crosshair link, analytics, propagation) works on it. See [rt-objects.md](rt-objects.md). **💾 Save as DICOM…** writes the active field out as a Deformable Spatial Registration. The IOD applies its grid between a pre- and a post-deformation matrix; both are written as the identity and the grid carries the whole mapping, `T(p) − p`. ## Propagating structures Once aligned, contours drawn on one workspace can be carried to the other - see [propagation.md](propagation.md). ## Transform simulator (registration QA) The *Simulation* module section applies an **exactly known** transform - rigid motion (translation + Euler rotation about the volume centre) plus an optional local Gaussian deformation (amplitude vector + σ, centred at the crosshair) - to a loaded workspace and generates the result into the other slot: the CT is resampled through the inverse transform; structure contours, dose grids and plan isocentres are carried along. The applied parameters stay displayed as ground truth. Any workspace, original or simulated, can then be exported as DICOM (see [export-and-tools.md](export-and-tools.md)). ## Accuracy verification `tests/registration.rs` registers analytically known transforms on a synthetic phantom, with the same tolerances for every engine: * **elastix rigid** recovers a known rotation + translation to ≈ 0.6 mm (asserted 1.5 mm); the inverse round-trips to 10⁻⁶ mm; the six-DOF analysis reproduces it to 10⁻³ degrees, zero residual, unit Jacobian. * **elastix B-spline** and **plastimatch B-spline** each recover a 7 mm Gaussian-bump deformation to ≈ 0.3 mm (asserted 3 mm), with no folding. * **plastimatch mutual information** recovers the same bump between images with *inverted* soft-tissue contrast, where mean squares has no minimum at the truth at all. * **landmark warp** lands on every landmark to 10⁻⁴ mm with all three kernels; the thin-plate spline reproduces a global shift everywhere, including far outside the landmark hull, and the Wendland kernel leaves points beyond its radius at exactly zero. * **local registration** recovers a displacement applied inside one blob and leaves every probe outside the region at exactly zero displacement; a refinement on top of a global result changes nothing outside its region. * **the vector field** reproduces the transform it was sampled from to < 0.05 mm. Unit tests check the Parzen window and its derivative against finite differences, the bending energy's gradient against a central difference and its vanishing on an affine field, the Procrustes fit against a reflection, region dilation by an exact margin, and the local lattice's coverage of its region. ## Notes Deformable results are intensity-driven: displacements inside large uniform regions are interpolated from the control lattice rather than measured - the Jacobian and per-structure displacements tell the difference. Mean squares assumes comparable intensities (CT-CT); for CT-MR use the plastimatch engine with mutual information, or place landmarks.