qsm_core/

lib.rs

//! # QSM-Core
//!
//! A Rust library for Quantitative Susceptibility Mapping (QSM) of the brain.
//!
//! QSM-Core reconstructs magnetic susceptibility maps from MRI phase data:
//! brain extraction, phase unwrapping, background field removal, dipole
//! inversion, and susceptibility source separation.
//!
//! ## Which API should I use?
//!
//! The crate offers two entry points at different levels:
//!
//! - **[`pipeline`] — the high-level API.** Describe a scan with [`ScanMetadata`]
//!   and a [`QsmPipelineConfig`], then call the `run_*` stage functions. This is
//!   the easiest way to go from phase data to a susceptibility map and the
//!   recommended starting point.
//! - **Algorithm building blocks — the low-level API.** Each algorithm
//!   ([`bgremove::vsharp()`], [`inversion::tv_admm()`], …) is a plain function taking
//!   a [`Grid`], a `*Params` struct, and (for iterative methods) a progress
//!   callback. Use these when you want to wire stages together yourself.
//!
//! ```no_run
//! use qsm_core::{Grid, bet, unwrap, bgremove, inversion};
//! use qsm_core::bet::BetParams;
//! use qsm_core::bgremove::VsharpParams;
//! use qsm_core::inversion::TvParams;
//!
//! # fn run(phase: &[f64], magnitude: &[f64]) {
//! let grid = Grid::new(128, 128, 64, 1.0, 1.0, 1.0);
//! let bdir = (0.0, 0.0, 1.0);
//!
//! let mask = bet::run_bet(magnitude, &grid, &BetParams::default(), bet::FslLeftRightOrder::Radiological, |_, _| {});
//! let unwrapped = unwrap::laplacian_unwrap(phase, &mask, &grid, unwrap::LaplacianSolver::Dct);
//! let (local, eroded) = bgremove::vsharp(&unwrapped, &mask, &grid, &VsharpParams::default(), |_, _| {}).unwrap();
//! let chi = inversion::tv_admm(&local, &eroded, &grid, bdir, &TvParams::default(), |_, _| {});
//! # let _ = chi;
//! # }
//! ```
//!
//! ## Modules
//!
//! **High-level pipeline**
//! - [`pipeline`] — config-driven full reconstruction
//!
//! **Algorithm building blocks**
//! - [`bet`] — brain extraction (BET)
//! - [`unwrap`] — phase unwrapping (ROMEO, Laplacian, best path, PRELUDE)
//! - [`bgremove`] — background field removal (V-SHARP, SHARP, RESHARP, PDF, iSMV, LBV, HARPERELLA)
//! - [`inversion`] — dipole inversion (TKD, TSVD, Tikhonov, TV, NLTV, RTS, MEDI incl. MEDI-mSMV, iLSQR, STAR-QSM, TGV)
//! - [`separation`] — paramagnetic/diamagnetic source separation (χ-separation, R2\*-QSM, WaveSep)
//! - [`swi`] — susceptibility weighted imaging (CLEAR-SWI, SMWI)
//! - [`complex`] — real/imaginary input to magnitude/phase
//! - [`fieldmap`] — multi-echo phase combination and B0 field mapping
//! - [`r2star`] — R2\*/T2\* mapping (ARLO)
//! - [`mask`] — mask thresholding and morphology
//! - [`homogeneity`] — receive-field bias correction
//!
//! **Core types & I/O**
//! - [`Grid`] — 3D volume descriptor shared by every algorithm
//! - [`io`] — NIfTI read/write
//!
//! ## Feature Flags
//!
//! - **`parallel`** — enables [Rayon](https://docs.rs/rayon)-based multi-threading for FFT and iterative solvers
//! - **`simd`** — enables SIMD acceleration via the [`wide`](https://docs.rs/wide) crate
//!
//! ## Algorithms
//!
//! | Stage | Methods |
//! |-------|---------|
//! | Brain extraction | BET |
//! | Phase unwrapping | ROMEO, Laplacian, best path (3D-SRNCP), PRELUDE |
//! | Background removal | V-SHARP, SHARP, RESHARP, PDF, iSMV, LBV, SDF |
//! | Dipole inversion | TKD, TSVD, Tikhonov, TV-ADMM, NLTV, RTS, MEDI, TGV, iLSQR, STAR-QSM |
//! | Combined unwrap+BFR | HARPERELLA, iHARPERELLA |
//! | SWI | CLEAR-SWI |
//! | Separation | Chi-separation (Shin 2021 iLSQR-initialized, MEDI-based) |
//! | Input | Real/imaginary to magnitude/phase (slice-chopping correction) |
//! | Multi-echo | MCPC-3D-S, R2\*/T2\* (ARLO), bias correction |
//! | Utilities | Frangi vesselness, surface curvature, Otsu thresholding, QSMART |

// ============================================================================
// Internal plumbing — public so advanced/downstream code can reach it, but
// hidden from the documented API surface. No stability guarantees.
// ============================================================================
#[macro_use]
#[doc(hidden)]
pub mod par;
#[doc(hidden)]
pub mod grid;
#[doc(hidden)]
pub mod fft;
#[doc(hidden)]
pub mod priority_queue;
#[doc(hidden)]
pub mod region_grow;
#[doc(hidden)]
pub mod kernels;
#[doc(hidden)]
pub mod solvers;
#[doc(hidden)]
pub mod utils;

// ============================================================================
// Core types
// ============================================================================
// `SliceGapError` is deliberately not re-exported here: it is reached as
// `qsm_core::grid::SliceGapError`, matching `models::onnx::OnnxError` and
// `models::download::DownloadError`. Each error type stays with the concern that
// produces it rather than accumulating at the crate root. The one shared by every
// algorithm entry point, for input it cannot process, is `qsm_core::error::InputError`.
pub use grid::Grid;
pub use utils::echo_time::EchoTime;

/// [`InputError`](error::InputError): what an algorithm returns instead of panicking when it is
/// handed input it cannot process.
pub mod error;

/// Physical constants (gyromagnetic ratio).
pub mod constants;

// ============================================================================
// Deep-learning model registry and weight management
// ============================================================================
/// Registry of deep-learning QSM models and their downloadable ONNX weights.
///
/// Query the catalog with [`models::all_models`] / [`models::find_model`]. A host can
/// supply weights from memory with [`models::provide_weights`] (the only source in WASM);
/// native hosts can also use `$QSM_MODEL_DIR`, the cache, or fetch on use with the `download`
/// feature. The `onnx` feature runs them via the pure-Rust `tract` engine (native and WASM).
pub mod models;

// ============================================================================
// High-level pipeline
// ============================================================================
pub mod pipeline;

// ============================================================================
// Algorithm building blocks
// ============================================================================
pub mod bet;
pub mod segment;
pub mod unwrap;
pub mod bgremove;
pub mod inversion;
pub mod separation;
pub mod swi;

/// Real/imaginary input: conversion to magnitude and phase.
///
/// [`complex_to_mag_phase`](complex::complex_to_mag_phase) also detects and undoes
/// alternate-slice π "chopping" (a reconstruction artefact common in GE real/imaginary
/// exports), deciding once for all echoes, and can negate the phase for vendors (GE) whose
/// real/imaginary phase runs opposite to the Siemens/Philips convention.
pub mod complex {
    pub use crate::utils::complex::{
        complex_to_mag_phase, slice_chopping_coherence, ChopCorrection, ComplexToMagPhaseParams, MagPhaseEchoes,
        SliceChopping, CHOP_AMBIGUOUS_COHERENCE,
    };
}

/// Multi-echo phase combination and B0 field mapping.
///
/// Building blocks for turning multi-echo wrapped phase into a B0 field map:
/// phase-offset removal (MCPC-3D-S / ASPIRE), weighted B0 estimation,
/// linear multi-echo fitting, and MEDI's complex fit (`Fit_ppm_complex`).
pub mod fieldmap {
    pub use crate::utils::multi_echo::{
        phase_offset_removal, mcpc3ds_combine, calculate_b0_weighted, multi_echo_linear_fit,
        bipolar_correction, field_to_hz, nan_box_smooth_3d, nan_box_smooth_3d_phase, gaussian_box_sizes,
        PhaseOffsetParams, PhaseOffsetRemoval, CoilCombinationResult, LinearFitParams, LinearFitResult, B0WeightType,
        multi_echo_robust_fit, RobustFitParams, RobustFitResult, RobustLoss, EchoQuality,
        fit_ppm_complex, fit_ppm_complex_noise, FitPpmComplexParams, FitPpmComplexResult,
        echo_spacing_problem, FIT_PPM_COMPLEX_SPACING_TOLERANCE_S,
    };
}

/// R2\*/T2\* mapping from multi-echo magnitude (ARLO).
pub mod r2star {
    pub use crate::utils::r2star::{
        arlo, arlo_echo_times_problem, t2star_from_r2star,
    };
}

/// R2/T2 mapping from multi-echo spin-echo (EPG), and R2' = R2* − R2.
///
/// EPG-based fitting models imperfect refocusing (B1 < 1) so it removes the
/// stimulated-echo bias that a mono-exponential fit suffers. [`r2prime`] combines
/// the spin-echo R2 with a gradient-echo R2* (from [`r2star`]) for chi-separation,
/// and returns the region R2 was measured over alongside the map — a spin-echo
/// acquisition is often a slab rather than the whole head, and R2' is undefined
/// outside it.
///
/// With no spin-echo acquisition there is no R2 to subtract; [`r2primenet`] predicts
/// R2' from the GRE-derived R2* with a trained network (`onnx` feature), which is what
/// lets the R2'-consuming [`separation`] methods run on GRE-only data.
pub mod relaxometry {
    pub use crate::utils::epg::{
        epg_cpmg_echoes, r2_epg, r2prime, R2EpgParams, R2PrimeMap, EPG_TE_TOLERANCE_S,
    };
    // The parameters and limits are always exported, so hosts built without `onnx` can carry
    // R2PRIMEnet settings and check a patch; only the inference functions need the feature.
    pub use crate::utils::r2primenet::{
        R2PrimeNetNorm, R2PrimeNetParams, AUTHORS_PATCH, WASM_PATCH, B0_TOLERANCE_T, TRAINED_B0_T,
        POOL_MULTIPLE,
    };
    #[cfg(feature = "onnx")]
    pub use crate::utils::r2primenet::{r2primenet, r2primenet_from_magnitude};
}

/// MP-PCA denoising for multi-volume data (e.g. multi-echo magnitude).
///
/// Random-matrix-theory denoising (Veraart 2016) that removes noise along the
/// volume dimension while preserving spatial edges. Applying [`mppca_denoise`]
/// to multi-echo magnitude before R2*/R2 fitting reduces relaxation-rate
/// variance without blurring structure.
pub mod denoise {
    pub use crate::utils::denoise::mppca_denoise;
}

/// Gibbs-ringing removal by subvoxel shifts for k-space-truncated images: Bautista et al.'s 3-D
/// method for 3-D acquisitions (the default), Kellner et al.'s 2-D method for 2-D multi-slice
/// acquisitions (MRtrix3's implementation), and DIPY's 2-D implementation as a parity option.
/// See [`GibbsMode`](crate::utils::gibbs::GibbsMode). Complementary to [`denoise`]: unringing
/// removes truncation ringing, MP-PCA removes random noise. Recommended before R2*/R2 fitting.
pub mod unring {
    pub use crate::utils::gibbs::{gibbs_unring, gibbs_unring_masked, gibbs_unring_volume, GibbsMode};
    pub use crate::utils::gibbs_mrtrix::{
        gibbs_unring_2d, gibbs_unring_2d_complex, gibbs_unring_3d, gibbs_unring_3d_complex, MrtrixUnringParams,
    };
}

/// Brain-mask thresholding and morphology.
///
/// Mask generation ([`otsu_threshold`](crate::mask::otsu_threshold)) plus
/// morphological operations (erode/dilate/close/fill-holes, signal-gated erosion) and sphere
/// masks.
pub mod mask {
    pub use crate::utils::mask::{
        create_sphere_mask, apply_mask_zero, scanner_masked_phase, ScannerMaskedPhase,
        SCANNER_FILL_MIN_FRACTION,
    };
    pub use crate::utils::threshold::otsu_threshold;
    pub use crate::utils::morphology::{
        erode_mask, dilate_mask, morphological_close, fill_holes,
        signal_gated_erosion, SignalErosionParams,
    };
}

/// Receive-field (B1−) bias correction for magnitude images.
pub mod homogeneity {
    pub use crate::utils::bias_correction::{
        makehomogeneous, get_sensitivity, HomogeneityParams,
    };
}

// ============================================================================
// I/O
// ============================================================================
/// Cropping reconstruction to the region that carries signal, and putting the answer back.
///
/// FFT-based stages cost `O(N log N)` in the whole grid, not in the brain. See [`crop`] for why
/// the box is also rounded up to FFT-friendly sizes, and for the wrap-around caveat.
pub mod crop;

/// Scan geometry from the NIfTI affine: B0 direction, obliquity, and resampling to an axial grid.
///
/// The dipole kernel lives in the voxel grid, so an oblique acquisition must either supply the
/// true B0 direction or be resampled to a cardinal-aligned grid. Wrapped phase has to be
/// resampled in the complex domain — see [`geometry::resample_complex_to_axial`].
pub mod geometry;

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

/// Retrospective motion correction: co-register a series of repeats, runs or echoes onto one of
/// them, in the complex domain, and re-derive B0 for each.
///
/// The image-domain half of motion correction — what a library consuming reconstructed NIfTI can
/// do, as opposed to the navigator-driven prospective methods that need raw k-space. Two things
/// it exists to get right: wrapped phase only ever moves through the complex domain
/// ([`motion::MotionSeries::apply_complex`]), and a series that rotated has one B0 direction per
/// volume ([`motion::VolumeMotion::bdir`]), not one for the series.
pub mod motion;

pub mod io;