qsm_core/pipeline/two_pass.rs
1//! Two-pass artefact reduction.
2//!
3//! Strong susceptibility sources (air–tissue interfaces, haemorrhage, implants) leave the phase
4//! unrecoverable around them, and a phase-quality mask drops those voxels — leaving holes. Filling
5//! the holes and reconstructing once puts unreliable field values back into the inversion, and the
6//! streaking that follows spreads well beyond the hole. Dropping them instead leaves nothing to
7//! report where the source actually is.
8//!
9//! Two-pass does both and keeps the better half of each: reconstruct once on the mask with its
10//! holes intact ("reliable"), once on the mask the user configured ("filled"), then take the
11//! reliable reconstruction wherever it is defined and the filled one everywhere else.
12//!
13//! Stewart AW, Robinson SD, O'Brien K, et al. QSMxT: Robust masking and artifact reduction for
14//! quantitative susceptibility mapping. Magn Reson Med. 2022;87(3):1289-1300.
15//! <https://doi.org/10.1002/mrm.29048>
16//!
17//! The orchestration — running the two passes, caching them, naming their outputs — belongs to
18//! the consumer, since each drives the pipeline stages itself. What lives here is the part both
19//! consumers must agree on: the default reliable-mask recipe, and the combination itself.
20
21use super::config::{
22 MaskOp, MaskSection, MaskThresholdMethod, MaskingInput, PipelineError,
23};
24
25/// The reliable-pass mask recipe used when the consumer does not supply one.
26///
27/// An Otsu threshold on the ROMEO phase-quality map, and nothing else. The refinements the
28/// `robust-threshold` recipe applies — dilate, fill-holes, erode — are deliberately absent: the
29/// holes *are* the signal here, and filling them would collapse this pass onto the other one.
30pub fn default_reliable_sections() -> Vec<MaskSection> {
31 vec![MaskSection {
32 input: MaskingInput::PhaseQuality,
33 generator: MaskOp::Threshold {
34 method: MaskThresholdMethod::Otsu,
35 value: None,
36 },
37 refinements: vec![],
38 }]
39}
40
41/// Restrict the reliable mask to the region being reconstructed.
42///
43/// The reliable mask decides where the reliable pass's values are *preferred*, so anywhere it
44/// reaches past the main mask it would promote a reconstruction over a region the caller excluded
45/// outright. Being a subset of the reconstruction domain is a property of the method rather than
46/// a setting, so this is applied whatever recipe the two masks came from.
47pub fn restrict_reliable_mask(reliable: &[u8], main: &[u8]) -> Result<Vec<u8>, PipelineError> {
48 if reliable.len() != main.len() {
49 return Err(PipelineError::DimensionMismatch {
50 expected: main.len(),
51 got: reliable.len(),
52 });
53 }
54 Ok(reliable.iter().zip(main).map(|(&r, &m)| r & m).collect())
55}
56
57/// Combine the two passes: the reliable reconstruction where it is defined, the filled one
58/// elsewhere.
59///
60/// `support` is where `reliable` carries a reconstruction — not the mask the pass was *asked*
61/// for. Background removal erodes, so the two differ by a rim of voxels in which `reliable` is
62/// zero; selecting on the requested mask would keep those zeros and ring the holes with a seam of
63/// empty voxels. Passing the support the inversion actually produced makes the seam impossible.
64pub fn combine_two_pass(
65 reliable: &[f64],
66 filled: &[f64],
67 support: &[u8],
68) -> Result<Vec<f64>, PipelineError> {
69 if reliable.len() != filled.len() {
70 return Err(PipelineError::DimensionMismatch {
71 expected: reliable.len(),
72 got: filled.len(),
73 });
74 }
75 if support.len() != reliable.len() {
76 return Err(PipelineError::DimensionMismatch {
77 expected: reliable.len(),
78 got: support.len(),
79 });
80 }
81
82 Ok(reliable
83 .iter()
84 .zip(filled)
85 .zip(support)
86 .map(|((&r, &f), &s)| if s != 0 { r } else { f })
87 .collect())
88}
89
90#[cfg(test)]
91mod tests {
92 use super::*;
93
94 #[test]
95 fn default_recipe_leaves_the_holes_alone() {
96 let sections = default_reliable_sections();
97 assert_eq!(sections.len(), 1);
98 assert_eq!(sections[0].input, MaskingInput::PhaseQuality);
99 assert!(matches!(
100 sections[0].generator,
101 MaskOp::Threshold { method: MaskThresholdMethod::Otsu, value: None }
102 ));
103 // Filling or closing here would defeat the purpose of the pass.
104 assert!(sections[0].refinements.is_empty());
105 }
106
107 #[test]
108 fn combine_picks_each_pass_where_it_belongs() {
109 let reliable = [1.0, 2.0, 3.0, 4.0];
110 let filled = [10.0, 20.0, 30.0, 40.0];
111 let support = [1u8, 0, 1, 0];
112 let out = combine_two_pass(&reliable, &filled, &support).unwrap();
113 assert_eq!(out, vec![1.0, 20.0, 3.0, 40.0]);
114 }
115
116 /// A zero inside the reliable pass's support is a real reconstructed value, not a gap: the
117 /// filled pass must not overwrite it.
118 #[test]
119 fn a_reconstructed_zero_is_kept() {
120 let out = combine_two_pass(&[0.0], &[9.0], &[1]).unwrap();
121 assert_eq!(out, vec![0.0]);
122 }
123
124 /// The eroded rim belongs to the filled pass — that is the point of selecting on the support
125 /// rather than on the mask the pass was asked for.
126 #[test]
127 fn the_eroded_rim_comes_from_the_filled_pass() {
128 // Mask asked for [1,1,1]; background removal eroded it to [0,1,0].
129 let reliable = [0.0, 5.0, 0.0];
130 let filled = [7.0, 8.0, 9.0];
131 let out = combine_two_pass(&reliable, &filled, &[0, 1, 0]).unwrap();
132 assert_eq!(out, vec![7.0, 5.0, 9.0]);
133 }
134
135 #[test]
136 fn combine_rejects_mismatched_lengths() {
137 assert!(combine_two_pass(&[1.0, 2.0], &[1.0], &[1, 1]).is_err());
138 assert!(combine_two_pass(&[1.0, 2.0], &[1.0, 2.0], &[1]).is_err());
139 }
140
141 #[test]
142 fn restrict_intersects() {
143 let out = restrict_reliable_mask(&[1, 1, 0, 1], &[1, 0, 1, 0]).unwrap();
144 assert_eq!(out, vec![1, 0, 0, 0]);
145 }
146
147 #[test]
148 fn restrict_rejects_mismatched_lengths() {
149 assert!(restrict_reliable_mask(&[1, 1], &[1]).is_err());
150 }
151}