Skip to main content

qsm_core/
lib.rs

1//! # QSM-Core
2//!
3//! A Rust library for Quantitative Susceptibility Mapping (QSM) of the brain.
4//!
5//! QSM-Core reconstructs magnetic susceptibility maps from MRI phase data:
6//! brain extraction, phase unwrapping, background field removal, dipole
7//! inversion, and susceptibility source separation.
8//!
9//! ## Which API should I use?
10//!
11//! The crate offers two entry points at different levels:
12//!
13//! - **[`pipeline`] — the high-level API.** Describe a scan with [`ScanMetadata`]
14//!   and a [`QsmPipelineConfig`], then call the `run_*` stage functions. This is
15//!   the easiest way to go from phase data to a susceptibility map and the
16//!   recommended starting point.
17//! - **Algorithm building blocks — the low-level API.** Each algorithm
18//!   ([`bgremove::vsharp()`], [`inversion::tv_admm()`], …) is a plain function taking
19//!   a [`Grid`], a `*Params` struct, and (for iterative methods) a progress
20//!   callback. Use these when you want to wire stages together yourself.
21//!
22//! ```no_run
23//! use qsm_core::{Grid, bet, unwrap, bgremove, inversion};
24//! use qsm_core::bet::BetParams;
25//! use qsm_core::bgremove::VsharpParams;
26//! use qsm_core::inversion::TvParams;
27//!
28//! # fn run(phase: &[f64], magnitude: &[f64]) {
29//! let grid = Grid::new(128, 128, 64, 1.0, 1.0, 1.0);
30//! let bdir = (0.0, 0.0, 1.0);
31//!
32//! let mask = bet::run_bet(magnitude, &grid, &BetParams::default(), |_, _| {});
33//! let unwrapped = unwrap::laplacian_unwrap(phase, &mask, &grid);
34//! let (local, eroded) = bgremove::vsharp(&unwrapped, &mask, &grid, &VsharpParams::default(), |_, _| {});
35//! let chi = inversion::tv_admm(&local, &eroded, &grid, bdir, &TvParams::default(), |_, _| {});
36//! # let _ = chi;
37//! # }
38//! ```
39//!
40//! ## Modules
41//!
42//! **High-level pipeline**
43//! - [`pipeline`] — config-driven full reconstruction
44//!
45//! **Algorithm building blocks**
46//! - [`bet`] — brain extraction (BET)
47//! - [`unwrap`] — phase unwrapping (ROMEO, Laplacian)
48//! - [`bgremove`] — background field removal (V-SHARP, SHARP, RESHARP, PDF, iSMV, mSMV, LBV, HARPERELLA)
49//! - [`inversion`] — dipole inversion (TKD, TSVD, Tikhonov, TV, NLTV, RTS, MEDI, iLSQR, TGV)
50//! - [`separation`] — paramagnetic/diamagnetic source separation (χ-separation, R2\*-QSM, WaveSep)
51//! - [`swi`] — susceptibility weighted imaging (CLEAR-SWI)
52//! - [`fieldmap`] — multi-echo phase combination and B0 field mapping
53//! - [`r2star`] — R2\*/T2\* mapping (ARLO)
54//! - [`mask`] — mask thresholding and morphology
55//! - [`homogeneity`] — receive-field bias correction
56//!
57//! **Core types & I/O**
58//! - [`Grid`] — 3D volume descriptor shared by every algorithm
59//! - [`io`] — NIfTI read/write
60//!
61//! ## Feature Flags
62//!
63//! - **`parallel`** — enables [Rayon](https://docs.rs/rayon)-based multi-threading for FFT and iterative solvers
64//! - **`simd`** — enables SIMD acceleration via the [`wide`](https://docs.rs/wide) crate
65//!
66//! ## Algorithms
67//!
68//! | Stage | Methods |
69//! |-------|---------|
70//! | Brain extraction | BET |
71//! | Phase unwrapping | ROMEO, Laplacian |
72//! | Background removal | V-SHARP, SHARP, RESHARP, PDF, iSMV, LBV, SDF |
73//! | Dipole inversion | TKD, TSVD, Tikhonov, TV-ADMM, NLTV, RTS, MEDI, TGV, iLSQR |
74//! | Combined unwrap+BFR | HARPERELLA, iHARPERELLA |
75//! | SWI | CLEAR-SWI |
76//! | Separation | Chi-separation (Shin 2021 iLSQR-initialized, MEDI-based) |
77//! | Multi-echo | MCPC-3D-S, R2\*/T2\* (ARLO), bias correction |
78//! | Utilities | Frangi vesselness, surface curvature, Otsu thresholding, QSMART |
79
80// ============================================================================
81// Internal plumbing — public so advanced/downstream code can reach it, but
82// hidden from the documented API surface. No stability guarantees.
83// ============================================================================
84#[macro_use]
85#[doc(hidden)]
86pub mod par;
87#[doc(hidden)]
88pub mod grid;
89#[doc(hidden)]
90pub mod fft;
91#[doc(hidden)]
92pub mod priority_queue;
93#[doc(hidden)]
94pub mod region_grow;
95#[doc(hidden)]
96pub mod kernels;
97#[doc(hidden)]
98pub mod solvers;
99#[doc(hidden)]
100pub mod utils;
101
102// ============================================================================
103// Core types
104// ============================================================================
105pub use grid::Grid;
106
107// ============================================================================
108// Deep-learning model registry and weight management
109// ============================================================================
110/// Registry of deep-learning QSM models and their downloadable ONNX weights.
111///
112/// Query the catalog with [`models::all_models`] / [`models::find_model`]. Native
113/// hosts fetch weights on use with the `download` feature; the `onnx` feature runs
114/// them via the pure-Rust `tract` engine from a byte buffer (native and WASM).
115pub mod models;
116
117// ============================================================================
118// High-level pipeline
119// ============================================================================
120pub mod pipeline;
121
122// ============================================================================
123// Algorithm building blocks
124// ============================================================================
125pub mod bet;
126pub mod unwrap;
127pub mod bgremove;
128pub mod inversion;
129pub mod separation;
130pub mod swi;
131
132/// Multi-echo phase combination and B0 field mapping.
133///
134/// Building blocks for turning multi-echo wrapped phase into a B0 field map:
135/// phase-offset removal (MCPC-3D-S / ASPIRE), weighted B0 estimation, and
136/// linear multi-echo fitting.
137pub mod fieldmap {
138    pub use crate::utils::multi_echo::{
139        phase_offset_removal, mcpc3ds_combine, calculate_b0_weighted, multi_echo_linear_fit,
140        bipolar_correction, field_to_hz, nan_box_smooth_3d, nan_box_smooth_3d_phase, gaussian_box_sizes,
141        PhaseOffsetParams, CoilCombinationResult, LinearFitParams, LinearFitResult, B0WeightType,
142    };
143}
144
145/// R2\*/T2\* mapping from multi-echo magnitude (ARLO).
146pub mod r2star {
147    pub use crate::utils::r2star::{
148        r2star_arlo, t2star_from_r2star, use_arlo,
149    };
150}
151
152/// R2/T2 mapping from multi-echo spin-echo (EPG), and R2' = R2* − R2.
153///
154/// EPG-based fitting models imperfect refocusing (B1 < 1) so it removes the
155/// stimulated-echo bias that a mono-exponential fit suffers. [`r2prime`] combines
156/// the spin-echo R2 with a gradient-echo R2* (from [`r2star`]) for chi-separation.
157pub mod relaxometry {
158    pub use crate::utils::epg::{
159        epg_cpmg_echoes, r2_epg, r2prime, R2EpgParams,
160    };
161}
162
163/// MP-PCA denoising for multi-volume data (e.g. multi-echo magnitude).
164///
165/// Random-matrix-theory denoising (Veraart 2016) that removes noise along the
166/// volume dimension while preserving spatial edges. Applying [`mppca_denoise`]
167/// to multi-echo magnitude before R2*/R2 fitting reduces relaxation-rate
168/// variance without blurring structure.
169pub mod denoise {
170    pub use crate::utils::denoise::mppca_denoise;
171}
172
173/// Gibbs-ringing removal (Kellner 2016 subvoxel shifts) for k-space-truncated
174/// images. Complementary to [`denoise`]: unringing removes truncation ringing,
175/// MP-PCA removes random noise. Recommended before R2*/R2 fitting.
176pub mod unring {
177    pub use crate::utils::gibbs::{gibbs_unring, gibbs_unring_masked, gibbs_unring_volume};
178}
179
180/// Brain-mask thresholding and morphology.
181///
182/// Mask generation ([`otsu_threshold`](crate::mask::otsu_threshold)) plus
183/// morphological operations (erode/dilate/close/fill-holes) and sphere masks.
184pub mod mask {
185    pub use crate::utils::mask::{
186        create_sphere_mask, apply_mask_zero, erode_mask, dilate_mask,
187    };
188    pub use crate::utils::threshold::otsu_threshold;
189    pub use crate::utils::curvature::morphological_close;
190    pub use crate::utils::bias_correction::fill_holes;
191}
192
193/// Receive-field (B1−) bias correction for magnitude images.
194pub mod homogeneity {
195    pub use crate::utils::bias_correction::{
196        makehomogeneous, get_sensitivity, HomogeneityParams,
197    };
198}
199
200// ============================================================================
201// I/O
202// ============================================================================
203/// Scan geometry from the NIfTI affine: B0 direction, obliquity, and resampling to an axial grid.
204///
205/// The dipole kernel lives in the voxel grid, so an oblique acquisition must either supply the
206/// true B0 direction or be resampled to a cardinal-aligned grid. Wrapped phase has to be
207/// resampled in the complex domain — see [`geometry::resample_complex_to_axial`].
208/// Cropping reconstruction to the region that carries signal, and putting the answer back.
209///
210/// FFT-based stages cost `O(N log N)` in the whole grid, not in the brain. See [`crop`] for why
211/// the box is also rounded up to FFT-friendly sizes, and for the wrap-around caveat.
212pub mod crop;
213
214pub mod geometry;
215
216pub mod io;