Skip to main content

qsm_core/pipeline/
config.rs

1//! Pipeline configuration types
2//!
3//! Defines the configuration structs and enums for the QSM pipeline runner.
4//! These are pure algorithm config types (no serde). Consumers that need
5//! serialization (e.g. qsmxt-config) provide their own serde wrappers
6//! and convert to these types.
7
8use crate::bgremove::{
9    IsmvParams, LbvParams, PdfParams, ResharpParams, SharpParams, SdfParams, VsharpParams,
10    HarperellaParams, MsmvParams,
11};
12use crate::inversion::{
13    IlsqrParams, MediParams, NltvParams, RtsParams, TgvParams, TikhonovParams, TkdParams, TvParams,
14    NdiParams, FansiParams, L1QsmParams, WhQsmParams, HdQsmParams, TfiParams, AmpPeParams,
15};
16use crate::separation::{
17    ChiSepIlsqrParams, ChiSepParams, DecomposeParams, HcChisepParams, R2starQsmParams,
18    WaveSepParams,
19};
20use crate::unwrap::romeo::RomeoParams;
21use crate::utils::multi_echo::{B0WeightType, LinearFitParams};
22use crate::utils::QsmartParams;
23
24
25// =========================================================================
26// Selection enums
27// =========================================================================
28
29/// Phase unwrapping algorithm
30#[derive(Clone, Copy, Debug, PartialEq)]
31pub enum UnwrappingAlgorithm {
32    Romeo,
33    Laplacian,
34}
35
36/// Whether an algorithm is valid when B0 does not lie along the voxel `+z` axis.
37///
38/// The dipole relationship is direction-dependent, and the FFT that implements it lives in the
39/// voxel grid, so an oblique acquisition has to be handled deliberately. There are three cases,
40/// and the difference matters because getting it wrong is silent — the reconstruction completes
41/// and the numbers are simply wrong.
42#[derive(Debug, Clone, Copy, PartialEq, Eq)]
43pub enum OrientationSupport {
44    /// The B0 direction is an explicit parameter — the dipole kernel is built from it — so the
45    /// method reconstructs correctly on the grid the data was acquired on, whatever its
46    /// orientation. See [`crate::geometry::b0_direction_from_affine`].
47    Arbitrary,
48    /// The method never uses B0. The SMV-family background removals rest on the spherical mean
49    /// value property of harmonic fields, which is rotation-invariant.
50    NotApplicable,
51    /// Assumes B0 along `+z` and offers no way to say otherwise. Deep-learning methods learned
52    /// the dipole relationship from axially-acquired training data, so there is no direction to
53    /// rotate; oblique data must be resampled to a cardinal grid first
54    /// ([`crate::geometry::resample_complex_to_axial`]).
55    ///
56    /// Unrolled networks with a physics data-consistency term (LPCNN, MoDL-QSM, NeXtQSM) build
57    /// that term from the true direction, so they are partly corrected — but their learned prior
58    /// is still axial, so they belong here.
59    AxialOnly,
60}
61
62impl OrientationSupport {
63    /// Whether oblique data must be resampled before this method can be trusted.
64    pub fn requires_axial(&self) -> bool {
65        matches!(self, OrientationSupport::AxialOnly)
66    }
67}
68
69/// Background field removal algorithm
70#[derive(Clone, Copy, Debug, PartialEq)]
71pub enum BgRemovalAlgorithm {
72    Vsharp,
73    Pdf,
74    Lbv,
75    Ismv,
76    Sharp,
77    Resharp,
78    Harperella,
79    Iharperella,
80    /// BFRnet deep-learning background removal (requires the `onnx` feature and
81    /// the `bfrnet` model weights; see [`crate::models`]).
82    Bfrnet,
83}
84
85/// Dipole inversion algorithm
86#[cfg_attr(feature = "introspection", derive(serde::Serialize))]
87#[derive(Clone, Copy, Debug, PartialEq)]
88pub enum InversionAlgorithm {
89    Tkd,
90    Tsvd,
91    Tikhonov,
92    Tv,
93    Rts,
94    Nltv,
95    Medi,
96    /// Preconditioned Total Field Inversion (single-step, total field).
97    Tfi,
98    Ilsqr,
99    Tgv,
100    Qsmart,
101    /// Nonlinear Dipole Inversion (FANSI ndi.m).
102    Ndi,
103    /// FANSI nonlinear Total Variation (nlTV).
104    Fansi,
105    /// FANSI nonlinear Total Generalized Variation (nlTGV).
106    FansiTgv,
107    /// L1 data-fidelity QSM (FANSI nlL1TV / PI-QSM).
108    L1qsm,
109    /// Weak-Harmonic QSM (FANSI WH_nlTV).
110    Whqsm,
111    /// Hybrid two-stage L1→L2 QSM (HD-QSM).
112    Hdqsm,
113    /// Approximate Message Passing with built-in Parameter Estimation (AMP-PE).
114    AmpPe,
115    /// xQSM deep-learning dipole inversion (requires the `onnx` feature and the
116    /// `xqsm` model weights; see [`crate::models`]).
117    Xqsm,
118    /// QSMnet deep-learning dipole inversion (requires the `onnx` feature and the
119    /// `qsmnet` model weights; see [`crate::models`]).
120    Qsmnet,
121    /// QSMnet+ deep-learning dipole inversion (susceptibility-scaling augmented;
122    /// requires the `onnx` feature and the `qsmnet-plus` weights).
123    QsmnetPlus,
124    /// AutoQSM single-step reconstruction (requires the `onnx` feature and the
125    /// `autoqsm` weights). NOTE: takes the **total** field — it does its own
126    /// background removal, so the `local_field_ppm` argument should be the total field.
127    Autoqsm,
128    /// QSMGAN deep-learning dipole inversion (local field → χ; `onnx` + `qsmgan` weights).
129    Qsmgan,
130    /// IR2QSM unrolled deep-learning dipole inversion (local field → χ; `onnx` + `ir2qsm` weights).
131    Ir2qsm,
132    /// LPCNN learned-proximal deep-learning dipole inversion (local field → χ; `onnx` + `lpcnn` weights).
133    Lpcnn,
134    /// MoDL-QSM model-based deep-learning dipole inversion (local field → χ33/STI component;
135    /// `onnx` + `modl-qsm` weights).
136    ModlQsm,
137    /// NeXtQSM single-step reconstruction (requires `onnx` + the two `nextqsm` weight files).
138    /// NOTE: takes the **total** field — it does its own background removal.
139    Nextqsm,
140    /// iQSM single-step reconstruction from wrapped **phase** (end-to-end: joint unwrap +
141    /// background removal + inversion). Not a dipole-inversion-stage option — use
142    /// [`super::run_iqsm`]. Rejected by `run_dipole_inversion`.
143    Iqsm,
144    /// iQSM+ single-step reconstruction from wrapped **phase** with orientation-adaptive
145    /// feature editing (uses the B0 direction). End-to-end; use [`super::run_iqsm_plus`].
146    IqsmPlus,
147}
148
149/// B0 estimation method
150#[derive(Clone, Copy, Debug, PartialEq)]
151pub enum B0EstimationMethod {
152    WeightedAvg,
153    LinearFit,
154}
155
156/// QSM referencing method
157#[derive(Clone, Copy, Debug, PartialEq)]
158pub enum QsmReference {
159    Mean,
160    None,
161}
162
163// =========================================================================
164// Masking types
165// =========================================================================
166
167/// Input data source for mask generation
168#[derive(Clone, Copy, Debug, PartialEq)]
169pub enum MaskingInput {
170    MagnitudeFirst,
171    Magnitude,
172    MagnitudeLast,
173    PhaseQuality,
174}
175
176/// Threshold method for mask generation
177#[derive(Clone, Copy, Debug, PartialEq)]
178pub enum MaskThresholdMethod {
179    Otsu,
180    Fixed,
181    Percentile,
182}
183
184/// A single mask operation (generator or refinement)
185#[derive(Clone, Debug, PartialEq)]
186pub enum MaskOp {
187    Threshold {
188        method: MaskThresholdMethod,
189        value: Option<f64>,
190    },
191    Bet {
192        fractional_intensity: f64,
193    },
194    Erode {
195        iterations: usize,
196    },
197    Dilate {
198        iterations: usize,
199    },
200    Close {
201        radius: usize,
202    },
203    FillHoles {
204        max_size: usize,
205    },
206    GaussianSmooth {
207        sigma_mm: f64,
208    },
209    /// Signal-gated erosion: peel only low-signal boundary voxels (skull-base / sinus dropout)
210    /// down to a depth cap. See [`crate::utils::signal_gated_erosion`].
211    ///
212    /// Gates on the **magnitude** image, never on the section's input: it divides out a
213    /// receive-coil bias estimate and compares against the in-mask median, which only means
214    /// anything for magnitude. Errors if no magnitude is supplied, so a phase-quality mask input
215    /// can still be refined with it as long as the magnitude is available.
216    SignalErode(crate::utils::SignalErosionParams),
217    /// HD-BET deep-learning brain extraction from the **magnitude** (a generator, like `Bet`),
218    /// whatever the section's input is. Requires the `onnx` feature and the `hd-bet` model
219    /// weights; see [`crate::models`].
220    HdBet(crate::bet::HdBetParams),
221}
222
223impl MaskOp {
224    /// Registry id of the deep-learning model this op runs, if any.
225    pub fn dl_model_id(&self) -> Option<&'static str> {
226        match self {
227            Self::HdBet(_) => Some("hd-bet"),
228            Self::Threshold { .. } | Self::Bet { .. } | Self::Erode { .. } | Self::Dilate { .. }
229            | Self::Close { .. } | Self::FillHoles { .. } | Self::GaussianSmooth { .. }
230            | Self::SignalErode(_) => None,
231        }
232    }
233}
234
235/// A mask section: input source + generator + refinements
236#[derive(Clone, Debug, PartialEq)]
237pub struct MaskSection {
238    pub input: MaskingInput,
239    pub generator: MaskOp,
240    pub refinements: Vec<MaskOp>,
241}
242
243impl MaskSection {
244    /// Get all operations (generator + refinements) in order
245    pub fn all_ops(&self) -> Vec<MaskOp> {
246        let mut ops = vec![self.generator.clone()];
247        ops.extend(self.refinements.iter().cloned());
248        ops
249    }
250}
251
252// =========================================================================
253// Per-stage config structs
254// =========================================================================
255
256/// Masking configuration
257#[derive(Clone, Debug)]
258pub struct MaskingConfig {
259    pub inhomogeneity_correction: bool,
260    pub homogeneity_sigma_mm: f64,
261    pub homogeneity_nbox: usize,
262    pub sections: Vec<MaskSection>,
263}
264
265impl Default for MaskingConfig {
266    fn default() -> Self {
267        Self {
268            inhomogeneity_correction: true,
269            homogeneity_sigma_mm: 7.0,
270            homogeneity_nbox: 15,
271            sections: vec![MaskSection {
272                input: MaskingInput::PhaseQuality,
273                generator: MaskOp::Threshold {
274                    method: MaskThresholdMethod::Otsu,
275                    value: None,
276                },
277                refinements: vec![
278                    MaskOp::Dilate { iterations: 1 },
279                    MaskOp::FillHoles { max_size: 0 },
280                    MaskOp::Erode { iterations: 1 },
281                ],
282            }],
283        }
284    }
285}
286
287/// Field mapping configuration
288#[derive(Clone, Debug)]
289pub struct FieldMappingConfig {
290    pub unwrapping_algorithm: UnwrappingAlgorithm,
291    pub phase_offset_removal: bool,
292    pub phase_offset_sigma: [f64; 3],
293    pub bipolar_correction: bool,
294    pub b0_estimation: B0EstimationMethod,
295    pub b0_weight_type: B0WeightType,
296    pub romeo_params: RomeoParams,
297    pub linear_fit_params: LinearFitParams,
298}
299
300impl Default for FieldMappingConfig {
301    fn default() -> Self {
302        Self {
303            unwrapping_algorithm: UnwrappingAlgorithm::Romeo,
304            phase_offset_removal: true,
305            phase_offset_sigma: [10.0, 10.0, 5.0],
306            bipolar_correction: false,
307            b0_estimation: B0EstimationMethod::WeightedAvg,
308            b0_weight_type: B0WeightType::PhaseSNR,
309            romeo_params: RomeoParams::default(),
310            linear_fit_params: LinearFitParams::default(),
311        }
312    }
313}
314
315/// Background removal configuration
316#[derive(Clone, Debug)]
317pub struct BgRemovalConfig {
318    pub algorithm: BgRemovalAlgorithm,
319    pub vsharp: VsharpParams,
320    pub pdf: PdfParams,
321    pub lbv: LbvParams,
322    pub ismv: IsmvParams,
323    pub sharp: SharpParams,
324    pub resharp: ResharpParams,
325    pub harperella: HarperellaParams,
326    pub sdf: SdfParams,
327    /// mSMV refinement parameters (`b0`/`te` are overridden from scan metadata by
328    /// the dispatcher). Used by the `msmv_refine` post-step.
329    pub msmv: MsmvParams,
330    /// Apply mSMV boundary-shadow refinement after the primary BFR (Roberts 2024).
331    /// mSMV is a refinement, not a standalone primary remover, so it is exposed
332    /// only as this post-step (redundant after `Ismv`, which is already SMV-based).
333    pub msmv_refine: bool,
334}
335
336impl Default for BgRemovalConfig {
337    fn default() -> Self {
338        Self {
339            algorithm: BgRemovalAlgorithm::Vsharp,
340            vsharp: VsharpParams::default(),
341            pdf: PdfParams::default(),
342            lbv: LbvParams::default(),
343            ismv: IsmvParams::default(),
344            sharp: SharpParams::default(),
345            resharp: ResharpParams::default(),
346            harperella: HarperellaParams::default(),
347            sdf: SdfParams::default(),
348            msmv: MsmvParams::default(),
349            msmv_refine: false,
350        }
351    }
352}
353
354/// Dipole inversion configuration
355#[derive(Clone, Debug)]
356pub struct InversionConfig {
357    pub algorithm: InversionAlgorithm,
358    pub tkd: TkdParams,
359    pub tsvd: TkdParams,
360    pub tikhonov: TikhonovParams,
361    pub tv: TvParams,
362    pub rts: RtsParams,
363    pub nltv: NltvParams,
364    pub medi: MediParams,
365    pub tfi: TfiParams,
366    pub ilsqr: IlsqrParams,
367    pub tgv: TgvParams,
368    pub qsmart: QsmartParams,
369    pub ndi: NdiParams,
370    /// Shared by the `Fansi` (nlTV) and `FansiTgv` (nlTGV) algorithms; the
371    /// dispatcher sets `is_tgv` from the selected algorithm.
372    pub fansi: FansiParams,
373    pub l1qsm: L1QsmParams,
374    pub whqsm: WhQsmParams,
375    pub hdqsm: HdQsmParams,
376    /// AMP-PE (`b0` is overridden from scan metadata by the dispatcher).
377    pub amp_pe: AmpPeParams,
378    /// Overlap-tiling for the deep-learning inversions, as `(core, halo)` in voxels. `None`
379    /// (default) runs the net whole-volume; `Some` runs it patch-by-patch (bounded memory) via the
380    /// `*_tiled` variants — an approximation, mainly for memory-constrained targets (e.g. WASM).
381    /// Stored as a plain tuple so this config type stays available without the `onnx` feature;
382    /// the dispatcher converts it to a [`crate::inversion::TileConfig`]. Ignored by the
383    /// natively-patch-based nets (autoqsm/qsmgan) and by non-DL algorithms.
384    pub tile: Option<(usize, usize)>,
385}
386
387impl Default for InversionConfig {
388    fn default() -> Self {
389        Self {
390            algorithm: InversionAlgorithm::Rts,
391            tkd: TkdParams::default(),
392            tsvd: TkdParams::default(),
393            tikhonov: TikhonovParams::default(),
394            tv: TvParams::default(),
395            rts: RtsParams::default(),
396            nltv: NltvParams::default(),
397            medi: MediParams::default(),
398            tfi: TfiParams::default(),
399            ilsqr: IlsqrParams::default(),
400            tgv: TgvParams::default(),
401            qsmart: QsmartParams::default(),
402            ndi: NdiParams::default(),
403            fansi: FansiParams::default(),
404            l1qsm: L1QsmParams::default(),
405            whqsm: WhQsmParams::default(),
406            hdqsm: HdQsmParams::default(),
407            amp_pe: AmpPeParams::default(),
408            tile: None,
409        }
410    }
411}
412
413/// Susceptibility source-separation algorithm (χ+ / χ−).
414#[cfg_attr(feature = "introspection", derive(serde::Serialize))]
415#[derive(Clone, Copy, Debug, PartialEq)]
416pub enum SeparationAlgorithm {
417    /// Shin 2021 projected-CG, QSM-initialized (local field + R2' + magnitude + QSM).
418    ChiSepIlsqr,
419    /// MEDI-based Gauss-Newton (local field + R2' + magnitude).
420    ChiSepMedi,
421    /// Closed-form from a QSM + R2* (Dimov 2022; R2* fit from magnitude if absent).
422    R2starQsm,
423    /// Wavelet-L1 proximal-gradient from a QSM + R2' (Fang 2023).
424    WaveSep,
425    /// Signal-domain per-voxel fit from a QSM + multi-echo magnitude (Chen 2021).
426    Decompose,
427    /// Hollow-cylinder fit from a QSM + R2' + multi-echo magnitude (Wharton & Bowtell).
428    HcChisep,
429    /// SUSEP-Net deep-learning separation from QSM + R2' + local field (requires
430    /// the `onnx` feature and the `susep-net` weights; see [`crate::models`]).
431    SusepNet,
432    /// χ-sepnet (SNU-LIST) deep-learning separation from local field + QSM + R2'
433    /// (requires the `onnx` feature and the `chi-sepnet` weights; see [`crate::models`]).
434    ChiSepNet,
435}
436
437impl BgRemovalAlgorithm {
438    /// How this method behaves when B0 is not along voxel `+z`.
439    ///
440    /// Only PDF projects onto dipole fields and therefore needs the direction; every other
441    /// method here is harmonic/SMV-based and is direction-independent. BFRnet is a network, but
442    /// background removal has no dipole orientation to get wrong in the way inversion does — it
443    /// separates harmonic from non-harmonic content, so it is treated as direction-independent.
444    pub fn orientation_support(&self) -> OrientationSupport {
445        match self {
446            BgRemovalAlgorithm::Pdf => OrientationSupport::Arbitrary,
447            BgRemovalAlgorithm::Vsharp
448            | BgRemovalAlgorithm::Lbv
449            | BgRemovalAlgorithm::Ismv
450            | BgRemovalAlgorithm::Sharp
451            | BgRemovalAlgorithm::Resharp
452            | BgRemovalAlgorithm::Harperella
453            | BgRemovalAlgorithm::Iharperella
454            | BgRemovalAlgorithm::Bfrnet => OrientationSupport::NotApplicable,
455        }
456    }
457}
458
459impl InversionAlgorithm {
460    /// How this method behaves when B0 is not along voxel `+z`.
461    ///
462    /// Every classical inversion builds its dipole kernel from the supplied direction, directly
463    /// or through the shared ADMM/FANSI spectral setup, so all of them handle oblique data on
464    /// the acquired grid. The deep-learning methods take no direction at all.
465    pub fn orientation_support(&self) -> OrientationSupport {
466        match self {
467            InversionAlgorithm::Tkd
468            | InversionAlgorithm::Tsvd
469            | InversionAlgorithm::Tikhonov
470            | InversionAlgorithm::Tv
471            | InversionAlgorithm::Rts
472            | InversionAlgorithm::Nltv
473            | InversionAlgorithm::Medi
474            | InversionAlgorithm::Tfi
475            | InversionAlgorithm::Ilsqr
476            | InversionAlgorithm::Tgv
477            | InversionAlgorithm::Qsmart
478            | InversionAlgorithm::Ndi
479            | InversionAlgorithm::Fansi
480            | InversionAlgorithm::FansiTgv
481            | InversionAlgorithm::L1qsm
482            | InversionAlgorithm::Whqsm
483            | InversionAlgorithm::Hdqsm
484            | InversionAlgorithm::AmpPe => OrientationSupport::Arbitrary,
485            InversionAlgorithm::Xqsm
486            | InversionAlgorithm::Qsmnet
487            | InversionAlgorithm::QsmnetPlus
488            | InversionAlgorithm::Autoqsm
489            | InversionAlgorithm::Qsmgan
490            | InversionAlgorithm::Ir2qsm
491            | InversionAlgorithm::Lpcnn
492            | InversionAlgorithm::ModlQsm
493            | InversionAlgorithm::Nextqsm
494            | InversionAlgorithm::Iqsm
495            | InversionAlgorithm::IqsmPlus => OrientationSupport::AxialOnly,
496        }
497    }
498}
499
500impl SeparationAlgorithm {
501    /// How this method behaves when B0 is not along voxel `+z`.
502    ///
503    /// The model-based separations carry the direction through their field terms. The rest
504    /// either consume an already-reconstructed χ map and R2\* (so orientation was settled
505    /// upstream) or are networks trained on axial data.
506    pub fn orientation_support(&self) -> OrientationSupport {
507        match self {
508            SeparationAlgorithm::ChiSepIlsqr | SeparationAlgorithm::ChiSepMedi => {
509                OrientationSupport::Arbitrary
510            }
511            SeparationAlgorithm::R2starQsm
512            | SeparationAlgorithm::WaveSep
513            | SeparationAlgorithm::Decompose
514            | SeparationAlgorithm::HcChisep => OrientationSupport::NotApplicable,
515            SeparationAlgorithm::SusepNet | SeparationAlgorithm::ChiSepNet => {
516                OrientationSupport::AxialOnly
517            }
518        }
519    }
520}
521
522// ─── Deep-learning model-registry mapping ───
523//
524// These map each stage enum's deep-learning variants to their [`crate::models`]
525// registry id. Exhaustive matches: adding a variant forces a decision here, and the
526// `registry_models_are_pipeline_wired` test asserts every registry model of a pipeline
527// stage is claimed by some variant — so a model added to the registry but not wired
528// into an enum fails the build's tests (registry↔pipeline drift guard).
529
530impl BgRemovalAlgorithm {
531    /// Model-registry id for deep-learning variants, else `None` (classical methods).
532    pub fn dl_model_id(self) -> Option<&'static str> {
533        match self {
534            Self::Bfrnet => Some("bfrnet"),
535            Self::Vsharp | Self::Pdf | Self::Lbv | Self::Ismv | Self::Sharp
536            | Self::Resharp | Self::Harperella | Self::Iharperella => None,
537        }
538    }
539    /// Every variant, for exhaustiveness in tests/tools.
540    pub const VARIANTS: &'static [Self] = &[
541        Self::Vsharp, Self::Pdf, Self::Lbv, Self::Ismv, Self::Sharp,
542        Self::Resharp, Self::Harperella, Self::Iharperella, Self::Bfrnet,
543    ];
544}
545
546impl InversionAlgorithm {
547    /// Model-registry id for deep-learning variants, else `None` (classical methods).
548    pub fn dl_model_id(self) -> Option<&'static str> {
549        match self {
550            Self::Xqsm => Some("xqsm"),
551            Self::Qsmnet => Some("qsmnet"),
552            Self::QsmnetPlus => Some("qsmnet-plus"),
553            Self::Autoqsm => Some("autoqsm"),
554            Self::Qsmgan => Some("qsmgan"),
555            Self::Ir2qsm => Some("ir2qsm"),
556            Self::Lpcnn => Some("lpcnn"),
557            Self::ModlQsm => Some("modl-qsm"),
558            Self::Nextqsm => Some("nextqsm"),
559            Self::Iqsm => Some("iqsm"),
560            Self::IqsmPlus => Some("iqsm-plus"),
561            Self::Tkd | Self::Tsvd | Self::Tikhonov | Self::Tv | Self::Rts | Self::Nltv
562            | Self::Medi | Self::Tfi | Self::Ilsqr | Self::Tgv | Self::Qsmart | Self::Ndi
563            | Self::Fansi | Self::FansiTgv | Self::L1qsm | Self::Whqsm | Self::Hdqsm
564            | Self::AmpPe => None,
565        }
566    }
567    pub const VARIANTS: &'static [Self] = &[
568        Self::Tkd, Self::Tsvd, Self::Tikhonov, Self::Tv, Self::Rts, Self::Nltv, Self::Medi,
569        Self::Tfi, Self::Ilsqr, Self::Tgv, Self::Qsmart, Self::Ndi, Self::Fansi, Self::FansiTgv,
570        Self::L1qsm, Self::Whqsm, Self::Hdqsm, Self::AmpPe, Self::Xqsm, Self::Qsmnet,
571        Self::QsmnetPlus, Self::Autoqsm, Self::Qsmgan, Self::Ir2qsm, Self::Lpcnn, Self::ModlQsm,
572        Self::Nextqsm, Self::Iqsm, Self::IqsmPlus,
573    ];
574}
575
576impl SeparationAlgorithm {
577    /// Model-registry id for deep-learning variants, else `None` (classical methods).
578    pub fn dl_model_id(self) -> Option<&'static str> {
579        match self {
580            Self::SusepNet => Some("susep-net"),
581            Self::ChiSepNet => Some("chi-sepnet"),
582            Self::ChiSepIlsqr | Self::ChiSepMedi | Self::R2starQsm | Self::WaveSep
583            | Self::Decompose | Self::HcChisep => None,
584        }
585    }
586    pub const VARIANTS: &'static [Self] = &[
587        Self::ChiSepIlsqr, Self::ChiSepMedi, Self::R2starQsm, Self::WaveSep,
588        Self::Decompose, Self::HcChisep, Self::SusepNet, Self::ChiSepNet,
589    ];
590}
591
592/// Configuration for the χ-separation stage.
593#[derive(Clone, Debug)]
594pub struct SeparationConfig {
595    pub algorithm: SeparationAlgorithm,
596    /// `cf` on the chi-sep params is overridden from scan metadata by the dispatcher.
597    pub chi_sep_ilsqr: ChiSepIlsqrParams,
598    pub chi_sep_medi: ChiSepParams,
599    /// `b0` overridden from scan metadata by the dispatcher.
600    pub r2star_qsm: R2starQsmParams,
601    pub wavesep: WaveSepParams,
602    /// `b0` overridden from scan metadata by the dispatcher.
603    pub decompose: DecomposeParams,
604    /// `b0` overridden from scan metadata; `se_echo_times` supplies the SE echoes.
605    pub hc_chisep: HcChisepParams,
606}
607
608impl Default for SeparationConfig {
609    fn default() -> Self {
610        Self {
611            algorithm: SeparationAlgorithm::ChiSepIlsqr,
612            chi_sep_ilsqr: ChiSepIlsqrParams::default(),
613            chi_sep_medi: ChiSepParams::default(),
614            r2star_qsm: R2starQsmParams::default(),
615            wavesep: WaveSepParams::default(),
616            decompose: DecomposeParams::default(),
617            hc_chisep: HcChisepParams::default(),
618        }
619    }
620}
621
622// =========================================================================
623// Top-level pipeline config
624// =========================================================================
625
626/// Complete QSM pipeline configuration.
627///
628/// Contains all per-stage configs. Consumers call individual stage functions
629/// (e.g. `run_field_mapping`, `run_bg_removal`) passing the relevant section.
630/// Masking config is used when the consumer needs to generate a mask.
631#[derive(Clone, Debug)]
632pub struct QsmPipelineConfig {
633    pub masking: MaskingConfig,
634    pub field_mapping: FieldMappingConfig,
635    pub bg_removal: BgRemovalConfig,
636    pub inversion: InversionConfig,
637    pub reference: QsmReference,
638}
639
640impl Default for QsmPipelineConfig {
641    fn default() -> Self {
642        Self {
643            masking: MaskingConfig::default(),
644            field_mapping: FieldMappingConfig::default(),
645            bg_removal: BgRemovalConfig::default(),
646            inversion: InversionConfig::default(),
647            reference: QsmReference::Mean,
648        }
649    }
650}
651
652// =========================================================================
653// Scan metadata and stage result types
654// =========================================================================
655
656/// Metadata about the scan
657#[derive(Clone, Debug)]
658pub struct ScanMetadata {
659    /// Volume dimensions (nx, ny, nz)
660    pub dims: (usize, usize, usize),
661    /// Voxel size in mm (vsx, vsy, vsz)
662    pub voxel_size: (f64, f64, f64),
663    /// Echo times in seconds
664    pub echo_times: Vec<f64>,
665    /// Main field strength in Tesla
666    pub field_strength: f64,
667    /// B0 direction as unit vector in voxel coordinates
668    pub b0_direction: (f64, f64, f64),
669}
670
671impl ScanMetadata {
672    /// Get a Grid from this metadata's dimensions and voxel sizes.
673    #[inline]
674    pub fn grid(&self) -> crate::Grid {
675        crate::Grid {
676            dims: self.dims,
677            voxel_size: self.voxel_size,
678        }
679    }
680}
681
682/// Results from field mapping stage
683pub struct FieldMappingResult {
684    /// B0 field map in ppm
685    pub b0_field_ppm: Vec<f64>,
686    /// Phase offset map (if phase offset removal was used)
687    pub phase_offset: Option<Vec<f64>>,
688}
689
690/// Results from background removal stage
691pub struct BgRemovalResult {
692    /// Local field in ppm
693    pub local_field_ppm: Vec<f64>,
694    /// Eroded mask
695    pub eroded_mask: Vec<u8>,
696}
697
698/// Pipeline error type
699#[derive(Debug)]
700pub enum PipelineError {
701    /// Invalid configuration
702    InvalidConfig(String),
703    /// Invalid input data
704    InvalidInput(String),
705    /// Algorithm failure
706    AlgorithmError(String),
707    /// Dimension mismatch
708    DimensionMismatch { expected: usize, got: usize },
709}
710
711impl std::fmt::Display for PipelineError {
712    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
713        match self {
714            Self::InvalidConfig(msg) => write!(f, "invalid config: {}", msg),
715            Self::InvalidInput(msg) => write!(f, "invalid input: {}", msg),
716            Self::AlgorithmError(msg) => write!(f, "algorithm error: {}", msg),
717            Self::DimensionMismatch { expected, got } => {
718                write!(f, "dimension mismatch: expected {}, got {}", expected, got)
719            }
720        }
721    }
722}
723
724impl std::error::Error for PipelineError {}
725
726#[cfg(test)]
727mod orientation_tests {
728    use super::*;
729
730    /// Every classical inversion builds its kernel from the supplied direction, so it must claim
731    /// Arbitrary. If a new one is added without wiring bdir through, this is where it shows up.
732    #[test]
733    fn classical_inversions_take_a_direction() {
734        for a in [
735            InversionAlgorithm::Tkd, InversionAlgorithm::Tsvd, InversionAlgorithm::Tikhonov,
736            InversionAlgorithm::Tv, InversionAlgorithm::Rts, InversionAlgorithm::Nltv,
737            InversionAlgorithm::Medi, InversionAlgorithm::Tfi, InversionAlgorithm::Ilsqr,
738            InversionAlgorithm::Tgv, InversionAlgorithm::Qsmart, InversionAlgorithm::Ndi,
739            InversionAlgorithm::Fansi, InversionAlgorithm::FansiTgv, InversionAlgorithm::L1qsm,
740            InversionAlgorithm::Whqsm, InversionAlgorithm::Hdqsm, InversionAlgorithm::AmpPe,
741        ] {
742            assert_eq!(a.orientation_support(), OrientationSupport::Arbitrary, "{a:?}");
743            assert!(!a.orientation_support().requires_axial(), "{a:?}");
744        }
745    }
746
747    /// Networks have no direction input, so oblique data must be resampled for them.
748    #[test]
749    fn learned_inversions_need_axial_data() {
750        for a in [
751            InversionAlgorithm::Xqsm, InversionAlgorithm::Qsmnet, InversionAlgorithm::QsmnetPlus,
752            InversionAlgorithm::Autoqsm, InversionAlgorithm::Qsmgan, InversionAlgorithm::Ir2qsm,
753            InversionAlgorithm::Lpcnn, InversionAlgorithm::ModlQsm, InversionAlgorithm::Nextqsm,
754            InversionAlgorithm::Iqsm, InversionAlgorithm::IqsmPlus,
755        ] {
756            assert!(a.orientation_support().requires_axial(), "{a:?}");
757        }
758    }
759
760    /// SMV-family background removal is harmonic and rotation-invariant; only PDF cares.
761    #[test]
762    fn background_removal_is_direction_independent_except_pdf() {
763        assert_eq!(BgRemovalAlgorithm::Pdf.orientation_support(), OrientationSupport::Arbitrary);
764        for a in [
765            BgRemovalAlgorithm::Vsharp, BgRemovalAlgorithm::Sharp, BgRemovalAlgorithm::Resharp,
766            BgRemovalAlgorithm::Ismv, BgRemovalAlgorithm::Lbv, BgRemovalAlgorithm::Harperella,
767            BgRemovalAlgorithm::Iharperella, BgRemovalAlgorithm::Bfrnet,
768        ] {
769            assert_eq!(a.orientation_support(), OrientationSupport::NotApplicable, "{a:?}");
770            assert!(!a.orientation_support().requires_axial(), "{a:?}");
771        }
772    }
773
774    /// No background removal or separation method should ever force a resample on its own.
775    #[test]
776    fn only_learned_methods_force_a_resample() {
777        for a in [SeparationAlgorithm::ChiSepIlsqr, SeparationAlgorithm::ChiSepMedi] {
778            assert_eq!(a.orientation_support(), OrientationSupport::Arbitrary, "{a:?}");
779        }
780        for a in [SeparationAlgorithm::SusepNet, SeparationAlgorithm::ChiSepNet] {
781            assert!(a.orientation_support().requires_axial(), "{a:?}");
782        }
783        for a in [
784            SeparationAlgorithm::R2starQsm, SeparationAlgorithm::WaveSep,
785            SeparationAlgorithm::Decompose, SeparationAlgorithm::HcChisep,
786        ] {
787            assert!(!a.orientation_support().requires_axial(), "{a:?}");
788        }
789    }
790}
791
792#[cfg(test)]
793mod tests {
794    use super::*;
795    use crate::models::{all_models, ModelStage};
796
797    /// Registry↔pipeline drift guard: every deep-learning model whose `stage` maps to a
798    /// pipeline stage enum must be claimed by some enum variant's `dl_model_id()`. A model
799    /// added to the registry but not wired into an enum (so no pipeline consumer could
800    /// select it) fails here. `PhaseToField` (iQFM) is intentionally exempt — it's exposed
801    /// as the standalone `run_iqfm` field-preparation building block, not a stage enum.
802    #[test]
803    fn registry_models_are_pipeline_wired() {
804        let inv: Vec<&str> = InversionAlgorithm::VARIANTS.iter().filter_map(|v| v.dl_model_id()).collect();
805        let bfr: Vec<&str> = BgRemovalAlgorithm::VARIANTS.iter().filter_map(|v| v.dl_model_id()).collect();
806        let sep: Vec<&str> = SeparationAlgorithm::VARIANTS.iter().filter_map(|v| v.dl_model_id()).collect();
807        for m in all_models() {
808            let (wired, enum_name) = match m.stage {
809                ModelStage::BackgroundRemoval => (bfr.contains(&m.id), "BgRemovalAlgorithm"),
810                ModelStage::DipoleInversion | ModelStage::SingleStep => (inv.contains(&m.id), "InversionAlgorithm"),
811                ModelStage::ChiSeparation => (sep.contains(&m.id), "SeparationAlgorithm"),
812                ModelStage::BrainExtraction => {
813                    (MaskOp::HdBet(Default::default()).dl_model_id() == Some(m.id), "MaskOp")
814                }
815                // iQFM: standalone run_iqfm, no stage enum.
816                ModelStage::PhaseToField => continue,
817            };
818            assert!(wired, "registry model '{}' (stage {:?}) has no {} variant — wire it or add a dl_model_id mapping", m.id, m.stage, enum_name);
819        }
820    }
821
822    #[test]
823    fn test_default_config() {
824        let config = QsmPipelineConfig::default();
825        assert_eq!(config.field_mapping.unwrapping_algorithm, UnwrappingAlgorithm::Romeo);
826        assert_eq!(config.bg_removal.algorithm, BgRemovalAlgorithm::Vsharp);
827        assert_eq!(config.inversion.algorithm, InversionAlgorithm::Rts);
828        assert_eq!(config.reference, QsmReference::Mean);
829        assert!(config.field_mapping.phase_offset_removal);
830        assert!(!config.field_mapping.bipolar_correction);
831    }
832
833    #[test]
834    fn test_default_masking_config() {
835        let config = MaskingConfig::default();
836        assert_eq!(config.sections.len(), 1);
837        assert_eq!(config.sections[0].input, MaskingInput::PhaseQuality);
838        assert_eq!(config.sections[0].refinements.len(), 3);
839    }
840
841    #[test]
842    fn test_mask_section_all_ops() {
843        let section = MaskSection {
844            input: MaskingInput::Magnitude,
845            generator: MaskOp::Threshold { method: MaskThresholdMethod::Otsu, value: None },
846            refinements: vec![MaskOp::Erode { iterations: 1 }, MaskOp::Dilate { iterations: 2 }],
847        };
848        let ops = section.all_ops();
849        assert_eq!(ops.len(), 3);
850        assert!(matches!(ops[0], MaskOp::Threshold { .. }));
851        assert!(matches!(ops[1], MaskOp::Erode { iterations: 1 }));
852        assert!(matches!(ops[2], MaskOp::Dilate { iterations: 2 }));
853    }
854}