RITK: Medical Image Processing and Registration

RITK is a Rust toolkit for loading medical images, preserving their physical geometry, applying ITK-style filters, evaluating similarity, and producing classical or differentiable registrations. The book is organized as a user-facing workflow:

1. load a volume through ritk-io;
2. inspect shape, spacing, origin, and direction;
3. preprocess intensities and spatial structure;
4. choose a same-modality or multi-modal metric;
5. resample and validate the result; and
6. write the output without losing the spatial contract.

The examples use real RITK APIs and committed fixtures where a dataset is needed. Synthetic examples are deterministic, small, and input-sensitive so they can generate figures in CI.

First runnable workflow

use coeus_core::SequentialBackend;
use ritk_filter::IntensityWindowingFilter;
use ritk_io::read_image_native;

let backend = SequentialBackend;
let input = read_image_native("volume.nii.gz")?;
let windowed = IntensityWindowingFilter::new(-160.0, 240.0, 0.0, 1.0)
    .apply_native(&input, &backend)?;
println!("shape = {:?}, spacing = {:?}", windowed.shape(), windowed.spacing());

The NativeImage alias is Image<f32, SequentialBackend, 3>. RITK stores voxels in [depth, row, column] order while origin, spacing, and direction describe the physical frame. Filters generally preserve that metadata; resampling and transforms are the deliberate exceptions.

Workspace map

Build the book and examples

# Run tests with the repository's bounded native-test runner
cargo nextest run -p ritk-filter
cargo nextest run -p ritk-segmentation

# Build docs
mdbook build docs/book
mdbook test docs/book
cargo build -p ritk-filter --examples
cargo build -p ritk-io --examples
cargo build -p ritk-codecs --examples
cargo build -p ritk-analyze --examples
cargo build -p ritk-registration --examples
cargo build -p ritk-segmentation --examples
cargo build -p ritk-diffusion --examples

For native tests use the repository’s configured cargo nextest command. Doctests use cargo test –doc. A figure is valid only after its generating example succeeds and the rendered artifact has been inspected.

How to read the chapters

• Part I explains format boundaries and spatial-axis conventions.

• Part II builds filtering and segmentation pipelines from intensity, spatial, morphology, diffusion, and seeded-label operations.

• Part III develops diffusion-MRI acquisition, Q-ball ODFs, and tractography.

• Part IV covers metrics, transforms, classical and differentiable registration, and post-registration validation.

• Part V explains backend dispatch, zero-copy boundaries, and measurements.

• Part VI maps the public RITK surface onto Coeus, Leto, and Moirai.

• Provider migration summary

• Architecture decision records

DICOM Format Boundary

Single source of truth for DICOM file parsing and pixel-frame decode.

Ownership

ritk-io::format::dicom owns the DICOM Part 10 file parser and pixel frame decoder. ritk-dicom provides the backend trait implementations.

Boundary Surface

• DicomParseBackend: parses a Part 10 file into a backend-owned object.
• PixelDecodeBackend: decodes one frame from a backend-owned object using DecodeFrameRequest.
• DicomBackend: combines parse and decode without dynamic dispatch.

Spatial Contract

DICOM file-axis [x,y,z] maps to RITK [depth,row,col] via spatial.rs. Physical-space metadata (origin, spacing, direction) is preserved through the boundary.

Codec Ownership

ritk-codecs owns JPEG, JPEG-LS, JPEG 2000, RLE, PackBits, and native pixel primitive implementations. Native-owned JPEG syntaxes route exclusively through NativeCodecBackend.

Pixel Precision

Frame requests keep DICOM BitsAllocated and BitsStored separate. BitsAllocated selects the byte container; BitsStored identifies the meaningful magnitude and sign bits. Encapsulated JPEG, JPEG-LS, and JPEG 2000 headers must agree with BitsStored before samples reach the modality transform. Lossless signed JPEG uses the codestream precision’s sign bit. Lossy DCT JPEG with signed DICOM metadata is rejected because the transform output does not define a signed stored-sample representation.

Image readers reject a missing or malformed BitsStored attribute instead of inferring it from BitsAllocated. When HighBit is present, the DICOM object boundary requires it to equal BitsStored minus one, which establishes the right-justified sample layout used by native decoding.

decode_stored_pixel_frame returns typed stored integers before modality rescale; display decoding remains a separate rescaled f32 path.

Pixel output

Series writers encode unsigned 16-bit samples. Multi-frame output uses unsigned 8-bit samples for baseline JPEG and unsigned 16-bit samples for native, JPEG-LS, JPEG 2000, lossless JPEG, and RLE transfer syntaxes. BitsAllocated and BitsStored equal that sample width, HighBit equals the width minus one, and PixelRepresentation is zero. These attributes describe the encoded samples as required by DICOM PS3.5 section 8.1.1; compressed fragment lengths describe the codestream, not the decoded sample width. Source bit attributes are checked for consistency but never copied over the output representation.

RT Dose emits unsigned 32-bit samples after checking the supplied positive dose scaling and each rounded sample’s range. SEG derives one-bit or eight-bit storage from BINARY or FRACTIONAL segmentation, rejecting a contradictory declared width. BINARY pixels follow the least-significant-bit-first stream in PS3.5 D.1: frames share boundary bytes, and padding follows the complete value. The former MSB-first, per-frame-padded reader and writer contradicted this wire contract; both now use the specified packing.

Each series slice has its own linear rescale; a multi-frame object has one rescale for the entire volume. Constant input encodes as zero with the input constant as intercept. Nonconstant input maps its finite range onto the full unsigned sample range with nearest-integer rounding. Non-finite samples, unrepresentable ranges, invalid dimensions, malformed source bit attributes, and invalid spatial values return a DicomWriteError cause through the existing writer result. Decimal String components fit the 16-byte wire limit. Spacing is positive and direction cosines are orthonormal within a bound derived from one rounding of each component to single precision followed by the double-precision dot product. No orientation normalization occurs silently.

Writers prepare and serialize all output before creating a directory or replacing any file. This preflight retains the serialized volume in memory and preserves existing output on input or serialization failure. Filesystem failures during persistence can still leave partial output; this is not an atomic transaction across a series.

Diffusion metadata

read_dicom_gradient_scheme_from_file reads one classic single-frame volume, while read_dicom_gradient_scheme_from_files accepts one representative file per volume in explicit acquisition order. The reader uses only the standard top-level Diffusion b-value (0018,9087) and Diffusion Gradient Orientation (0018,9089) attributes defined by DICOM PS3.3 C.8.13.5.9. It validates finite s/mm² values and three finite direction components, then constructs a physically typed LPS GradientScheme. It does not infer volume grouping from a directory or guess private vendor tags; enhanced functional groups require a separate sequence-aware reader. For an unweighted frame, a finite zero b-value with no orientation is mapped to the required zero vector; nonzero weighting still requires an orientation.

Invariant

Every DICOM loader must reject before constructing Image<B,3> when the object declares SamplesPerPixel ≠ 1.

Example: DICOM to NIfTI Conversion

Analyze 7.5 to NIfTI converter example.

Source

crates/ritk-io/examples/dicom_to_nifti.rs

Description

This example loads an Analyze 7.5 image pair and saves it as a NIfTI file. It exercises the ritk-io DICOM reader and NIfTI writer boundaries.

Usage

cargo run --example dicom_to_nifti -- <input_dir> <output_nifti>

Verification

• Reads a scalar Analyze 7.5 image
• Writes a NIfTI file with correct affine metadata
• Round-trips through ritk-image boundary types

Example: DICOM Dump Utility

Dump DICOM file metadata and pixel values.

Source

crates/ritk-io/examples/dump_dicom.rs

Description

Utility to inspect DICOM file metadata, pixel value ranges, and spatial transformation data. Exercises the ritk-dicom parser boundary.

Usage

cargo run --example dump_dicom -- <dicom_file> [--json]

Verification

• Parses DICOM Part 10 file
• Extracts patient/study/series metadata
• Reports pixel value statistics (min, max, mean)

NIfTI Format Boundary

ritk-nifti is RITK’s native single-source-of-truth implementation for the single-file NIfTI boundary. It reads NIfTI-1 and NIfTI-2 .nii files, detects gzip-wrapped .nii.gz input, and writes either header version explicitly.

Ownership

ritk-nifti owns the NIfTI file reader and writer. ritk-io::format::nifti is a facade re-export. Analyze 7.5 .hdr/.img pairs belong to ritk-analyze; they are not interpreted as NIfTI by this crate.

The image convenience APIs support:

• three-dimensional f32 scalar images;
• four-dimensional f32 acquisition series;
• three-dimensional u32 label maps;
• NIfTI sform and qform spatial metadata; and
• NIfTI-1 and NIfTI-2 single-file streams, compressed or uncompressed.

NiftiDocument is the lossless transport surface. It retains the complete uncompressed single-file stream, including unprojected header fields, the extension indicator and blocks, and exact sample bits. Writing to .nii reproduces those bytes; writing to .nii.gz changes only the gzip framing. The gzip reader drains the stream through its checksum trailer before acceptance.

NiftiDocument also parses all ten fixed-width ritk-codecs::SampleType representations. It validates each standard datatype code against bitpix and the declared payload length, then retains the sample bytes exactly, including wide integers, signed zero, and NaN payloads. Transcoding between .nii and .nii.gz changes the compression framing only. The f32 image and u32 label convenience readers retain their narrower conversion contracts; use NiftiDocument for byte-preserving transport.

The typed header view exposes both transform codes and classifies the active forms as qform-only, sform-only, compatible in handedness, or conflicting in handedness. Both forms remain in the document when they differ: qform can describe scanner coordinates while sform describes a standard space. A handedness conflict is reported explicitly because the NIfTI-1 specification states that operations on such an image are unspecified. See the official NIfTI-1 FAQ, questions 19 and 21.

transcode_nifti_document performs .nii/.nii.gz framing conversion inside RITK. It validates and, for gzip output, finishes compression before opening the destination, so format or compression errors leave an existing output unchanged.

Converting a stored series

Format adapters that produce ritk-image-io::StoredSeries can construct a NIfTI document without decoding samples into f32 or creating an intermediate file. The conversion boundary checks the target’s declared capabilities and cross-volume constraints before returning the document:

let document = NiftiDocument::from_stored_series(
    "nrrd",
    &stored_series,
    NiftiVersion::Two,
    source_metadata_losses,
)?;
document.write("converted.nii.gz")?;

The source identifier names the adapter that produced the series. Pass every format-specific metadata field that the adapter cannot preserve in the shared model as a scoped FormatMetadataLoss; unresolved losses prevent document construction. This input describes source-side loss only. The NIfTI target also rejects unsupported coordinate maps, sample types, axes, calibration, and incompatible grids before a destination exists.

The ten fixed-width scalar sample representations retain their bit patterns; payload bytes are normalized to NIfTI’s little-endian order. Linear calibration is represented by NIfTI’s global slope and intercept without changing stored samples. Per-frame calibration is accepted only when every frame in each volume uses the same transform, and every volume must share that transform. NIfTI-1 writes spatial and scaling fields as 32-bit floats and rejects an affine that becomes singular after narrowing; NIfTI-2 writes them as 64-bit floats. Both write physical geometry in RAS millimeters with the RITK [depth,row,col] to NIfTI [x,y,z] axis mapping.

A SingleVolume series is rank 3. List and Unspecified acquisition axes are rank 4 even when they contain one volume, so the axis meaning survives conversion. Diffusion axes are rejected because NIfTI’s scalar fourth axis cannot preserve their gradient metadata through this adapter. DICOM/NIfTI adapter integration is a separate conversion step: the source reader must first create the shared stored series and report source metadata losses.

Spatial Contract

NIfTI file-axis RAS maps to RITK [depth, row, col] through the format boundary. The reader constructs each image directly as [nz, ny, nx] from X-fastest NIfTI raw bytes; the writer emits RITK ZYX flat data in that file order.

Affine Conversion

• Reader: file affine columns [x,y,z] become internal metadata columns [depth,row,col] = [z,y,x].
• Writer: sform columns are emitted as [internal_col, internal_row, internal_depth].

RITK physical metadata uses LPS coordinates, while NIfTI affines use RAS. The boundary performs the LPS/RAS sign conversion; callers must not pre-flip their images.

Acquisition Series

A repeated acquisition is represented by a NIfTI rank-4 image:

dim[0] = 4
dim[1..=3] = [nx, ny, nz]
dim[4] = number of volumes

The fourth axis can represent diffusion gradient directions, functional timepoints, or another repeated measurement. NIfTI stores this axis slowest, so the complete voxels for volume 0 are followed by volume 1, then volume 2, and so on. read_nifti_series preserves that acquisition order and returns one Image<f32, B, 3> per volume.

Every volume shares one shape and one physical grid because a NIfTI series carries one spatial transform. The series writers therefore reject an empty series or any volume whose shape, origin, spacing, or direction differs from volume 0. This prevents one header from silently describing only part of the written data.

Rank behavior

The single-volume and series APIs are deliberately asymmetric:

Returning volume 0 through the single-image API would discard the rest of an acquisition while reporting success. Conversely, a rank-3 image is a valid series of one and needs no caller-side rank branch.

Writing follows the same canonical representation. A one-image slice passed to write_nifti_series is written as rank 3 and remains readable by read_nifti. Two or more images are written as rank 4 with their count in dim[4].

Public API

use coeus_core::SequentialBackend;
use ritk_nifti::{
    read_nifti_series, write_nifti2_series, write_nifti_series,
};

let backend = SequentialBackend;

// All images must have the same shape and physical metadata.
write_nifti_series("diffusion.nii.gz", &volumes, &backend)?;
let decoded = read_nifti_series("diffusion.nii.gz", &backend)?;

// Select NIfTI-2 explicitly when its wider header fields are required.
write_nifti2_series("diffusion-nifti2.nii", &decoded, &backend)?;

read_nifti_series_from_bytes provides the same decoding contract for an in-memory .nii or .nii.gz payload.

Validation and Failure Semantics

The boundary validates dimensions, voxel byte ranges, datatype widths, spatial metadata, and gzip expansion limits before constructing images. A declared payload that ends after only some volumes is an error; the reader does not return the complete prefix as a partial series.

The writer reports the zero-based position and mismatched grid property when a volume cannot share the first volume’s spatial transform. Choose the single-volume API for one image and the series API whenever the input may contain repeated acquisitions.

Invariant

NIfTI parser/writer dependency changes stay behind ritk-nifti; callers in ritk-io, CLI, and viewer code consume the same authoritative API.

NRRD Format Boundary

Single source of truth for NRRD file I/O.

Ownership

ritk-nrrd owns the NRRD file reader and writer. ritk-io::format::nrrd is a facade re-export.

read_nrrd_header parses the bounded header without decoding the payload. Its result exposes canonical standard fields, the declared format version, comment lines, an effective custom key/value map, and source-ordered decoded key/value records, including repeated keys. The effective map follows NRRD’s rule that the last value for a repeated key is the stored value; the record list retains earlier values for inspection during format conversion. The parser limits headers to 16 MiB and 65,536 retained entries across standard fields, comments, and key/value records. These header records do not imply that payload conversion preserves every format’s metadata.

Field identifiers are case-insensitive. Teem’s equivalent spellings (byteskip, lineskip, datafile, axismins, axismaxs, centerings, blocksize, oldmin, oldmax, and sampleunits) map to canonical lowercase field names before duplicate checks. Empty comment strings are ignored; retained comment lines keep their leading # and source order. NRRD header lines use ASCII; non-ASCII bytes return a typed parse error. These delimiter, key, field, alias, and comment rules follow the Teem NRRD format definition, Sections 1.2, 1.6, 5, and 6.

Spatial Contract

NRRD file-axis [x,y,z] maps to RITK [depth,row,col] via crates/ritk-nrrd/src/spatial.rs. The reader constructs the tensor directly as [nz,ny,nx] from X-fastest NRRD raw bytes; the writer emits RITK ZYX flat data directly.

Direction Vectors

• Reader: NRRD space directions vectors [x,y,z] become internal metadata columns [depth,row,col] = [z,y,x].
• Writer: NRRD space directions are generated from internal columns [col,row,depth].

Every finite positive spacing retains its direction, including values below 1e-9. Stored-volume construction rejects an axis component whose direction-times-spacing product overflows or underflows to zero, before a format writer can emit unrepresentable geometry.

A rank-2 array can be embedded in either a two- or three-dimensional world. For a named patient space such as LPS or RAS, each direction vector and the origin use three world coordinates even though the array has two axes. RITK derives the missing slice direction from the plane normal and promotes the array to [1,Y,X] with a one-millimeter through-plane step. When no named space is declared, rank-2 vectors retain the two-component interpretation and are promoted into the XY plane. Component count follows the world-space dimension, independently of array rank, as specified by Teem’s NRRD space fields.

NRRD’s named patient spaces are normalized to RITK’s LPS millimeter coordinates. RAS reverses the first two physical components; LAS reverses the anterior/posterior component. Supported space units (mm, cm, m, um, and nm) scale both origins and direction vectors before spacing and orientation are derived. Unknown spaces, anonymous coordinate frames, and unrecognized units fail instead of being relabeled as patient coordinates. These mappings follow Section 4 of the Teem NRRD format specification.

Invariant

NRRD parser/writer dependency changes stay behind ritk-nrrd; callers in ritk-io, CLI, and viewer code consume the same authoritative API.

Diffusion gradient metadata

read_nrrd_gradient_scheme implements the NA-MIC DWI convention. One nominal DWMRI_b-value is combined with each DWMRI_gradient_XXXX squared norm to recover the per-volume effective b-value. The measurement frame maps gradient coordinates to world coordinates; RAS world coordinates are converted once to RITK LPS. When the field is absent, the reader uses the identity frame. A nonzero gradient therefore uses scanner coordinates directly unless an explicit frame is present; an all-zero baseline also needs no frame. See Teem’s NRRD specification, section 4 for the coordinate- frame contract. Every encoded nonzero weighting is preserved, including values below the scanner-input baseline threshold used by ritk-diffusion-scheme constructors. Missing indices, non-finite values, DWMRI_NEX, and B-matrix encodings fail explicitly rather than being guessed.

Document construction and retained metadata

NrrdDocument::new(series, comments, records) consumes a StoredSeries, a Vec<String> of comments, and a Vec<(String, String)> of custom key/value records. It constructs a document without an intermediate file. The borrowed series(), comments(), and records() accessors expose its samples and retained metadata. Construction uses the same validation as the writer: identity calibration, supported series axes, physical geometry, matching volume shapes, sample types, geometry and coordinate maps, and serialized header byte and entry limits. Generated header entries count toward the 16 MiB and 65,536-entry limits.

read_nrrd_document(path, budget) checks retained metadata and the regenerated header limits before decoding samples. Its header and attached payload remain bound to one opened source file. The resulting document keeps the stored sample representation and values through StoredSeries.

Retained order

Comments retain their leading #, text, and relative order. Non-empty forms such as # note and ##note are supported. The header parser drops comment lines containing only # and spaces; the constructor rejects those forms, non-ASCII text, embedded carriage returns or newlines, and caller-supplied copies of # Complete NRRD file written by ritk. The document reader filters that generated banner, and each write emits it once.

Custom records are decoded key/value strings, kept in their supplied or source order, including repeated keys. For example, source:=scanner followed by source:=reviewed remains two records in records(), while the header’s effective map reports reviewed. Writing escapes backslashes again. Keys must be non-empty ASCII, must not start with # or contain :=, and keys and values must not contain carriage returns or newlines. Values must also be ASCII.

The writer groups retained comments before custom records, after generated fields. It preserves order within each list, rather than their original interleaving with one another or with standard fields. Whitespace, escaping, field spelling, and the original header bytes are not a document round-trip contract.

Regenerated fields and typed loss

Serialization derives type, dimension, sizes, kinds, space, space units, space directions, and space origin from the typed series. Output uses patient LPS millimeter coordinates, encoding: raw, endian: little, and an attached payload. Acquisition volumes use the trailing contiguous axis. Coordinate-map records and diffusion metadata (modality:=DWMRI, DWMRI_b-value, gradient records, and the identity measurement frame) are regenerated where applicable. Output uses NRRD0005 for diffusion and NRRD0004 otherwise; source version and encoding are not retained as formatting choices. Callers supply the corresponding typed values, rather than copies of generated records in the custom-record list.

The current document reader rejects standard fields outside that generated structural set, including content, labels, data file, spacings, and skip fields. The constructor likewise rejects standard field names passed as custom records. These failures return NrrdDocumentError::UnsupportedField; the error’s field may describe a metadata category rather than the specific source field. Reserved coordinate-map or DWMRI records, malformed or repeated generated records on input, and unsupported DWMRI extensions are also rejected. read_nrrd_header remains available to inspect those fields. The stored-sample reader’s support for a payload form, such as detached data, does not imply that the document boundary can retain its header semantics.

Validation before destination mutation

write_nrrd_document builds the bounded output header and validates the whole series before File::create can create or truncate the destination. Metadata conflicts, unsupported fields, header limits, calibration, and series mismatch errors therefore leave an existing destination unchanged. Construction performs these checks too; writing rechecks them before opening the destination. This is a validation guarantee: a later I/O error during header or payload writing, or flushing, can leave partial output.

Stored samples and format conversion

NrrdDocument combines a validated StoredSeries with retained comments and custom records without an intermediate file. Writing derives structural fields from the samples, validates first, and leaves an existing destination unchanged on rejection.

The document boundary retains only the generated structural fields and the custom records it can emit unchanged. It rejects standard fields such as content, labels, data file, and spacings, detached payload references, and malformed generated records rather than dropping them. It also rejects empty # and ## comments, non-ASCII comments, and the writer’s generated banner. The header parser retains non-empty forms such as # note and ##note; the document boundary preserves those values unchanged. The document constructor runs the same header, series, geometry, calibration, and entry-limit checks as its writer, so construction cannot produce a value that serialization would reject.

Use read_nrrd_stored when NRRD is an input to a format conversion. It returns ritk_image_io::StoredVolume, retaining the element type, each stored value, the spatial metadata, the coordinate map, and the calibration state. The ordinary read_nrrd API remains the compute path and converts values to f32.

The stored reader supports signed and unsigned 8-, 16-, 32-, and 64-bit integers plus IEEE 754 32- and 64-bit floats. It accepts the aliases listed in the NRRD type table and both endian markers. NRRD’s encoding field is required. RITK reads raw, ascii (text, txt), and gzip (gz). Binary multi-byte samples require an explicit endian marker; one-byte samples and ASCII samples do not, as specified in Section 5 of the Teem NRRD format. ASCII values are whitespace-delimited and each token is limited to 128 bytes to bound parser scratch space. Readers consume the declared array payload and ignore following bytes, as permitted by the NRRD encoding rules in Section 5. Binary reads and writes retain signed zero and NaN payload bits. ASCII values are parsed from text, so they do not carry binary NaN payload bits. Element type names follow the Teem NRRD type table and payload rules, including the standard int8 through uint64 aliases. RITK also accepts the historical char alias as signed 8-bit for existing files, although Teem does not list bare char as a NRRD type descriptor.

read_nrrd_stored_series returns one stored volume per acquisition entry. It preserves acquisition order for both leading interleaved axes and trailing contiguous axes. write_nrrd_stored_series writes the trailing contiguous layout; the existing compute-image series writer keeps its leading NA-MIC layout.

Diffusion series output uses the NRRD0005 magic and writes an explicit identity measurement frame because its gradient vectors are already in LPS coordinates. Other stored series retain NRRD0004. There is no implicit diffusion default: modality:=DWMRI, a declared acquisition axis, one nominal DWMRI_b-value, and contiguous gradient entries are required before a diffusion series is constructed. The reader accepts the standard axismins, axismaxs, and centerings aliases and canonicalizes them before checking duplicate fields. Per-axis physical units without space directions or spacings are rejected because they cannot determine sample geometry. Stored scalar or list data with a measurement frame but no DWMRI scheme is rejected before payload reading because the stored-volume model cannot retain that frame. Axis support bounds and cell/node centering are likewise rejected instead of being discarded.

The custom RITK coordinate-map field serializes Cartesian, curvilinear, phased-array, and per-slice transform maps. Slice-series transforms store nine row-major direction components and three translation coordinates per depth slice, with round-trip decimal formatting.

NRRD does not have a standard modality-calibration field. The stored writer therefore accepts only value-identity calibration and returns NrrdStoredWriteError::UnsupportedCalibration before it creates an output file when the volume carries a non-identity linear or lookup-table transform. Convert to a target format that can encode the calibration, or apply it explicitly before choosing a format whose contract stores only pixel values.

Detached data files must be relative to the NRRD header’s directory. Absolute paths and parent traversal are rejected at the read boundary; raw and decompressed payloads are bounded by the declared shape and element type. Detached data currently accepts one file name. line skip and byte skip fields are parsed before payload decoding; positive byte skips on gzip data apply after decompression, while byte skip: -1 is accepted only for raw data.

Stored reads accept an ImageReadBudget that caps encoded payload bytes, decoded sample bytes, gzip-expanded payload bytes including discarded byte-skip data, and the number of volumes in an acquisition series. The default limits are 1 GiB for each byte count and 65,536 volumes. The adapter checks declared counts before allocating decoded sample storage. The stored writer validates format semantics before opening the destination and streams samples through a buffered file writer without making a second volume-sized encoded payload.

MetaImage Format Boundary

MetaImage is a lightweight lossless boundary for header-driven volume interchange. The native facade accepts both single-file MHA and header-plus-raw MHD inputs and preserves validated shape, spacing, origin, and direction.

let image = ritk_io::read_image_native("volume.mha")?;
ritk_io::write_image_native("volume_copy.mha", &image)?;

Use MHA when one self-contained file is preferable. Use MHD with a raw payload when an external pipeline already expects separate header and voxel files. A round trip must compare shape and physical metadata in addition to voxel values. Readers construct Coeus-backed images directly on the selected native backend; writers extract host data only at the format boundary.

Example Summary

Analyze 7.5 Format Boundary

Analyze 7.5 stores one logical dataset in two files with the same stem:

brain.hdr   348-byte binary header
brain.img   contiguous binary voxel payload

The original format description is the Analyze 7.5 File Format specification, especially the header organization and image-dimension definitions on pages 3–8. The header gives dimensions, scalar type, bit depth, spacing, and payload offset. The payload contains no delimiters: its expected length follows from the dimension product and bits per voxel.

RITK’s supported subset

RITK deliberately reads a narrow, unambiguous subset:

Big-endian files, four-dimensional series, complex values, RGB values, and negative per-image offsets are rejected with contextual errors. Rejecting an unsupported variant is preferable to decoding it with the wrong byte order or geometry.

The extensions are not sufficient format identification. A paired NIfTI-1 dataset uses ni1\0 at header bytes 344–347 and commonly includes the four-byte extension indicator after the shared 348-byte header. RITK reports that case as paired NIfTI instead of decoding NIfTI affine fields as Analyze history fields. Use the .nii single-file form with RITK’s NIfTI reader. The NIfTI-1 FAQ documents the shared pair layout and extension indicator; the SimpleITK I/O list documents NiftiImageIO as the standard handler for these extensions.

The public API is backend-bound but format behavior is shared:

use coeus_core::SequentialBackend;
use ritk_analyze::{read_analyze, write_analyze};

let backend = SequentialBackend;
let image = read_analyze("brain.hdr", &backend)?;
write_analyze("copy.hdr", &image, &backend)?;

Passing brain.img to read_analyze is equivalent; the reader derives both sibling paths from the stem.

Axis order and byte order

The file describes dimensions as [x, y, z]. X varies fastest in the payload:

file_index(x, y, z) = x + nx·y + nx·ny·z

RITK stores a three-dimensional image with shape [z, y, x] and the same X-fastest flat order:

ritk_index(z, y, x) = z·ny·nx + y·nx + x

The flat sequences are identical, so reading does not transpose or copy an intermediate volume. The header spacing does require a semantic reorder: file [sx, sy, sz] becomes RITK tensor-axis spacing [sz, sy, sx].

Scalar conversion and scaling

datatype selects the stored scalar and bitpix must match it. RITK checks that pair before calculating payload bytes. Each stored value is converted to f32; integer values outside binary32’s exact integer range and finite f64 values can round under that explicit output contract.

The historical funused1 field is used by several Analyze-derived writers as an intensity scale. RITK applies a finite nonzero scale after scalar decoding; zero means one. This convention is not uniform across every Analyze variant, so a foreign pipeline should verify representative values rather than infer calibration from the filename.

Spatial metadata limits

Analyze pixdim[1..3] carries voxel spacing. Non-finite values are rejected; legacy zero or negative spacing is normalized to unit spacing for compatibility.

That normalization deliberately differs from the NIfTI reader, which reports a non-positive pixdim as a scoped spatial loss. NIfTI has sform/qform alternatives and can still describe the volume when pixdim is unusable; Analyze has no other spatial field, so a loss would leave the volume with no geometry at all. Normalizing is the only outcome that returns an image, so the divergence is a property of the formats rather than an inconsistency between the readers.

The original header defines originator as ten history bytes, not a complete world-space transform. RITK’s writer uses the common five-i16 convention and stores rounded voxel coordinates. A physical origin therefore round-trips only to the nearest voxel. The format has no direction matrix, so the RITK reader returns identity direction and the writer rejects a non-identity direction before creating either file rather than writing a dataset whose geometry silently disagrees with the source. Analyze files cannot establish scanner-space orientation as precisely as NIfTI, DICOM, NRRD, or MGH.

The NIfTI-1 rationale explains why NIfTI retained the 348-byte pair layout while adding explicit coordinate-system semantics. Because both formats can use .hdr and .img, select the reader from an authoritative source; extensions alone do not prove which format is present.

Bounded decoding and writing

Before allocating output, the reader validates signed dimensions, checked voxel and byte products, datatype/bit-depth agreement, finite metadata, offset, and exact file length. It then fallibly reserves the final Vec<f32> and streams conversion through an 8 KiB fixed buffer. Peak decoder-owned storage is therefore the returned f32 volume plus constant scratch, not the complete encoded payload plus the returned volume.

The writer validates dimensions, value count, checked byte size, finite spacing representable in the header’s f32 fields, and origin voxel coordinates within the format’s i16 range before creating either file. Header spacing can round from RITK’s f64 metadata to f32. The writer streams little-endian voxel bytes through an 8 KiB buffer and publishes the header after the payload completes; it does not construct a second volume-sized byte vector.

Failure behavior

Reading returns no partial image when any contract fails. Errors identify:

• invalid or unsupported dimensionality;
• unsupported endianness or scalar type;
• paired NIfTI presented through the ambiguous .hdr/.img extensions;
• mismatched datatype and bitpix;
• non-finite spacing, scaling, or offset;
• offset and byte-count overflow;
• truncated or trailing payload data;
• allocation, seek, read, or image-construction failure.

Writing rejects zero or oversized dimensions, storage/shape disagreement, non-identity direction, and invalid spatial metadata before creating output. I/O failures remain visible with their source chain.

Checked conversion matrix

Each row maps one Analyze pair property to the observable RITK behavior, the representable-or-lossy verdict, and the value test that pins it. Tests are identified as path::test relative to repos/ritk; run them with cargo test -p ritk-analyze (and -p ritk-io for the dispatch row).

A round trip alone cannot prove the axis order, because a self-consistent transposition would pass it. analyze_writer_emits_pixdim_in_file_axis_order therefore reads the raw header bytes at offsets 80, 84, and 88 and pins pixdim[1..3] = [3.75, 2.5, 1.25] for core spacing [1.25, 2.5, 3.75].

Next

The Analyze round-trip example uses a shared display scale and a separate absolute-difference panel so visual similarity is not mistaken for proof of equality.

Example: Analyze 7.5 Round Trip

This example writes and reads a deterministic 4 × 64 × 64 f32 volume through the public Analyze API. The phantom contains concentric structures, an off-center marker, a z-dependent diagonal, and a per-slice intensity offset. Those features make axis reversal, slice duplication, and value corruption visible and detectable.

How to read the figure

The first two panels use the same display range [0, 1100]. They should look identical because the writer stores f32 and the reader returns those values without quantization. Similar appearance alone is not the oracle.

The third panel is the explicit comparison:

• dark blue means the displayed source and decoded voxel bits are equal;
• red would identify any mismatch;
• max |decoded − source| = 0 is calculated from the actual slices.

The executable also compares every voxel in all four slices before generating the SVG. A wrong shape, geometry field, slice order, or value stops the example.

Public API

use coeus_core::SequentialBackend;
use ritk_analyze::{read_analyze, write_analyze};

let backend = SequentialBackend;
write_analyze("phantom.hdr", &source, &backend)?;
let decoded = read_analyze("phantom.hdr", &backend)?;

The .hdr file is 348 bytes. For shape [4, 64, 64], the .img payload is 4 × 64 × 64 × 4 = 65,536 bytes. The final panel reports those measured file sizes and the geometry that the format can reproduce under RITK’s convention.

Direction is identity because Analyze 7.5 has no full direction matrix. The chosen origin is an exact integer multiple of spacing on each physical axis; otherwise the originator convention would round it to the nearest voxel.

After producing the figure, the executable appends one byte to the temporary payload and confirms that the reader rejects the trailing data. The malformed pair never becomes an image.

Run it

From the repository root:

cargo run -p ritk-analyze --example book_analyze -- \
  docs/book/figures/analyze_roundtrip.svg

The Analyze pair is temporary. Only the deterministic SVG is retained. The complete source is crates/ritk-analyze/examples/book_analyze.rs.

MGH and MGZ Format Boundary

MGH is FreeSurfer’s native volume format. MGZ stores the same byte stream inside gzip compression. RITK exposes both through the ritk-mgh crate:

• read_mgh reads one-volume .mgh, .mgz, and .mgh.gz files, with ASCII case-insensitive suffix matching, and rejects multi-frame input;
• read_mgh_series reads every frame of an acquisition or time series;
• write_mgh writes the representation selected by the path extension;
• MghReader is a stateless adapter whose read method receives the backend;
• MghWriter retains a backend for repeated operations.

FreeSurfer describes MGH as an internal format whose field contract is defined by its reader and writer, and its tools report dimensions, frame count, voxel size, orientation, and voxel-to-RAS transforms. See the FreeSurfer MGH format note, mri_info reference, and mri_convert reference.

On-disk organization

An uncompressed RITK-readable file has two regions:

byte 0                                                        end of file
┌──────────────────────── 284 bytes ────────────────────────┬──────────────┐
│ version, dimensions, frames, type, DOF, RAS geometry, pad │ voxel bytes  │
└───────────────────────────────────────────────────────────┴──────────────┘
                                                               x → y → z

All numeric header fields and voxel scalars are big-endian. The first voxel axis varies fastest, followed by the second and third axes. RITK maps that ordering to the crate’s [z, y, x] image shape without transposing the stored voxel sequence.

The reader accepts four on-disk scalar types:

The writer emits MRI_FLOAT. Float input therefore round-trips bit for bit, including signed zero and finite values. Integer input is converted according to Rust’s integer-to-f32 conversion; integers outside the exact binary32 integer range can round.

Frames are a dimensional contract

The MGH header can describe several consecutive frames with common geometry. A diffusion series or time series may therefore have nframes > 1. RITK’s current MGH API returns Image<_, _, 3>, which can represent one volume but cannot represent a fourth frame axis.

The single-volume read_mgh entry point consequently requires nframes == 1. It rejects a multi-frame file and names its declared frame count. Use read_mgh_series when every frame of an acquisition or time series is required. Returning only frame zero would be more dangerous than rejecting the file: the operation would report success after silently discarding the rest of an acquisition.

This is a type-boundary decision, not a limitation of gzip or scalar decoding. The series API exposes that additional dimension as an ordered collection of 3-D images with shared geometry.

RAS geometry

When goodRASFlag == 1, MGH stores:

• voxel spacing D = diag(d_x, d_y, d_z);
• three direction-cosine columns in Mdc;
• the physical RAS coordinate c_ras of the volume center.

RITK stores the physical coordinate of voxel index zero. For dimensions (width, height, depth), define

h = [(width - 1)/2, (height - 1)/2, (depth - 1)/2]ᵀ

Then the header’s RAS corner is

origin_ras = c_ras - Mdc · D · h

The reader moves that point into the stored model’s LPS frame by negating its x and y components, and negates the same two components of every direction column. Writing applies both inverses; the flip is its own inverse:

c_ras = flip(origin + Mdc · D · h),   flip(v) = [-v_x, -v_y, v_z]

The header orders spacing and Mdc by the x, y, z voxel axes, while RITK orders image axes [depth, row, col] = [z, y, x]. The reader therefore reverses both: RITK spacing[0] is the header’s d_z, and RITK direction column 0 is the header’s Mdc[:, 2]. This is the same file-axis to internal-axis reversal NIfTI, NRRD, and MetaImage apply (docs/architecture.md §7–§9).

The frame conversion matches ritk-nifti, which negates the same two rows of its RAS affine (docs/architecture.md §21). An MGH volume and the same volume read from NIfTI therefore agree on both axis order and handedness: one physical point has one coordinate through either reader.

Mdc is stored transposed — FreeSurfer’s header descriptor names the field the transpose of the direction-cosine matrix — so its three consecutive triples are the columns consumed above. No transpose is applied on read or write.

When the RAS flag is absent, the reader uses zero origin, unit spacing, and identity direction. Applications that require scanner-space agreement should inspect geometry before combining volumes; equal array dimensions alone do not imply equal physical space.

Bounded streaming decode

The reader validates version, dimensions, frame count, scalar type, and geometry before constructing an image. It then converts the payload through a fixed 16 KiB input scratch buffer directly into the final Vec<f32>.

The output allocation grows only after the corresponding input bytes have been read. This matters for untrusted files: a header can declare a large volume, but a truncated payload cannot force the reader to commit the complete decoded allocation before proving that data exists. Multiplication of dimensions and byte counts uses checked arithmetic, and allocation failure is returned as an error.

For a 256 × 256 × 256 MRI_FLOAT volume, the decoded image itself is 64 MiB. The former whole-payload path additionally retained another 64 MiB encoded buffer while converting it. The streaming path retains approximately the decoded output plus 16 KiB of input scratch. The output vector can have unused geometric capacity, and a reallocation can temporarily involve both its old and new allocations. This is an allocation model, not a process-RSS claim: allocator, backend, gzip, and image-construction state still contribute to observed resident memory.

On the committed 128 × 128 × 64 public-reader benchmark, streaming also reduced median read time by 11.7% for MGH and 15.9% for MGZ on the development host. The benchmark includes file open, optional decompression, endian conversion, and image construction; it does not isolate disk hardware or claim the same percentage for every host.

Failure behavior

The reader returns errors for:

• unsupported format versions or scalar type codes;
• non-positive or overflowing dimensions;
• frame counts other than one;
• truncated headers or voxel payloads;
• invalid spatial metadata or image construction;
• allocation failure.

Truncation errors name the first voxel whose complete encoded value could not be confirmed. MGZ decompression errors retain their gzip context.

Next

The round-trip example writes both representations, checks their voxel and geometry contracts, and makes the exact reconstruction visible with an absolute-difference panel.

Example: MGH and MGZ Round Trip

This example creates one deterministic 15 × 96 × 96 MR-like volume with non-default origin, spacing, and direction. It writes the volume as both MGH and MGZ, reads both files through the public API, and checks:

1. shape equality;
2. bit-exact f32 voxel equality;
3. origin, spacing, and direction equality;
4. rejection of a complete two-frame MGH input.

The central-slice images are expected to look identical. The fourth panel is the decisive comparison: it displays the per-voxel absolute error, uses red for any nonzero mismatch, and reports the maximum error numerically.

In this deterministic input, gzip reduces MGZ to a small fraction of the MGH size because the synthetic tissues contain large repeated regions. That ratio is an example result, not a general medical-image compression guarantee.

Public API

The core workflow is:

use coeus_core::SequentialBackend;
use ritk_mgh::{read_mgh, write_mgh};

let backend = SequentialBackend;
write_mgh(&source, "brain.mgh", &backend)?;
write_mgh(&source, "brain.mgz", &backend)?;

let decoded_mgh = read_mgh("brain.mgh", &backend)?;
let decoded_mgz = read_mgh("brain.mgz", &backend)?;

Compression is selected by .mgz or .mgh.gz, matched without ASCII case sensitivity; the caller does not select a second codec API. The uncompressed .mgh path is useful when another tool requires direct field access, while .mgz is the normal storage-efficient representation.

read_mgh deliberately rejects a file whose header declares more than one frame. Use read_mgh_series to preserve every volume of an acquisition or time series; the single-volume API never returns only frame zero.

Why the explicit difference panel matters

Placing source and reconstruction side by side is necessary but insufficient: a few changed voxels, a small intensity bias, or a display-window mismatch can be invisible at page scale. The runnable example compares f32::to_bits() for every voxel before producing the figure. The image panel then communicates the same result:

error(x, y) = max(
    |MGH(x, y) - source(x, y)|,
    |MGZ(x, y) - source(x, y)|
)

The all-blue error panel means every value is zero. Any nonzero value would be red, and the figure generator fails before writing the SVG.

Geometry check

The example uses:

origin    = [12.5, -18.25, 32.0] mm
spacing   = [0.75, 0.75, 1.5] mm
direction = [[0, -1, 0],
             [1,  0, 0],
             [0,  0, 1]]

These binary-exact values exercise the RAS-center conversion without requiring a plotting tolerance. The example fails if either decoded image changes any geometry component.

Run it

From the repository root:

cargo run -p ritk-mgh --example book_mgh -- \
  docs/book/figures/mgh_roundtrip.svg

The program uses temporary files for MGH, MGZ, and the malformed multi-frame fixture. Only the deterministic SVG is retained. The complete runnable source is crates/ritk-mgh/examples/book_mgh.rs.

MINC2 Format Boundary

MINC2 stores medical images in an HDF5 hierarchy. The format separates voxel arrays from named dimensions and supporting scan metadata. The MINC2 format reference defines the hierarchy and attributes; the later MINC2 design paper explains its HDF5 organization, coordinate model, scaling, chunking, and multiresolution design.

Hierarchy and RITK’s current profile

A full MINC2 namespace can contain:

/minc-2.0
├── dimensions
│   ├── xspace
│   ├── yspace
│   └── zspace
├── image
│   └── 0
│       ├── image
│       ├── image-min
│       └── image-max
└── info

The standard’s minimal-file rule requires the image and its dimensions; the three principal groups form the recommended framework. RITK currently reads one three-dimensional contiguous image dataset and writes one contiguous little-endian f32 dataset. The reader supports global and first-spatial-axis per-slice integer scaling. It does the same for boolean storage, treating false and true as stored codes 0 and 1. It does not yet read chunked/compressed datasets, expose arbitrary metadata under info, or write multiresolution levels. Those cases return an error where they are detectable.

This restricted profile is useful for exact RITK round trips. It is not a claim that every MINC2 variant is supported. Validate representative foreign files before using the reader in a quantitative or clinical pipeline.

Dimensions and physical coordinates

MINC names spatial axes xspace, yspace, and zspace. The dimorder attribute on the image dataset gives the slowest-to-fastest array order. The dimension attribute specification defines:

RITK’s three-dimensional shape is [z, y, x]; x remains the fastest-changing voxel index. With ordered dimension records (d_0,d_1,d_2), the physical point for image index (i) is

[ p(i) = o + D,\operatorname{diag}(s),i, ]

where (o_k=\text{start}(d_k)), (s_k=|\text{step}(d_k)|), and column (k) of (D) is the corresponding direction cosine with the step sign absorbed. The writer requires finite origins, positive finite spacing, and orthonormal direction columns before creating a file.

Scalar conversion

The reader accepts contiguous HDF5 boolean, signed integer, unsigned integer, f32, and f64 payloads and returns f32. The current writer accepts Image<f32, B, 3> and preserves the little-endian IEEE-754 bits exactly.

For integer and boolean image data, the MINC pixel-conversion specification maps each stored value (v) to a real intensity (r):

[ r = r_{\min} + (v-v_{\min})\frac{r_{\max}-r_{\min}}{v_{\max}-v_{\min}}. ]

Here ([v_{\min},v_{\max}]) is the image dataset’s valid_range, whose endpoint order is insignificant. If it is absent, RITK uses the complete stored integer range. The real range ([r_{\min},r_{\max}]) comes from image-min and image-max. Those datasets may be scalar for one global mapping or contain one value per slice along the first spatial image axis. If both are absent, the MINC default real range is [0, 1]. Equal valid_range endpoints are malformed because the conversion denominator is zero, so the reader rejects them. For boolean storage, the default valid range is [0, 1]; false and true then follow the same mapping.

Values outside valid_range denote missing or uninitialized data. RITK’s Image<f32, B, 3> has no missing-value mask, so the reader returns a contextual error naming the first invalid voxel instead of silently inventing a value. Floating-point image datasets bypass image-min/image-max scaling, as the MINC conversion contract requires. Integer values beyond binary64’s exact integer range and finite f64 image values can still round under the explicit f32 output contract.

Bounded reading and writing

The reader validates exact dataset and dimension shapes and calculates element, slice, and byte counts with checked arithmetic. It reads at most 8 KiB of raw voxel bytes at a time and decodes directly into the returned f32 volume. A hostile shape that claims more voxels than the file backs therefore fails on the first unbacked chunk without reserving a second volume-sized byte buffer. Scaling is applied while each chunk enters the final output.

Before file creation, the writer checks the voxel product, the format’s i32 dimension limit, payload bytes, storage length, and physical metadata. It then converts at most 2,048 voxels at a time into one stack-backed 8 KiB little-endian scratch buffer. Writer-owned scratch is therefore constant with volume size; it no longer duplicates the complete volume as bytes.

use coeus_core::SequentialBackend;
use ritk_minc::{read_minc, write_minc};

let backend = SequentialBackend;
let image = read_minc("brain.mnc", &backend)?;
write_minc(&image, "copy.mnc", &backend)?;

Failure behavior

Reading reports malformed HDF5 structure, absent image or dimensions, invalid dimension attributes, unsupported storage layout or scalar type, byte-count overflow, truncated voxel data, malformed scaling ranges, incomplete image-min/image-max pairs, out-of-range stored values, and image-construction failure. Writing reports empty or unrepresentable axes, storage/shape disagreement, invalid physical metadata, allocation failure, and positioned-write failure. Logical preflight errors occur before the output path is created.

The MINC2 round-trip example makes the visual and numerical comparison explicit.

Example: MINC2 Round Trip

This example constructs a deterministic 9 × 64 × 64 MR-like phantom with an off-center lesion, asymmetric marker, slice-dependent texture, anisotropic spacing, nonzero origin, and a non-identity direction matrix. It writes and reads the volume through write_minc and read_minc.

How to read the figure

The source and decoded panels use the same display range [0, 1200]. They are expected to look identical, but appearance is not the correctness test. The third panel is a direct bitwise comparison:

• dark blue means the source and decoded f32 bits match;
• red would mark a corrupted voxel;
• the displayed mismatch count and maximum absolute error are computed from the actual central slices.

The executable also compares every voxel and every geometry field over the complete volume before writing the SVG. A wrong shape, axis order, origin, spacing, direction, or value stops generation.

The lower panel is intentionally different. It writes a second, foreign-style i16 MINC2 fixture with stored codes [0, 25, 50, 100] in each of two slices. Both slices use valid_range = [0, 100], but their real ranges differ:

• slice 0 maps to [-1000, 1000];
• slice 1 maps to [0, 200].

The displayed decoded rows come from read_minc. They expose the mapping as numbers so the scaling effect cannot be hidden by automatic image contrast.

What crosses the format boundary

The last panel reports the measured file size and the actual source geometry. The payload is one x-fastest contiguous little-endian f32 dataset. Shape, voxel bits, origin, spacing, and direction round-trip exactly in RITK’s current profile. The lower panel separately exercises per-slice quantitative integer scaling. Compression, chunking, patient/study metadata, integer-valued writing through the public writer, and multiresolution images remain outside this example.

The writer converts the 36,864 voxels in 2,048-value chunks. Its conversion scratch remains 8 KiB rather than growing to another 147,456-byte volume-sized buffer.

Run it

From the repository root:

cargo run -p ritk-minc --example book_minc -- \
  docs/book/figures/minc_roundtrip.svg

The .mnc file is temporary. Only the deterministic SVG is retained. The complete source is crates/ritk-minc/examples/book_minc.rs.

MRtrix .mif Format Boundary

.mif is MRtrix3’s container format: a plain-text header terminated by a bare END line, followed directly by raw binary voxel data in the same file (inline), or in a sidecar .mif.dat named by the header’s file: key (detached). RITK exposes both layouts through the ritk-mif crate:

• read_mif reads one-volume .mif and rejects multi-frame input;
• read_mif_series reads every frame of an acquisition or diffusion series;
• write_mif writes one volume; write_mif_series writes an interleaved multi-frame file;
• MifReader is a thin adapter whose read method receives the backend;
• MifWriter retains a backend for repeated operations.

MRtrix3 documents the format in its image format reference.

On-disk organization

byte 0                                                        end of file
┌────────────── text header, `END`-terminated ──────────────┬──────────────┐
│ mrtrix image: version 3.0, dim, vox, layout, datatype, …  │ voxel bytes  │
└───────────────────────────────────────────────────────────┴──────────────┘
                                                               x → y → z → frame

Header values are single-line (dim: 128 128 60) or multi-line blocks (transform: followed by four rows). A trailing backslash continues a line; # starts a comment.

The reader accepts float32, float64, int32, uint32, int16, uint16, int8, and uint8. The writer emits float32. A BE suffix selects big-endian; the default is little-endian. Every axis must span at least one voxel — a zero extent is rejected, because zero voxels would also make a truncated payload satisfy the length check.

Spatial convention

RITK tensors are [depth, row, col] = [z, y, x], while the .mif file names [x, y, z]. The raw payload is already in RITK’s flat order — X varies fastest in both — so no data permutation is needed; only the metadata order differs:

• transform is decomposed into RITK origin, spacing, and direction; the writer reorders its columns from internal [depth, row, col] to file [x, y, z].
• A header with no transform is axis-aligned identity, and its vox: triple — file [x, y, z] — is reversed into [Δdepth, Δrow, Δcol]. Reading the three components index-for-index would transpose the spacing against the dim line, which the same file decodes as [nz, ny, nx].

A round-trip test cannot catch that transposition: the writer always emits a transform, so the transform-less branch is reached only by files this crate did not write, and a writer/reader that agree on a permutation round-trip exactly. The axis order is therefore pinned by a hand-built-header oracle in ritk-mif’s reader tests, not by the round-trip suite.

ritk-io route

ritk-io::format::mif re-exports the codec and defines the local ImageReader/ImageWriter adapters, so ImageFormat::Mif and read_image_native/write_image_native reach .mif without a consumer depending on ritk-mif directly.

TIFF and BigTIFF Format Boundary

TIFF is a tagged raster container rather than one fixed pixel encoding. A file starts with a byte-order and version header, then links one or more Image File Directories (IFDs). Each IFD describes one raster through tags such as width, height, sample representation, photometric interpretation, compression, and strip or tile locations. The TIFF 6.0 specification is the format reference.

Classic TIFF stores file offsets in 32 bits. BigTIFF changes the header version and directory layout to use 64-bit offsets while retaining the same image-data model. The BigTIFF design documents the header and IFD differences. RITK reads both forms through the same API. Its public writer emits classic TIFF.

RITK’s volume model

TIFF defines a sequence of images, not a medical three-dimensional coordinate system. RITK gives that sequence one explicit interpretation:

IFD page 1    IFD page 2                     IFD page N
    │             │                               │
    ▼             ▼                               ▼
  z = 0         z = 1             ...         z = N - 1

each page: rows outside, columns inside  →  Image shape [z, y, x]

Every page must have the same width, height, and accepted color model. Page order is the linked IFD order. RITK does not sort pages by a tag or filename.

The public surfaces are:

Palette, grayscale-alpha, RGBA, CMYK, and YCbCr pages are rejected by these loaders. The grayscale loader also rejects an RGB page that appears later in an otherwise grayscale stack; a matching sample count cannot substitute for a matching color contract.

Scalar conversion

The grayscale reader accepts decoded unsigned integers, signed integers, and IEEE f32 or f64 pages. Every value becomes f32 in the returned image. The conversion is exact for every u8, u16, i8, and i16 value. Larger integers can round outside binary32’s exact integer range, and finite f64 values can round or overflow to infinity. The RGB reader applies the same conversion independently to its three explicit channels.

The writer stores Gray32Float, so an f32 voxel stack round-trips bit for bit. It does not quantize to an integer TIFF representation.

Physical-space metadata

The current RITK TIFF API does not encode or interpret a complete medical origin, three-axis spacing, and direction matrix. Reading therefore assigns:

origin    = [0, 0, 0]
spacing   = [1, 1, 1]
direction = identity

Writing preserves pixels and page order but omits the source image’s RITK geometry. This is an intentional boundary, not a registration result. A caller must restore geometry from an authoritative companion source before combining the volume with scanner-space DICOM, NIfTI, MGH, or another physical image.

Bounded decoding and memory

RITK validates each page’s color model and checked width × height × channels sample count before asking the decoder for its raster. Every decoded page is then appended directly into the final Vec<f32> after a fallible reserve:

• an f32 first page becomes the final allocation directly;
• later f32 pages append to that allocation;
• integer and f64 pages convert from the decoder-owned page straight into final storage, without constructing a second page-sized Vec<f32>.

The locked tiff 0.9.1 decoder also applies per-image, intermediate-buffer, and IFD-value limits through its Limits contract. Those dependency limits bound one decode operation; the complete output still scales with the number of accepted IFD pages. RITK detects allocation failure while growing that output and returns an error.

Failure behavior

Reading returns an error for:

• invalid TIFF/BigTIFF headers or directory structures;
• zero or overflowing page dimensions and RGB sample counts;
• unsupported color models;
• inconsistent page dimensions or a later page with another color model;
• decoded sample counts that disagree with declared geometry;
• decoder-limit, decompression, I/O, allocation, or image-construction failure.

The operation returns no partial image. The error identifies the page when the failure occurs.

Next

The multi-page round-trip example proves pixel and page-order equality with a difference panel, then makes the separate geometry boundary visible instead of implying that TIFF preserves scanner coordinates.

Example: Multi-page TIFF Round Trip

This example constructs a deterministic 5 × 64 × 64 volume. Slice intensity and a diagonal marker both vary with z, so reversing or duplicating IFD pages would change visible data and fail the exact voxel comparison. The source uses non-default medical geometry to exercise the separate metadata contract.

The first two panels intentionally use the same display range. They should look identical. The third panel is the decisive comparison: blue means every displayed voxel is bit-identical, red would mark any mismatch, and the numeric maximum must remain zero. The final panel is expected to differ because the current TIFF writer and reader do not carry RITK’s full physical geometry.

Public API

use coeus_core::SequentialBackend;
use ritk_tiff::{read_tiff, write_tiff};

let backend = SequentialBackend;
write_tiff(&source, "stack.tiff", &backend)?;
let decoded = read_tiff("stack.tiff", &backend)?;

The decoded shape is [number_of_IFDs, rows, columns]. In this example:

page 1 → z = 0
page 3 → z = 2
page 5 → z = 4

The program compares f32::to_bits() at every voxel before writing the figure. A page-order, value, or shape defect therefore stops the example instead of producing a plausible image.

Geometry restoration is explicit

The source uses origin [10, -20, 30] mm, spacing [0.5, 0.75, 2] mm, and a 90-degree axial rotation. The decoded TIFF correctly has zero origin, unit spacing, and identity direction under the current API contract. Applications must replace those defaults from a trusted source before physical resampling or registration.

Run it

From the repository root:

cargo run -p ritk-tiff --example book_tiff -- \
  docs/book/figures/tiff_roundtrip.svg

The TIFF itself is temporary. Only the deterministic SVG is retained. The complete runnable source is crates/ritk-tiff/examples/book_tiff.rs.

PNG Format Boundary

PNG is a practical import and export format for screenshots, microscopy slices, QA artifacts, and simple test data. The native readers decode a single image or stack a lexically ordered directory into a leading depth axis; the native writer emits a single [1, rows, cols] slice as 8-bit grayscale.

use coeus_core::SequentialBackend;
use ritk_io::format::png::native::{PngReader, PngSeriesReader, PngWriter};
use ritk_io::{ImageReader, ImageWriter};

let slice = ImageReader::read(&PngReader::new(SequentialBackend), "slice.png")?;
let volume = ImageReader::read(&PngSeriesReader::new(SequentialBackend), "slices")?;
assert_eq!(slice.shape()[0], 1);
assert!(volume.shape()[0] >= 1);

// One slice in, one file out. A volume is rejected, not truncated.
ImageWriter::write(&PngWriter, "out.png", &slice)?;

Series stacking is lexical, so zero-pad slice names when numeric ordering is required. PNG does not carry the same medical frame metadata as a volumetric format; assign or validate spacing and direction before registration.

Writing is windowed and lossy

PNG’s grayscale format is 8-bit unsigned, so writing an Image<f32> requires choosing a stored range. The writer maps the image’s own [min, max] onto [0, 255] and records nothing about it — PNG has no tag to record it in. A reader gets back the rank ordering and the shape, not the original scale. Callers who need the source values in the file’s own units must pre-window the image themselves. write_png requires nz == 1; use write_png_volume to emit a [depth, rows, cols] image as a directory of slices.

Example Summary

JPEG Format Boundary

JPEG is a visualization and interchange boundary, not a quantitatively exact archival format. It is appropriate for previews, overlays, reports, and lightweight exports after metric-sensitive computation has completed.

The native facade infers the format from the path:

let preview = ritk_io::read_image_native("aligned_preview.jpg")?;
ritk_io::write_image_native("aligned_preview_copy.jpg", &preview)?;

JPEG quantization changes intensities, so a round trip must use a bounded error or perceptual check rather than bitwise equality. Keep NIfTI, NRRD, or MetaImage as the quantitative source for registration and metric evaluation. After decoding, the result is an ordinary RITK image and can enter the same filter pipeline as any other input.

Quality-75 reconstruction oracle

The writer uses the Annex K luminance table and the quality scaling implemented by the locked jpeg-encoder 0.7.1 provider. At quality 75 the scale is 200 - 2 * 75 = 50, and each quantizer is Q = clamp((K * 50 + 50) / 100, 1, 255). The executable oracle checks the complete emitted table and the coefficients needed by four analytical 8 x 8 blocks: Q(0,0) = 8, Q(0,4) = 12, Q(4,0) = 9, and Q(4,4) = 34.

A constant centered level a has DCT coefficient F(0,0) = 8a. The fourth orthonormal basis has signs +--++--+; horizontal, vertical, and product blocks therefore have F(0,4), F(4,0), or F(4,4) = 8a. The coefficients 8a = 96, 72, and 136 are exact multiples of those quantizers, so forward quantization and inverse DCT reconstruct every integer sample exactly: the derived error bound is zero. Tests also require exact agreement with the independent decoder already used by the crate. A second oracle applies the Annex A forward DCT, emitted quality-75 quantizers, and inverse DCT to the rounded source gradient. Its exact decoded samples preserve source-sensitive coverage across the gradient’s odd-frequency coefficients. A changed quality table fails the table oracle, while an inverse-DCT normalization, covered basis, or level-shift mutation fails the analytical pixels.

The file reader maps the JPEG frame’s full encoded range to 0–255 with nearest integer rounding before grayscale conversion or RGB image construction. This display policy applies equally to low-precision lossless and 12-bit DCT input. DICOM decoding uses the clinical path instead: it preserves full integer sample values and requires the frame precision to match BitsStored.

Example Summary

JPEG-LS Native Codec

RITK implements grayscale JPEG-LS in Rust. The encoder and decoder do not call CharLS, GDCM, or another C/C++ codec through FFI. This keeps malformed marker handling, dimension validation, entropy state, and reconstructed samples inside Rust’s checked ownership boundary.

JPEG-LS is a predictive codec designed for low complexity. It is distinct from the block-DCT JPEG family and from wavelet-based JPEG 2000. The normative baseline is ITU-T T.87, equivalent to ISO/IEC 14495-1. DICOM PS3.5 section 8.2.3 encapsulates JPEG-LS streams and assigns transfer syntax UIDs:

• 1.2.840.10008.1.2.4.80 for lossless coding;
• 1.2.840.10008.1.2.4.81 for near-lossless coding.

For MONOCHROME1 and MONOCHROME2, DICOM permits one sample per pixel, 8- or 16-bit allocation, and 2–16 stored bits. RITK’s decoder accepts that 2–16-bit precision range. The current public encoder accepts unsigned single-component 8–16-bit samples.

Prediction and local gradients

For the sample at row y, column x, JPEG-LS uses four causal neighbors:

c  b  d       previous reconstructed row
a  x          current reconstructed row

The edge-detecting predictor is

          min(a, b),        if c >= max(a, b)
Px =      max(a, b),        if c <= min(a, b)
          a + b - c,        otherwise

The first row and column use the boundary rules defined by T.87 Annex A. Three local gradients,

g1 = d - b
g2 = b - c
g3 = c - a

are quantized into a context. Flat neighborhoods enter run mode; all other neighborhoods enter regular mode.

Regular mode

Regular mode predicts x, computes the prediction error, maps its signed value to a nonnegative integer, and writes a Golomb-Rice code. Context statistics adapt the Rice parameter and a small bias correction as samples are processed. The decoder performs the inverse mapping and reconstructs the same causal sample before moving to the next column.

For lossless coding, the quantized error equals the exact prediction error:

NEAR = 0  =>  reconstructed[x] = source[x]

For near-lossless coding, the error is quantized in steps of 2 × NEAR + 1. The reconstruction contract is the exact per-sample bound

abs(reconstructed[x] - source[x]) <= NEAR

This is not an average quality target. Every decoded sample must satisfy it.

Run mode

When all three quantized gradients are zero, the neighborhood is locally flat. Run mode emits adaptive run-length segments using the standard J table. When a run ends before the row boundary, one run-interruption sample is coded against the appropriate neighbor context.

Run mode explains why a smooth medical image can compress well without a transform: long predictable regions require few coded bits, while boundaries return to regular mode.

Encoding through the public API

use ritk_codecs::jpeg_ls::encoder::encode_grayscale_jpeg_ls;

let pixels = [0, 128, 512, 4095];
let lossless = encode_grayscale_jpeg_ls(&pixels, 2, 2, 12, 0)?;
let near_lossless = encode_grayscale_jpeg_ls(&pixels, 2, 2, 12, 3)?;

The encoder validates all external metadata before allocating entropy state:

• rows and columns are nonzero and fit the 16-bit SOF55 fields;
• the checked product equals pixels.len();
• precision is in the implemented 8–16-bit encoder range;
• every sample is no greater than 2^precision - 1; and
• NEAR fits the SOS field and is no greater than MAXVAL / 2.

Invalid input returns JpegLsEncodeError. RITK does not clamp geometry, truncate NEAR, or emit a partial stream.

Bounded reconstruction memory

Prediction reads only the previous row and the reconstructed prefix of the current row. RITK therefore shares one ReconstructionRows implementation between encoder and decoder:

previous row: cols × i32
current row:  cols × i32
total:        2 × cols × 4 bytes

The previous implementation retained (rows + 1) × cols values. At 512 × 512, reconstruction scratch decreases from 1,050,624 bytes to 4,096 bytes. The new bound is O(cols) and does not grow with image height. The decoded output still requires one value per pixel because the public API returns the complete image.

DICOM boundary

DICOM image metadata and the JPEG-LS interchange header must agree on rows, columns, sample count, and precision. RITK validates that geometry before scan allocation. One marker pass checks every declared segment bound and returns the exact entropy slice, avoiding a second scan through compressed bytes. It supports the single-component, non-interleaved profile used by the grayscale native path. Unsupported color/interleaved streams, mapping tables, restart intervals, non-preset LSE records, and nonzero point transforms return errors; they do not fall back to a foreign codec.

Verification

RITK checks the codec at complementary levels:

1. exact lossless round trips over uniform, gradient, edge, run-interruption, 8-, 12-, and 16-bit fixtures;
2. randomized lossless equality and near-lossless error-bound properties;
3. complete DICOM encapsulation and decode round trips;
4. malformed geometry, precision, NEAR, marker lengths, interleave, mapping table, and restart-interval rejection, plus bounded arbitrary marker bytes;
5. a structural assertion that the reconstruction workspace is exactly two rows; and
6. unchanged 512 × 512 Criterion encode/decode workloads.

Continue with the worked JPEG-LS example to see why the source and near-lossless anatomy panels can look almost identical and how the magnified error panel exposes the actual changes.

Example: Lossless and Near-Lossless JPEG-LS

This example creates one deterministic 96 × 96, 12-bit medical-style phantom. It encodes and decodes the same image twice through RITK’s public native API: lossless NEAR = 0 and near-lossless NEAR = 3.

The first three panels deliberately use the same [0, 4095] display range. That makes their anatomy directly comparable, but an error of at most three gray levels is difficult to see across a 12-bit range. The fourth panel solves that visibility problem: it displays abs(near_lossless - source) on [0, 3]. It is an error map, not another contrast-stretched anatomy image.

The metrics below the panels are computed from the arrays that are rendered:

• the lossless mismatch count is an exact equality oracle and must be zero;
• each encoded size is the actual returned stream length;
• changed-sample count reports how many near-lossless values differ; and
• maximum absolute error must be no greater than the declared NEAR = 3.

Source and command

Source: crates/ritk-codecs/examples/book_jpeg_ls.rs

cargo run -p ritk-codecs --example book_jpeg_ls -- \
  docs/book/figures/jpeg_ls_codec.svg

The example fails instead of writing a misleading figure if lossless reconstruction differs at any sample, near-lossless error exceeds NEAR, the near-lossless path produces a degenerate all-zero error map, malformed input does not return PixelCountMismatch, or any rendered array disagrees with the declared geometry.

Adapt the workflow

let stream = encode_grayscale_jpeg_ls(
    &samples,
    rows,
    columns,
    bits_stored,
    near,
)?;

let reconstructed = decode_jpeg_ls_fragment(
    &stream,
    PixelLayout {
        rows: usize::try_from(rows)?,
        cols: usize::try_from(columns)?,
        samples_per_pixel: 1,
        bits_allocated: if bits_stored <= 8 { 8 } else { 16 },
        pixel_representation: PixelSignedness::Unsigned,
        rescale_slope: 1.0,
        rescale_intercept: 0.0,
    },
)?;

Use near = 0 when every stored integer must be preserved. A nonzero value permits a bounded per-sample change; whether that is clinically acceptable is a dataset and workflow policy outside the codec.

JPEG 2000 Native Codec

RITK implements the JPEG 2000 codestream path in Rust. It does not call OpenJPEG, openjp2, or another C codec through FFI. This matters at a medical image boundary: malformed dimensions, precision, or packet lengths remain Rust errors instead of crossing an unsafe foreign interface.

Decoder contract

The native decoder currently accepts a complete grayscale codestream with one component, LRCP progression, no progression changes, no multiple-component transform, inline packet headers, one tile-part per tile, and unit component sampling. Packet coding uses default precincts, 64 × 64 nominal code-blocks, default code-block style, no SOP/EPH markers, and at most 32-bit component precision, matching the decoder’s signed coefficient representation. These constraints match the encoder and the grayscale DICOM workflow shown in this chapter. The decoder preflights main and tile headers plus complete tile coverage before allocating output. Color or other multi-component streams, chroma subsampling, progression overrides, packed packet headers, multi-part tiles, non-default packet coding, and MCT are typed errors. Returning an error is necessary here: replaying the same packet cursor independently for each component can produce plausible-looking but duplicated channels.

The codestream boundary is exact. Every marker segment must include a valid length and remain inside its tile-part or codestream boundary, SOD is found by structural marker parsing rather than a byte-pattern search, every tile declared by SIZ must appear, and EOC must terminate the stream. Every expected LRCP packet header must also be present. Psot = 0 extends only the final tile-part to the validated terminal EOC; marker-looking bytes inside a length-delimited COM payload are never treated as boundaries. Arbitrary bytes after EOC are rejected. The one permitted exception is a single zero byte when an odd-length codestream must be padded to an even DICOM Fragment Item length, as required by DICOM PS3.5 Section 8.2.

An included code-block must declare a nonzero packet-body length. Tier-1 also checks the pass count against the code-block’s bit-plane budget and tracks MQ terminal fill consumption. JPEG 2000 arithmetic termination permits bounded look-ahead into an artificial 0xFF 0xFF marker; consuming more than the two terminal reads used by OpenJPEG’s predictable-termination check is treated as truncation. A decoded first tile followed by a truncated marker tail or entropy body therefore returns an error instead of a partially populated image whose missing coefficients or voxels appear as valid zeros.

JPEG 2000 is a transform codec, not the older block-DCT JPEG format. The Part 1 pipeline is:

integer samples
    │
    ├─ unsigned DC level shift
    ▼
discrete wavelet transform (5/3 or 9/7)
    ▼
quantization (none for reversible 5/3)
    ▼
EBCOT code-block coding and tier-2 packets
    ▼
SOC | SIZ | COD | QCD | tile-part | EOC

The normative codestream and transform definitions are in ITU-T T.800 (11/2015), identical to ISO/IEC 15444-1, especially Annexes A, B, D, E, F, and G. DICOM embeds a raw J2K codestream in encapsulated Pixel Data rather than a JP2 file container. DICOM PS3.5 Section 8.2.4 defines that encoding boundary, while Section A.4.4 maps the reversible path to transfer syntax UID 1.2.840.10008.1.2.4.90 and the reversible or irreversible path to .91.

Reversible 5/3 transform

The reversible transform uses integer lifting. For an interleaved one- dimensional signal, the forward predict and update steps are

d[n] = x[2n+1] - floor((x[2n] + x[2n+2]) / 2)
s[n] = x[2n]   + floor((d[n-1] + d[n] + 2) / 4)

Whole-sample symmetric extension supplies samples at image boundaries. RITK applies the one-dimensional transform by columns and rows at each resolution, placing LL, HL, LH, and HH bands in Mallat layout. The inverse lifting steps are exact integer inverses. With no quantization, decoded samples therefore equal the source samples bit-for-bit.

Irreversible 9/7 transform

The irreversible path uses the floating-point 9/7 lifting transform and dead-zone scalar quantization. QuantizationStep controls the coefficient interval width Δ: larger steps map more nearby coefficients to the same integer index, usually reducing packet bytes while increasing reconstruction error. JPEG 2000 stores Δ as a five-bit exponent and eleven-bit mantissa in QCD. RITK rounds the requested positive finite step to the nearest representable QCD value for each subband, then uses that exact represented value for both coefficient quantization and metadata.

This is scalar quality control, not target-rate control. The same Δ can produce different byte lengths for different anatomy, noise, precision, and wavelet depth. Treat codestream size and PSNR as measured outputs. Target-byte or target-bitrate encoding requires multiple coding-pass truncation points and rate-distortion optimization, which this one-layer encoder does not claim.

The worked example shows unit and coarse (Δ = 32) reconstructions side by side, prints each actual byte count and PSNR, and renders the coarse absolute error with its own magnified black→red→yellow scale. That separate palette and scale make changes visible even where the shared 12-bit anatomy display range hides them.

Encoding through the public API

use ritk_codecs::jpeg_2000::encoder::{
    encode_grayscale_j2k, Jpeg2000Encoding,
};
use ritk_codecs::PixelSignedness;

let pixels = [0, 128, 512, 4095];
let codestream = encode_grayscale_j2k(
    &pixels,
    2,
    2,
    12,
    PixelSignedness::Unsigned,
    Jpeg2000Encoding::Lossless {
        decomposition_levels: 1,
    },
)?;

The encoder validates before constructing the transform buffer:

• rows and columns are nonzero and their product fits usize;
• pixels.len() equals rows × columns;
• precision is in the current native encoder range, 1–16 bits;
• decomposition depth does not exceed the geometry’s meaningful resolution count; and
• every sample fits the declared signed or unsigned range.

Lossy construction validates the scalar step separately:

use ritk_codecs::jpeg_2000::encoder::{Jpeg2000Encoding, QuantizationStep};

let encoding = Jpeg2000Encoding::Lossy {
    decomposition_levels: 3,
    quantization_step: QuantizationStep::new(8.0)?,
};

Zero, negative, NaN, and infinite steps fail at QuantizationStep::new. A finite step whose QCD exponent is unavailable for any requested subband fails before wavelet allocation.

These are errors, not panic conditions. The current i32 entropy path intentionally rejects DICOM’s wider legal precision range rather than truncating values or risking lifting overflow. Color components, chroma subsampling, multiple tiles, custom precincts, target-rate control, and JP2 containers are also outside the current encoder contract.

Memory and packet structure

Unsigned JPEG 2000 components are centered by subtracting 2^(precision-1). RITK applies that offset while constructing the Mallat transform buffer. It does not first allocate a second, complete shifted image. For a reversible image of N samples, this removes 4N bytes of peak full-image storage and one pass that writes those bytes. Code-block coefficient vectors remain bounded by the 64 × 64 code-block geometry. Irreversible encode similarly retains one f32 Mallat volume and quantizes one code block at a time; it no longer materializes a second complete i32 coefficient volume.

The encoder emits one image-wide tile, one quality layer, LRCP progression, one precinct per resolution, and 64 × 64 nominal code-blocks. Before output allocation, the decoder validates SIZ/SOT geometry, tile-header marker extents, single-part tile coverage, and EOC. Packet decode then requires each LRCP packet header and bounds its declared body length before consuming packet data.

Interoperability evidence

RITK tests the native paths at three levels:

1. analytical round trips require zero error for reversible 5/3;
2. quantization tests require each requested step to round to one QCD pair and require packet coefficients to use the reconstructed QCD value;
3. captured OpenJPEG 2.5.4 codestreams exercise independent producer/consumer interoperability without loading OpenJPEG at test time; and
4. malformed geometry, unsupported component/progression overrides, truncated markers, zero-length or exhausted entropy bodies, incomplete packets, multi-part declarations, missing tiles, invalid Psot = 0 boundaries, trailing payload, and missing EOC require typed rejection instead of a panic or partial image.

Continue with the worked codec example to read the reconstruction and error panels.

Example: JPEG 2000 Quality and Size

This example creates one deterministic 96 × 96, 12-bit medical-style phantom. It encodes and decodes that same image three times through RITK’s public native API: three-level reversible 5/3, three-level irreversible 9/7 with Δ = 1, and the same irreversible path with Δ = 32.

All four anatomy panels use the same [0, 4095] display range. The reversible panel must match the source both visually and numerically. Comparing Δ = 1 with Δ = 32 isolates scalar quantization: transform depth and source stay fixed. The fifth panel is not another contrast-stretched anatomy image; it is abs(coarse - source), displayed on [0, max error] with a black→red→yellow heat map. Its independent scale, distinct palette, and printed maximum error expose where coarse quantization removes detail.

The metrics below the panels are computed from the same arrays that are rendered:

• the 5/3 mismatch count is an exact equality oracle and must be zero;
• each codestream size is the actual returned byte length;
• both 9/7 maximum errors are measured per sample;
• PSNR uses the declared 12-bit peak value, 4095; and
• the executable fails if Δ = 32 does not reduce bytes and increase error relative to Δ = 1 on this fixed phantom.

Source and command

Source: crates/ritk-codecs/examples/book_jpeg_2000.rs

cargo run -p ritk-codecs --example book_jpeg_2000 -- \
  docs/book/figures/jpeg_2000_codec.svg

The example fails rather than writing a misleading figure if reversible reconstruction differs at any sample, either lossy error panel is degenerate, the coarse size/error tradeoff is reversed, or rendered arrays do not match the declared geometry.

Adapt the workflow

For an unsigned grayscale image:

use ritk_codecs::jpeg_2000::encoder::{encode_grayscale_j2k, Jpeg2000Encoding};
use ritk_codecs::{decode_jpeg2000_fragment, PixelLayout, PixelSignedness};

let codestream = encode_grayscale_j2k(
    &samples,
    rows,
    columns,
    bits_stored,
    PixelSignedness::Unsigned,
    Jpeg2000Encoding::Lossless {
        decomposition_levels,
    },
)?;

let reconstructed = decode_jpeg2000_fragment(
    &codestream,
    PixelLayout {
        rows: usize::try_from(rows)?,
        cols: usize::try_from(columns)?,
        samples_per_pixel: 1,
        bits_allocated: u16::try_from(bits_stored)?,
        pixel_representation: PixelSignedness::Unsigned,
        rescale_slope: 1.0,
        rescale_intercept: 0.0,
    },
)?;

Use Jpeg2000Encoding::Lossless for the lossless DICOM transfer syntax. For an irreversible stream, construct QuantizationStep::new(delta) and pass it in Jpeg2000Encoding::Lossy. Increase delta only after measuring reconstruction error on the intended image class. The value controls scalar quantization, not a promised file size or bitrate.

VTK Format Boundary

The VTK boundary supports scalar image handoff and surface-oriented output. For scalar image paths, use the same native facade as the other lossless formats:

let image = ritk_io::read_image_native("volume.vtk")?;
ritk_io::write_image_native("volume_copy.vtk", &image)?;

Surface output uses explicit mesh and polydata writer APIs exported by ritk-io, including write_mesh_as_vtk and the OBJ, PLY, STL, VTP, and glTF helpers. A surface writer does not infer image spacing from a raw point list; construct the mesh in the intended physical frame before exporting it.

Direction-aware volume handoff

ritk-snap::LoadedVolume uses [depth, row, column] axes and keeps its direction cosine matrix beside the decoded samples. The VTK boundary converts that contract to a zero-copy [ritk_vtk::VtkImageVolume] in VTK order [x, y, z] = [column, row, depth]:

let volume = ritk_vtk::VtkImageVolume::try_from(&loaded_volume)?;
assert_eq!(volume.dimensions(), [columns, rows, depth]);
let samples = volume.scalars();

The carrier preserves origin, spacing, direction, channel count and x-fastest interleaved samples while sharing the source allocation. Call to_vtk_image_data only at an existing serializer or filter boundary that requires VTK’s owned Vec<f32> attribute arrays; that operation is explicit and retains the direction matrix on the carrier.

Example Summary

Descriptive Statistics and Histograms

An image is both a spatial field and a population of intensity samples. Descriptive statistics answer population questions before an algorithm uses spatial relationships:

• What range of values exists?
• Where is the center of the distribution?
• How broad is the distribution?
• Does a foreground mask select a materially different population?

These questions are useful for quality control, choosing display windows, normalization, and checking whether a mask contains the anatomy it claims to contain. They do not explain where a value occurs; position extrema, label statistics, and spatial filters answer different questions.

The population contract

For finite samples (x_1,\ldots,x_N), RITK computes

[ \mu = \frac{1}{N}\sum_{i=1}^{N} x_i ]

and

[ \sigma = \sqrt{\frac{\sum_{i=1}^{N}(x_i-\mu)^2}{N-\mathrm{ddof}}}. ]

Use ddof = 0 for population standard deviation and ddof = 1 for sample standard deviation. This is the same divisor convention documented by NumPy std. RITK rejects ddof >= N; returning zero, infinity, or NaN would hide that no valid divisor remains.

Mean and variance use two f64 accumulation passes even though image samples are f32. The wider accumulator prevents a long CT-scale sum from losing individual voxel contributions once the running total’s f32 spacing exceeds their magnitude. Results remain f32 at the public image boundary.

Quartile convention

RITK reports three discrete order statistics:

[ Q_1 = x_{(\lfloor N/4\rfloor)},\quad Q_2 = x_{(\lfloor N/2\rfloor)},\quad Q_3 = x_{(\lfloor 3N/4\rfloor)}. ]

Indices are zero-based after ordering the samples. These are observed sample values, not interpolated estimates. NumPy exposes multiple percentile methods; its default linear method can therefore differ between adjacent samples. RITK’s floor-rank convention is stable and explicit.

The implementation does not fully sort the image. It copies a borrowed slice once, then isolates the three ranks in place with progressively smaller quickselect suffixes. The average complexity is (O(N)), and the workspace is one f32 buffer of length (N). Masked statistics collect foreground values once and reuse that allocation as the percentile workspace.

Full image and masked foreground

use ritk_statistics::{
    compute_statistics_from_slice,
    masked_statistics_from_slices,
};

let image = [12.0, 15.0, 82.0, 91.0, 154.0, 160.0];
let mask  = [ 0.0,  0.0,  1.0,  1.0,   1.0,   1.0];

let full = compute_statistics_from_slice(&image, 0)?;
let foreground = masked_statistics_from_slices(&image, &mask, 0)?;

assert!(foreground.mean > full.mean);
Ok::<(), ritk_statistics::StatisticsError>(())

A mask value greater than 0.5 selects a foreground sample. Image and mask buffers must have identical lengths. An empty foreground is an error; RITK does not substitute full-image statistics because that changes the requested population.

Histograms

A histogram divides an explicit finite range [min, max] into equal-width bins. Every bin except the last is half-open. The last includes max:

[ B_k = [e_k,e_{k+1}),\qquad B_{K-1} = [e_{K-1},e_K]. ]

This matches the edge convention documented by NumPy histogram. Finite values outside the range are excluded. NaN and infinite values are rejected rather than discarded or assigned to bin zero. Histogram::bin_width() returns f64, matching the wide bin-coordinate arithmetic even when subtracting two finite f32 bounds would overflow in f32.

use ritk_statistics::histogram_from_slice;

let values = [0.0, 1.0, 2.0, 3.0, 4.0];
let histogram = histogram_from_slice(&values, 0.0, 4.0, 4)?;

assert_eq!(histogram.counts, vec![1, 1, 1, 2]);
assert_eq!(histogram.total(), values.len());
Ok::<(), ritk_statistics::StatisticsError>(())

Histogram allocation is fallible. A hostile or accidental bins value that cannot reserve its count buffer returns HistogramAllocationFailed instead of relying on an unchecked allocation.

Failure handling

StatisticsError distinguishes failures that require different caller responses:

Python bindings map these invalid values to ValueError, the Python category for an argument of the correct type whose value violates the operation’s contract.

The worked example renders a source field, foreground mask, normalized histogram comparison, and numeric cross-checks from the same computed arrays.

Example: Full-image and Masked Distributions

This example creates a deterministic intensity field with dark background, two foreground tissue groups, and a bright lesion. It then computes full-image and masked-foreground statistics through RITK’s public slice APIs.

The panels answer three separate questions:

1. Intensity field — what values are present and where?
2. Foreground mask — exactly which samples define the masked population?
3. Distribution comparison — how do the selected populations differ on one intensity axis?

Blue bars are the share of all image samples in each bin. Orange bars are the share of masked samples in the same bins. Each population is normalized by its own sample count, so bar heights compare distribution shape rather than mask size. Solid lines mark means, dotted lines mark medians, and the narrow bands at the base mark interquartile ranges.

The numeric table makes the change measurable. The full-image median remains in the dark background peak, while the masked median lies in foreground tissue. This is why an empty-mask fallback would be misleading: full-image and masked statistics are not interchangeable views of one population.

Source and command

Source: crates/ritk-statistics/examples/book_descriptive_statistics.rs

cargo run -p ritk-statistics --example book_descriptive_statistics -- \
  docs/book/figures/descriptive_statistics.svg

The example fails before writing the figure if:

• RITK’s min, max, mean, standard deviation, or floor-rank quartiles disagree with an independent full-sort reference;
• either histogram total differs from its source population;
• the image or mask violates the finite-input contract; or
• the phantom does not produce a large enough mean and median separation to make the comparison visually meaningful.

Adapt the workflow

extern crate ritk_statistics;

use ritk_statistics::{
    compute_statistics_from_slice, histogram_from_slice, masked_statistics_from_slices,
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let image_values = [0.0_f32, 40.0, 80.0, 120.0];
    let mask_values = [0.0_f32, 1.0, 1.0, 0.0];

    let full = compute_statistics_from_slice(&image_values, 0)?;
    let masked = masked_statistics_from_slices(&image_values, &mask_values, 0)?;
    let full_histogram = histogram_from_slice(&image_values, 0.0, 120.0, 4)?;

    assert_eq!(full.min, 0.0);
    assert_eq!(full.max, 120.0);
    assert_eq!(masked.min, 40.0);
    assert_eq!(masked.max, 80.0);
    assert_eq!(full_histogram.total(), image_values.len());
    Ok(())
}

Use identical histogram edges when comparing populations. Different ranges or bin counts can make the same data appear different and invalidate a visual comparison.

Intensity Transformations

Pointwise intensity remapping filters: windowing, rescaling, thresholding, sigmoid, arithmetic operations, and histogram equalization.

Design

All intensity filters operate on flat host buffers via substrate-agnostic pure functions. Each filter follows the extract → compute → reconstruct sequence through ritk-image::Image boundary. No Coeus tensor is constructed in the hot path.

Filter Families

• Windowing: IntensityWindowingFilter — maps input range to output range
• Rescaling: RescaleIntensityFilter — affine remap with saturation clipping
• Thresholding: ThresholdImageFilter — binary threshold (inside/outside)
• Sigmoid: SigmoidImageFilter — logistic intensity remap
• Arithmetic: Add, Subtract, Multiply, Divide — unary/binary ops
• Equalization: HistogramEqualizationFilter, AdaptiveHistogramEqualizationFilter
• Clamp/Shift-Scale: ClampImageFilter, ShiftScaleImageFilter

The windowing and rescaling example uses the RIRE CT fixture to show why a tissue-specific HU window is not equivalent to a global output rescale. Its histogram and labeled panels are the reference visual for this chapter.

Spatial Intensity Non-Uniformity

N4BiasFieldCorrectionFilter estimates a smooth multiplicative MRI bias field in log-intensity space with histogram sharpening and multi-resolution cubic B-spline fitting. It is a spatial correction rather than a pointwise remap; the N4 book example shows the source slice, the corrected slice, the relative change map, and the estimated field on the RIRE MR fixture.

Verification

Each filter is differentially tested against its Coeus-generic counterpart via assert_coeus_matches_coeus. Analytical oracles verify endpoint mapping and saturation behavior.

Example: Windowing and Rescaling

CT voxels carry Hounsfield units (HU), but the useful display range depends on the tissue being inspected. Windowing is a saturating affine map:

[ W(I) = \operatorname{clamp}(I, a, b), \qquad O(I) = \frac{W(I)-a}{b-a}(o_{max}-o_{min}) + o_{min}. ]

Values below a map to o_min, values above b map to o_max, and values inside the window retain their relative order. Rescaling is different: it uses the global image minimum and maximum as a and b, so it changes the representation range without selecting a tissue-specific interval.

The runnable example applies both native filters to the real RIRE Patient 001 CT volume. The generated figure is deliberately labeled and uses one fixed display convention per panel:

1. Input CT uses a fixed [-1000, 1000] HU display window.
2. Soft-tissue window applies [-160, 240] HU and emits [0, 1].
3. Lung window applies [-1000, 400] HU and emits [0, 1].
4. Global rescale uses the observed CT extrema and emits [0, 255].
5. Input distribution plots the same axial slice and marks both HU windows.
6. Filter contract states the saturation and output-range behavior next to the image panels.

The histogram prevents a common interpretation error: the soft-tissue and lung panels are not alternative contrast adjustments applied after rendering. They are different numeric maps of the same CT slice, and the marked intervals show which input HU values are expanded or saturated.

Source and command

Source: crates/ritk-io/examples/book_windowing_rescale.rs

cargo run -p ritk-io --example book_windowing_rescale -- \
  docs/book/figures/windowing_rescale.svg

The source image, filtered outputs, and histogram all come from the in-tree RIRE fixture. The example uses the Coeus-native path and does not write a second copy of the image through a separate display library.

Verification

• Each filter preserves the CT shape and physical metadata through the native image boundary.
• Windowed values stay inside [0, 1]; the global rescale stays inside [0, 255].
• The figure is regenerated from the actual source and filter outputs.
• Analytical endpoint and saturation behavior is covered by the filter’s native tests; the figure checks the real-data visual behavior separately.

Example: A Complete Processing Pipeline

This example is the reference workflow for the filter chapters. It creates one deterministic three-dimensional scalar image, applies several real RITK native filters, and renders every result with an explicit display range.

Read the panels left-to-right, then continue on the next row:

1. the smooth scalar input with local perturbations;
2. a bounded sigmoid remap;
3. threshold suppression;
4. physical gradient magnitude;
5. binary opening;
6. binary closing;
7. grayscale opening;
8. Perona–Malik diffusion;
9. the absolute diffusion change;
10. curvature flow;
11. the absolute curvature change; and
12. the contract summary.

The input is a smooth phantom with two deterministic local perturbations, making the diffusion change visible as spatial structures rather than a directional stripe or fine-grain noise pattern. It also contains a small isolated foreground island and a small background hole, replicated through the three slices. With the radius-one cubic structuring element used below, opening removes the island while closing fills the hole, so the two binary panels are intentionally not interchangeable. The input, sigmoid, grayscale, and diffusion panels share [0, 1]. The binary panels use the same range because their values are exactly 0 or 1. Gradient and change panels use their own data-derived upper bounds, which are printed under the panel. The phantom has three identical depth slices and the figure shows the center slice, so radius-one 3-D morphology and diffusion operate on a genuine volume rather than a depth-one edge case. This avoids the common documentation error where every output is independently contrast-stretched and therefore appears unchanged. The binary panel subtitles also report the number of foreground voxels removed by opening and background voxels filled by closing.

Source and command

Source: crates/ritk-filter/examples/book_processing_pipeline.rs

cargo run -p ritk-filter --example book_processing_pipeline -- \
  docs/book/figures/processing_pipeline.svg

The example uses Image::from_flat_on with [depth, row, column] shape and unit physical spacing. Every filter receives the same image boundary and the native SequentialBackend; no plotting library supplies the processed values.

Reuse the individual stages

The core calls are:

let remapped = SigmoidImageFilter::new(0.42, 0.10, 0.0, 1.0)
    .apply_native(&input, &backend)?;
let mask = BinaryThresholdImageFilter::new(0.62, 1.0, 1.0, 0.0)
    .apply_native(&remapped, &backend)?;
let edges = GradientMagnitudeFilter::unit().apply_native(&input)?;
let denoised = DiffusionConfig {
    num_iterations: 12,
    time_step: 0.0625,
    conductance: 0.08,
    ..DiffusionConfig::default()
}
.apply_native(&input, &backend)?;

Binary opening and closing are explicit compositions in this example: erode → dilate and dilate → erode. That makes the structuring-element semantics visible and lets a caller substitute another radius or foreground value without changing the image contract.

Verification

• The example checks image construction and every native filter result through Result; no stage is replaced by a display-only approximation.
• All panels are generated from the same source and use fixed, labeled scales.
• The change panel is computed from the input and diffusion arrays, so it exposes whether the filter actually changed the data.
• Filter-specific algebraic and native-path tests remain the authoritative correctness checks; the SVG is the visual check for structure and readability.

Example: Thresholding

Thresholding changes voxel values but not the image frame. RITK exposes a retaining threshold and a binary indicator:

• ThresholdImageFilter::outside(lower, upper, outside) keeps values in the closed interval and replaces values outside it.
• BinaryThresholdImageFilter::new(lower, upper, foreground, background) writes the foreground value inside the interval and background elsewhere.

The complete pipeline figure shows both behavior classes. The threshold panel is the suppression stage; the binary morphology panels consume a binary indicator made from the same sigmoid output.

let retained = ThresholdImageFilter::outside(0.0, 0.58, 0.0)
    .apply_native(&sigmoid, &backend)?;
let mask = BinaryThresholdImageFilter::new(0.62, 1.0, 1.0, 0.0)
    .apply_native(&sigmoid, &backend)?;
assert_eq!(retained.shape(), sigmoid.shape());
assert_eq!(mask.shape(), sigmoid.shape());

The interval is inclusive for the binary filter. That boundary matters for label maps and should be tested with values exactly at both endpoints. Values outside the interval are not clipped into the interval; they become the configured background.

Source and verification

Source: crates/ritk-filter/examples/book_processing_pipeline.rs

cargo run -p ritk-filter --example book_processing_pipeline -- \
  docs/book/figures/processing_pipeline.svg

The native filter tests cover endpoint inclusion, outside suppression, shape preservation, and native/generic parity. The figure verifies that the foreground topology is visible and that thresholding is not an independently contrast-stretched display artifact.

Example: Sigmoid and Arithmetic

The sigmoid maps an input intensity through

S(I) = (M - m) / (1 + exp(-(I - alpha) / beta)) + m.

Alpha is the midpoint and beta controls the transition width. RITK keeps the output range explicit, which makes a preprocessing contract reproducible across displays and downstream metrics.

let normalized = SigmoidImageFilter::new(0.42, 0.10, 0.0, 1.0)
    .apply_native(&input, &backend)?;

Pointwise arithmetic filters use the same image boundary. Add, subtract, multiply, and divide are binary image operations: both operands must have compatible shapes and metadata. Use them to combine a corrected image with a mask or to apply a calibrated scale after the nonlinear remap. Do not use a display-only conversion as a substitute for a numeric arithmetic stage.

The pipeline figure makes the sigmoid effect readable by keeping the input and sigmoid panels on the same [0, 1] display range. The output is visibly different because the remap is data-derived, not because the renderer chooses a new contrast window.

Source and verification

Source: crates/ritk-filter/examples/book_processing_pipeline.rs

cargo run -p ritk-filter --example book_processing_pipeline -- \
  docs/book/figures/processing_pipeline.svg

The sigmoid tests cover midpoint, monotonicity, bounded output, and the degenerate beta step contract. Binary arithmetic tests cover shape mismatch and value-semantic results.

Example: N4 Bias-Field Correction

MRI intensity is not spatially uniform even when the underlying tissue is. The N4 model writes the measured image as

[ I(x) = S(x) B(x), \qquad \log I(x) = \log S(x) + \log B(x), ]

where (S) is the tissue signal and (B) is a slowly varying multiplicative bias field. RITK estimates the log-bias with histogram sharpening and a multi-resolution cubic B-spline fit, then returns (\exp(\log I - \widehat b)). The implementation follows Tustison et al., N4ITK: Improved N3 Bias Correction, IEEE TMI 29(6), 1310–1320 (2010), which is also cited by the filter’s Rustdoc.

The figure uses the real RIRE Patient 001 T1 volume:

The four panels use a two-by-two layout:

1. the source MR slice;
2. the native N4-corrected slice;
3. the estimated field’s percentage deviation from its foreground median, using blue for lower and red for higher multiplicative bias; and
4. the estimated multiplicative field (I / I_{corrected}).

The first two panels intentionally use one shared robust intensity window formed from both slices’ 2nd–98th percentile ranges. Their labels repeat that window so the comparison cannot be mistaken for independently normalized images. Because N4 correction is spatially smooth and can remain subtle in anatomical views, the third panel removes the field’s global median scale and amplifies its spatial deviation. It is masked outside the bright foreground and clipped only to its 98th-percentile absolute deviation for display; the underlying correction is not altered.

The third and fourth panels are computed from the same input and corrected output. They are diagnostics of the correction estimated by N4, not hand-authored illustrations.

Source and command

Source: crates/ritk-registration/examples/book_n4_bias.rs

cargo run -p ritk-registration --example book_n4_bias -- \
  docs/book/figures/n4_bias_correction.svg

The example reads test_data/registration/rire/training_001_mr_T1.mha, runs the Coeus-native N4BiasFieldCorrectionFilter, and preserves the volume shape and physical metadata. It uses two fitting levels and twelve iterations per level so the documentation example remains bounded while exercising the real multi-resolution path. The full filter defaults remain available through N4Config::default() for production workloads.

Verification

• The source and corrected arrays have identical voxel counts and geometry.
• The corrected output contains finite positive intensities.
• The figure is regenerated by the runnable Rust example from the RIRE input.
• Filter correctness is covered by the package’s analytical and native-path N4 tests; this figure checks the real-data visual behavior separately.

Spatial Filtering

Convolution-based spatial filters: Gaussian smoothing, gradient magnitude, Canny edge detection, and separable gradient computation.

Design

All spatial filters operate on flat z-major host buffers via substrate-agnostic pure functions. The gradient magnitude and Canny detectors share a common host core that implements gradient computation, non-maximum suppression, and hysteresis. Public apply_native methods preserve the ritk-image::Image geometry boundary while executing those kernels without constructing a second tensor representation.

Gaussian Smoothing

GaussianFilter builds one normalized sampled kernel per axis from physical sigma and voxel spacing, then applies separable zero-padded convolution. The native path and the generic path share the kernel and host-core contracts.

Gradient Magnitude

Computes the magnitude of the gradient via central differences, then optionally smooths with a discrete Gaussian kernel.

Canny Edge Detection

Four stages form the binary edge map:

1. Gaussian smoothing.
2. Gradient magnitude and direction.
3. Non-maximum suppression along the continuous gradient direction.
4. Hysteresis thresholding.

The complete filter gallery uses the real public pipeline on a deterministic phantom and writes the figure below:

Run it from the repository root with:

cargo run -p ritk-filter --example book_filter_gallery -- \
  docs/book/figures/filter_gallery.svg

Verification

The filter example is source-linked from Gaussian Smoothing and Canny Edge Detection. Package tests provide the value-semantic and differential coverage; the figure is a visual smoke check of the same public calls.

Example: Gaussian Smoothing

This example is the Gaussian stage of the shared filter gallery. It constructs a deterministic scalar phantom, wraps it in ritk_image::Image<f32, _, 3>, and applies GaussianFilter::apply_native with a physical sigma of 2.0.

The figure compares the input with the smoothed image and the Canny edge map that consumes the same image boundary:

Source and command

Source: crates/ritk-filter/examples/book_filter_gallery.rs

cargo run -p ritk-filter --example book_filter_gallery -- \
  docs/book/figures/filter_gallery.svg

The image shape and spatial metadata are unchanged by the filter. The package tests cover native/generic parity; the committed SVG is the visual output of this example’s figure contract. Rerun the command after changing the example to regenerate the artifact from the computed values.

Example: Gradient Magnitude

For a scalar image I on a regular grid, RITK estimates physical derivatives with central differences and reports

|grad I| = sqrt((dI/dz)^2 + (dI/dy)^2 + (dI/dx)^2).

Spacing is part of the calculation. A two-voxel intensity change represents a different physical slope when the voxel spacing changes.

let gradient = GradientMagnitudeFilter::new(*input.spacing())
    .apply_native(&input)?;

GradientMagnitudeFilter::unit() is useful only when the image has unit spacing. For physical measurements, pass the image spacing explicitly. The native implementation uses a zero-flux boundary stencil and preserves the input geometry.

Source and verification

Source: crates/ritk-filter/examples/book_processing_pipeline.rs

cargo run -p ritk-filter --example book_processing_pipeline -- \
  docs/book/figures/processing_pipeline.svg

The figure uses one data-derived upper bound for the gradient panel and prints that bound below the image. This prevents a gradient map from looking identical to the input merely because both were independently normalized. Filter tests cover constant fields, linear fields with known spacing, boundary behavior, and native/generic parity.

Example: Canny Edge Detection

The filter gallery runs the complete Canny pipeline on the same deterministic phantom as the Gaussian example. CannyEdgeDetector::apply_native performs Gaussian pre-smoothing, central-difference gradients, continuous-direction non-maximum suppression, and 26-connected hysteresis.

Source and command

Source: crates/ritk-filter/examples/book_filter_gallery.rs

cargo run -p ritk-filter --example book_filter_gallery -- \
  docs/book/figures/filter_gallery.svg

The output is binary (1.0 at retained edges, 0.0 elsewhere) and preserves the input geometry. Threshold behavior is covered by the package’s value semantic tests; the SVG makes the connected edge structure inspectable.

Morphological Operations

Binary and grayscale morphological filters: erosion, dilation, opening, closing, top-hat transforms, hit-or-miss, and morphological reconstruction.

Design

All morphological operations operate on flat host buffers via substrate-agnostic pure functions. Binary operations use safe-border padding (zero padding at boundaries); grayscale operations preserve spatial metadata.

Binary Morphology

• BinaryErodeFilter, BinaryDilateFilter — flat rectangular structuring elements with one voxel radius per axis
• BinaryFillholeFilter — flood-fill based hole filling
• BinaryMorphologicalClosing, BinaryMorphologicalOpening — dilation followed by erosion (closing) or erosion followed by dilation (opening)

For anisotropic volumes, a voxel radius does not represent one physical distance. voxel_radii_for_physical_radius converts a millimetre radius to independent [z, y, x] radii using floor(radius / spacing_i). Pass the result to with_axis_radii; for 3.0 mm slices and 0.4 mm in-plane spacing, a 3.2 mm support becomes [1, 8, 8] rather than [8, 8, 8]. The latter would span 24 mm through-plane and can corrupt registration masks even though its axial view looks plausible.

Grayscale Morphology

• GrayscaleErosion, GrayscaleDilation — min/max filter with safe-border padding
• GrayscaleClosingFilter, GrayscaleOpeningFilter — extensive and anti-extensive operations, respectively
• WhiteTopHatFilter, BlackTopHatFilter — residue of opening/closing
• GrayscaleGradientFilter — erosion minus dilation

Reconstruction

• GrayscaleGeodesicErosion, GrayscaleGeodesicDilation — geodesic operations using reconstruction

Verification

Each operation is differentially tested against its Coeus-generic counterpart via assert_coeus_matches_coeus. Analytical oracles verify:

• Radius 0 is identity
• All-foreground/all-background edge cases
• Opening is anti-extensive; closing is extensive away from the finite-image boundary where the documented zero-background padding applies

The binary opening and closing example constructs a mask with both removable foreground artifacts and fillable background holes, then renders explicit red removal and green fill maps from the native filter outputs.

Example: Binary Opening and Closing

Binary opening and closing are dual compositions with visibly different purposes:

• open_B(A) = (A ⊖ B) ⊕ B removes structures that cannot contain the structuring element;
• close_B(A) = (A ⊕ B) ⊖ B fills gaps that cannot contain the reflected structuring element.

Opening is anti-extensive: it removes foreground structures that cannot contain the structuring element. Closing is extensive: it fills background gaps that cannot contain the reflected element.

The runnable example uses one deterministic 3-D mask containing both defect classes. The isolated speck and one-voxel spur give opening something to remove; the small internal holes give closing something to fill. The red and green change maps make those distinct effects explicit:

Source and command

Source: crates/ritk-filter/examples/book_binary_morphology.rs

cargo run -p ritk-filter --example book_binary_morphology -- \
  docs/book/figures/binary_morphology.svg

The example executes BinaryMorphologicalOpening::apply_native and BinaryMorphologicalClosing::apply_native with radius one. It fails unless:

• opening preserves geometry and every displayed interior output voxel is less than or equal to its input voxel;
• closing preserves geometry and every displayed interior output voxel is greater than or equal to its input voxel;
• opening removes at least one foreground voxel;
• closing fills at least one background voxel; and
• the opening and closing outputs differ.

The contract is checked on the center slice because RITK’s documented zero-background boundary condition can remove foreground at the outermost volume planes during the erosion stage of closing.

The complete processing pipeline shows the same operations in a longer filter chain. Tests additionally cover radius-zero identity, border behavior, foreground values, and native/generic parity.

Example: Grayscale Morphology

Grayscale opening is a local minimum followed by a local maximum. It removes bright protrusions smaller than the structuring element while retaining the image geometry. Grayscale closing performs the dual maximum-then-minimum operation and fills small dark holes.

let opened = GrayscaleOpeningFilter::new(2)
    .apply_native(&sigmoid, &backend)?;
let closed = GrayscaleClosingFilter::new(2)
    .apply_native(&sigmoid, &backend)?;

The radius is measured in voxels, not physical units. Choose it from the feature size you intend to remove and account for anisotropic spacing before using the filter on a clinical volume. Both operations use replicate padding at the boundary, so their behavior differs from binary erosion’s background-outside policy.

Source and verification

Source: crates/ritk-filter/examples/book_processing_pipeline.rs

cargo run -p ritk-filter --example book_processing_pipeline -- \
  docs/book/figures/processing_pipeline.svg

The figure renders the opening stage on [0, 1]; the closing call above follows the same image contract and can be inspected with the same display helper. The opening result therefore shows local intensity removal directly. Tests cover radius zero identity, monotonicity, idempotence, safe-border behavior, and native/generic parity.

Diffusion Filtering

PDE-based image smoothing: Perona-Malik, curvature flow, curvature diffusion, and gradient anisotropic diffusion.

Design

All diffusion filters operate on flat host buffers via substrate-agnostic pure functions. Each filter implements a host core shared by the Coeus-generic apply path and the native apply_native path. No Coeus tensor is constructed in the hot path.

Perona-Malik

Edge-preserving smoothing via the Perona-Malik PDE:

dI/dt = div(g(|∇I|) ∇I)

where g is the edge-stopping function. Two variants: g1 (exponential) and g2 (inverse quadratic).

Curvature Flow

dI/dt = ∇·(∇I/|∇I|) |∇I|

Smooths along edges while preserving across-edge structure.

Curvature Diffusion

∂I/∂t = |∇I| ∇·(c(|∇I|) ∇I/|∇I|)

RITK implements the Modified Curvature Diffusion Equation: each intensity level set evolves by conductance-modulated mean curvature. Unlike heat-equation diffusion, the gradient-normalized flux and edge-stopping conductance preserve strong boundaries while smoothing along their geometry.

Gradient Anisotropic Diffusion

dI/dt = ∇·(g(|∇I|) ∇I)

where g is a function of gradient magnitude. Smoothing is anisotropic: strong along edges, weak across edges.

Verification

Each filter is differentially tested against its Coeus-generic counterpart via assert_coeus_matches_coeus.

Example: Perona–Malik Diffusion

Perona–Malik diffusion evolves

dI/dt = div(c(|grad I|) grad I),

where the conductance c is large in homogeneous regions and small across strong gradients. RITK provides exponential and quadratic conductance strategies and an explicit Euler configuration.

let config = DiffusionConfig {
    num_iterations: 12,
    time_step: 0.0625,
    conductance: 0.08,
    function: ConductanceFunction::Exponential,
};
let diffused = config.apply_native(&input, &backend)?;

For a three-dimensional unit grid, the explicit time step must satisfy the stability bound dt <= 1/6. The example uses 1/16. Conductance is expressed in the same intensity-gradient units as the image; it is not a display contrast parameter.

The figure uses a shared [0, 1] scale for input and output and adds an absolute change panel. That diagnostic is required when visual contrast alone cannot show whether a denoiser changed the data.

Source and verification

Source: crates/ritk-filter/examples/book_processing_pipeline.rs

cargo run -p ritk-filter --example book_processing_pipeline -- \
  docs/book/figures/processing_pipeline.svg

Tests cover constant-field invariance, conductance selection, stability-sized steps, metadata preservation, and native/generic parity.

Example: Curvature Flow

Curvature flow evolves a scalar field according to mean curvature. Unlike ordinary isotropic diffusion, the update depends on local surface geometry and is intended to smooth small-scale boundary roughness while retaining larger features for longer.

let filtered = CurvatureFlowImageFilter::new(CurvatureFlowConfig {
    num_iterations: 5,
    time_step: 0.0625,
})
.apply_native(&input, &backend)?;

The explicit Euler stability bound for a three-dimensional unit grid is dt <= 1/6. RITK uses double-precision stencil arithmetic internally for the curvature numerator and writes the result back to the image scalar type. Spacing is read from the image metadata.

The complete pipeline source runs five steps and renders the result beside an absolute input-to-output change panel. The change panel is not a second contrast-adjusted image; it is a direct difference diagnostic.

Source and verification

Source: crates/ritk-filter/examples/book_processing_pipeline.rs

cargo run -p ritk-filter --example book_processing_pipeline -- \
  docs/book/figures/processing_pipeline.svg

Curvature-flow tests cover constant images, stability-sized steps, boundary stencils, iteration behavior, and native/generic parity.

Seeded Segmentation

Segmentation assigns a discrete anatomical or material label to each voxel. RITK provides automatic thresholding, connected components, region growing, clustering, watershed, level-set evolution, and ensemble consensus methods in ritk-segmentation. This chapter starts with GrowCut because it makes the interaction between image evidence and user-provided constraints explicit.

GrowCut state and update

GrowCut is a synchronous cellular automaton. Voxel (p) carries an intensity (I_p), label (L_p), and confidence (C_p\in[0,1]). Seed voxels start with their supplied label and unit confidence; unseeded voxels start at label zero and confidence zero. For each face-connected neighbor (q), RITK computes

[ g(p,q)=1-\frac{|I_p-I_q|}{I_{\max}-I_{\min}},\qquad A(q,p)=C_q,g(p,q). ]

If (A(q,p)>C_p), the neighbor wins the local competition and both state variables update:

[ (L_p,C_p)\leftarrow(L_q,A(q,p)). ]

All voxels read the previous iteration and write the next one, so traversal order does not change the result. The process stops only when the complete ((L,C)) state is stable or the configured iteration limit is reached. A same-label confidence increase matters: it can strengthen a later attack at a class boundary even though the current label did not change.

The worked GrowCut figure makes this update visible over time. Orange and cyan are competing labels, while the underlying light/dark image is the evidence controlling attack strength. Pixels that retain that grayscale image are still unlabeled. The figure evaluates (g) on each side of its known boundary so the final circle follows from the update rule rather than from the appearance of the last panel alone.

This follows the state and transition rule in Vezhnevets and Konouchine, “GrowCut — Interactive Multi-Label N-D Image Segmentation by Cellular Automata,” §2.1. The paper establishes algorithm structure; RITK’s executable tests and the analytical example establish this implementation’s behavior.

Public native API

use coeus_core::SequentialBackend;
use ritk_segmentation::GrowCutFilter;

let backend = SequentialBackend;
let labels = GrowCutFilter::new(200)
    .apply_native(&image, &seeds, &backend)?;

image and seeds must have identical three-dimensional geometry. Positive integer-valued seed voxels are immutable labels; zero denotes an unlabeled voxel. apply_native preserves the input shape, origin, spacing, and direction and reports shape or storage failures through Result.

Choosing seeds

• Place at least one seed inside each class that the result must contain.
• Keep seeds away from uncertain boundaries.
• Add corrective seeds when different tissues have overlapping intensities.
• Treat the iteration limit as a resource bound, not as a convergence test.

The two-level phantom in the worked GrowCut example has an exact solution. Real CT and MR volumes do not: intensity overlap, partial-volume voxels, artifacts, and seed placement all affect clinical validity. Validate real segmentations against independent labels and report overlap and boundary metrics.

Example: GrowCut from Sparse Seeds

The runnable example constructs a deterministic two-tissue image with a circular target, places a (3\times3) background seed and a (3\times3) target seed, and executes GrowCutFilter::apply_native.

Read the four numbered panels from left to right:

• orange marks the background label and cyan marks the circular target label;
• the first panel overlays only the two (3\times3) seed regions on the grayscale input;
• the next two panels are actual GrowCutFilter outputs after 8 and 40 synchronous sweeps, not illustrative drawings;
• the last panel is the converged RITK output. It reports foreground Dice and the exact label-error count against the analytical circle.

At each sweep, labeled voxels try to transfer their label to face-connected neighbors. A transfer is strong when the two intensities are similar. In this phantom, equal-intensity neighbors have (g=1), so both colored fronts spread without losing confidence. At the circular boundary the intensity difference equals the complete image range, giving (g=0); neither label can attack across it. The calculation under the panels uses these exact phantom values.

Pixels that retain the grayscale input in the two middle panels remain undecided at that iteration. They are not a third output class.

Source and command

Source: crates/ritk-segmentation/examples/book_growcut.rs

cargo run -p ritk-segmentation --example book_growcut -- \
  docs/book/figures/growcut.svg

The example fails unless GrowCut:

• preserves shape, origin, spacing, and direction;
• labels every voxel exactly as the analytical circle specifies;
• returns zero false labels; and
• produces foreground Dice (=1).

The truth image remains an independent oracle generated from the circle equation, even though it is no longer a standalone panel. The analytical phantom verifies label propagation, competition at a high-contrast boundary, and spatial metadata. It does not establish clinical accuracy on heterogeneous anatomy.

Registration Metrics

Metric selection follows the intensity relationship, not the optimizer.

The differentiable path exposes native MSE, NCC, LNCC, and NGF evaluators. The classical path complements them with mutual information, packed MIND-SSC, and sealed translation metrics. Both paths consume the RITK image geometry contract; the classical path performs an explicit conversion to Leto arrays at its numeric boundary.

Packed MIND-SSC

MIND describes anatomy by comparing nearby patches within each image instead of comparing CT and MR intensities directly. RITK’s SSC pattern uses the 12 cross-axis pairs among six axial neighbours. A 3×3×3 patch produces 12 summed squared patch distances; local mean distance normalizes contrast, and the minimum-distance subtraction makes the largest response one. Five unary bits encode six levels (0–5) and pack all 12 components into the low 60 bits of one u64. XOR/popcount is therefore the exact L1 distance between quantized descriptors.

MindSscFixedPrep stores packed descriptors only for a deterministic fixed set of complete-support centers. The default cap is 8,192. It follows the Hoeffding sample requirement of 6,623 for a bounded [0,1] mean loss with ±0.02 error at 99% confidence under uniform random sampling, rounded up to the next power of two. RITK fixes the hash seed and stratifies in physical space for reproducible coverage. That statistical design explains the cap; one fixed realization does not constitute population or clinical validation.

During each pose evaluation, RITK maps every required fixed patch point into moving physical space and trilinearly samples it. Support inside ITK’s half-voxel field uses replicate-edge interpolation; support outside is zero background and remains in the fixed denominator. This preserves a pose-invariant objective while keeping persistent memory proportional to selected centers, not image voxels.

Heinrich et al. define MIND in section 3 of the 2012 paper (https://doi.org/10.1016/j.media.2012.05.008) and the compact 12-component SSC pattern in section 3.1 of the 2013 paper (https://doi.org/10.1007/978-3-642-40811-3_24).

Evaluate the identity transform before optimization. A post-registration metric is meaningful only when the identity value, transform convention, sampled frame, and optimizer state are recorded beside it.

Example Summary

Optimization and Registration

Metrics only become registration algorithms once ritk couples them to an optimizer. On the differentiable side, ritk-registration::metric::autodiff supplies a reusable gradient_descent driver configured by GradientDescentConfig, rebuilding a transform from trainable Coeus Vars each iteration and stepping parameters with the native SGD helper. On the classical side, ImageRegistration uses ClassicalConfig and MutualInformationMetric to run deterministic CPU optimization over rigid or affine parameters. This chapter covers that seam: parameterization, iteration budgets, tolerances, step sizes, and why optimizer behavior must be read together with the chosen similarity metric.

Atlas integration is split but coherent. Coeus provides the autodiff graph, tensor execution, and backend flexibility for learning-style registration loops, while Leto supports the classical numeric path where predictable CPU array behavior is preferred. RITK keeps those implementation details behind a common image boundary so callers can reason about transforms, convergence, and output geometry without rewriting file I/O or preprocessing around each optimizer family.

Bounded rigid capture and refinement

search_rigid_pose searches a six-degree-of-freedom rigid transform in physical millimetres. It starts from a validated RigidSearchAnchor, performs four coarse-to-fine coordinate-descent levels, and polishes the capture objective with bounded Nelder–Mead. A second structural objective defaults to a half-range of one terminal capture cell. with_structural_half_range_cells accepts a NonZeroU8 when a consumer needs to test a wider local basin. The configured radius scales both the structural bounds and initial simplex, but the original global rigid bounds remain authoritative. The effective interval is their intersection. When capture ends on a global boundary, each simplex edge points toward the side with available room rather than constructing an invalid outward vertex. This separation supports multimodal registration where normalized mutual information (NMI) finds the broad basin and MIND-SSC or Normalized Gradient Fields (NGF) resolves local soft-tissue structure without enabling an unbounded second search. The tested Rustdoc on RigidSearchConfig::with_structural_half_range_cells is the copyable API example.

The capture schedule and all simplex operations preserve finite proposals and their bounded-objective rejection semantics. If an otherwise valid finite resolution multiplied by a schedule factor would overflow, only that non-finite proposal is replaced by the finite endpoint in its direction before an objective sees it. A separate finite-transform check converts overflow from centroid and residual composition into a typed numerical failure before metric evaluation.

RigidSearchResult::capture_saturated and structural_saturated report when an optimum touches its permitted boundary. Structural saturation covers either the configured local half-range or a tighter global bound. Saturation is a quality-control signal: it means the configured search region may have clipped the optimum, not that the transform is invalid. The caller still owns image sampling, fixed-domain support, overlap acceptance, and the final choice between the capture and structural candidates.

Objective errors propagate through the search result. A malformed or empty metric sample set therefore cannot be converted into a plausible pose.

When centroids alone do not place the orientation inside this bounded search, fit_symmetric_trimmed_rigid estimates the complete anchor from direction- specific FixedToMovingCorrespondence and MovingToFixedCorrespondence physical-space matches. It fits each direction independently, retains the half with smallest Euclidean residual under a deterministic LTS fit, inverts the reverse result, and combines the two fixed-to-moving transforms with a log-Euclidean mean. Exact endpoint pairs with conflicting directions are discarded because selecting either would make the directional schedules depend on image-role ordering. The fit requires a strict majority supporting one non-collinear rigid consensus; exactly two equal half-population consensuses are not identifiable. See the robust rigid capture initializer for the executable known-transform and outlier case.

Example Summary

Diffusion Gradient Schemes

ritk-diffusion-scheme is RITK’s single-source-of-truth for validated diffusion-MRI acquisition metadata. A diffusion series is meaningful only when each volume carries a physically typed weighting, a unit gradient direction, and a declared coordinate frame. This crate owns that format-neutral contract; format crates convert DICOM, NRRD, MRtrix, or companion-file metadata at their trust boundaries.

Why a separate crate?

Diffusion model fitting consumes three inputs — a 4-D image, a gradient scheme, and (for constrained methods) a response function. If the scheme is just two f64 arrays, each consumer re-implements the same validation: non-finite checks, zero/unit-vector contracts, b-value thresholding, shell grouping, frame tracking, and reorientation. A single validated type eliminates this duplication and prevents the single most common defect class in this domain: a rigid-motion correction that does not rotate the gradient table.

Diffusion weighting

DiffusionWeighting stores the diffusion sensitization factor with dimension time per area. Internal storage uses canonical SI seconds per square meter through the Aequitas quantity system. The scanner-facing constructor and accessor use the MRI convention seconds per square millimeter:

use ritk_diffusion_scheme::DiffusionWeighting;

let b1000 = DiffusionWeighting::from_seconds_per_square_millimeter(1_000.0)?;
assert_eq!(b1000.seconds_per_square_millimeter(), 1_000.0);

Construction rejects negative, NaN, or infinite values, and catches overflow when converting to canonical SI. The is_unweighted method checks whether the value is exactly zero in the scanned representation.

Gradient directions and frames

A GradientDirection pairs a weighting with a validated direction vector. The invariant depends on the weighting:

• An unweighted (b = 0) entry requires an exact zero vector.
• A weighted entry requires a finite unit vector within 1e-6 Euclidean norm.

Both are enforced at construction:

use ritk_diffusion_scheme::{GradientDirection, DiffusionWeighting};
use ritk_spatial::Vector;

let b0_weighting = DiffusionWeighting::from_seconds_per_square_millimeter(0.0)?;
let b1000_weighting = DiffusionWeighting::from_seconds_per_square_millimeter(1_000.0)?;

// b0 entry: zero vector.
let b0 = GradientDirection::new(b0_weighting, Vector::new([0.0, 0.0, 0.0]))?;

// Weighted entry: unit vector.
let dir = GradientDirection::new(b1000_weighting, Vector::new([1.0, 0.0, 0.0]))?;

// This fails — weighted volume with non-unit vector.
assert!(GradientDirection::new(
    b1000_weighting,
    Vector::new([0.5, 0.5, 0.5]),
).is_err());

GradientFrame is the coordinate frame of every direction in a scheme:

A direction is numerically valid only within its declared frame. Applying an ImageAxis direction in LPS or vice versa produces a silently wrong tensor field. The frame is pinned at construction and persistent: every format codec states which frame it produces, and the reader converts once.

The GradientScheme type

GradientScheme is an ordered collection of validated gradient directions with one declared frame. Construction is infallible once all entries are validated:

use ritk_diffusion_scheme::{GradientScheme, GradientFrame};
use ritk_spatial::Vector;

let scheme = GradientScheme::from_seconds_per_square_millimeter(
    vec![
        (0.0, Vector::new([0.0, 0.0, 0.0])),
        (1_000.0, Vector::new([1.0, 0.0, 0.0])),
        (1_000.0, Vector::new([0.0, 1.0, 0.0])),
        (2_000.0, Vector::new([0.0, 0.0, 1.0])),
    ],
    GradientFrame::ImageAxis,
)?;

assert_eq!(scheme.len(), 4);

B0 thresholding

from_seconds_per_square_millimeter applies a default 50 s/mm² threshold. Scanner baseline acquisitions at or below this value are canonicalized to exact zero weighting and direction, because their gradient orientation is either absent or not physically meaningful. Values above the threshold must carry a unit direction.

Query

• directions() — all entries in acquisition order.
• frame() — the declared GradientFrame.
• b0_indices(threshold) and dwi_indices(threshold) — split the scheme at a given weighting threshold.
• shells() — unique nonzero shell weightings sorted in ascending order. Multi-shell acquisitions (e.g. b = 1 000, 2 000, 3 000) produce the corresponding list.

Gradient reorientation

When a rigid or affine correction is applied to a diffusion series, the gradient directions must rotate with the image. A motion correction that does not rotate the gradient table produces a silently wrong tensor field — this is the single most common defect class in diffusion MRI.

reorient applies one proper orthonormal rotation to every weighted direction. b = 0 volumes pass through unchanged because their placeholder vector carries no physical meaning:

// 90° rotation about the z-axis.
let rotation = [
    [0.0, -1.0, 0.0],
    [1.0,  0.0, 0.0],
    [0.0,  0.0, 1.0],
];
let reoriented = scheme.reorient(rotation)?;

reorient_per_volume is the motion-correction path: each volume gets its own rotation matrix from the per-volume registration. It requires exactly one rotation per volume and rejects any matrix that is non-finite, not orthonormal within 1e-9, or has a determinant other than +1.

Format codecs

Every codec states which frame it produces and converts once. Any transform applied to the series reorients the directions with it.

FSL .bval / .bvec

FSL stores gradient metadata as companion text files. The .bval file is a single line of space-separated s/mm² values. The .bvec file is three lines — one per spatial component (x, y, z) — of space-separated unit direction components.

Directions are in image-axis coordinates (GradientFrame::ImageAxis).

use ritk_diffusion_scheme::fsl::{read_fsl_scheme, write_fsl_scheme};

let bval = "0 1000 1000";
let bvec = "0 1 0.5\n0 0 0.866\n0 0 0";
let scheme = read_fsl_scheme(bval, bvec)?;

let (bval_out, bvec_out) = write_fsl_scheme(&scheme);

parse_fsl_bval and parse_fsl_bvec provide standalone parsing when the caller already has validated b-values or directions separately.

MRtrix .mif DW_scheme

MRtrix embeds the gradient table in the .mif container header as a DW_scheme key. The value is a matrix with four columns per volume: gradient x, y, z (unit direction, or [0,0,0] for b = 0) and b-value in s/mm².

Directions are in image-axis coordinates (GradientFrame::ImageAxis).

use ritk_diffusion_scheme::mrtrix::{read_mrtrix_scheme, write_mrtrix_scheme};

let header_text = "\
DW_scheme: 3,4
0,0,0,0
1,0,0,1000
0.5,0.866,0,1000
";
let scheme = read_mrtrix_scheme(header_text)?;
let dw_block = write_mrtrix_scheme(&scheme);

The reader validates the four-column structure, rejects missing or malformed dimension declarations, and checks that every component is parseable and finite. The writer emits the DW_scheme: N,4 header block ready for inclusion in a .mif file.

NA-MIC DWI NRRD

The NA-MIC DWI convention stores gradient metadata as NRRD key-value pairs. One nominal DWMRI_b-value is combined with per-volume DWMRI_gradient_NNNN direction keys whose four-digit index starts at 0000.

Directions are read through the NRRD measurement frame and returned in RITK physical LPS coordinates (GradientFrame::Lps).

use ritk_nrrd::read_nrrd_gradient_scheme;

let scheme = read_nrrd_gradient_scheme("dwi.nrrd")?;

The reader validates consecutive index ordering, rejects non-contiguous or out-of-range indices, and handles DWMRI_NEX and B-matrix encodings by failing explicitly rather than guessing.

DICOM

DICOM diffusion metadata lives in two standard top-level tags:

• Diffusion b-value (0018,9087) — s/mm²
• Diffusion Gradient Orientation (0018,9089) — unit direction in the patient coordinate system (GradientFrame::Lps)

use ritk_dicom::diffusion::read_dicom_gradient_scheme_from_files;

let scheme = read_dicom_gradient_scheme_from_files(&dicom_paths)?;

When standard tags are absent, the reader attempts extraction from vendor private blocks. Siemens CSA headers ((0029,1020) series or (0029,1010) image level) carry B_value and DiffusionGradientDirection in a binary SV10 table. GE and Philips private groups are recognised as present but return Ok(None) until a vendor-specific decoder is added.

The reader does not infer volume grouping from a directory or guess enhanced functional groups; those require a separate sequence-aware parser. For an unweighted frame, a finite zero b-value with no orientation is mapped to the required zero vector.

Cross-codec verification

ADR 0036 verification condition 8 requires that every codec round-trip recovers the same typed scheme, and that all four codecs agree on one dataset expressed in all four conventions. The cross-codec test in crates/ritk-diffusion-scheme/tests/cross_codec.rs writes and reads a single-shell scheme through all four codecs and asserts identity:

FSL round-trip    → scheme_fsl
MRtrix round-trip → scheme_mrtrix
NRRD round-trip   → scheme_nrrd
DICOM round-trip  → scheme_dicom

assert scheme_fsl == scheme_mrtrix == scheme_nrrd == scheme_dicom

A multi-dataset variant repeats the check against the ds002087 (single-shell b = 1 000) and ds004666 (multi-shell) OpenNeuro datasets when their test data has been downloaded.

Error types

GradientSchemeError is a #[non_exhaustive] enum covering every validation failure with an acquisition index where applicable:

Every error carries enough context to identify the exact volume or metadata field that failed, supporting actionable diagnostics at the format boundary.

Invariant

ritk-diffusion-scheme is the only crate that constructs GradientScheme. Format crates convert raw bytes to validated entries at their boundaries; consumers receive an already-validated scheme. This prevents diffusion models from silently operating on unvalidated or frame-mismatched metadata.

Diffusion Models

ritk-diffusion owns orientation models derived from voxel signals. Every model consumes a validated GradientScheme and one voxel’s signal vector; every model produces a typed result that carries its coordinate frame forward. This chapter documents all five implemented models.

Model-to-solver map

The model type owns only the estimation and the derived metrics. Every solver is dispatched to its owner (Leto for linear algebra, Coeus for nonlinear optimisation, Apollo for the SH basis) — RITK adds no local solver, no local gradient, and no local SH implementation.

Relationship to the gradient scheme

Every model in this crate assumes a validated GradientScheme. The Diffusion Gradient Schemes chapter establishes the guarantees this chapter depends on:

1. Every direction is validated (b = 0 → zero vector, weighted → unit vector).
2. Every direction has one declared frame (ImageAxis or Lps).
3. Gradient reorientation is a scheme-level operation — it must be called before the scheme enters any model fitter.
4. Shell grouping is provided by the scheme (b0_indices, dwi_indices, shells).

No model in this chapter re-validates direction norms or re-derives shell structure. The scheme chapter owns format-codec concerns; this chapter owns the transformation from validated scheme + signals → fitted field.

DTI — Diffusion Tensor Imaging

DTI is the baseline diffusion model. It assumes a 3-D Gaussian displacement distribution characterised by a symmetric 3×3 diffusion tensor \(D\). The normalised signal at gradient direction \(\mathbf g\) and b-value \(b\) is

\[ \frac{S}{S_0} = \exp\!\left(-b\,\mathbf g^{\!T} D\,\mathbf g\right). \]

Taking the log gives a linear system in the six unique tensor elements:

\[ \ln(S/S_0) = -b \cdot [g_x^2, g_y^2, g_z^2, 2g_x g_y, 2g_x g_z, 2g_y g_z] \cdot \mathbf d \]

where \(\mathbf d = [D_{xx}, D_{yy}, D_{zz}, D_{xy}, D_{xz}, D_{yz}]^{\!T}\). RITK assembles the design matrix over all weighted acquisitions and solves via leto_ops::solve_least_squares. Fractional anisotropy, mean diffusivity, and the principal eigenvector are derived from the fitted tensor through a closed-form 3×3 symmetric eigendecomposition, per ADR 0036 decision 2.

Configuration

DtiConfig carries one parameter: the b0_threshold separating reference and weighted volumes. The default is 50 s/mm², the same threshold used by GradientScheme::from_seconds_per_square_millimeter.

use ritk_diffusion::dti::{estimate_dti, DtiConfig};

let config = DtiConfig::default();
let tensor = estimate_dti(&scheme, &signals, config)?;

Output

DiffusionTensor provides:

Error conditions

DtiError covers signal/scheme length mismatch, non-finite signals, missing b0 or weighted volumes, fewer than six DWI directions (Underdetermined), invalid baseline or normalised signals, a failed linear solve, and non-positive eigenvalues in the fitted tensor.

DKI — Diffusion Kurtosis Imaging

DKI extends DTI with the kurtosis tensor \(W\), a fourth-order symmetric tensor that captures non-Gaussian diffusion. The normalised signal is

\[ \frac{S}{S_0}= \exp\!\left(-b\,D(\mathbf g)

• \frac{b^2}{6}\,\mathrm{MD}^2\,W(\mathbf g)\right) \]

where \(D(\mathbf g)=\mathbf g^{\!T}D\mathbf g\) is the apparent diffusion coefficient, \(\mathrm{MD} = \operatorname{tr}(D)/3\), and \(W(\mathbf g)\) is the full contraction of the kurtosis tensor with the gradient direction.

The model has 21 parameters (6 for D, 15 for W). Estimation proceeds in two stages:

1. A log-linear DTI fit supplies the initial D guess and baseline S₀.
2. Levenberg-Marquardt refines D and W simultaneously, using the analytic Jacobian of the DKI residual.

The solver requires at least 21 DWI directions distributed across multiple shells — a single-shell scheme cannot disambiguate the quadratic D term from the quartic W term.

Configuration

KtiConfig carries a b0_threshold and a LevenbergMarquardtConfig. The default tolerances are appropriate for synthetic data; in-vivo fitting typically requires relaxed convergence criteria.

use ritk_diffusion::dki::{estimate_dki, KtiConfig};

let config = KtiConfig::default();
let kt = estimate_dki(&scheme, &signals, &config)?;

Kurtosis metrics

The apparent kurtosis coefficient is \(K(\mathbf g) = (\mathrm{MD}/D(\mathbf g))^2 \cdot W(\mathbf g)\).

Output

DiffusionKurtosisTensor provides elements_d() (six D elements), elements_w() (fifteen W elements), eigenvalues(), fa(), md(), mk(), ak(), rk(), principal_eigenvector(), baseline_signal(), residual_norm(), converged(), iterations(), gradient_norm(), frame(), predict_signal(), quadratic_form(), and kurtosis_at_direction().

Error conditions

KtiError covers signal/scheme length mismatch, non-finite signals, missing b0/weighted volumes, invalid baseline or normalised signals, a failed DTI initial fit, and Levenberg-Marquardt solver failures.

ODF — Analytical Q-ball

Q-ball imaging estimates an orientation distribution function (ODF) without assuming a single diffusion tensor. RITK implements the analytical construction of Descoteaux et al. (2007): the normalised signal is expanded in Apollo’s real, even-degree spherical harmonic basis, and the Funk–Radon transform converts signal coefficients to ODF coefficients.

Signal model

\[ E(\mathbf g)=\sum_{l=0,2,4,\dots}^{l_{\max}} \sum_{m=-l}^{l} c_{lm}Y_{lm}(\mathbf g) \]

Only even degrees appear because diffusion is antipodally symmetric. The design matrix is evaluated at the scheme’s scattered gradient directions by apollo_sht::RealSphericalHarmonicBasis::design_matrix.

Regularised least squares

The coefficients are estimated by solving

\[ \hat{\mathbf c} = (B^{\!T} B + \lambda L)^{-1} B^{\!T} \mathbf e \]

where \(\lambda \ge 0\) is the Laplace–Beltrami penalty weight and \(L\) has diagonal entries \(l(l+1)\) for each degree-\(l\) coefficient. The solve routes through leto_ops::solve_least_squares.

Funk–Radon transform

The fitted signal coefficients become Q-ball ODF coefficients through

\[ \psi_{lm}=2\pi P_l(0)c_{lm} \]

where \(P_l(0)\) is the Legendre polynomial of degree \(l\) at zero. This produces an ODF, not a fibre orientation density — the two must not be interpreted interchangeably.

Single-shell constraint

Analytical Q-ball fits one q-space shell. OdfConfig carries a shell_tolerance (maximum absolute b-value difference among weighted volumes). A multi-shell acquisition returns MixedShells instead of entering a physically invalid angular fit.

Configuration

use ritk_diffusion::odf::{estimate_odf, OdfConfig};
use ritk_diffusion_scheme::DiffusionWeighting;

let config = OdfConfig::new(
    4,                                    // l_max
    0.006,                                // regularization
    DiffusionWeighting::from_seconds_per_square_millimeter(50.0)?,
    DiffusionWeighting::from_seconds_per_square_millimeter(0.0)?,
)?;
let odf = estimate_odf(&scheme, &signals, config)?;

Output

OdField provides coefficients(), l_max(), baseline_signal(), normalized_signal_residual(), frame(), evaluate(theta, phi), evaluate_at_direction([x, y, z]), and evaluate_on_grid(theta_samples, phi_samples).

Error conditions

OdfError covers signal/scheme length mismatch, non-finite signals, missing b0/weighted volumes, mixed shells, underdetermined basis, invalid baseline, invalid regularisation, invalid spherical grid, invalid evaluation directions, basis errors (propagated from Apollo), and failed solves.

CSD — Constrained Spherical Deconvolution

CSD deconvolves the diffusion signal with an axially symmetric response function to recover the fibre orientation distribution (fODF). In the spherical harmonic basis the convolution becomes a diagonal rescaling of each degree block:

\[ s_{lm} = \sqrt{\frac{4\pi}{2l+1}} \cdot r_l \cdot f_{lm} \]

where \(r_l\) are the rotational harmonics of the single-fibre response. The deconvolution matrix is \(B_{\text{resp}} = B \cdot \operatorname{diag}(\kappa_l)\) with \(\kappa_l = 4\pi/(2l+1) \cdot r_l\). The fODF coefficients are recovered by solving

\[ \min \|B_{\text{resp}} \mathbf f - \mathbf S/S_0\|_2 \quad \text{subject to} \quad \mathbf f \ge 0 \]

through leto_ops::nnls (Lawson–Hanson active-set NNLS). The result is a non-negative fODF whose peaks correspond to fibre directions.

Response function

ResponseFunction carries rotational harmonics \(r_0, r_2, r_4, \dots\) with \(r_0 = 1.0\) (normalised to unit baseline). The convenience constructor from_tensor(b_value, ad, rd, l_max) computes them from an axially symmetric diffusion tensor by numerical projection onto the Legendre polynomials.

Configuration

use ritk_diffusion::csd::{estimate_fod, CsdConfig, ResponseFunction};
use leto_ops::NnlsConfig;

let response = ResponseFunction::from_tensor(3_000.0, 0.0017, 0.0003, 8)?;
let config = CsdConfig::new(8, b0_threshold, NnlsConfig::default())?;
let fod = estimate_fod(&scheme, &signals, &response, &config)?;

Output

FodField provides coefficients() (guaranteed non-negative), l_max(), baseline_signal(), residual_norm(), nnls_iterations(), nnls_converged(), frame(), evaluate(theta, phi), evaluate_at_direction([x, y, z]), evaluate_on_grid(theta, phi), and find_peaks(grid_theta, grid_phi, relative_threshold).

Peak extraction

find_peaks samples the fODF on a dense equiangular grid, locates every sample exceeding all eight neighbours and a configurable relative-amplitude floor, then deduplicates peaks within 5° of each other (to handle pole-adjacent duplicates from the φ-periodic grid). Peaks are returned sorted by descending amplitude.

Volume type

FodVolume stores a 3-D grid of fODF coefficients and supports trilinear interpolation for sub-voxel direction queries during whole-brain tractography via direction_at(point, grid_theta, grid_phi, relative_threshold).

Error conditions

CsdError covers signal/scheme length mismatch, non-finite signals, missing b0/weighted volumes, invalid baseline, underdetermined system, response degree too low for the requested basis, invalid \(r_0\), invalid evaluation/grid parameters, Apollo basis errors, NNLS failures, and volume construction validation (shape, spacing, origin, coefficient count).

Comparison with ODF

NODDI — Neurite Orientation Dispersion and Density Imaging

NODDI (Zhang et al., 2012) is a 3-compartment tissue model that separates water diffusion into intra-neurite (restricted), extra-neurite (hindered), and CSF (free) pools:

\[ \frac{S}{S_0} = (1 - f_{\text{iso}})\left[f_{\text{intra}}A_{\text{ic}}

• (1 - f_{\text{intra}})A_{\text{ec}}\right]
• f_{\text{iso}}A_{\text{iso}} \]

• \(A_{\text{ic}}\) — Watson-averaged stick signal evaluated by Monte Carlo quadrature over 300 quasi-uniform directions on the sphere.
• \(A_{\text{ec}} = \exp(-b \cdot d_{\text{ec}})\) — hindered extra-cellular compartment.
• \(A_{\text{iso}} = \exp(-b \cdot d_{\text{iso}})\) — free-water CSF compartment.

Biophysical constants are fixed: \(d_{\parallel} = 1.7 \times 10^{-3}\) mm²/s, \(d_{\text{ec}} = 0.8 \times 10^{-3}\) mm²/s, \(d_{\text{iso}} = 3.0 \times 10^{-3}\) mm²/s.

Parameters and metrics

ODI = 0 means perfectly aligned sticks (κ → ∞); ODI = 1 means isotropic dispersion (κ → 0). The extra-cellular fraction is computed as \(f_{\text{ec}} = (1 - f_{\text{iso}})(1 - f_{\text{intra}})\).

Solver

Fitting uses Levenberg-Marquardt with five parameters. The Jacobian has two analytic columns (\(f_{\text{intra}}\), \(f_{\text{iso}}\)) and three finite-difference columns (ODI, θ, φ). The initial guess uses DTI’s PEV for the fibre direction and conservative defaults for the volume fractions.

Configuration

use ritk_diffusion::noddi::{estimate_noddi, NoddiConfig};

let config = NoddiConfig::default();
let fit = estimate_noddi(&scheme, &signals, &config)?;

Output

NoddiFit provides ndi() / f_intra(), f_iso(), odi(), f_extra(), principal_direction(), baseline_signal(), residual_norm(), converged(), iterations(), gradient_norm(), frame(), and predict_signal(direction, b_value).

Volume type

NoddiVolume stores a 3-D grid of principal directions and supports nearest-neighbour spatial lookup for tractography via direction_at(point). Unlike FodVolume, no peak extraction is needed — NODDI intrinsically yields a single fibre orientation per voxel.

Error conditions

NoddiError covers signal/scheme length mismatch, non-finite signals, missing b0/weighted volumes, invalid baseline, solver failures (including DTI initial-fit failures), and volume construction validation.

From model fields to streamlines

Every model in this crate produces a type that can feed a direction field for tractography. The direction-field helpers live in ritk-tractography (see the tractography chapter):

The Deterministic Streamline Tractography chapter documents the integration algorithm that consumes these direction fields. For a complete DTI volume, the reusable ritk-tractography surface also selects FA-threshold seeds and applies the volume’s mask and anisotropy floor through DtiTractographyConfig. The signal-to-streamlines example closes the loop end-to-end: known tensor → synthetic signals → model fit → direction field → streamlines.

Diffusion MRI Acquisition and Q-ball ODFs

Diffusion-weighted MRI measures how water self-diffusion varies with gradient direction and strength. RITK separates three concerns:

This chapter focuses on the physics that relates gradient pulses to signal attenuation — the bridge between the validated metadata in the scheme chapter and the model-fitting algorithms documented in the Diffusion Models chapter.

From gradient pulse to b-value

For the idealized Stejskal–Tanner pulsed-gradient experiment, the diffusion weighting is

b = γ² G² δ² (Δ − δ/3)

where γ is the gyromagnetic ratio, G the gradient amplitude, δ the pulse duration, and Δ the pulse separation. The measured signal for a Gaussian diffusion tensor D is

S(b, g) = S₀ exp(−b gᵀ D g)

Thus (b) has dimensions of time per area, conventionally s/mm²; it is not a diffusivity. DiffusionWeighting stores the canonical SI value (s/m²) internally through the Aequitas quantity system and exposes the scanner-facing s/mm² value at the API boundary. See the scheme chapter for the type-level contract and Stejskal and Tanner, “Spin Diffusion Measurements”.

Relationship to the gradient scheme

The Diffusion Gradient Schemes chapter establishes four guarantees that every model fitting depends on: validated directions, unambiguous coordinate frames, scheme-level reorientation, and shell grouping. The Diffusion Models chapter documents how all five implemented models (DTI, DKI, ODF, CSD, NODDI) consume these guarantees — this chapter focuses on the physics shared by all of them.

From scheme to signal to model

The pipeline established by the three chapters in this part is:

1. Validated scheme — the Diffusion Gradient Schemes chapter ensures every direction carries a declared frame, a physically typed weighting, and a validated unit-vector contract.
2. Orientation model — the Diffusion Models chapter documents the transformation from validated scheme + voxel signals → fitted field. Five models are implemented: DTI (log-linear tensor), DKI (nonlinear kurtosis), ODF (analytical Q-ball), CSD (non-negative fODF), and NODDI (3-compartment biophysical).
3. Streamline geometry — the Deterministic Streamline Tractography chapter integrates the fitted orientation field into Gaia polyline curves.

use ritk_diffusion_scheme::{GradientFrame, GradientScheme};
use ritk_spatial::Vector;

let scheme = GradientScheme::from_seconds_per_square_millimeter(
    vec![
        (0.0, Vector::new([0.0, 0.0, 0.0])),
        (1_000.0, Vector::new([1.0, 0.0, 0.0])),
        (1_000.0, Vector::new([0.0, 1.0, 0.0])),
    ],
    GradientFrame::ImageAxis,
)?;
assert_eq!(scheme.len(), 3);
Ok::<(), ritk_diffusion_scheme::GradientSchemeError>(())

Format coordinate contracts

• FSL bvec directions are image-axis coordinates. The three rows are x, y, and z components; the companion bval sequence supplies s/mm² values. DiffusionMaps retains this ImageAxis frame. When a fitted field is placed in DtiVolume, the volume validates that frame and converts the external [column, row, depth] direction order to RITK’s image-index [depth, row, column] order exactly once. LPS maps are rejected at this boundary because physical-to-index conversion requires image geometry that DtiVolume does not own.
• Standard DICOM Diffusion Gradient Orientation (0018,9089) is in the patient coordinate system and RITK labels it LPS. Diffusion b-value (0018,9087) is s/mm². The current reader accepts top-level attributes from classic one-instance-per-volume inputs; it does not guess private vendor fields or enhanced functional groups. See DICOM PS3.3 C.8.13.5.9.
• NA-MIC DWI NRRD stores one nominal DWMRI_b-value. A gradient’s squared norm relative to the largest gradient norm scales that volume’s effective b-value. The measurement frame maps the gradient into world coordinates, after which RAS is converted once to RITK’s LPS physical frame. See the NA-MIC DWI NRRD convention.

These distinctions prevent a numerically valid direction from being applied in the wrong physical frame.

Analytical Q-ball estimation

Q-ball imaging estimates an orientation distribution function (ODF) without assuming a single diffusion tensor. RITK first normalizes weighted signals by the mean b0 signal and fits an even-degree real spherical-harmonic expansion:

E(g) = Σₗ,ₘ cₗₘ Yₗₘ(g)

The weighted samples must occupy one q-space shell. OdfConfig makes the allowed absolute b-value tolerance explicit; an off-shell acquisition returns MixedShells instead of entering a physically invalid angular fit. The fitted OdField retains the scheme frame, so evaluation directions and reported peaks keep an unambiguous coordinate meaning.

The least-squares system is augmented with a Laplace–Beltrami penalty. For coefficient degree l, the diagonal penalty factor is l(l + 1), so the configured nonnegative weight suppresses high-degree oscillation. The Funk–Radon transform then acts degree by degree:

ψₗₘ = 2π Pₗ(0) cₗₘ

This is the analytical Q-ball construction of Descoteaux et al., “Regularized, Fast, and Robust Analytical Q-Ball Imaging”. It produces an ODF, not a constrained-spherical-deconvolution fiber orientation density; those quantities must not be interpreted interchangeably.

OdfConfig validates harmonic degree, regularization, the b0 threshold, and the single-shell tolerance. Estimation fails on non-finite signals, non-finite values created by baseline normalization or the numerical solve, missing b0 or weighted samples, nonpositive baseline, an underdetermined basis, or a failed solve. Evaluation also rejects overflow created by finite coefficients and finite basis values. OdField stores coefficients contiguously and reuses its basis metadata for repeated evaluation. evaluate_on_grid returns one flat row-major allocation rather than a pointer-chasing grid of rows.

The signal-to-streamlines example uses a known tensor to verify that the recovered antipodal ODF peak matches its analytical x axis.

The signal-to-streamlines example closes the loop: known tensor → synthetic signals → model fit → direction field → streamlines.

Creating and Validating Deterministic Tractography

ritk-tractography integrates a local orientation field into curves, producing Gaia polyline geometry. It does not directly reconstruct axons and it does not resolve the biological cause of a local orientation. The current algorithm is a deterministic Euler baseline — explicit stepping with direction continuity, turn-angle gating, and step-count bounding — designed for reproducible examples and baseline algorithms.

This chapter shows how to create a tractogram through the public API, validate it against an analytical field, and export it for comparison with reference toolchains. It also documents the integration algorithm and the direction-field helpers that bridge diffusion models to tractography.

Integration rule

For a physical seed point \(\mathbf x_0\), step size \(h > 0\), and local unit orientation \(\mathbf v(\mathbf x)\), one forward step is

\[ \mathbf x_{k+1} = \mathbf x_k + h \, \mathbf v(\mathbf x_k). \]

At each sample the integrator validates the returned direction: non-finite components or a norm that deviates from unity beyond \(10^{-6}\) produce a typed InvalidDirection error rather than a silently dropped streamline.

Sign continuity

Diffusion orientations are antipodally symmetric: \(\mathbf v\) and \(-\mathbf v\) describe the same physical axis. An ODF peak extractor, spherical harmonic evaluator, or eigendecomposition may return either sign arbitrarily. At each step the integrator computes the dot product between the current direction and the returned direction. When the dot product is negative, the direction is flipped so the streamline continues forward rather than reversing:

if current_direction.dot(&next_direction) < 0.0 {
    next_direction = -next_direction;
}

This prevents arbitrary sign choices from making a streamline reverse direction mid-integration.

Bidirectional tracking

TrackingDirection::Bidirectional integrates both signs from the seed — one half along +v and one along -v — reverses the backward half, and joins the two halves with the seed appearing exactly once. The joined polyline runs from the backward termination point through the seed to the forward termination point.

TrackingDirection::Forward integrates only along the initial direction returned at the seed. The backward termination reason is None in the output.

Configuration

TractographyConfig carries four validated parameters:

Construction validates every parameter and returns InvalidStepSize, InvalidMaxSteps, or InvalidTurnLimit for out-of-range values:

use ritk_tractography::{TractographyConfig, TrackingDirection};

let config = TractographyConfig::new(0.5, 1_000, 45.0, TrackingDirection::Bidirectional)?;

The turn-angle cosine limit is precomputed once at construction: cos(max_turn_degrees.to_radians()). A turn-angle check fails when the dot product between consecutive (sign-aligned) directions falls below this threshold.

Termination is part of the result

Every half-streamline records exactly one TerminationReason:

The distinction between FieldBoundary and TurningAngle matters for interpretation. FieldBoundary means that the caller-defined direction field returned no orientation at the proposal. That boundary may come from an image extent, tissue mask, anisotropy threshold, or model-fit failure; it is not by itself evidence of an anatomical endpoint. TurningAngle means that the next valid orientation exceeded the configured limit and may reflect noise, a crossing-fibre region, or genuine curvature rejected by the configuration.

A seed at which the field returns None at the first query is a normal untrackable seed — it produces no streamline and is not an error.

Direction fields

A direction field is any closure Fn(&Point<3>) -> Option<Vector<3>>. ritk-tractography provides five pre-built helpers that bridge diffusion models to the integration algorithm, and the caller can supply a custom closure for ad-hoc or synthetic fields.

Single-voxel fields

dti_pev_direction_field and fod_peak_direction_field extract one direction from a single-voxel fit and return the same direction at every query point. They are useful for bootstrapping and unit testing:

use ritk_tractography::{euler_tractography, dti_pev_direction_field};

let tensor = ritk_diffusion::dti::estimate_dti(&scheme, &signals, config)?;
let result = euler_tractography(&seeds, track_config, dti_pev_direction_field(&tensor))?;

dti_pev_direction_field returns None when the fitted tensor has near-zero FA (\(< 10^{-10}\)), i.e. it is isotropic and the PEV is not trackable.

fod_peak_direction_field pre-extracts the strongest fODF peak via FodField::find_peaks(50, 100, 0.1) and returns None when no peak meets the 10% relative-amplitude threshold:

use ritk_tractography::fod_peak_direction_field;

let result = euler_tractography(&seeds, track_config, fod_peak_direction_field(&fod))?;

Whole-brain fields

fod_volume_direction_field and noddi_direction_field perform spatial neighbourhood lookups at each integration step — they query a 3-D volume of pre-computed model fits rather than a single voxel.

fod_volume_direction_field trilinearly interpolates fODF coefficients from the surrounding 2×2×2 voxel neighbourhood at each step, then extracts the strongest peak via a 50×100 spherical-grid search with a 10% relative-amplitude threshold:

use ritk_tractography::fod_volume_direction_field;

let result = euler_tractography(&seeds, track_config, fod_volume_direction_field(&volume))?;

noddi_direction_field uses nearest-neighbour spatial lookup — NODDI intrinsically yields a single fibre orientation per voxel, so no peak extraction is needed:

use ritk_tractography::noddi_direction_field;

let result = euler_tractography(&seeds, track_config, noddi_direction_field(&volume))?;

Custom fields

Any closure satisfying Fn(&Point<3>) -> Option<Vector<3>> is a valid direction field. The closure defines the tracking domain: return None outside the domain boundary and a unit Vector<3> inside it:

let result = euler_tractography(&seeds, config, |point| {
    let [x, y, z] = point.to_array();
    // Track only within a 10 mm radius sphere.
    let r_sq = x * x + y * y + z * z;
    if r_sq > 100.0 {
        return None;
    }
    // Constant horizontal field — every trackable point yields +x.
    Some(Vector::new([1.0, 0.0, 0.0]))
})?;

The integrator validates unit-norm at every sample; a non-unit direction is a typed error, not silently corrected.

Reusable DTI-volume pipeline

DtiVolume already owns the validated image-index grid, fitted-voxel mask, and tracking anisotropy floor. ritk-tractography owns the remaining policy in DtiTractographyConfig: inclusive FA-threshold seed selection, an optional evenly strided seed cap, and Euler integration through the volume’s direction field. This keeps downstream applications from reimplementing the policy at a CLI or application boundary.

use ritk_tractography::{
    DtiTractographyConfig, TrackingDirection, TractographyConfig,
    dti_volume_seed_points_with_mask, dti_volume_tractography_with_mask,
};

let tracking = TractographyConfig::new(
    0.5,
    1_000,
    60.0,
    TrackingDirection::Bidirectional,
)?;
let policy = DtiTractographyConfig::new(0.25, 10_000, tracking)?;
let seed_mask = vec![true; volume.maps().len()];
let seeds = dti_volume_seed_points_with_mask(
    &volume,
    policy.seed_anisotropy(),
    policy.max_seeds(),
    Some(&seed_mask),
)?;
let tracks = dti_volume_tractography_with_mask(&volume, policy, Some(&seed_mask))?;

The seed threshold is inclusive and the volume mask is authoritative, so an unfitted voxel cannot become a seed at threshold zero. max_seeds == 0 selects all qualifying voxels; a nonzero cap uses a stride through the qualifying storage order rather than truncating at the first cap-sized prefix. The orchestration function returns NoSeeds with the threshold and eligible fitted FA peak instead of producing an apparently successful empty tractogram. The mask-aware entry points require one boolean per DTI voxel in [depth, row, column] order and apply it only to starting voxels; passing None preserves the unmasked entry points’ behavior.

The reusable DTI-volume example fits a known two-regime synthetic volume, verifies the seed and streamline counts, and generates the figure used on that page. It verifies software behavior only; it does not claim anatomical or clinical validity.

Create and validate a tractogram

A smooth tractogram is not a correctness oracle. The strongest local check starts from a direction field whose domain and expected path behavior are known before tracking. RITK’s runnable book example uses the curved analytical bundle below so creation and validation exercise the same public API without using the implementation as its own oracle.

1. Define a physical field and seeds

The closure returns a unit tangent inside the bundle and None outside it. This makes the dashed bundle boundary in the figure an executable stopping oracle rather than a visual annotation:

use ritk_spatial::{Point, Vector};

fn bundle_center(x: f64) -> f64 {
    20.0 + 5.0 * (x / 8.0).sin()
}

fn direction_field(point: &Point<3>) -> Option<Vector<3>> {
    let [x, y, z] = point.to_array();
    if !(2.0..=38.0).contains(&x)
        || (y - bundle_center(x)).abs() > 4.0
        || z.abs() > 0.5
    {
        return None;
    }

    let slope = 0.625 * (x / 8.0).cos();
    let norm = (1.0 + slope * slope).sqrt();
    Some(Vector::new([1.0 / norm, slope / norm, 0.0]))
}

let seed_x = 20.0;
let seeds = [-2.5, -1.25, 0.0, 1.25, 2.5]
    .map(|offset| Point::new([seed_x, bundle_center(seed_x) + offset, 0.0]));

The seeds, field, and output use physical coordinates. A voxel-index field must be converted with the image geometry before export; use TractographyResult::map_points when tracking occurred in index space.

2. Create streamlines

Construct a validated configuration, then integrate the seeds. The selected step size, turn limit, and field domain are part of the experiment and must be reported with any downstream result:

use ritk_tractography::{
    TrackingDirection, TractographyConfig, euler_tractography,
};

let config = TractographyConfig::new(
    0.35,
    160,
    20.0,
    TrackingDirection::Bidirectional,
)?;
let result = euler_tractography(&seeds, config, direction_field)?;

3. Validate value semantics before rendering

The analytical example requires one output per seed, boundary termination in both directions, and containment of every emitted point. These assertions fail if tracking silently drops a valid seed, records the wrong stopping cause, or appends the first out-of-domain proposal:

use ritk_tractography::TerminationReason;

assert_eq!(result.seeds_attempted(), seeds.len());
assert_eq!(result.streamlines_generated(), seeds.len());

for streamline in result.streamlines() {
    assert_eq!(
        streamline.forward_termination(),
        TerminationReason::FieldBoundary,
    );
    assert_eq!(
        streamline.backward_termination(),
        Some(TerminationReason::FieldBoundary),
    );

    assert!(streamline.geometry().points().iter().all(|point| {
        direction_field(&Point::new([point.x, point.y, point.z])).is_some()
    }));
}

The runnable source also verifies the independently known diffusion axis before rendering. Run it from the repository root:

cargo run --locked -p ritk-diffusion --example book_diffusion_tractography -- \
  docs/book/figures/diffusion_tractography.svg

The command writes the checked SVG plus .trk, .tck, and .tsf artifacts. See Signal to Streamlines for the signal model, angular oracle, exported scalar layout, and interpretation of each panel.

Validation ladder

Each check supports a different claim. Passing a lower level does not imply a higher one:

The first four levels are covered by the deterministic book example and crate tests. Exporting a file does not complete the interoperability level, and the real-subject figure demonstrates execution on scanner data rather than ground truth. Diffusion tractography can produce coherent false-positive bundles even when local orientations are accurate; Maier-Hein et al.’s ground-truth challenge therefore cautions against treating tractography alone as anatomical evidence (Nature Communications 8, 1349, 2017, Results: “Tractograms contained more invalid than valid bundles”).

Streamline export

TractographyResult provides export methods for all three tractogram interchange formats. Tractography points are in physical millimetre coordinates (the native coordinate system for all three formats).

.trk — DSI Studio / TrackVis

to_trk writes the header with an identity voxel-to-RAS affine. Callers that need anatomical space should set dim and voxel_size to match the reference image:

let trk = result.to_trk(
    [128, 128, 60],   // dim
    [1.5, 1.5, 2.0],  // voxel_size (mm)
);
// trk can be written via ritk_trk::write_trk(&trk, path)?;

to_trk_header accepts an optional custom vox_to_ras affine:

let trk = result.to_trk_header(dim, voxel_size, Some(vox_to_ras_affine));

to_trk_with_scalars attaches per-point scalar values (e.g. FA, MD) for DSI Studio colour-coding. Each inner Box<[f32]> must contain n_points × n_scalars values in row-major (per-point scalar stride) order:

let trk = result.to_trk_with_scalars(
    dim, voxel_size,
    &["FA", "MD"],
    fa_md_scalars,
);

.tck — MRtrix3

to_tck writes with a default header (Float32LE datatype, no transform):

let tck = result.to_tck();
// tck can be written via ritk_tck::write_tck(&tck, path)?;

to_tck_header accepts optional mrtrix_version, comments, and a 4×4 voxel-to-scanner transform:

let tck = result.to_tck_header(
    Some("3.0.4".into()),
    Some("RITK Euler tractography".into()),
    Some(transform_matrix),
);

Per-point scalars are not stored natively in the .tck format; use the MRtrix3 tckmap weights sidecar (ritk_tck::write_tck_weights) when scalar export is needed.

.trx — Tractography Reference eXchange

to_trx writes with a default header ("float32" dtype):

let trx = result.to_trx();
// trx can be written as a directory via ritk_trx::write_dir(&trx, path)?;

to_trx_with_dpv attaches per-vertex data arrays (DPV) such as FA and MD. The caller must populate the header’s dpv map with TrxArrayDef entries declaring each array’s dtype and component count, and provide the corresponding raw encoded byte buffers:

use std::collections::HashMap;
let mut trx = result.to_trx_with_dpv(HashMap::new());
trx.header.dpv.insert(
    "FA".into(),
    ritk_trx::TrxArrayDef { dtype: "float32".into(), n_components: 1 },
);
trx.dpv_data.insert("FA".into(), fa_bytes);

Format comparison

Error types

TractographyError covers every failure mode with seed and step indices for localisation:

A seed at which the field returns None, or whose integration produces fewer than two points, is an expected untrackable seed — it produces no line and is not an error.

Output types

Streamline pairs Gaia polyline geometry with termination diagnostics:

TractographyResult aggregates the output of one call to euler_tractography:

What the current algorithm establishes

The implementation and tests establish deterministic geometry, sign continuity, bounded memory growth, explicit stopping reasons, and rejection of invalid field samples. They do not establish anatomical validity, clinical utility, crossing-fibre resolution, uncertainty, or invariance to acquisition and preprocessing choices. Apply the validation ladder above and state the highest level for which an independent oracle was actually run.

The signal-to-streamlines example closes the loop end-to-end: known tensor → synthetic signals → model fit → direction field → streamlines. Each visual element — seeds, orientation glyphs, domain boundaries, generated trajectories — has a defined meaning.

The human tractography and connectomics example then runs the same public tracking boundary on a checksummed 160-volume human HARDI acquisition. It adds exact DWI/parcellation alignment, endpoint accounting, an 84-region connectivity matrix, and an explicit separation between internal software validation and biological validation.

Example: Signal to Streamlines

This runnable example connects two independently verifiable stages: a single-fiber diffusion-tensor signal is converted to an analytical Q-ball ODF, then a separate curved unit-vector field is integrated into streamlines.

Read the numbered panels from left to right:

1. Blue points are the 48 actual unit gradient directions supplied to GradientScheme, all at (b=1500) s/mm². Two additional b0 measurements establish (S_0).
2. Orange points show the generated tensor signal against alignment with the known fiber axis. Stronger alignment causes greater attenuation. The blue polar shape is the ODF evaluated from RITK’s fitted coefficients; the dashed red line is the independent analytical axis. The reported angular error is computed before the SVG is written.
3. Gray line segments are local direction samples, orange points are seeds, dashed curves bound the trackable field, and blue curves are Gaia polylines returned by euler_tractography. No curve contains the first proposal outside the dashed boundary.

The ODF panel and tractography panel intentionally use separate fields. This prevents the figure from implying that one voxel’s ODF is enough to form a whole tract: tractography requires a spatial field of local orientations.

Source and command

Source: crates/ritk-diffusion/examples/book_diffusion_tractography.rs

cargo run -p ritk-diffusion --example book_diffusion_tractography -- \
  docs/book/figures/diffusion_tractography.svg

The example fails unless:

• a deterministic one-degree full-sphere search places the analytical Q-ball peak within two degrees of the known antipodal x axis;
• all five seeds produce streamlines;
• every emitted point remains inside the analytical vector-field domain; and
• the figure can be written without an unbounded or infallible allocation assumption in the library path.

This synthetic result verifies the stated numerical and geometry contracts. It is not a claim of tractography accuracy on patient data.

Scalar export

The example also estimates a diffusion tensor from the same signals to compute FA and MD, then exports the tractography result in all three interchange formats:

• .trk — with per-point FA and MD scalars for DSI Studio colour-coding via to_trk_with_scalars.
• .tck — plain streamline geometry via to_tck_header, with provenance metadata (MRtrix version and a comment).
• .tck weights sidecar — a FA scalar file (.tsf) via write_tck_weights, compatible with MRtrix3 tckmap and mrview -tck_weights_in.

let dti_tensor = estimate_dti(&scheme, &signals, DtiConfig::default())?;
let fa = dti_tensor.fa();
let md = dti_tensor.md();

// .trk with FA + MD interleaved per point.
let fa_md_scalars: Vec<Box<[f32]>> = result
    .streamlines()
    .iter()
    .map(|s| {
        let n = s.geometry().len();
        let mut scalars = Vec::with_capacity(n * 2);
        for _ in 0..n {
            scalars.push(fa as f32);
            scalars.push(md as f32);
        }
        scalars.into_boxed_slice()
    })
    .collect();
let trk = result.to_trk_with_scalars(
    [128, 128, 1], [1.0, 1.0, 1.0], &["FA", "MD"], fa_md_scalars,
);
trk.write(&mut file)?;

// .tck with provenance metadata.
let tck = result.to_tck_header(
    Some("3.0.4".into()),
    Some("RITK book example: synthetic tensor → Q-ball → Euler tractography".into()),
    None,
);
tck.write(&mut tck_file)?;

// .tck weights sidecar — one FA scalar per point.
let fa_only: Vec<Box<[f32]>> = result
    .streamlines()
    .iter()
    .map(|s| {
        let n = s.geometry().len();
        (0..n).map(|_| fa as f32).collect::<Vec<_>>().into_boxed_slice()
    })
    .collect();
ritk_tck::write_tck_weights(&fa_only, ritk_tck::TckDatatype::Float32LE, &mut tsf_file)?;

All three files are written alongside the SVG figure:

The .trk file can be opened in DSI Studio or TrackVis; the .tck file and .tsf sidecar can be loaded in MRtrix3’s mrview with -tck_weights_in.

Example: Reusable DTI-Volume Tracking

This example exercises the complete reusable DTI tractography boundary. It fits a known prolate tensor for eight image-axis voxels and a low-anisotropy tensor for four voxels, places those maps in a DtiVolume, selects seeds by an inclusive FA threshold, and tracks through the validated volume.

The first panel plots the fitted FA for every voxel. Blue bars are the values returned by DiffusionMaps::fractional_anisotropy; the red line is the seed threshold. The orange points are the seeds returned by dti_volume_seed_points, not hand-selected coordinates. The second panel plots the streamlines returned by dti_volume_tractography; its boundary marks the first low-anisotropy voxel where DtiVolume stops the field.

The example asserts before writing the SVG that:

• the fitted high- and low-anisotropy regimes lie on the expected sides of the seed and tracking thresholds;
• the tracking result attempted and generated exactly the computed seed count;
• every generated line terminates at a field boundary or the configured step limit; and
• the SVG contains exactly one bar, seed marker, and streamline primitive for each corresponding computed value.

Source and command

Source: crates/ritk-tractography/examples/book_dti_volume_tractography.rs

cargo run --locked -p ritk-tractography --example book_dti_volume_tractography -- \
  docs/book/figures/dti_volume_tractography.svg

The output is deterministic because the gradient scheme, tensor values, and rendering coordinates are fixed. This is a code and numerical-contract verification example. It is not anatomical or clinical validation.

Human Tractography and Connectomics

The analytical signal-to-streamlines example can test recovery against a known fibre axis. Human diffusion MRI has no equivalent voxelwise ground truth. This workflow therefore makes a narrower claim: RITK can reproducibly fit a complete public human acquisition, create a whole-brain tractogram, assign endpoints to an aligned anatomical parcellation, and expose the validation evidence and residual limitations with the result.

Data and provenance

The input is the Stanford HARDI dataset: 81 × 106 × 76 human diffusion MRI voxels, 150 directions at (b=2000) s/mm², 10 (b=0) volumes, and a DWI-aligned reduced FreeSurfer parcellation. The repository record states that the diffusion series was motion-corrected to the mean (b=0) image but was not eddy-current corrected.

The data use the Open Data Commons PDDL 1.0 public-domain dedication. Stanford’s record also prohibits attempts to identify participants or infringe their privacy. This example performs aggregate method demonstration only.

Imaging bytes are gitignored. The download script resolves each Stanford object directly and verifies the five MD5 digests published with the dataset fetcher:

bash test_data/diffusion/download.sh
cargo run --release -p ritk-diffusion --example book_brain_tractography

The executable exits without writing when the data are absent. A successful run replaces both the figure above and the complete connectivity matrix JSON.

Creation workflow

The example executes one deterministic path:

1. Parse all 160 FSL gradient entries. Six-decimal unit vectors are normalised only when their norm error is within the derived (10^{-5}) rounding envelope; larger errors remain invalid input.
2. Read the DWI and reduced parcellation through RITK’s native NIfTI boundary. Their shape, origin, spacing, and direction matrix must match exactly.
3. Fit a diffusion tensor at every valid voxel. Seed image labels 1 and 2 (cerebral white matter and corpus callosum) at fractional anisotropy (\mathrm{FA}\geq0.25).
4. Track bidirectionally with 0.5-voxel Euler steps, sign-invariant trilinear orientation interpolation, a 60-degree turn limit, and a continuation floor of (\mathrm{FA}=0.15). The fitted maps retain GradientFrame; DtiVolume validates the FSL ImageAxis frame and performs the single (column,row,depth) to [depth,row,column] permutation at its image-grid boundary before integration. No example-local adapter is involved.
5. Convert the tracker’s [depth,row,column] coordinates to parcellation [x,y,z], then count each endpoint pair in the undirected regional matrix. White-matter labels are background for endpoint assignment; only grey-matter labels actually present in the image become graph nodes.
6. Transform each streamline through the reference image geometry before reporting its physical length in millimetres.

These are explicit method choices, not tuned claims of anatomical optimality. The source is book_brain_tractography.rs.

Computed result

The committed artifacts are generated from the checksummed acquisition:

The strongest non-self edges are descriptive outputs of this tracking and endpoint policy:

The complete JSON preserves all 84 labels, the fully symmetric 84 × 84 weight array — both (i, j) and (j, i) carry the weight — and the streamline accounting: how many were supplied, how many produced an inter-region edge, how many stayed within one region, and how many could not be assigned. A table entry is a streamline count, not an axon count, connection probability, or effect size.

Reading the figure

1. Directionally encoded FA. Brightness is fractional anisotropy. Hue is the absolute principal-eigenvector direction in anatomical axes: red left-right, green anterior-posterior, blue superior-inferior. The callosal arc, bilateral projection anatomy, dark ventricles, and darker cortical rim are visible. Colour channels are derived from the NIfTI direction matrix instead of assuming storage axes.

2. Whole-brain streamlines. The middle panel projects a bounded sample of the three-dimensional tractogram over the same axial plane. Each polyline’s colour encodes its endpoint-to-endpoint direction. All 9,737 streamlines contribute to the reported measures and connectome even though at most 900 are drawn.

3. Endpoint connectome. The symmetric heat map uses (\log(1+w)) colour scaling so both weak and strong nonzero edges remain visible. The upper-left and lower-right blocks are within-hemisphere edges; the off-diagonal blocks are inter-hemisphere edges. Empty cells are zero, not missing observations.

Validation evidence

The workflow checks claims at distinct levels:

These checks establish reproducibility and internal consistency. The DIPY Stanford connectivity example uses the same aligned DWI/label relationship and endpoint-count interpretation, which supplies an independent convention check rather than a numeric oracle.

Limits of the result

• The input is one subject and retains acquisition artefacts: no eddy-current or susceptibility correction is applied here.
• A tensor model uses one orientation per voxel even though the HARDI acquisition can support multi-fibre modelling. Crossing fibres therefore remain unresolved.
• Tracking is deterministic and is not anatomically constrained. The 24.1 % exact endpoint-assignment rate exposes that many tracks terminate in white matter or background under this policy; the example does not dilate labels or search nearby cortex to inflate the count.
• Raw streamline counts depend on seeding, stopping, length, geometry, and parcellation. No SIFT-class correction, null model, population normalisation, or test-retest analysis is applied.
• A coherent tractogram can contain false-positive bundles. The international tractography challenge found this even among state-of-the-art pipelines (Maier-Hein et al., Nature Communications 8, 1349, 2017).

Consequently this is method and software validation, not anatomical ground truth or clinical validation. Stronger biological evidence requires an independent reference such as physical phantoms, tracer or histological data, and a prospectively specified population protocol.

Anatomical Parcellation

A parcellation is a label volume: every voxel carries the identifier of the anatomical region it belongs to, with the reserved label 0 meaning background. It is the shared vocabulary between whatever produces one — an atlas registered onto a subject, a segmentation, a surface annotation rasterised into a volume — and whatever consumes one: connectome construction, region-wise statistics, targeting.

ritk-parcellation owns the type and depends on nothing but spatial primitives. ritk-registration produces one by warping a labelled atlas onto a subject.

Geometry: why the direction cosines are not optional

A label volume is only useful if a physical point can be mapped to the voxel containing it. ParcellationGrid carries the standard medical-imaging affine:

\[ p = \mathbf{o} + R \, (\mathbf{s} \odot \mathbf{i}) \]

where \(\mathbf{i}\) is the continuous voxel index, \(\mathbf{s}\) the per-axis spacing, \(R\) the direction cosine matrix whose columns are the physical directions of the index axes, and \(\mathbf{o}\) the physical position of voxel \((0,0,0)\)’s centre.

Reducing this to \(\mathbf{o} + \mathbf{s} \odot \mathbf{i}\) — treating the index axes as parallel to the physical ones — is correct only for an axis-aligned volume. Almost no acquired volume is: a scanner stores an oblique slice stack, and NIfTI’s qform/sform, NRRD’s space directions, and DICOM’s image orientation all exist to record that obliquity.

Dropping it does not fail loudly. It returns the label of a different region, displaced by up to the volume’s extent times the sine of the obliquity. Connectome edges built on such lookups are wrong in a way nothing downstream catches, because every label involved is a legitimate label.

use ritk_parcellation::{Parcellation, ParcellationGrid};

// The general case: shape, spacing, origin, and row-major direction cosines.
let grid = ParcellationGrid::new(
    [256, 256, 128],           // [nx, ny, nz]
    [1.0, 1.0, 1.5],           // spacing in mm, by spatial axis
    [-127.5, -127.5, -95.0],   // physical position of voxel (0, 0, 0)
    direction_cosines,         // [f64; 9], row-major
)?;

// The axis-aligned special case.
let grid = ParcellationGrid::axis_aligned([256, 256, 128], [1.0, 1.0, 1.5], [0.0; 3])?;

Labels are stored z-major — index [ix, iy, iz] at offset iz·ny·nx + iy·nx + ix. That is this crate’s contract, not a universal one: NIfTI, NRRD, and MetaImage store volumes in this order, but a format that does not converts at the I/O boundary rather than here.

Construction and lookup

Parcellation::new rejects a label array that does not cover the grid, and one in which every voxel is background: a parcellation with no regions cannot answer any question asked of it, so it fails where it is built rather than returning empty results at every call site.

let parcellation = Parcellation::new(labels, grid, region_names)?;

// Nearest-neighbour lookup at a physical point.
let label = parcellation.label_at(&point);

label_at distinguishes two answers that are often conflated:

Outside the field of view is not the same claim as inside it but unassigned, so they are not collapsed. A caller counting discarded streamlines needs both.

Region statistics

region_statistics() walks the volume once and returns per-region size, position, and extent — one pass rather than one pass per region, which for a hundred-region atlas over a million voxels is the difference between a million reads and a hundred million.

Two of these are load-bearing beyond description:

Volume normalises a connectome. A raw streamline count between two regions grows with how big those regions are, because a larger region presents a larger surface for streamlines to terminate on. Comparing a raw count across regions of different size, or across subjects whose regions differ in size, therefore compares anatomy to arithmetic.

Centroid turns a region into a target point. Note that for a concave or disconnected region the centroid can lie outside the region itself — a C-shaped gyrus has its centre of mass in the gap. That is a property of centroids, not a defect; a point guaranteed to be inside the region is a different query.

Coarsening an atlas

A fine-grained atlas is collapsed to a coarser one through remap_labels, which is a total mapping: a label the caller does not mention keeps its value, and returning 0 removes a region. retain_regions keeps a named subset — the usual preparation for a cortex-only connectome.

// Merge the two hemispheres' matching parcels.
let bilateral = parcellation.remap_labels(|label| label % 1000, names)?;

// Keep only the cortical parcels.
let cortex = parcellation.retain_regions(&cortical_labels)?;

Nearest labelled voxel

Assigning a streamline endpoint to a region by reading the label directly under it discards most of a tractogram. Streamlines are tracked through white matter and stop where the orientation field stops being coherent, at or just short of the grey-matter boundary — while a cortical parcellation labels the grey matter and leaves the white matter background. The endpoint lands in a region-less voxel and the streamline is dropped despite ending exactly where it should.

NearestLabelSearch recovers those. The offsets within a radius are enumerated and sorted once, so a per-streamline loop does not re-derive the same neighbourhood a million times.

The walk does not simply return the first labelled voxel it meets. Offsets are ordered by the distance between voxel centres, while what matters is the distance from the point, which sits anywhere inside its voxel — the two differ by up to half a voxel diagonal. For an endpoint near a parcel boundary that is exactly the difference between the parcel it is in and the one across the border, so every candidate is scored and the nearest kept, with an early exit once no remaining offset can improve on the best found.

use ritk_parcellation::NearestLabelSearch;

let search = NearestLabelSearch::new(parcellation.grid(), 2.0)?;
if let Some(found) = search.find(&parcellation, &endpoint) {
    // found.label, found.index, found.distance
}

The radius trades recall against specificity, and the trade is not free in either direction. Too small and the tractogram is decimated. Too large and endpoints reach past their own gyrus into a neighbouring parcel across a sulcus — anatomically adjacent, functionally unconnected — and manufacture edges no fibre supports. Cortical thickness is 2–4 mm and sulcal walls approach within a millimetre of each other, so a few millimetres is where the two failure modes balance.

The search measures distance, not connectivity: it cannot distinguish “the nearest parcel across the sulcus” from “the parcel this fibre actually entered”. That is why the assignment reports the distance it required.

Full-brain parcellation from an atlas

ritk-registration parcellates a subject by deforming an already-labelled brain onto it and carrying the labels along.

atlas intensity ──► register to subject ──► deformation
atlas labels    ──────────────────────────────┴──► warp ──► Parcellation

use ritk_registration::{
    AtlasParcellationConfig, LabelledAtlas, parcellate_with_atlas_set,
};

let result = parcellate_with_atlas_set(
    &subject_t1,
    &atlases,
    &AtlasParcellationConfig::default(),
)?;

let parcellation = result.parcellation;   // on the subject's own grid
let agreement = result.agreement;         // per-voxel, in [0, 1]
let quality = result.registration_quality; // final CC per atlas

Running it without writing Rust

The same pipeline is a command and a Python call, so a parcellation can be produced by whichever surface the rest of the analysis lives on.

ritk parcellate atlas --subject T1.nii.gz \
                      --atlas-intensity a1.nii.gz --atlas-labels a1_dseg.nii.gz \
                      --atlas-intensity a2.nii.gz --atlas-labels a2_dseg.nii.gz \
                      --output dseg.nii.gz --agreement agreement.nii.gz

import ritk

result = ritk.registration.parcellate_with_atlases(
    subject, atlas_intensities, atlas_labels, fusion="majority"
)
parcellation = result.parcellation          # feeds build_connectivity_matrix
agreement = result.agreement                # [Z, Y, X], in [0, 1]

Every atlas must already lie on the subject’s grid. All three surfaces reject a mismatch rather than resampling it: a registration recovers a deformation, never a resampling, so an atlas of the wrong size quietly accepted would produce labels for a different brain.

Labels are warped, never interpolated

Label values are identifiers, not measurements. Region 17 and region 19 do not average to region 18, and any interpolation that produced 18 would invent an anatomical claim out of two unrelated ones — silently, since 18 is a valid label. Every resampling in the pipeline is nearest-neighbour for that reason. The cost is a jagged boundary at the voxel scale, which is the honest representation of what a label map can say.

Why one atlas is usually not enough

A single atlas transfers not only its anatomy but its idiosyncrasies, and wherever the registration is locally wrong, the labels are locally wrong with no signal that they are. Registering several independently labelled brains and fusing their votes is the standard remedy: an error must be shared by a majority of the atlases to survive, and disagreement becomes measurable rather than invisible.

Majority voting treats every atlas as equally trustworthy everywhere — right when the atlases are interchangeable, wrong when some registered better than others in a particular region. Joint label fusion lets a well-registered atlas outvote a poorly registered one in the region where that is true, at the cost of a patch comparison and a small dense solve per voxel.

The returned agreement map is the per-voxel confidence. Low agreement marks where the result is a coin toss between neighbouring parcels — usually the boundaries, which is exactly where the answer matters most for a connectome, since that is where streamlines end.

Checking that it worked

A parcellation always looks plausible, so the question is not whether one came back but whether it is right. The atlas parcellation example synthesises a subject whose correct parcellation is known, deforms three atlases onto it — one of them deliberately mislabelled — and reports Dice against the truth alongside the agreement map. It is also where the choice between the two fusion rules becomes a measurement rather than a preference.

What atlas propagation cannot do

It transfers a predefined parcellation. It cannot discover a region the atlas does not contain, and it cannot represent anatomy the registration could not reach — a resected cavity, a large lesion, or a malformation has no counterpart in a healthy atlas, and the labels warped over it are meaningless rather than absent. The agreement map is the only signal of that, and it is a weak one when every atlas is equally wrong.

FreeSurfer formats

The freesurfer module reads and writes the FreeSurfer surface family:

use ritk_parcellation::freesurfer::{lut, SurfaceAnnotation};

let table = lut::read(lut_file)?;                // names and colours per label
let annotation = SurfaceAnnotation::read(file)?; // per-vertex structure indices
let names = lut::region_names(annotation.color_table());

An annotation identifies each vertex’s structure by colour: the file stores red + green·2⁸ + blue·2¹⁶ per vertex, and the reader resolves it against the embedded table to a structure index. A value no entry has, or two entries with one colour, is rejected rather than guessed. The lookup table’s fourth column is transparency, not opacity — FreeSurfer derives alpha = 255 - t, and the reader stores that alpha in the table’s RGBA colour. Both tables are the segmentation editor’s own LabelTable, so a FreeSurfer parcellation opens with its names and colours intact.

Every reader treats its input as hostile: counts are bounded, storage grows only as real input backs it, and every failure is a typed FreeSurferError.

A SurfaceAnnotation labels vertices of a mesh, not voxels. Converting one to a volumetric Parcellation needs the geometry those vertices belong to — Surface reads the binary triangle format — and a rasterisation of the cortical ribbon, which rasterise_ribbon performs.

use ritk_parcellation::freesurfer::{Surface, rasterise_ribbon};

let white = Surface::read(std::fs::File::open("lh.white")?)?.translated(c_ras);
let pial = Surface::read(std::fs::File::open("lh.pial")?)?.translated(c_ras);

let (parcellation, report) = rasterise_ribbon(&white, &pial, &annotation, &grid, 16)?;
println!("{} of {} columns filled", report.columns - report.unfilled_columns, report.columns);

The frame, which is where this goes wrong

FreeSurfer stores surfaces in surface RAS (tkrRAS), not the scanner frame a volume carries. The two differ by a translation: surface RAS puts the origin at the centre of the conformed 256³ volume, scanner RAS puts it where the scanner did. The offset is the volume’s c_ras, typically tens of millimetres — enough to place a cortical ribbon outside the brain entirely, without ever failing.

The reader returns coordinates as stored and does not guess, because the surface file does not contain the offset. Surface::translated applies it once known. Rasterising with the wrong frame lands no column inside the volume, which is rejected rather than returned as an empty parcellation.

What the rasterisation does and does not do

The two surfaces share a vertex numbering: vertex i of lh.white and vertex i of lh.pial are the inner and outer end of the same cortical column. Walking that segment and stamping the vertex’s label into every voxel it crosses fills the ribbon one column at a time — the same approach as FreeSurfer’s mri_surf2vol --fill-ribbon.

It fills the ribbon; it does not tessellate it. A voxel no column crosses stays background even if it lies geometrically inside the ribbon, which happens where the mesh is coarse relative to the voxel grid. steps closes gaps along a column and does nothing for gaps between columns, so a coarse mesh on a fine grid is a limitation of the input rather than of the setting. RibbonReport::unfilled_columns is the signal that the two disagree in resolution, and contested_voxels counts where two parcels wanted the same voxel — concentrated at boundaries and inside sulcal folds, which is exactly where a connectome’s endpoints land.

lut::region_names is separately usable as the region_names of a volumetric parcellation from any source.

Next

The Connectome Construction and Graph Measures chapter builds the region graph from a parcellation and a tractogram.

Example: Atlas Parcellation of a Synthetic Subject

Atlas propagation is easy to run and hard to check. A parcellation is a volume of integers that always looks plausible, so a pipeline that has silently mis-registered returns something a reader cannot distinguish from a correct result. This example removes that problem by synthesising the subject: the correct parcellation is known exactly, so every number below is an overlap against a truth rather than a claim that the pipeline returned something.

cargo run --release -p ritk-registration --example book_parcellation

What is synthesised

The subject is a 16 × 20 × 24 volume at 2.0 × 1.5 × 1.0 mm containing three structures of different size, position, and aspect ratio, each at its own intensity. The shape and the spacing are unequal on every axis on purpose: a cubic volume on an isotropic grid cannot fail an axis-order test, because reversing the axes leaves every array the same length and every voxel the same size, so a transposition would cancel itself out and the example would pass while the geometry was wrong.

Each atlas is that anatomy displaced by a different whole-voxel shift. That is the situation atlas propagation exists for — the anatomy corresponds but the coordinates do not — and it is what the registration has to undo before the labels mean anything. The displacement is whole voxels because a label volume cannot be interpolated: shifting by half a voxel would require inventing labels at the boundary, which is the one thing the chapter says never to do.

The third atlas is mislabelled: its anatomy matches the others exactly, and only its label volume swaps two structures. The distinction matters. Swapping its intensities as well would make it a different brain, which a registration is entitled to match poorly, and the fusion would then be resolving a disagreement about anatomy. Mislabelling is the case fusion actually exists for: every atlas fits the subject equally well and only the naming is in dispute.

What the run reports

atlases, before any registration:
  0: shifted by [1, 0, -1] voxels, Dice vs truth [0.69, 0.69, 0.73]
  1: shifted by [-1, 1, 1] voxels, Dice vs truth [0.57, 0.51, 0.58]
  2: shifted by [0, -1, 2] voxels, Dice vs truth [0.60, 0.00, 0.00]

majority vote:
  Dice vs truth [1.00, 1.00, 1.00]
  3 regions, mean agreement 0.964
  final cross-correlation per atlas [0.9642, 0.9254, 0.9001]

mean agreement inside each structure, majority vote:
  structure 1: 0.952
  structure 2: 0.667
  structure 3: 0.667

Three things in that output are worth reading carefully.

The atlases start between 0.51 and 0.73 Dice and the fused result reaches 1.00. That gap is the registration doing the work; without it the labels would arrive at the displaced positions and the parcellation would be the average of three wrong answers.

The mislabelled atlas scores 0.00 on the two structures it swapped, and the fused result is still 1.00 on both. That is majority voting surviving one dissenting atlas, which is the reason for using more than one.

The agreement is 0.952 on the structure all three atlases named correctly and exactly 0.667 on the two the third atlas renamed — two of three. The agreement map is a measured vote share, not a decoration, and the figure’s third panel shows it as the grey blocks over the two disputed structures against the white unanimous background. The thin grey strips elsewhere are the parcel boundaries, where the residual registration error is: precisely where a connectome’s streamline endpoints land.

Where weighting has nothing to work with

The example also runs both fusion rules over atlases with none mislabelled. Nothing is then in dispute except each atlas’s own registration error:

control — interchangeable atlases, none mislabelled:
  majority vote: Dice vs truth [1.00, 1.00, 1.00]
  joint fusion:  Dice vs truth [0.88, 0.95, 0.91]

Voting wins, and the reason is the assumption each rule makes rather than a fault in either. Joint label fusion buys its advantage by trusting whichever atlas matches the subject best in each neighbourhood — which requires the atlases to differ in local registration quality and for the intensities to reveal that difference. Here all three are the same anatomy under different whole-voxel shifts, so they fit equally well everywhere, the weights come out near-equal, and following the locally best one means inheriting that one atlas’s warp error. Voting instead averages three independent errors, which is exactly what averaging is good at.

Read that as a statement about the fixture, not a ranking. On real cohorts the atlases are different brains that register with genuinely different local quality, which is the situation joint fusion was designed for and this synthetic one deliberately is not. What the control establishes is the boundary of the claim in the chapter: majority voting is right when the atlases are interchangeable, and this is what “interchangeable” costs the alternative.

What the example does not show

It does not show whether the pipeline works on a brain. Synthetic blocks have sharp intensity boundaries and no noise, no partial volume, no lesion, and no anatomical variability, so the registration’s job here is far easier than on real data. What the example establishes is that the pipeline is wired correctly end to end and that its reported agreement means what it says — necessary conditions, not sufficient ones. Validation against a labelled cohort is a different tier of evidence, and this is not it.

GIFTI Surface Interchange

FreeSurfer’s own surface files (lh.white, lh.curv, lh.aparc.annot) are read by ritk-parcellation::freesurfer. GIFTI is the format other tools use to exchange the same content, and ritk-gifti reads and writes it. The reference is the GIFTI Surface Data Format 1.0 specification (14 January 2011).

Layout

A GIFTI file is XML: a root GIFTI element holding optional file metadata, an optional label table, and one or more DataArray elements. Each data array declares:

The file-type extensions only name conventions: a .surf.gii holds one point set and one triangle array, a .func.gii or .shape.gii one float per vertex, and a .label.gii a label table plus integer keys per vertex.

use ritk_gifti::GiftiImage;

let surface = GiftiImage::read(std::fs::File::open("lh.pial.surf.gii")?)?.surface()?;
println!("{} vertices, {} triangles", surface.vertices.len(), surface.triangles.len());

The encoding named gzip is zlib

GZipBase64Binary is base64 over a zlib stream, not a gzip file: the specification says the data “is compressed using ZLIB” (section 5.0), and nibabel inflates it with zlib.decompress. A reader that expects the gzip header rejects every conforming file.

What the reader checks

The reader treats every document as hostile input:

• Element nesting follows the DTD. A Data outside a DataArray, or a second root element, is an error. Elements the DTD does not define are skipped.
• A shape has one to six axes and at most 10⁸ values, checked before any payload is decoded.
• Each payload must decode to exactly the declared number of values. A zlib stream is inflated no further than one byte past that length, so a stream that expands beyond its shape is rejected without being inflated.
• NumberOfDataArrays must match the arrays present, and label colours must lie in [0, 1].

Column-major arrays are kept as stored, and DataArray::row_major transposes them. The writer always emits row-major, little-endian data, the forms every reader accepts.

What it does not do

ExternalFileBinary data is reported as unsupported. Resolving it needs the directory of the XML file, and a reader of bytes does not have that path. CIFTI, which embeds GIFTI-like arrays in NIfTI-2 containers, is out of scope.

Connectome Construction and Graph Measures

A connectome is the graph whose nodes are anatomical regions and whose edges summarise the streamlines running between them. ritk-connectome turns a tractogram plus a Parcellation into that graph, and then measures it.

Parcellation (labels) ──┐
                        ├──► ConnectivityMatrix ──► graph measures
Streamlines (polylines) ┘

The crate sits at the end of the diffusion-MRI pipeline documented in this part: the Gradient Schemes chapter validates acquisition metadata, Diffusion Models fits orientation models, Tractography produces streamlines, the Parcellation chapter supplies the regions, and this chapter builds the graph.

What a connectome edge does and does not mean

An edge weight is a property of the tractogram, not of the brain. Streamline counts are known not to be proportional to axonal density: tracking systematically favours short, straight, high-anisotropy paths and systematically loses long, curved, or crossing ones.

Two of the weightings below exist to remove the leading geometric parts of that dependence — pathway length and region size — and none of them turns a count into a measurement of connection strength. Comparison across subjects processed identically is the defensible use; absolute interpretation of a single weight is not.

Construction

Two decisions separate a tractogram from a connectome, and both change the answer.

use ritk_connectome::{
    ConnectomeConfig, EdgeWeighting, EndpointAssignment, build_connectivity_matrix,
};

let config = ConnectomeConfig::new()
    .with_assignment(EndpointAssignment::RadialSearch { radius_mm: 2.0 })
    .with_weighting(EdgeWeighting::InverseNodeVolume);

let matrix = build_connectivity_matrix(&parcellation, &streamlines, &config)?;

The command line and Python reach the same builder:

ritk tract connectome --tractogram tracks.tck --labels dseg.nii.gz \
                      --output matrix.json --measures measures.json

matrix = ritk.connectome.build_connectivity_matrix(
    parcellation, streamlines, assignment_radius=2.0
)

1. Which region does an endpoint belong to?

Terminal is exact, and exactly what discards streamlines terminating in white matter: tracking stops at the grey/white boundary while a cortical parcellation labels only grey matter, so rejection rates above half are ordinary. It is correct when the parcellation covers the whole brain including white matter, or when the streamlines were tracked to terminate inside grey matter by construction.

RadialSearch searches the endpoint’s own voxel first, so it can only add assignments that Terminal would have dropped — it never moves an endpoint that already sat inside a region. The radius trade-off is set out in the Parcellation chapter.

2. What does an edge weigh?

Pathway length. A streamline is reconstructed step by step, and the probability of surviving to its endpoint falls with the number of steps. Long-range connections are therefore systematically under-counted relative to short ones, by roughly the pathway length. InverseLength divides it back out.

Region size. A large region presents a large surface for streamlines to terminate on, so its edges are heavier for reasons of geometry rather than anatomy. This confounds comparison between regions of different size and between subjects whose regions differ in size. InverseNodeVolume divides by the summed node volumes.

Mean length is not a connectivity measure at all: an edge of weight 60 means the streamlines joining those regions averaged 60 mm, whatever their number. Useful paired with a count matrix, not as a substitute.

Accounting

Every streamline lands in exactly one bucket, and the buckets sum to the whole:

The three outcomes partition the input exactly: assigned + intra_region + unassigned = total. There is no bucket for a malformed streamline, because a Polyline cannot hold fewer than two finite points — its constructor refuses to build one.

This is kept alongside the matrix because a connectome is not interpretable without it. A matrix built from a tractogram of which four fifths were discarded is a different claim from one built from a tractogram of which a twentieth were, and nothing in the weights distinguishes the two.

let accounting = matrix.accounting();
println!("{:.1}% of streamlines produced an edge", 100.0 * accounting.assigned_fraction());

The human tractography and connectomics example runs this over a real diffusion dataset, and the atlas parcellation example produces the parcellation such a run consumes.

The connectivity matrix

ConnectivityMatrix stores a dense, fully symmetric \(n \times n\) weight matrix: both \((i,j)\) and \((j,i)\) carry the weight. A triangular layout would halve the memory — for a whole-brain atlas a few hundred kilobytes either way — at the cost of an index-ordering branch inside every measure that walks a row, and the graph algorithms walk rows constantly.

Self-connections sit on the diagonal. They are recorded, because a tractogram’s intra-region streamlines are real, but excluded from degree, density, and every path-based measure, where a self-loop is not a connection between two nodes.

Graph measures

let measures = matrix.measures();

println!("global efficiency  {:.3}", measures.global_efficiency());
println!("mean clustering    {:.3}", GraphMeasures::mean(measures.clustering()));
println!("communities        {}",    measures.communities().count());
println!("modularity         {:.3}", measures.communities().modularity());

They are computed together rather than one at a time because they share the all-pairs shortest-path solution, which dominates the cost.

The one convention everything rests on

A connectome edge weight is a measure of connection strength: a heavier edge means the regions are more closely linked. Every shortest-path measure needs the opposite — a distance, where a stronger link is a shorter step. The two are related by inversion,

\[ \ell(i, j) = 1 / w(i, j) \]

and the crate applies that inversion at the single point where the graph is converted for path finding. Getting it backwards is the classic error in weighted network analysis, and it does not fail loudly: it produces a path length that is largest exactly where the connection is strongest, so every derived measure inverts its meaning while remaining a plausible number.

Segregation

The binary clustering coefficient throws away everything the weights say: a triangle closed by three heavy edges counts the same as one closed by three negligible ones. The Onnela form replaces the triangle count with the sum of the geometric means of the three normalised weights. The geometric mean is the choice that makes the measure behave — it is zero if any of the three edges is absent, so a missing edge closes no triangle, and it reduces exactly to the binary coefficient when every present edge carries the maximum weight.

Community detection maximises modularity, the excess of within-community weight over a strength-preserving null model:

\[ Q = \frac{1}{2m} \sum_{ij} \left[ w_{ij} - \frac{k_i k_j}{2m} \right] \delta(c_i, c_j) \]

The implementation visits nodes in index order rather than at random. The published Louvain algorithm randomises, which makes the result vary between runs on the same input — unacceptable for a measure that has to be compared between subjects or reproduced from a paper. The cost is a partition that may differ from the best one a randomised search would find, and the modularity reported is always the modularity of the partition actually returned.

Modularity maximisation cannot resolve communities smaller than roughly \(\sqrt{2m}\) in total weight; below that, merging two genuinely separate communities increases \(Q\). This is a property of the objective, not of the search, so no amount of optimisation effort avoids it.

Integration

A real connectome is often not connected — an isolated region, or a parcellation whose regions the tractogram never reached. That makes the characteristic path length infinite, since some pair has no path. Two responses are in use and both are wrong in one direction: averaging over only the reachable pairs understates the cost of disconnection, and dropping the measure loses it entirely.

Global efficiency avoids the problem by averaging the reciprocal distance, where an unreachable pair contributes exactly zero. It is therefore reported unconditionally, and is the better-behaved summary. characteristic_path_length() is reported over reachable pairs only, and must be read together with reachable_pair_fraction(): a graph in fragments can show a short characteristic path precisely because the long paths do not exist.

Centrality and hubs

Betweenness identifies hubs whose position matters rather than whose degree is large: a node of modest degree bridging two otherwise separate modules carries enormous traffic, and degree does not see it. It is computed by Brandes’ algorithm, which never enumerates a path.

The rich-club coefficient asks whether the high-degree nodes are preferentially wired to each other. \(\Phi(k)\) rises with \(k\) in any graph, because high-degree nodes have more edges and so are likelier to be connected by chance alone — so a rising raw curve is not by itself evidence of a rich club. The measure is the ratio against a null model that keeps every degree and rewires everything else:

\[ \Phi_{\text{norm}}(k) = \Phi(k) / \langle \Phi_{\text{random}}(k) \rangle \]

normalised_rich_club computes it. The ensemble is built by repeated double-edge swaps, which preserve every node’s degree exactly — so the club membership at each threshold is identical in every sample, and only the edges among the club change. That is what makes the ratio a statement about wiring rather than about degree.

use ritk_connectome::measures::rich_club::{RandomisationConfig, normalised_rich_club};

let (levels, report) = normalised_rich_club(&matrix, RandomisationConfig::new(1000, 42))?;
for level in &levels {
    if let Some(ratio) = level.ratio {
        println!(
            "k={} ratio {ratio:.2} (random spread {:.2})",
            level.observed.degree, level.random_deviation
        );
    }
}
println!("rewiring acceptance {:.1}%", 100.0 * report.acceptance());

Ensemble size, swaps per edge, and the seed are explicit because each is a study-design choice with no defensible default. The seed is fixed rather than drawn, so a reported ratio can be reproduced.

Read the ratio with the acceptance fraction. A ratio near one means either that the wiring is unremarkable or that the degree sequence never allowed it to be otherwise — four nodes of degree five among eight degree-one leaves has exactly six slots for edges between them, every pair, so every graph with that sequence has a complete club and one is the true answer. A low acceptance fraction means the graph was too constrained to rewire and the ensemble never left where it started.

Each level also reports mean_weight, the weighted companion: a club can be topologically complete while its edges are individually weak.

Persistence

let json = matrix.to_json()?;
let restored = ConnectivityMatrix::from_json(&json)?;

References

• Rubinov, M. & Sporns, O. (2010). Complex network measures of brain connectivity: Uses and interpretations. NeuroImage 52(3):1059–1069.
• Jones, D. K., Knösche, T. R. & Turner, R. (2013). White matter integrity, fiber count, and other fallacies. NeuroImage 73:239–254.
• Blondel, V. D. et al. (2008). Fast unfolding of communities in large networks. J. Stat. Mech. 2008:P10008.
• Brandes, U. (2001). A faster algorithm for betweenness centrality. J. Math. Sociol. 25(2):163–177.
• Onnela, J.-P. et al. (2005). Intensity and coherence of motifs in weighted complex networks. Phys. Rev. E 71:065103.
• Colizza, V. et al. (2006). Detecting rich-club ordering in complex networks. Nature Physics 2:110–115.

.trk — DSI Studio / TrackVis Tractogram Format

ritk-trk reads and writes the DSI Studio / TrackVis .trk binary tractogram format. The format stores a fixed 1000-byte header followed by per-streamline data blocks, with an embedded voxel-to-RAS affine that maps on-disk voxel coordinates to physical space.

Format layout

A .trk file consists of:

1. A 1000-byte binary header starting with the 6-byte magic TRACK\0. All multi-byte integer fields are little-endian.
2. For each streamline declared by n_count:
   • A 4-byte n_points count (i32 LE)
   • n_points × (3 + n_scalars) float32 values (position x, y, z followed by per-point scalars)
   • n_properties float32 values for per-streamline properties

Coordinate system

Streamlines in memory are Gaia Polyline<f64> values in physical RAS+mm coordinates. On read, the header affine vox_to_ras (a 4×4 row-major matrix) is applied to every voxel coordinate triple to produce physical coordinates. On write, the inverse affine converts back to voxel space.

Reading

TrkTractogram::read accepts any impl Read and returns the parsed header, streamlines, per-point scalars, and per-streamline properties:

use std::fs::File;
use ritk_trk::TrkTractogram;

let mut file = File::open("tracks.trk")?;
let trk = TrkTractogram::read(&mut file)?;
println!("{} streamlines, {} scalars per point",
    trk.streamlines.len(),
    trk.header.n_scalars,
);

Validation

The reader checks:

• Magic bytes (TRACK\0)
• Header size equals 1000
• Streamline count is non-negative and bounded (≤ 100M)
• Per-streamline point count is non-negative and bounded (≤ 100K)
• Every coordinate is finite after affine application
• Every polyline is valid Gaia geometry

Writing

TrkTractogram::write serialises the header and streamlines, applying the inverse of vox_to_ras to convert physical coordinates back to voxel space:

let mut file = File::create("output.trk")?;
trk.write(&mut file)?;

Header fields

TrkHeader is the fixed 1000-byte header. Key fields:

The default header uses an identity vox_to_ras, RAS voxel order, version 2, and hdr_size 1000.

Per-point scalars and properties

When n_scalars > 0, per-point scalar values are stored inline after the position triplet. The scalar_name field is a space-separated list of names (e.g. "FA MD"). Scalars are returned as Vec<Box<[f32]>> — one flat array per streamline in n_points × n_scalars row-major order.

When n_properties > 0, per-streamline properties follow the streamline data block. Properties are returned as Vec<Box<[f32]>> — one array per streamline.

Error types

TrkError covers every failure mode:

Usage in the diffusion pipeline

ritk-tractography provides TractographyResult::to_trk, to_trk_header, and to_trk_with_scalars that produce TrkTractogram values directly from tractography output. See the Tractography chapter for the export API.

.tck — MRtrix3 Tractogram Format

ritk-tck reads and writes the MRtrix3 .tck tractogram format. The format uses a human-readable text header of key: value pairs terminated by the line END, followed by binary streamline data with NaN-delimited streamlines and an Inf barrier at end-of-file.

Format layout

A .tck file consists of:

1. A text header section:
   • The first line is mrtrix tracks
   • Key-value pairs in the form key: value, one per line
   • The line END terminates the header
2. Binary streamline data after END:
   • Streamline points encoded as 3-tuples of the declared datatype
   • Streamlines separated by a NaN triplet
   • End-of-file signalled by an Inf triplet (or physical EOF)

Coordinate system

Streamlines are stored in scanner-space millimetre coordinates — the native .tck coordinate system. The optional transform header entry maps voxel space to scanner space as a 4×4 row-major matrix, stored in TckHeader::transform for consumers that need to convert.

Datatypes

TckDatatype supports four binary encodings:

The default is Float32LE. The datatype is parsed from the datatype: header line and used for both reading and writing.

Reading

TckTractogram::read accepts any impl Read and returns the parsed header and streamlines:

use std::fs::File;
use ritk_tck::TckTractogram;

let file = File::open("tracks.tck")?;
let tck = TckTractogram::read(file)?;
println!("{} streamlines, datatype {:?}",
    tck.streamlines.len(),
    tck.header.datatype,
);

Validation

The reader checks:

• The first header line is mrtrix tracks
• The datatype header value is recognised
• Header lines parse as key: value pairs
• Every coordinate is finite (mixed NaN/Inf in a non-delimiter position is an error)
• Every polyline is valid Gaia geometry
• The transform value (if present) contains exactly 16 floats

The reader gracefully handles EOF mid-streamline (flushing any partial streamline) and skips empty header lines.

Writing

TckTractogram::write serialises the header and streamlines, using the declared datatype for all binary encoding:

let mut file = File::create("output.tck")?;
tck.write(&mut file)?;

The writer emits all well-known header keys in canonical order (datatype, mrtrix_version, file, comments, count, total_count, transform) before any remaining raw fields, then END, then the binary data.

Header fields

TckHeader carries parsed header metadata with typed accessors alongside a raw fields: HashMap<String, String> for unknown keys:

Weights sidecar

The .tck format does not store per-point scalars inline. MRtrix3 uses a separate weights file with the same binary layout: one scalar per point instead of three, with the same NaN delimiter / Inf barrier convention.

write_tck_weights writes a weights sidecar file:

use ritk_tck::{write_tck_weights, TckDatatype};

let fa_values: Vec<Box<[f32]>> = /* one per streamline, one f32 per point */;
let mut file = File::create("fa.tsf")?;
write_tck_weights(&fa_values, TckDatatype::Float32LE, &mut file)?;

read_tck_weights reads a weights file back into per-streamline scalar arrays:

use ritk_tck::read_tck_weights;

let file = File::open("fa.tsf")?;
let scalars = read_tck_weights(file)?;

Error types

TckError covers every failure mode:

Usage in the diffusion pipeline

ritk-tractography provides TractographyResult::to_tck and to_tck_header that produce TckTractogram values directly from tractography output, with optional mrtrix_version, comments, and transform parameters. Use write_tck_weights to export per-point scalars (FA, MD) alongside the tractogram. See the Tractography chapter.

.trx — Tractography Reference eXchange Format

ritk-trx reads and writes the TRX (Tractography Reference eXchange) format. TRX stores tractography data as a directory containing a JSON header and raw binary arrays, with explicit support for per-vertex data (DPV), per-streamline data (DPS), and hierarchical grouping.

Format layout

A TRX directory contains:

mytracks.trx/
├── header.json      — JSON metadata
├── positions.raw    — flat array of point coordinates (nb_points × 3)
├── offsets.raw      — cumulative point indices (nb_streamlines + 1)
└── dpv/             — per-vertex data arrays (optional)
    ├── FA.raw
    └── MD.raw

Positions and offsets

Points from all streamlines are flattened into a single positions.raw array of nb_points × 3 values. The offsets.raw array has nb_streamlines + 1 entries; streamline i occupies positions offsets[i]..offsets[i+1]. The sentinel offsets[nb_streamlines] must equal nb_points.

This contiguous layout eliminates per-streamline count headers and enables zero-copy slicing of the position array.

Data types

The dtype field in header.json declares the binary encoding for positions. Supported values are "float32" (4-byte LE f32) and "float64" (8-byte LE f64). Offset arrays are always uint64 LE. DPV/DPS arrays declare their own dtype in the header’s dpv/dps maps.

Coordinate system

Streamlines are stored in physical millimetre coordinates — the native TRX coordinate system. The optional reference field in the header carries NIfTI reference metadata: an affine (16 f64 values), voxel dimensions, and voxel sizes, enabling conversion to/from voxel space when a reference image is available.

Reading

TrkTractogram::read_dir reads from a directory path, expecting header.json, positions.raw, and offsets.raw. If the header declares DPV arrays, dpv/<name>.raw files are also read:

use ritk_trx::TrxTractogram;

let trx = TrxTractogram::read_dir("mytracks.trx")?;
println!("{} streamlines, {} total points",
    trx.streamlines.len(),
    trx.header.nb_points,
);

from_raw and from_raw_with_dpv read from in-memory buffers, enabling round-trip fidelity through encode/decode cycles.

Validation

The reader checks:

• positions.raw length matches nb_points × 3 for the declared dtype
• offsets.raw has exactly nb_streamlines + 1 entries
• The sentinel offset equals nb_points
• Every offset is monotonic and within bounds
• Every coordinate is finite
• Every polyline is valid Gaia geometry

Writing

TrkTractogram::write_dir creates the directory structure, writes header.json, positions.raw, offsets.raw, and any DPV data files:

trx.write_dir("output.trx")?;

to_raw returns the header and raw byte buffers for in-memory encoding:

let (header, positions, offsets, dpv_data) = trx.to_raw()?;

Header fields

TrxHeader is a Serde-serialisable JSON structure:

DPV and DPS arrays

TrxArrayDef declares one data array:

DPV data maps key names to raw byte buffers in TrxTractogram::dpv_data. Each buffer must match the length implied by nb_points × n_components × sizeof(dtype).

NIfTI reference

TrxReference carries optional NIfTI image metadata:

Groups

TrxGroup defines a named collection of streamline indices, enabling hierarchical organisation of tractography output (e.g. by bundle):

{
  "name": "CST_left",
  "indices": [0, 1, 2, 3, 4]
}

Error types

TrxError covers every failure mode:

Usage in the diffusion pipeline

ritk-tractography provides TractographyResult::to_trx and to_trx_with_dpv that produce TrxTractogram values with per-vertex data arrays (e.g. FA, MD). The caller populates the header’s dpv map with TrxArrayDef entries and provides the corresponding raw byte buffers. See the Tractography chapter.

Comparison with .trk and .tck

Classical Registration

The classical registration stack in ritk collects deterministic, non-ML alignment methods under one API. ritk-registration::classical exposes landmark-based rigid registration via Kabsch SVD, translation search with sealed metric types, mutual-information driven rigid and affine optimization, and temporal synchronization helpers for paired acquisition streams. These algorithms are deliberately CPU-centric and geometry-conscious: the code assumes that fixed and moving images have already crossed the ritk-image boundary with correct origin, spacing, and direction, then performs the actual solve in a form that is easy to audit and validate.

Rigid-search transforms use physical [z,y,x] millimetres. Native image world points use metadata [x,y,z]. Use rigid_physical_affine_to_native when a classical search result must enter a native resampler; it applies the axis permutation to both sides of the linear part and to translation. index_affine_to_physical is only for transforms whose input contract is explicitly index space.

Under Atlas, this chapter is where Coeus and Leto meet most visibly. Images arrive as Coeus-backed Image<f32, B, 3> values, then image_to_leto_volume converts contiguous voxel storage into Leto arrays for the classical engine. After optimization, leto_volume_to_image restores the result in the original physical frame so later filters, writers, and benchmarks still operate on the standard ritk image contract. That division keeps format handling, tensor execution, and classical numerics separated without fragmenting the user-facing workflow.

For local CT/MR soft-tissue refinement, MindSscFixedPrep can be evaluated by the same rigid-search objective closure. Prepare it once before search. Its selected centers and denominator remain fixed across poses, while moving patch values are sampled on demand. The optimizer bounds and search policy are unchanged.

Example Summary

Deformable registration has its own chapter because it solves for a spatially varying displacement field rather than a finite-dimensional rigid or affine transform. The deterministic Thirion Demons example is documented in Deformable Registration.

Example: Geometry Validation

Verify RITK’s NIfTI import and index-to-world mapping before registration.

Source

crates/ritk-registration/examples/geometry_check.rs

Description

The example accepts a CT and MR NIfTI path, prints each image’s shape, spacing, origin, direction, and several index-to-world samples. The output is the geometry contract to check before a registration metric or resampler combines the images.

Usage

cargo run -p ritk-registration --example geometry_check -- \
  /data/ct.nii.gz /data/mr.nii.gz

The example does not assume a particular dataset or filesystem layout. For the repository’s downloaded fixture, replace the two paths with the CT/MR files available in your local data directory.

Verification

• Loads both NIfTI files through ritk_io::format::nifti::native::NiftiReader
• Computes index-to-world for corner and interior voxels that fit each volume
• Makes axis order explicit: RITK stores the native volume as [z, y, x]
• Exercises the ritk-io::format::nifti spatial boundary without private paths

Example: Deep Learning Registration

Rigid registration using Coeus autodiff and optimization.

Source

crates/ritk-registration/examples/dl_registration.rs

Description

Demonstrates the Coeus-native registration pipeline: builds a differentiable registration graph, optimizes rigid transform parameters via gradient descent, and outputs aligned images. Uses coeus-autograd for autodiff and coeus-optim for the optimizer.

Usage

cargo run --example dl_registration -- <fixed> <moving> <output>

Verification

• Optimizes rigid transform (translation + rotation)
• Outputs aligned image with matching geometry
• Uses assert_coeus_matches_coeus differential test against reference

Example: Deep Learning Training

Train a Coeus-native registration model end-to-end.

Source

crates/ritk-registration/examples/dl_train.rs

Description

Full training loop for a Coeus-native registration network: builds the autodiff graph, runs forward/backward passes, and reports training metrics. Uses coeus-nn for the neural network layer and coeus-autograd for gradient computation.

Usage

cargo run --example dl_train -- <dataset_dir> <checkpoint_dir>

Verification

• Trains for specified number of epochs
• Reports loss convergence
• Saves model checkpoint in Coeus-native format

Temporal Signal Synchronization

Spatial registration answers “where is this anatomy?” Temporal synchronization answers “when did this event occur?” A dynamic MR sequence, PET time-activity curve, physiological trace, or gated reconstruction can be spatially aligned while its samples still refer to different moments.

RITK estimates one constant delay between two equally sampled finite signals. It does not resample unequal acquisition rates, model a delay that changes over time, or perform dynamic time warping. Those are different contracts.

Lag convention

Let (x_i) be the reference signal and (y_i) the moving signal. At integer lag (k), RITK evaluates only the index pairs that exist in both arrays:

[ r(k)= \frac{ \sum_i (x_i-\bar{x}k)(y{i+k}-\bar{y}_k) }{ \sqrt{ \sum_i(x_i-\bar{x}k)^2 \sum_i(y{i+k}-\bar{y}_k)^2 } }. ]

The overlap-specific means prevent a shorter edge overlap from being compared against full-signal means. Pearson normalization also makes the estimated lag invariant to a positive scale and offset applied to either signal.

A positive (k) means the moving signal is delayed. To align it, sample the moving signal at (i+k). The reported shift in seconds is

[ \Delta t = k,T_\mathrm{frame}. ]

This sign convention is visible in the worked example: the unaligned orange curve lags the reference, and the aligned curve is evaluated at reference_index + shift_frames.

From integer peak to fractional estimate

The integer lag with maximum normalized correlation is the discrete estimate. When correlations exist on both sides of an interior maximum, RITK fits a three-point parabola:

[ \delta = \frac{r(k-1)-r(k+1)} {2[r(k-1)-2r(k)+r(k+1)]}, \qquad \hat{k}=k+\delta. ]

The offset is accepted only for a finite concave peak and is bounded to the neighboring interval. A peak at the configured search boundary remains integer-valued. Parabolic refinement is an estimator; digitized correlation curves can bias it. Céspedes et al. compare parabolic and cosine fits and describe this limitation in Methods for Estimation of Subsample Time Delays of Digitized Echo Signals.

Correlation ties are deterministic: RITK chooses the smallest absolute lag, then the smaller signed lag. This prevents scheduling or iterator order from changing a result.

Validated configuration and typed failures

use leto::Array1;
use ritk_registration::{
    TemporalSync, TemporalSyncConfig, TemporalSyncStatus,
};

let reference = Array1::from_vec(
    [8],
    vec![0.0, 0.2, 0.9, 0.4, -0.3, -0.8, -0.2, 0.3],
)?;
let moving = Array1::from_vec(
    [8],
    vec![0.0, 0.0, 0.2, 0.9, 0.4, -0.3, -0.8, -0.2],
)?;

let config = TemporalSyncConfig::try_new(
    0.04, // seconds per sample
    3,    // search ±3 frames
    0.8,  // acceptance threshold
)?;
let result = TemporalSync::with_config(config)
    .synchronize(&reference, &moving)?;

assert!(result.shift_frames() > 0.0);
assert_eq!(result.status(), TemporalSyncStatus::Accepted);
Ok::<(), Box<dyn std::error::Error>>(())

TemporalSyncConfig::try_new rejects non-finite or non-positive frame spacing, a zero search range, and a minimum correlation outside [0, 1]. Synchronization then rejects:

• unequal signal lengths;
• fewer than three samples;
• the first NaN or infinite sample, with signal identity and index; and
• a zero-variance signal whose lag is not identifiable.

Two identical constant signals are not “perfectly synchronized.” They contain no timing feature, so every lag explains them equally well. RITK reports TemporalSyncError::UnidentifiableSignal instead of inventing a zero shift.

Acceptance is not numerical success

minimum_correlation classifies a measured estimate:

• TemporalSyncStatus::Accepted, or
• TemporalSyncStatus::BelowMinimumCorrelation.

A below-threshold result still exposes its shift and correlation. This lets a quality-control pipeline display or record the observation without confusing weak evidence with a numerical failure.

The correlation coefficient can be negative. The acceptance threshold is non-negative because a negatively correlated peak is not evidence that the signals share the same polarity and timing.

Residual metrics and units

After estimating (\hat{k}), RITK linearly interpolates (y) at (i+\hat{k}). Only coordinates inside the moving signal contribute:

[ \mathrm{RMS} = \sqrt{\frac{1}{N_\mathrm{overlap}} \sum_{i\in\mathrm{overlap}} [x_i-\tilde{y}(i+\hat{k})]^2}. ]

residual_rms and residual_max_abs are in signal amplitude units, not seconds. shift_seconds is the timing quantity. overlap_samples states the exact RMS denominator, so a large shift cannot make the error look smaller by dividing by samples that were never compared.

Search allocation and diagnostic profiles

synchronize scans lags once and retains only the current best peak and its neighbors. Search scratch is constant-size; it does not allocate full lag and correlation arrays. Call correlation_profile only when a graph or diagnostic record needs every lag:

let profile = synchronizer.correlation_profile(&reference, &moving)?;
for sample in &profile {
    if let Some(correlation) = sample.correlation() {
        println!("lag={} r={correlation}", sample.lag_frames());
    }
}
Ok::<(), ritk_registration::TemporalSyncError>(())

A profile sample is None when that particular overlap is locally constant. The normal peak selector ignores it. Full-signal zero variance remains a typed input error.

The direct method performs (O(NS)) work for (N) samples and search radius (S). It is appropriate when the physically plausible delay window is small. The committed temporal_sync benchmark covers 4,096 samples and a ±64-frame search without changing the workload between implementations.

Method reference

Xiao, Ding, and Hu describe lagged cross-correlation and maximum-lag selection for multimodal physiological signals in Section 2.3, Equation 1, and Algorithm 1. RITK pins its overlap, normalization, tie, sign, interpolation, and failure conventions above because a method name alone does not define those choices.

Example: Before and After Temporal Alignment

The same delayed synthetic acquisition drives every panel and every number in this figure. The example fails before writing the SVG if the estimated delay, acceptance status, correlation, overlap, or independently recomputed residuals contradict the displayed result.

Read the panels in order

1. Before — blue and orange peaks refer to the same waveform but occur at different sample indices. This is the error being corrected.
2. Search — normalized correlation is evaluated over the configured lag window. The green vertical line marks the fractional estimate; the orange dashed line is the minimum accepted correlation.
3. After — the delayed moving signal is sampled at reference_index + shift_frames. The dashed green trace should lie over the blue reference through the valid overlap.
4. Verify — the residual has its own expanded vertical scale. This prevents a visually overlapping pair from hiding systematic error.

The summary strip reports the known and estimated shifts separately, gives seconds only for the timing result, and labels RMS residuals as signal amplitude. “Before RMS” compares signals at the same uncorrected indices; “After RMS” uses only the valid interpolated overlap.

Source and command

Source: crates/ritk-registration/examples/book_temporal_sync.rs

cargo run -p ritk-registration --example book_temporal_sync -- \
  docs/book/figures/temporal_synchronization.svg

The example uses a 7.25-frame delay, 40 ms frame spacing, a ±20-frame search, and a minimum correlation of 0.85. Small deterministic acquisition noise keeps the example realistic without making the expected delay ambiguous.

Before rendering, it verifies that:

• fractional refinement is closer to 7.25 frames than integer quantization;
• peak correlation is at least 0.99 and the result is accepted;
• independently interpolated RMS, maximum residual, and overlap equal RITK’s public diagnostics; and
• temporal alignment reduces RMS residual by at least 85%.

Interpreting a real acquisition

A high peak is necessary but not sufficient. Inspect whether:

• the peak is isolated rather than one of many periodic aliases;
• it lies inside the search range rather than at its boundary;
• enough overlap remains after shifting;
• residual structure is noise-like rather than a repeated physiological feature; and
• frame spacing is the actual acquisition cadence.

If delay changes over the recording, one global shift is the wrong model. Split the acquisition only when the domain supports piecewise stationarity, or use a method designed for time-varying correspondence.

Multi-modal Registration

CT, MR, PET, dose, and derived maps do not share a simple intensity relationship. Multi-modal registration therefore needs a modality-robust objective, trustworthy physical metadata, and visual validation.

The minimum workflow is:

1. load both volumes and inspect their physical frames;
2. choose mutual information, MIND-SSC, LNCC, or NGF rather than assuming MSE is valid;
3. resample the moving image in the fixed frame;
4. compare identity and candidate transforms; and
5. retain a numerical change map and a labeled overlay.

Histogram estimation and overlap

MutualInformationMetric makes two decisions explicit. Fixed and moving modalities receive separate intensity ranges, so a CT Hounsfield-unit window is not interpreted as an MR signal range. HistogramEstimator::Discrete is the default when exact discrete-information identities matter. HistogramEstimator::MovingLinearPartialVolume keeps the fixed sample in one bin while distributing the resampled moving intensity between adjacent bins. The resulting objective changes continuously inside each moving-intensity bin, which is the useful direction for transform optimization.

Masked sample evaluation accepts borrowed slices and an optional fixed mask. For optimization, keep that mask and the accepted moving-field support fixed across candidate poses. Recomputing the support from each candidate lets a bad pose improve its score by cropping difficult anatomy. Report overlap separately and reject candidates that violate the application-derived support floor.

SpatiallyConditionedMutualInformationMetric additionally associates every sample with a fixed zero-based region label. It computes

[ H(F\mid R)=\sum_r p(r)H(F\mid r),\qquad H(F,M\mid R)=\sum_r p(r)H(F,M\mid r), ]

and applies the chosen NMI normalization to the conditional fixed, moving, and joint entropies. This preserves coarse location information that one global histogram discards. Region labels belong to the fixed frame and must not change with the candidate pose. The reusable workspace allocates its regional histograms once and clears them between evaluations. Toews and Wells derive this local-region formulation in section 3.2, equations 8–9 (https://doi.org/10.1007/978-3-642-02498-6_36).

For CT/MR soft tissue, a bounded two-stage search is preferable to a fitted weighted sum: NMI captures the multimodal basin, then a local structural metric refines anatomy. The structural half-range defaults to one terminal NMI cell; callers can select a larger nonzero number of cells while the original global rigid bounds continue to constrain every candidate. A global MIND-SSC or NGF search can lock onto remote anatomy, skull, or air boundaries.

Those bounds are refinement limits, not a substitute for orientation capture. For a larger initial mismatch, generate normalized-cross-correlation block matches in both directions, convert both sets to physical fixed-to-moving correspondences, and pass them to fit_symmetric_trimmed_rigid. The resulting full transform enters search_rigid_pose through RigidSearchAnchor; residual parameters remain subject to the same finite bounds. This follows the rigid subset of Modat et al.’s symmetric block-matching design (https://doi.org/10.1117/1.JMI.1.2.024003, sections 2.1–2.3).

MIND-SSC provides a complementary local structural objective when corresponding soft-tissue patches retain self-similarity despite modality-dependent intensity. Prepare fixed descriptors once at complete-support centers. For each candidate pose, only the six moving patch neighbourhoods at those centers are sampled. The default deterministic 8,192-center cap bounds memory and runtime; use caller-provided indices when an anatomical mask or validation protocol defines the exact population. The image field follows ITK’s half-voxel [-0.5,size-0.5) convention: support inside it uses replicate-edge trilinear interpolation, while support outside it is explicit zero background and stays in the fixed denominator.

The CT/MR fixture uses moving-linear partial-volume NMI plus a reference transform so the example is reproducible rather than optimizer-seed dependent. The coarse figure grid is intentionally not used to demonstrate spatial conditioning: its shallow axis leaves too little entropy per local region. The registered overlay is visibly different from identity because red and green fringes expose residual misalignment, while the MR change panel exposes where resampling actually changed the sampled values.

Example Summary

Example: Robust Rigid Capture Initializer

This example turns bidirectional block correspondences into a full fixed-to- moving rigid transform, then uses that transform as the zero residual of the bounded six-parameter search. It manufactures a known 20-degree rotation and translation, injects inconsistent correspondences below the 50% breakdown boundary, and checks the recovered transform against the analytical result.

The two correspondence slices have explicit directions. Forward entries store fixed-to-moving matches. Reverse entries store moving-to-fixed matches. RITK fits each direction independently with least-trimmed squares, inverts the reverse result, and combines the two fixed-to-moving transforms with a log-Euclidean mean. Conflicting unordered endpoint directions are removed from both schedules before fitting. Each direction retains its best half and the returned root-mean-square residual covers both retained sets.

RigidSearchAnchor validates that the initializer is finite, homogeneous, orthonormal, and proper. Search parameters then describe only the residual rotation and translation around that anchor; all existing global and local bounds remain active. A zero residual reproduces the supplied transform.

//! Fit a robust full rigid anchor and refine only its bounded residual pose.
#![expect(
    clippy::print_stdout,
    reason = "example reports the recovered pose and retained correspondence count"
)]

use anyhow::{ensure, Result};
use ritk_registration::{
    fit_symmetric_trimmed_rigid, search_rigid_pose, AffineTransform, FixedToMovingCorrespondence,
    MovingToFixedCorrespondence, RigidSearchAnchor, RigidSearchConfig,
};

fn transform_point(point: [f64; 3]) -> [f64; 3] {
    let angle = 20.0_f64.to_radians();
    let (sine, cosine) = angle.sin_cos();
    [
        cosine * point[0] - sine * point[1] + 4.0,
        sine * point[0] + cosine * point[1] - 3.0,
        point[2] + 2.0,
    ]
}

fn transform_residual(candidate: &AffineTransform, expected: &AffineTransform) -> f64 {
    candidate
        .as_array()
        .iter()
        .zip(expected.as_array())
        .map(|(actual, expected)| (actual - expected).powi(2))
        .sum::<f64>()
        .sqrt()
}

fn main() -> Result<()> {
    let fixed_points = [
        [-2.0, -1.0, 0.0],
        [3.0, -1.0, 1.0],
        [-1.0, 4.0, 2.0],
        [2.0, 3.0, -2.0],
        [-3.0, 1.0, 4.0],
        [4.0, 2.0, 3.0],
    ];
    let mut forward = Vec::new();
    let mut reverse = Vec::new();
    for fixed in fixed_points {
        let moving = transform_point(fixed);
        forward.push(FixedToMovingCorrespondence::try_new(fixed, moving)?);
        reverse.push(MovingToFixedCorrespondence::try_new(moving, fixed)?);
    }
    // Forty percent of each direction is deliberately inconsistent. The fit
    // retains the best half, so the correct 60% consensus remains identifiable.
    for index in 0_u32..4 {
        let offset = f64::from(index);
        let fixed = [50.0 + offset, -80.0, 20.0];
        let moving = [-70.0, 40.0 + offset, -30.0];
        forward.push(FixedToMovingCorrespondence::try_new(fixed, moving)?);
        reverse.push(MovingToFixedCorrespondence::try_new(moving, fixed)?);
    }

    let fitted = fit_symmetric_trimmed_rigid(&forward, &reverse)?;
    let expected = AffineTransform([
        20.0_f64.to_radians().cos(),
        -20.0_f64.to_radians().sin(),
        0.0,
        4.0,
        20.0_f64.to_radians().sin(),
        20.0_f64.to_radians().cos(),
        0.0,
        -3.0,
        0.0,
        0.0,
        1.0,
        2.0,
        0.0,
        0.0,
        0.0,
        1.0,
    ]);
    let gamma = 1_024.0 * f64::EPSILON / (1.0 - 1_024.0 * f64::EPSILON);
    ensure!(
        transform_residual(&fitted.transform, &expected) <= gamma * 4.0,
        "trimmed fit did not recover the manufactured rigid transform"
    );

    let anchor = RigidSearchAnchor::try_new(fitted.transform, [0.0; 3])?;
    let objective = |candidate: &AffineTransform| Ok(-transform_residual(candidate, &expected));
    let search = search_rigid_pose(
        anchor,
        RigidSearchConfig::try_new(4.0, 4.0, 0.5, 0.5, 128)?,
        objective,
        objective,
    )?;
    ensure!(
        transform_residual(&search.structural_transform, &expected) <= gamma * 4.0,
        "bounded residual search moved away from the recovered anchor"
    );

    println!(
        "retained {}/{} correspondences; inlier RMS {:.3e} mm; search saturated: {}",
        fitted.inlier_count,
        fitted.correspondence_count,
        fitted.inlier_rms_mm,
        search.capture_saturated || search.structural_saturated
    );
    Ok(())
}

Run the same source used by this page:

cargo run -p ritk-registration --example book_rigid_initializer

The example is synthetic by design: it verifies the estimator and transform composition independently of an image metric. A CT/MR application still has to generate correspondences from physical image grids, reject low-confidence matches, and validate the final pose with held-out landmarks or segmentations.

Example: CT/MR Mutual-Information Registration

Visual validation of rigid CT-to-MR registration using the repository’s RIRE Patient 001 fixture. The dataset’s fiducial transform is compared with identity sampling and rendered through RITK’s native physical resampler. The figure is explicitly labeled so the registration change can be read without guessing which panel is pre-registration or post-registration.

Source

crates/ritk-registration/examples/book_registration.rs

Description

The example first samples both modalities onto a common coarse physical grid, evaluates moving-linear partial-volume NMI before and after applying the RIRE fiducial transform, prepares packed MIND-SSC descriptors on at most 8,192 deterministically selected complete-support CT centers within a conventional [-200, 200] HU soft-tissue mask, and scores MR directly at those centers before and after the transform. It then resamples the original MR volume onto the full CT grid. Because this fixture supplies a geometric registration standard, the example uses that transform for the reproducible output rather than publishing a blind optimizer result that can stop in a wrong MI basin. The generated SVG contains four labeled panels:

1. CT windowed to [-1000, 1000] HU.
2. Identity CT/MR overlay.
3. Registered CT/MR overlay using the RIRE fiducial transform through the native resampler.
4. Absolute MR intensity change between the identity and registered samples.

In the overlay panels, red is CT and green is MR. Correctly coincident anatomy appears yellow; red or green fringes identify residual misalignment. Their subtitles report fixed-denominator MIND-SSC similarity. The score uses the classical physical [z,y,x] affine through RITK’s axis-convention bridge, then samples native [x,y,z] image geometry; the example does not reinterpret an already-physical affine as an index transform. The fourth panel is not a second registration result: it is a data-derived diagnostic showing where the transform changes the sampled MR values. Its subtitle reports the maximum and mean absolute change for the displayed slice.

Usage

cargo run -p ritk-registration --example book_registration -- \
  docs/book/figures/ct_mri_registration.svg

Verification

• Loads the in-tree CT/MR pair from test_data/registration/rire/.
• Requires normalized mutual information to improve from identity to the dataset transform.
• Reports fixed-denominator MIND-SSC at identity and at the dataset transform, together with its deterministic center count and exact retained heap payload. The metric is an independent diagnostic here; this fixture does not assert that one raw CT soft-tissue mask validates MIND-SSC across modalities.
• Produces the registered MR on the full CT physical grid.
• Uses one parsed RIRE transform and one native resampling call for the registered panel, avoiding a duplicate reference panel.
• Reports the maximum and mean absolute MR change between identity and registered sampling.

The source volumes and reference transform are distributed with the RIRE fixture. See the RIRE fixture provenance and reference transform for provenance, license, geometry, and the reference transform.

Deformable Registration

RITK’s classic Thirion Demons implementation estimates a dense displacement field from the fixed-image gradient and the moving-to-fixed intensity difference. Each iteration applies a bounded force update, optional fluid regularization, optional diffusion regularization, and a forward warp. The result carries the warped moving image, the three displacement components, the final MSE, and the performed iteration count.

The public API is intentionally slice-oriented:

use ritk_registration::demons::{DemonsConfig, ThirionDemonsRegistration};

let registration = ThirionDemonsRegistration::new(DemonsConfig::default());
let result = registration.register(
    &fixed,
    &moving,
    [depth, height, width],
    [spacing_z, spacing_y, spacing_x],
)?;

The dims order is [z, y, x], matching RITK’s contiguous volume contract. The returned disp_x uses the forward-warp convention implemented by the engine. A registration result is not accepted merely because the call returns Ok: deformable registration examples must compare the output against an identity baseline and inspect the displacement field.

Dataset-backed translation example

The runnable book example loads the RIRE Patient 001 T1 MR volume, extracts a bounded central crop, normalizes it, and creates the moving image with a known three-voxel translation. Using one MR acquisition on both sides is important: classic Demons minimizes an intensity-difference objective, so applying it directly to raw CT and MR intensities would be mathematically invalid.

The identity overlay below has separated red/green anatomy and yellow overlap. After Demons, those fringes approach coincidence and the measured mean-squared error falls from 0.049194 to 0.005062. The final panel renders the recovered x-displacement field on a signed scale; its positive interior direction agrees with RITK’s forward-warp convention.

Regenerate the figure from the real dataset and algorithm:

cargo run -p ritk-registration --example book_demons_registration -- \
  docs/book/figures/demons_registration.svg

The example rejects non-finite displacement values and fails unless the measured post-registration MSE is lower than the identity baseline. Its 6 × 48 × 48 crop and 35 iterations preserve the anatomical registration signal while keeping the already-built example runtime below one second.

For cross-modality registration, use CT/MR Mutual-Information Registration. That example uses the RIRE fiducial transform as its geometric reference and a mutual-information metric that does not assume matching CT and MR intensities.

Validation and Benchmarking

Validation and performance are separate acceptance axes. A faster run is not an accuracy result, and a lower metric is not proof of correct geometry.

Registration validation should combine:

• shape, spacing, origin, and direction checks;
• metric values before and after;
• overlap measures when labels exist;
• convergence state and iteration budget; and
• a labeled overlay plus an input-to-output change map.

Benchmarking should measure loading, preprocessing, metric evaluation, resampling, and optimizer loops separately. This makes a regression attributable to an algorithm rather than a hidden boundary copy.

Example Summary

The CT/MR example records both the identity-to-registered metric change and the visible resampling change so those claims remain distinguishable.

Example: Registration Validation

Validation is a separate stage from optimization. A registration can improve a metric while violating geometry, producing an empty resample, or moving a label boundary in the wrong direction. A useful report combines:

• shape and physical-frame checks;
• metric values before and after;
• overlap measures when labels exist; and
• convergence state and iteration budget.

The CT/MR registration example demonstrates the visual part of this contract: it renders identity and registered overlays with red/green fringes and a data-derived MR resampling-change map.

For a label-space check, use the statistics facade after resampling both label maps onto one grid:

let dice = ritk_statistics::dice_coefficient(&fixed_labels, &moving_labels)?;
let hausdorff = ritk_statistics::hausdorff_distance(&fixed_labels, &moving_labels)?;

The exact statistics signature depends on the label representation, so the validation layer must keep the conversion at the boundary and report the input shapes with the metric values. Never accept only an is_ok result: record the value, units, reference frame, and threshold used for the decision.

Source and verification

The visual source is crates/ritk-registration/examples/book_registration.rs.

cargo run -p ritk-registration --example book_registration -- \
  docs/book/figures/ct_mri_registration.svg

The registration example requires normalized mutual information to improve from identity to the dataset transform, preserves the full CT grid, and reports maximum and mean absolute MR resampling change. The registration and statistics test suites provide the numerical oracles; the figures are checked for correct labels, shared display conventions, and visible pre/post change.

Benchmarking

Benchmarking keeps algorithm cost separate from file I/O, image extraction, metadata reconstruction, and dispatch. Many filters operate on flat buffers, so a useful benchmark compares the public image call with the reusable buffer core and records where any boundary copy occurs.

Use release builds, fixed inputs, and a benchmark closure that returns the computed result. Record median and confidence intervals. A benchmark that only checks that a call succeeds cannot detect a dead computation or wrong value.

Registration measurements

Registration uses the same separation. Time image loading and preprocessing outside the timed region, then measure metric evaluation, resampling, and optimizer steps independently. Keep fixed and moving images and the optimizer configuration in the benchmark input. Report accuracy and runtime separately.

The repository’s runnable benchmark lane is the source of truth for command names and budgets. Build it in release mode and run the smallest case before collecting a full baseline.

Example Summary

Example: Gradient Recursive Gaussian Benchmark

Micro-benchmark: GradientRecursiveGaussian buffer-level vs per-Image passes.

Source

crates/ritk-filter/examples/bench_gradient_rg.rs

Description

Compares the optimized gradient_recursive_gaussian_components (one tensor extraction, component-wise smoothing) against the per-Image filter path. Both exercise the identical substrate-agnostic host core via ritk-image::Image boundary. No Burn tensor is constructed.

Usage

cargo run --release --example bench_gradient_rg

Verification

• Reports timing for both buffer-level and per-Image paths
• Uses assert_coeus_matches_coeus to verify bitwise-identical results
• Demonstrates the native host core performance vs wrapper overhead

Backend Dispatch

Backend dispatch in ritk is designed to be explicit in the type system and cheap at runtime. ritk-image re-exports Atlas backends such as SequentialBackend and MoiraiBackend, while ritk-io binds readers and writers to a chosen backend through types like PngReader<B>, MetaImageWriter<B>, and VtkReader<B>. Deeper in the stack, many filters expose a native host-core path and then reconstruct the same ritk-image::Image boundary, so callers get one public API while the compiler monomorphizes concrete backend implementations underneath.

That matters for Atlas because Coeus is not just a storage detail: it is the execution substrate that lets the same algorithm run deterministically on the sequential backend or in parallel on Moirai. The chapter therefore covers where dispatch is compile-time, where a host extraction is unavoidable, and how the unified read_image_native and write_image_native helpers keep format inference separate from compute choice.

Example Summary

Zero-copy I/O

Zero-copy in ritk means minimizing ownership churn at the boundary, not pretending that every file decode is literally allocation-free. The ritk-image::Image contract keeps voxel storage in a Coeus tensor, exposes checked accessors such as data_slice() for contiguous CPU views, and lets writers extract host data only when a format encoder truly needs it. Supporting crates follow the same pattern: ritk-morphology::StructuringElement::offsets() returns a borrowed slice, and many filter kernels operate on flat buffers so they do not repeatedly rebuild tensor structures in the hot path.

Atlas integration sharpens the distinction between cheap views and necessary copies. Staying inside the Coeus-backed image pipeline often means metadata moves without touching the voxels, but some seams necessarily materialize data: JPEG encode or decode crosses a compression boundary, and classical registration widens f32 image values into f64 Leto volumes before optimization. This chapter is about recognizing those intentional costs, keeping them localized, and preserving a single image object through as much of the workflow as possible.

Example Summary

Coeus Nonlinear Least-Squares Solver

coeus-optim owns the Levenberg-Marquardt (damped Gauss-Newton) solver used by the nonlinear diffusion models in ADR 0036 decision 2. DKI, NODDI, and IVIM fit their parameters through this solver; RITK adds no local nonlinear optimiser.

Why a separate solver contract?

The gradient-descent optimisers in coeus-optim (SGD, Adam, AdamW, RMSProp, AdaGrad) step on gradients already accumulated into parameters, which suits network training. A least-squares solver instead re-evaluates the model at trial points to decide whether a step is acceptable, and exploits the Gauss-Newton curvature approximation \(J^{\!T}\!J \approx H\) that a bare gradient does not expose. The two are different contracts, not two spellings of one.

LeastSquaresProblem trait

LeastSquaresProblem<T> is the contract between a model and the solver. It requires:

The Jacobian is row-major: entry \((i, j)\) at index \(i \cdot n + j\) is \(\partial r_i / \partial p_j\).

The solver requires \(m \ge n\) — fewer residuals than parameters makes \(J^{\!T}\!J\) singular by construction, and damping would mask that rather than solve it.

Domain errors

ProblemError::Domain signals that the model is undefined at the trial parameters (e.g. a negative diffusivity under a square root). The solver treats this as a rejected step — it increases damping to pull the next trial back toward the last accepted point — rather than a failure. ProblemError::Evaluation is for unrecoverable failures.

Levenberg-Marquardt algorithm

Each iteration solves the damped normal equations:

\[ (J^{\!T}\!J + \lambda \cdot \operatorname{diag}(J^{\!T}\!J)) \delta = -J^{\!T}r \]

Damping is scaled by \(\operatorname{diag}(J^{\!T}\!J)\) rather than the identity (Marquardt’s modification), making the step invariant to rescaling of individual parameters — a diffusion model mixing diffusivities near \(10^{-3}\) with signal amplitudes near \(10^3\) is exactly the badly-scaled case that motivates it.

The algorithm proceeds as follows:

1. Evaluate residuals and Jacobian at the current parameters.
2. Check gradient convergence: \(\|J^{\!T}r\|_\infty \le\) gradient_tolerance.
3. Build \(J^{\!T}\!J\) and solve the damped system via Cholesky decomposition (the matrix is symmetric positive-definite for \(\lambda > 0\)).
4. If Cholesky fails, increase \(\lambda\) and retry (up to a cap).
5. Accept the step when cost decreases; otherwise reject and increase \(\lambda\). On acceptance, decrease \(\lambda\).
6. Check step-tolerance and cost-tolerance convergence.

Usage

use coeus_optim::{
    LeastSquaresProblem, LeastSquaresReport, LevenbergMarquardtConfig,
    ProblemError, levenberg_marquardt,
};

struct MyModel { /* data */ }

impl LeastSquaresProblem<f64> for MyModel {
    fn residual_count(&self) -> usize { /* ... */ }
    fn parameter_count(&self) -> usize { /* ... */ }
    fn residuals(&self, p: &[f64], r: &mut [f64]) -> Result<(), ProblemError> { /* ... */ }
    fn jacobian(&self, p: &[f64], j: &mut [f64]) -> Result<(), ProblemError> { /* ... */ }
}

let problem = MyModel { /* ... */ };
let initial = vec![0.0; problem.parameter_count()];
let config = LevenbergMarquardtConfig::default();
let report: LeastSquaresReport<f64> = levenberg_marquardt(
    &problem, &initial, &config,
)?;

Configuration

LevenbergMarquardtConfig<T> carries seven tuning parameters:

The tolerance defaults are \(\sqrt{\varepsilon}\) for the working scalar type — the standard choice for a first-order criterion in floating point: a residual computed to relative accuracy \(\varepsilon\) cannot certify stationarity below roughly \(\sqrt{\varepsilon}\).

Report and termination

LeastSquaresReport<T> carries the result:

Termination has four variants:

termination.is_converged() is true for all variants except IterationLimit.

Error types

SolverError covers solver-level failures:

ProblemError::Domain is handled internally by the solver (treated as a rejected step) and never propagates to the caller.

First-order optimisers

coeus-optim also provides five gradient-descent optimisers — SGD, Adam, AdamW, RMSProp, AdaGrad — used by the Coeus network training stack. These implement the Optimizer<T, B> trait and are not used by the diffusion pipeline, which requires the residual/curvature contract of LeastSquaresProblem.

The DKI and NODDI models in ritk-diffusion each implement LeastSquaresProblem<f64> and call levenberg_marquardt directly. See the Diffusion Models chapter for the model-specific usage.

Leto Linear Algebra Operations

leto-ops owns the dense linear algebra, matrix decompositions, and constrained solvers used by the diffusion pipeline. DTI and ODF estimation route through solve_least_squares; CSD routes through nnls; DKI and NODDI’s damped normal equations route through cholesky_solve. RITK adds no local solve, no local decomposition, and no local eigendecomposition.

Least-squares solve

solve_least_squares solves the overdetermined linear system \(\min \|Ax - b\|_2\) via QR decomposition. It is the workhorse for linear model fitting:

use leto_ops::solve_least_squares;
use leto::{Array1, Array2};

let solution = solve_least_squares(&design.view(), &rhs.view())?;

The DTI log-linear fit assembles one design matrix row per DWI volume and calls this function once per voxel. The ODF spherical harmonic estimator assembles an augmented design matrix (with Laplace–Beltrami penalty rows) and calls the same function.

NNLS — Non-Negative Least Squares

nnls implements the Lawson–Hanson active-set algorithm: it iteratively adds variables to an active set, solves the unconstrained subproblem on the active set, and removes variables with negative coefficients. The result is a solution vector \(x \ge 0\) that minimises \(\|Ax - b\|_2\).

Configuration

NnlsConfig controls convergence:

Result

NnlsResult carries:

Usage

use leto_ops::{nnls, NnlsConfig};

let config = NnlsConfig::default();
let result = nnls(&design.view(), &rhs.view(), &config)?;
assert!(result.solution.iter().all(|&x| x >= 0.0));

CSD uses NNLS to enforce the non-negativity constraint on the fibre orientation distribution. The deconvolution matrix \(B_{\text{resp}}\) is assembled by rescaling the Apollo SH design matrix with the rotational harmonics of the response function. See the Diffusion Models chapter.

Cholesky decomposition and solve

cholesky_decompose factorises a symmetric positive-definite matrix as \(A = LL^{\!T}\). cholesky_solve solves \(Ax = b\) using the Cholesky factorisation:

use leto_ops::cholesky_solve;

let step = cholesky_solve(&damped_normal.view(), &rhs.view())?;

This is the inner solve inside each Levenberg-Marquardt iteration. The damped matrix \(J^{\!T}\!J + \lambda \cdot \operatorname{diag}(J^{\!T}\!J)\) is symmetric positive-definite for \(\lambda > 0\), so Cholesky is the correct factorisation. When it fails, damping is increased and the solve is retried — see the Coeus Solver chapter.

cholesky_inv and cholesky_det provide the inverse and determinant from the Cholesky factor.

QR decomposition

qr_decompose factorises \(A = QR\) into an orthogonal matrix \(Q\) and an upper-triangular matrix \(R\). solve_least_squares uses QR internally. col_piv_qr provides column-pivoted QR for rank-deficient problems:

use leto_ops::qr_decompose;

let (q, r) = qr_decompose(&matrix.view())?;

Pseudoinverse

pinv computes the Moore–Penrose pseudoinverse \(A^+\) via SVD:

use leto_ops::pinv;

let a_pinv = pinv(&matrix.view())?;

Eigendecomposition

symmetric_eigen_jacobi computes all eigenvalues and eigenvectors of a real symmetric matrix via the Jacobi algorithm. eigenvalues computes eigenvalues of a general square matrix:

use leto_ops::symmetric_eigen_jacobi;

let decomposition = symmetric_eigen_jacobi(&matrix.view())?;
// decomposition.eigenvalues, decomposition.eigenvectors

RITK’s DTI estimator uses its own analytic 3×3 symmetric eigendecomposition (via the cubic characteristic polynomial) because the 3×3 case admits a closed form. Larger problems would route through these Leto functions.

Matrix norms and properties

Singular value decomposition

svd_decompose computes the full SVD \(A = U \Sigma V^{\!T}\). Variants provide rank-revealing, tolerance-controlled, and bidiagonal-reduction paths:

use leto_ops::svd_decompose;

let svd = svd_decompose(&matrix.view())?;
// svd.u, svd.singular_values, svd.vt

Iterative solvers

leto-ops no longer ships Krylov recurrences: stack ownership of the iterative-solver layer moved to Athena (Atlas ADR 0033), which provides CG, BiCGSTAB, restarted GMRES, and damped LSQR over backend-generic operators. The diffusion pipeline documented here needs none of them — every solve routes through the direct dense layer in the table above (solve_least_squares, nnls, cholesky_solve).

Relationship to the diffusion pipeline

Every solve routes through Leto; RITK owns only the design-matrix assembly and post-processing.

Apollo Spherical Harmonic Basis

apollo-sht owns the spherical harmonic basis used by the ODF and CSD models in ritk-diffusion. The real, even-order, antipodally symmetric basis is an upstream provider capability per ADR 0036 decision 2. RITK assembles the design matrix and calls apollo_sht for basis evaluation; it implements no local spherical-harmonic path, no local associated-Legendre recurrence, and no local SH normalization.

Why a separate SH provider?

The complex spherical harmonic transform in apollo-sht predates the diffusion pipeline — it was built for spectral analysis on the sphere (Gauss-Legendre grids, forward/inverse transforms, Parseval verification). The real even-order basis is a specialization of that general machinery: same associated-Legendre recurrence, same Condon-Shortley convention, same normalization constants, different output convention (real, antipodally symmetric, only even degrees). Placing the real basis in apollo-sht alongside the complex SHT means:

• The recurrence kernel and normalization are verified once, not duplicated.
• The Gauss-Legendre quadrature that certifies SHT exactness also certifies the real basis’s orthonormality on the sphere.
• Every consumer (ODF, CSD, and future FOD models) shares one basis instance.

Real Basis Convention

The real, orthonormal convention follows the MRtrix3 formulation. From the complex orthonormal harmonics \(Y_l^m(\theta, \phi)\) with Condon-Shortley phase:

\[ \begin{aligned} R_l^0(\theta, \phi) &= \operatorname{Re}\big(Y_l^0(\theta, \phi)\big) \\ R_l^m(\theta, \phi) &= \sqrt{2}\, \operatorname{Re}\big(Y_l^m(\theta, \phi)\big) \quad (m > 0) \\ R_l^{-m}(\theta, \phi) &= \sqrt{2}\, \operatorname{Im}\big(Y_l^m(\theta, \phi)\big) \quad (m > 0) \end{aligned} \]

The basis is orthonormal: \(\int_{S^2} R_l^m R_{l’}^{m’}\,d\Omega = \delta_{l,l’}\delta_{m,m’}\). Only even degrees \(l = 0, 2, 4, \ldots, l_{\max}\) are included, because the diffusion signal is antipodally symmetric (\(S(\mathbf{g}) = S(-\mathbf{g})\)) and odd-degree harmonics vanish.

The coefficient count is \((l_{\max} + 1)(l_{\max} + 2) / 2\) — for \(l_{\max} = 8\) that is 45 coefficients.

RealSphericalHarmonicBasis

RealSphericalHarmonicBasis is the main type consumed by the diffusion pipeline. It stores degree-major, order-minor metadata over even degrees: l=0: m=0, l=2: m=-2..=2, l=4: m=-4..=4, and so on.

Construction

use apollo_sht::RealSphericalHarmonicBasis;

// Create a basis for even degrees 0, 2, 4, 6, 8
let basis = RealSphericalHarmonicBasis::new(8)?;

Construction validates:

• l_max is even (odd degrees have no place in an antipodally symmetric basis).
• l_max >= 2 (trivial basis not useful for orientation estimation).
• l_max <= 85 (maximum stable degree in binary64 — the normalization product \((2l)!\) stays finite until degree 85).

Accessors

Pointwise Evaluation

real_spherical_harmonic(degree, order, theta, phi) evaluates one basis function at spherical coordinates \(\theta \in [0, \pi]\), \(\phi \in [0, 2\pi)\):

use apollo_sht::real_spherical_harmonic;

// R_2^2 at the equator, φ = 0
let value = real_spherical_harmonic(2, 2, std::f64::consts::FRAC_PI_2, 0.0)?;

For Cartesian directions, use the convenience method on the basis:

// Evaluate all 45 basis functions at a unit direction
let direction = [1.0_f64 / 3.0_f64.sqrt(); 3];
let row = basis.evaluate_at_direction(&direction)?;
// row.len() == 45

evaluate_at_direction validates that the input vector has unit-length squared norm (within \(32\varepsilon\)) and converts to spherical coordinates.

To evaluate all basis functions at spherical angles:

let row = basis.evaluate(theta, phi)?;

Design Matrix

design_matrix(directions) builds the full \(N \times K\) design matrix \(B\) where \(B_{i,k} = R_{l(k)}^{m(k)}(\theta_i, \phi_i)\) for \(N\) gradient directions and \(K\) coefficients:

let directions: Vec<[f64; 3]> = scheme.directions.iter()
    .map(|d| d.unit_direction().to_array())
    .collect();
let design: leto::Array2<f64> = basis.design_matrix(&directions)?;
// design.shape() == [N, K]

The operation costs \(O(N \cdot K \cdot l_{\max})\). No temporary row vectors are allocated — each basis value is pushed directly into the row-major buffer. The result is a Leto Array2<f64> ready for linear solves.

Error Types

RealShError is a comprehensive, non-exhaustive error enum:

All errors are typed — no panics, no .unwrap() in the public API.

Underlying Complex Machinery

The real basis sits on top of apollo-sht’s complex spherical harmonic infrastructure:

The complex transform machinery is verified by:

• Theorem 1 — Bonnet’s recurrence converges for all \(|x| \leq 1\) (numerically stable upward recurrence at fixed order).
• Theorem 2 — The \(n\)-point GL rule integrates polynomials of degree \(\leq 2n-1\) exactly (Golub & Welsch, 1969).
• Theorem 3 — The complex harmonics are orthonormal on \(S^2\) (Driscoll & Healy, 1994).
• Theorem 4 — Forward-inverse SHT round-trips exactly for band-limited fields when the grid satisfies the Shannon-Nyquist condition.
• Parseval verification — Energy \(\sum |a_{lm}|^2\) is conserved through the transform (proptest-verified with random band-limited fields).

Relationship to the Diffusion Pipeline

The real SH basis is consumed by two models in ritk-diffusion:

ODF — Orientation Distribution Function

The analytical Q-ball ODF estimator:

1. Builds the design matrix \(B\) via basis.design_matrix(&scheme.directions).
2. Appends Laplace-Beltrami penalty rows: \(\lambda \cdot l^2(l+1)^2\) on the diagonal for each \((l, m)\).
3. Solves the augmented system via leto_ops::solve_least_squares.
4. Evaluates the ODF on a display grid: \(\Psi(\theta, \phi) = \sum c_{lm} R_l^m(\theta, \phi)\).

use apollo_sht::RealSphericalHarmonicBasis;
use leto_ops::solve_least_squares;

let basis = RealSphericalHarmonicBasis::new(config.l_max())?;
let design = basis.design_matrix(&directions)?;
// ... add Laplace-Beltrami rows ...
let coeffs = solve_least_squares(&augmented_design.view(), &augmented_rhs.view())?;

CSD — Constrained Spherical Deconvolution

The fibre ODF estimator:

1. Builds the SH design matrix as above.
2. Computes the response function’s rotational harmonics \(r_l = 2\pi \int_{-1}^{1} R(\cos\alpha)\,P_l(\cos\alpha)\,d(\cos\alpha)\).
3. Forms the convolution matrix \(B_{\text{resp}}\) by scaling each SH column by the corresponding \(r_l\).
4. Solves via leto_ops::nnls to enforce non-negative fODF coefficients.
5. Extracts peaks from the reconstructed fODF for tractography.

let basis = RealSphericalHarmonicBasis::new(config.l_max())?;
let design = basis.design_matrix(&directions)?;
let conv_matrix = build_deconvolution_matrix(&design, &response, &basis)?;
let result = nnls(&conv_matrix.view(), &signal.view(), &nnls_config)?;
// result.solution contains non-negative fODF coefficients

The design matrix is the same object in both models — only the subsequent linear system changes (Tikhonov-regularized least-squares for ODF, non-negative least-squares for CSD).

Boundary

Under ADR 0036 decision 2, the real SH basis and its design matrix belong to Apollo, never to RITK. A RITK-local associated-Legendre recurrence or SH normalization path is a boundary violation. The chapter in ritk_diffusion.md documents the consumption side; this chapter documents the provider contract that RITK calls into.

References

• Driscoll, J.R. & Healy, D.M. (1994). “Computing Fourier transforms and convolutions on the 2-sphere.” Advances in Applied Mathematics, 15(2), 202-250.
• Golub, G.H. & Welsch, J.H. (1969). “Calculation of Gauss quadrature rules.” Mathematics of Computation, 23(106), 221-230.
• MRtrix3 spherical harmonic formulation: https://mrtrix.readthedocs.io/en/dev/concepts/spherical_harmonics.html
• Tournier, J.D., Calamante, F., & Connelly, A. (2007). “Robust determination of the fibre orientation distribution in diffusion MRI.” NeuroImage, 35(4), 1459-1472.

Gaia Polyline Geometry

Gaia owns the Polyline type — the canonical open curve that ritk-tractography produces as streamline output and ritk-connectome consumes as the input streamlines for connectivity matrix construction. This is an upstream provider contract per ADR 0036 decision 2: streamlines are Gaia polyline geometry; RITK owns the integration policy that produces them and does not define a local curve or polyline type.

Why Gaia?

Gaia is the Atlas geometry owner — it already carries the Point3<T> type (via Leto), exact geometric predicates, axis-aligned bounding boxes, and mesh builders for CFD milifluidic simulation. Polyline geometry sits naturally in the same domain: ordered point sequences in three dimensions, with validation, measurement, and spatial queries. Placing Polyline here means:

• Tractography output is typed in a geometry crate that downstream consumers (visualisation, meshing, surface reconstruction) already depend on.
• The same Aabb<T> spatial query is available for streamlines and meshes.
• Format-crate readers (ritk-trk, ritk-tck, ritk-trx) return Vec<Polyline<f64>> — the interchange format boundary is clean.

A RITK-local line or curve type is a boundary violation under ADR 0036.

Polyline Construction

Polyline::new(points) accepts a Vec<Point3<T>> and validates:

• At least two points. A single-point polyline is geometrically a point, not a curve; a zero-point polyline is meaningless.
• Every coordinate is finite. NaN or infinity in any vertex coordinate is rejected with the offending index.

use gaia::Polyline;
use leto::geometry::Point3;

let line = Polyline::<f64>::new(vec![
    Point3::new(0.0, 0.0, 0.0),
    Point3::new(3.0, 4.0, 0.0),
])?;
assert_eq!(line.segment_count(), 1);
assert!((line.arc_length() - 5.0).abs() < f64::EPSILON);

The type parameter T: Scalar defaults to f64. Storage is frozen as a boxed slice — no retained overcapacity — so memory never exceeds the vertex count once validated.

Error Types

PolylineError has two non-exhaustive variants:

Accessors

let streamline: Polyline<f64> = /* from tractography */;
for (from, to) in streamline.segments() {
    let length = (to - from).norm();
    // ... per-segment logic
}

Measurements

Arc Length

arc_length() computes the sum of Euclidean distances between consecutive vertices — the physical length of the streamline in millimetres:

let length: f64 = streamline.arc_length();

This is the quantity used by tractography termination (minimum length check) and connectome edge-weighting.

Axis-Aligned Bounding Box

aabb() returns the Aabb<T> enclosing all vertices:

let bounds = streamline.aabb();
// bounds.min / bounds.max — Point3<T>
// bounds.extents() — (width, height, depth)
// bounds.center() — geometric centre
// bounds.volume() — enclosed box volume

The Gaia Aabb<T> type also supports point containment, AABB intersection, and union operations. Tractography uses it for spatial indexing; connectome verification uses it to check streamlines against the parcellation volume.

Streamline Output in ritk-tractography

ritk-tractography defines a Streamline struct that wraps a Polyline<f64> plus per-point scalars:

pub struct Streamline {
    pub geometry: Polyline<f64>,
    pub scalars: Option<Box<[f64]>>,
}

The integration step converts accumulated integration points into a Polyline through points_to_polyline(), which validates the point sequence and returns a TractographyError on construction failure.

Export methods map Streamline to Gaia polylines for format-crate writers:

// .trk export (via ritk_trk)
let streamlines: Vec<Polyline<f64>> = self.streamlines
    .iter().map(|s| s.geometry().clone()).collect();
let trk = TrkTractogram::from_polylines(streamlines, header)?;

// .tck export (via ritk_tck) — same pattern
// .trx export (via ritk_trx) — same pattern

Each export method rips the geometry out of the Streamline wrapper and hands a Vec<Polyline<f64>> to the format-crate writer, which stores vertex coordinates in the target format’s coordinate system and binary encoding.

Connectivity Matrix Input in ritk-connectome

ritk-connectome::build_connectivity_matrix takes &[Polyline<f64>] as its streamline input:

use ritk_connectome::build_connectivity_matrix;

let matrix = build_connectivity_matrix(&parcellation, &streamlines)?;

The function samples each streamline’s vertices against a volumetric parcellation, mapping each transition between parcellation regions to an edge in the connectivity matrix. The first() and last() accessors determine endpoint regions for ROI-to-ROI connectivity; all intermediate vertices determine the path.

Relationship to the Aabb Type

Polyline and Aabb compose naturally for spatial queries:

// Quick reject: does the streamline intersect the parcellation volume?
if !parcellation_bbox.intersects(&streamline.aabb()) {
    continue; // skip this streamline entirely
}

// Per-vertex containment test
for point in streamline.points() {
    if parcellation_bbox.contains_point(&point) {
        let label = parcellation.label_at_point(&point);
        // ...
    }
}

This spatial pipeline — AABB cull → per-vertex query — is the core of the connectome construction algorithm.

Boundary

Under ADR 0036 decision 2, the polyline type and its geometric predicates belong to Gaia, never to RITK. RITK owns the integration and termination policy that produces Polyline<f64> instances from a direction field; it does not define a local polyline, curve, or line type. The format-crate readers and the connectome builder consume Polyline<f64> as their canonical curve representation.

References

• ADR 0036 — Neuroimaging and MR Ownership
• Deterministic Streamline Tractography — the RITK integration side
• Connectome Construction and Graph Measures — the RITK consumption side
• .trk Format, .tck Format, .trx Format — interchange format chapters

Moirai Parallel Execution Backend

Moirai is the Atlas-native hybrid concurrency library that provides the parallel execution backend for batched voxel-wise fitting in the diffusion MRI pipeline. It blends synchronous data-parallel operators (the Rayon-replacement surface) with a lock-free work-stealing scheduler, a unified hybrid runtime that serves both async and parallel work, and compile-time execution policies that monomorphize to zero-cost concrete paths.

At the RITK level, Moirai is surfaced through the MoiraiBackend type in ritk-image and the broader Coeus backend traits. Compute-heavy diffusion kernels — DTI log-linear solves, DKI Levenberg-Marquardt iterations, CSD non-negative deconvolution — distribute across CPU cores through Moirai’s data-parallel operators. The important design point is that Moirai does not require a separate image type or parallel-only API: readers, writers, filters, and registration code are parameterized over a backend, so the same Image<f32, B, 3> contract can execute on the sequential backend for deterministic single-threaded work or on Moirai for throughput.

Architecture

Moirai has three layers, each a separate domain crate:

A .par() worker task runs on the same unified hybrid scheduler as async work, not a separate thread pool. The sync return shape is a property of the operation, not an isolation boundary: a parallel kernel can spawn or drive async work on the same runtime, and sync operators compose cleanly inside pure compute kernels without async contagion.

The scheduler uses lock-free work-stealing: each worker owns a private Chase-Lev deque and a SplitDeque that keeps the freshest work on a local stack (zero atomics) and spills only the oldest half to a shared deque for thieves. Batch stealing transfers up to 16 items at a time, reducing steal contention when spawn rate greatly exceeds steal rate.

Execution Policies

Every data-parallel operator is generic over an ExecutionPolicy — a zero-sized type-level marker that decides whether to parallelize. Because the decision is an associated function, generic code over P: ExecutionPolicy monomorphizes to one concrete path with no value passed and no dynamic dispatch.

The Adaptive policy is the everyday entry point. It parallelizes only at or above 1,024 elements, where dispatch and join overhead is amortized; below that threshold the sequential path eliminates all scheduling cost. Because Adaptive is itself a zero-sized type, .par() is a fully monomorphized, zero-cost abstraction — the parallel/sequential decision is a cheap inlined runtime check.

use moirai_parallel::{Adaptive, Parallel, Sequential, ExecutionPolicy};

// Force sequential (deterministic / nested)
for_each_with::<Sequential, _, _>(&data, |x| process(x));

// Force parallel (ignore threshold)
for_each_with::<Parallel, _, _>(&data, |x| process(x));

// Adaptive default (everyday use)
data.par().for_each(|x| process(x));

Extension Traits

The primary surface is two extension traits on slices:

ParallelSlice<T>

pub trait ParallelSlice<T> {
    /// Adaptive, auto-routing parallel view over `&[T]`.
    fn par(&self) -> ParRef<'_, T, Adaptive>;
}

ParRef provides read-only parallel combinators:

ParallelSliceMut<T>

pub trait ParallelSliceMut<T> {
    /// Adaptive, auto-routing mutable parallel view over `&mut [T]`.
    fn par_mut(&mut self) -> ParMut<'_, T, Adaptive>;
}

ParMut provides mutable in-place combinators:

use moirai_parallel::{ParallelSlice, ParallelSliceMut};

let v: Vec<u64> = (0..1_000_000).collect();

// Map-reduce over a read-only slice
let sum = v.par().map_reduce(0u64, |&x| x, |a, b| a + b);

// In-place mutation of every element
let mut doubled = v.clone();
doubled.par_mut().for_each(|x| *x *= 2);

Data-Parallel Free Functions

For rare cases that need to force a specific policy via turbofish, every operator is also available as a *_with::<P> free function. These are the low-level override and are what the extension traits delegate to.

Read-only operators

Mutable operators

Index-domain operators

use moirai_parallel::{for_each_with, map_collect_with, map_reduce_with, Sequential, Parallel};

// Force sequential: deterministic, no scheduling overhead
for_each_with::<Sequential, _, _>(&small_data, |x| debug_assert!(x > 0));

// Force parallel: even for small inputs (benchmarking, nested regions)
let results: Vec<_> = map_collect_with::<Parallel, _, _, _>(&data, expensive);

// Elementwise product across multiple slices (index-aligned)
let dot = reduce_index_with::<Adaptive>(
    n, 0.0_f64,
    |i| a[i] * b[i],
    |x, y| x + y,
);

Chunk Operators

For batched workloads — like per-voxel model fitting where each “chunk” is a contiguous slice of the voxel array — chunk operators amortize dispatch overhead by assigning each worker a range of consecutive elements. The diffusion pipeline typically uses for_each_chunk_mut_with_state to give each worker a reusable scratch buffer (design matrix, workspace arrays) across its assigned voxels.

use moirai_parallel::for_each_chunk_mut_with_state;

// Per-voxel DTI fitting: each worker gets a private design matrix
for_each_chunk_mut_with_state::<Adaptive, _, _, _, _>(
    &mut voxel_results,
    256,                             // chunk size (voxels per task)
    || VoxelWorkspace::new(),         // init per-worker scratch
    |ws, chunk| {
        for result in chunk {
            *result = fit_dti(ws, &signal[result.index]);
        }
    },
);

The multi-buffer variants (pair, triple, quad) are fused kernels that update several output arrays in one authoritative pass — for example, writing FA, MD, and PEV fields simultaneously — without allocating intermediate tuples.

Borrowing Scope

scope provides a Rayon-style region for spawning parallel sub-tasks that capture non-'static references. It blocks until every spawned job has completed, so borrowed data cannot escape.

use moirai_parallel::scope;

let data: Vec<u64> = (0..1000).collect();
let sum = std::sync::atomic::AtomicU64::new(0);

scope(|s| {
    s.spawn(|| sum.fetch_add(data.iter().sum::<u64>(), Ordering::Relaxed));
    s.spawn(|| sum.fetch_add(data.len() as u64, Ordering::Relaxed));
});
// Both tasks guaranteed complete; `data` is still alive.

Join

join runs two closures concurrently. The left closure is scheduled on a worker thread; the right closure runs on the caller. If the scheduler refuses the left closure (shutting down or admission queue full), it runs on the caller instead, so join gracefully degrades to sequential but never drops work.

use moirai_parallel::join;

let (a, b) = join(
    || expensive_first_half(),
    || expensive_second_half(),
);

For explicit policy selection, use join_with::<P>.

Work-Stealing Deques

The scheduler’s lock-free work-stealing deques live in moirai-scheduler:

The split-deque design is the key performance optimisation for the diffusion pipeline: an owner’s push/pop against its private stack costs zero atomics, and only work old enough to be worth stealing reaches the shared deque. A thief that finds the victim empty returns StealResult::Empty (look elsewhere) rather than spinning; a thief that loses a race returns StealResult::Retry (the victim may still hold work, so retrying the same deque is worthwhile).

Batch stealing (steal_batch) transfers half the available items (up to 16) in one atomic CAS, cutting steal contention by an order of magnitude when many small tasks are spawned rapidly — the diffusion pipeline’s voxel-chunk dispatch pattern.

Memory ordering follows the weak-memory-correct Chase-Lev formulation of Lê, Pop, Cohen & Nardelli (PPoPP 2013), with a Morrison-Afek fence-free fast path for pop on x86/x86_64 when the deque holds enough items that no steal can reach the popped slot. Old buffers freed by resize are retired to an epoch-reclamation list and freed only when no accessor is in-flight.

NUMA Awareness

On multi-socket systems, the scheduler discovers hardware topology via CpuTopology and uses AdaptiveBackoff (spin → yield → sleep) for NUMA-aware victim selection. A thief prefers victims on the same NUMA node to reduce cross-socket memory latency, and the backoff strategy ramps from a brief spin (good when work is about to appear) to a park (good when the system is fully idle).

Diffusion Pipeline Integration

Moirai is the execution substrate for every compute-heavy pass in the diffusion MRI pipeline. Each model’s fitting kernel dispatches across voxels through Moirai’s chunk operators:

The pattern is the same in every case:

1. A flat Vec<ModelResult> is allocated with one element per voxel (all elements MaybeUninit).
2. for_each_chunk_mut_with_state parcels the result array into chunks and assigns each to a worker.
3. Each worker initialises a scratch workspace once (init closure), then processes its assigned voxels sequentially within the chunk, reusing the workspace.
4. On completion, every slot is initialised and the Vec<MaybeUninit<T>> is safely reinterpreted as Vec<T>.

This is zero-copy (no intermediate per-voxel allocations), cache-friendly (consecutive voxels share the same signal array and scheme), and NUMA-aware (the scheduler prefers to assign chunks to workers on the same socket as the memory).

// Pattern for batched voxel-wise fitting (simplified)
fn fit_volume(
    signal: &[f64],
    scheme: &GradientScheme,
    config: &DtiConfig,
) -> Vec<DtiResult> {
    let n_voxels = signal.len() / scheme.len();
    let mut results: Vec<MaybeUninit<DtiResult>> = Vec::with_capacity(n_voxels);
    unsafe { results.set_len(n_voxels); }

    for_each_chunk_mut_with_state::<Adaptive, _, _, _, _>(
        &mut results,
        256,                           // voxels per task
        || DtiWorkspace::new(scheme),  // one design matrix per worker
        |ws, chunk| {
            for (i, slot) in chunk.iter_mut().enumerate() {
                let voxel_start = (ws.chunk_base + i) * scheme.len();
                slot.write(ws.fit(&signal[voxel_start..]));
            }
        },
    );

    // SAFETY: every slot written by the parallel region above.
    unsafe { std::mem::transmute(results) }
}

Performance Characteristics

Coeus Tensor Integration

Coeus provides the N-dimensional tensor type, backend abstraction, and automatic differentiation engine that RITK’s image processing pipeline is built on. The Image<T, B, D> type in ritk-core is a Coeus Tensor<T, B> with fixed rank D and spatial metadata; every image filter, registration transform, and diffusion model that operates on pixel data does so through the Coeus tensor/backend contract.

This chapter documents the integration surface — the Coeus types and traits that RITK consumes — rather than the provider internals, which are documented in the Coeus book.

Architecture

The Coeus tensor stack has three layers:

RITK consumes all three through ritk-core’s Image type and the diffusion pipeline’s model fitting.

Backend Abstraction

ComputeBackend

ComputeBackend is the device abstraction — the trait that lets the same Tensor<T, B> code run on different hardware:

pub trait ComputeBackend: Send + Sync + Clone + 'static {
    type DeviceBuffer<T: Scalar>: StorageMut<T>;

    fn name(&self) -> &'static str;
    fn num_threads(&self) -> usize;
    fn allocate<T: Scalar>(&self, len: usize) -> Self::DeviceBuffer<T>;
    fn fill<T: Scalar>(&self, dst: &mut Self::DeviceBuffer<T>, val: T);
    fn copy_to_device<T: Scalar>(&self, src: &[T], dst: &mut Self::DeviceBuffer<T>);
    fn copy_to_host<T: Scalar>(&self, src: &Self::DeviceBuffer<T>, dst: &mut [T]);
}

Two concrete backends ship in coeus-core:

MoiraiBackend is the default for Tensor<T> (when no backend type parameter is supplied). It routes parallel work through the same unified hybrid scheduler that serves both async and parallel tasks — see the Moirai Execution Backend chapter.

Backend

The Backend trait extends ComputeBackend with a parallel_for method:

pub unsafe trait Backend: ComputeBackend + Default {
    fn parallel_for<F>(&self, start: usize, end: usize, f: F)
    where
        F: Fn(usize) + Send + Sync + 'static;
}

The Backend trait is unsafe because implementors must guarantee that parallel_for returns only after every invocation of f has completed. CPU kernels use this synchronization to keep scoped borrows alive across dispatch — the diffusion pipeline’s voxel-chunk pattern depends on this.

Storage Hierarchy

The storage traits form a sealed hierarchy that backs Tensor<T, B>:

CpuStorage<T> is the standard CPU-side implementation backed by the Mnemosyne allocator. CowStorage<T> wraps any StorageMut<T> with copy-on-write semantics: make_unique() triggers a deep copy when the allocation is shared.

Tensor Type

Tensor<T, B = MoiraiBackend> is the generic N-dimensional array:

pub struct Tensor<T: Scalar, B: ComputeBackend = MoiraiBackend> {
    storage: B::DeviceBuffer<T>,
    layout: Layout,
    _backend: PhantomData<B>,
}

Construction

Every constructor has an _on variant that accepts an explicit backend reference, and a default-backend variant that uses B::default().

Accessors

Views

Zero-copy views share the underlying storage with the parent tensor:

Views carry their own Layout but reference the same storage allocation. Mutation triggers copy-on-write: the view’s make_unique() deep-copies the shared buffer first.

Backend Transfer

to_backend(&new_backend) copies tensor memory to a different device. When source and destination are the same type, the transfer is a zero-copy clone. When they differ, the data is read from the source backend to host memory, then written to the destination backend.

use coeus_tensor::Tensor;
use coeus_core::{SequentialBackend, MoiraiBackend};

let t: Tensor<f32, SequentialBackend> = Tensor::from_slice([2, 3], &[1.0; 6]);
let t_moirai: Tensor<f32, MoiraiBackend> = t.to_backend(&MoiraiBackend::new());

Image Integration in RITK

RITK’s Image<T, B, D> is a rank-fixed wrapper around Tensor<T, B> with spatial metadata (origin, spacing, direction cosine matrix). The relationship is:

// Simplified — the actual Image type is in ritk-core
pub struct Image<T: Scalar, B: ComputeBackend, const D: usize> {
    tensor: Tensor<T, B>,   // pixel data, Coeus-owned
    origin: [f64; 3],       // physical origin (mm, LPS)
    spacing: [f64; 3],      // voxel spacing (mm)
    direction: [[f64; 3]; 3], // direction cosine matrix
}

This means:

• Every image operation (pixel access, filtering, resampling) runs through the Coeus tensor/backend contract.
• Backend substitution is seamless: Image<f32, SequentialBackend, 3> for deterministic single-threaded work, Image<f32, MoiraiBackend, 3> for parallel throughput — same pixel access API.
• Zero-copy views (slice, transpose, permute) share pixel storage with the parent image via Coeus’s COW semantics.

Autodiff Integration

coeus-autograd provides Var<T, B> — a differentiable variable that wraps a Tensor<T, B> and tracks operations for gradient computation. The diffusion pipeline uses this for nonlinear model fitting (DKI, NODDI):

use coeus_autograd::Var;
use coeus_tensor::Tensor;

let x: Var<f64> = Var::from_tensor(Tensor::from_vec(vec![1.0, 2.0, 3.0]));
let y = x.sin() + x.cos();
y.backward();
// x.grad() contains the gradient dy/dx

The DKI Levenberg-Marquardt path uses coeus-autograd to compute the analytic Jacobian of the kurtosis signal model with respect to the 21 free parameters. The NODDI model combines analytic gradients for the intra-axonal compartment with finite-difference gradients for the Watson dispersion integral.

Relationship to the Diffusion Pipeline

The Coeus tensor is the universal data carrier: image pixels, signal vectors, design matrices, and Jacobian matrices all flow through Tensor<T, B>. RITK adds spatial metadata and domain-specific interpretation; Coeus provides the generic N-D array operations.

Boundary

Coeus owns tensor storage, backends, autodiff, and the Tensor<T, B> type. RITK owns the Image<T, B, D> wrapper with spatial metadata, format-crate I/O, and domain-specific image processing. A RITK-local tensor type or local gradient computation is a boundary violation — the Coeus tensor is the single source of truth for pixel data.

References

• Coeus Nonlinear Least-Squares Solver — the LM solver that consumes Coeus tensors for Jacobian evaluation
• Moirai Parallel Execution Backend — the default backend for Tensor<T>
• Backend Dispatch — RITK’s backend selection and dispatch patterns
• Diffusion Models — the models that consume Coeus tensors for design-matrix assembly and fitting

Leto Operations Integration

Leto is the Atlas-owned N-dimensional strided array library that provides the array storage, layout descriptors, zero-copy views, and storage backends used by every numeric operation in the RITK pipeline. The diffusion models assemble design matrices into Array2<f64>, signal vectors into Array1<f64>, and pass them to leto-ops solvers through the .view() method.

This chapter documents the integration surface — the Leto array types, views, and storage that RITK consumes — rather than the solver operations themselves, which are covered in the Leto Linear Algebra Operations provider chapter. The two chapters are paired: the provider chapter documents what Leto computes; this chapter documents how RITK holds and passes the data.

Architecture

Leto has a single domain crate (leto) that carries arrays, views, and storage, plus a companion operations crate (leto-ops) that carries the linear algebra and statistical kernels. The crates are separated so that consumers that only need array storage (like apollo-sht) can depend on leto without pulling in the decompositions and solvers.

RITK consumes leto directly for Array1/Array2 construction, and leto-ops for every linear solve in the diffusion pipeline.

Array Type

Array<T, S, N> is the generic N-dimensional strided array:

pub struct Array<T, S, const N: usize> {
    layout: Layout<N>,    // shape, strides, offset
    storage: S,           // heap-allocated or borrowed backing
    _marker: PhantomData<T>,
}

The type is rank-generic (const N: usize) and storage-generic (S), which means the same array API works for 1-D signals, 2-D matrices, and 3-D image volumes, backed by any storage that implements Storage<T>.

Type Aliases

Most RITK code uses rank-specific type aliases that fix S = VecStorage<T>:

use leto::{Array1, Array2};

// DTI: one design matrix per voxel (N_dwi x 6)
let design: Array2<f64> = Array2::zeros([n_dwi, 6]);

// CSD: non-negative fODF coefficients
let coefficients: Array1<f64> = result.solution;

Construction