qsm_core/pipeline/
bg_removal.rs

1//! Background field removal stage
2//!
3//! Dispatcher that calls the appropriate background removal algorithm
4//! based on configuration, including algorithm-specific parameter defaults
5//! (V-SHARP radii computation, LBV's scan parameters, iSMV radius).
6
7use super::config::*;
8
9/// Run background field removal on a total field map.
10///
11/// # Arguments
12/// * `field_ppm` - Total field map in ppm
13/// * `mask` - Binary brain mask
14/// * `metadata` - Scan metadata
15/// * `config` - Background removal configuration
16/// * `magnitude` - The magnitude of every echo, stacked with the echo slowest. PDF requires it
17///   (its data weight is the `Fit_ppm_complex` noise map); the other methods ignore it.
18/// * `progress` - Progress callback (current_iter, max_iter)
19///
20/// # Returns
21/// `BgRemovalResult` with local field in ppm and eroded mask
22pub fn run_bg_removal(
23    field_ppm: &[f64],
24    mask: &[u8],
25    metadata: &ScanMetadata,
26    config: &BgRemovalConfig,
27    magnitude: Option<&[f64]>,
28    progress: &mut dyn FnMut(usize, usize),
29) -> Result<BgRemovalResult, PipelineError> {
30    // Every background removal here convolves with a kernel defined on the voxel grid, so a
31    // volume with gaps between its slices is refused rather than approximated. See
32    // ScanMetadata::require_contiguous_slices for why this is unconditional.
33    metadata.require_contiguous_slices()?;
34    let grid = metadata.grid();
35    let n_voxels = grid.n_total();
36
37    if field_ppm.len() != n_voxels {
38        return Err(PipelineError::DimensionMismatch {
39            expected: n_voxels,
40            got: field_ppm.len(),
41        });
42    }
43
44    let (local_field, eroded_mask) = match config.algorithm {
45        BgRemovalAlgorithm::Vsharp => {
46            if let crate::bgremove::VsharpGrid::Isotropic { spacing_mm } = config.vsharp.working_grid(&grid) {
47                if !(spacing_mm.is_finite() && spacing_mm > 0.0) {
48                    return Err(PipelineError::InvalidConfig(format!(
49                        "V-SHARP isotropic working grid spacing must be finite and > 0 mm, got {spacing_mm}"
50                    )));
51                }
52            }
53            crate::bgremove::vsharp(
54                field_ppm, mask, &grid, &config.vsharp, progress,
55            )
56            .map_err(|e| PipelineError::InvalidConfig(e.0))?
57        }
58        BgRemovalAlgorithm::Pdf => {
59            // PDF weights the fit by 1/N_std, the noise of the field, which MEDI's pipeline takes
60            // from Fit_ppm_complex (radians over the echo spacing, from the echo magnitudes). PDF
61            // takes ppm like the others, and N_std must be in the field's unit (it sets the CG
62            // tolerance), so it is converted with the factor that converts the field. That is
63            // the echo spacing by the noise map's unit, not a choice of echo time: the noise is
64            // that of the phase accrued over one echo spacing, and PDF itself is linear in the
65            // field, so it has no radian conversion of its own to choose.
66            let te = &metadata.echo_times;
67            let mag = super::inversion::echo_magnitudes("PDF", magnitude, te, n_voxels)?;
68            let delta_te = if te.len() == 1 { te[0] } else { te[1] - te[0] };
69            let rad_per_ppm = 2.0 * std::f64::consts::PI * crate::constants::GAMMA_BAR_HZ_PER_T * metadata.field_strength * delta_te * 1e-6;
70            let n_std: Vec<f64> = crate::utils::fit_ppm_complex_noise(mag.data, mag.n_echoes)
71                .into_iter()
72                .map(|v| v / rad_per_ppm)
73                .collect();
74            let local = crate::bgremove::pdf(
75                field_ppm, &n_std, mask, &grid,
76                metadata.b0_direction, &config.pdf, progress,
77            )?;
78            (local, mask.to_vec())
79        }
80        BgRemovalAlgorithm::Lbv => {
81            // LBV's tol is absolute in MEDI's units (radians over the echo spacing by default),
82            // so it converts the field itself; field strength and echo times are scan
83            // parameters from metadata.
84            let params = crate::bgremove::LbvParams {
85                b0: metadata.field_strength,
86                echo_times: super::inversion::scan_echo_times("LBV", config.lbv.te, metadata)?,
87                ..config.lbv.clone()
88            };
89            crate::bgremove::lbv(field_ppm, mask, &grid, &params, progress)?
90        }
91        BgRemovalAlgorithm::Ismv => {
92            crate::bgremove::ismv(
93                field_ppm, mask, &grid, &config.ismv, progress,
94            )?
95        }
96        BgRemovalAlgorithm::Sharp => {
97            crate::bgremove::sharp(
98                field_ppm, mask, &grid, &config.sharp,
99            )
100        }
101        BgRemovalAlgorithm::Resharp => {
102            crate::bgremove::resharp(
103                field_ppm, mask, &grid,
104                &config.resharp,
105                progress,
106            )
107        }
108        // HARPERELLA / iHARPERELLA are integrated unwrap + BFR methods, but this stage has the
109        // unwrapped total field in ppm, so they take it as such (ordinary Laplacian). Wrapped
110        // phase input is `HarperellaInput::WrappedPhase`, called directly.
111        BgRemovalAlgorithm::Harperella => {
112            crate::bgremove::harperella(
113                crate::bgremove::HarperellaInput::FieldPpm(field_ppm), mask, &grid,
114                &config.harperella,
115                progress,
116            )?
117        }
118        BgRemovalAlgorithm::Iharperella => {
119            crate::bgremove::iharperella(
120                crate::bgremove::HarperellaInput::FieldPpm(field_ppm), mask, &grid,
121                &config.iharperella,
122                progress,
123            )?
124        }
125        BgRemovalAlgorithm::Bfrnet => {
126            let local = run_bfrnet(field_ppm, mask, &grid, config.tile, progress)?;
127            (local, mask.to_vec())
128        }
129    };
130
131    Ok(BgRemovalResult {
132        local_field_ppm: local_field,
133        eroded_mask,
134    })
135}
136
137/// Source the BFRnet weights and run inference, whole-volume or tiled (`tile`). Requires the
138/// `onnx` feature; weights resolve through [`crate::models`] (bytes supplied with
139/// [`crate::models::provide_weights`], then `$QSM_MODEL_DIR`/cache/download on native).
140#[cfg(feature = "onnx")]
141fn run_bfrnet(
142    field_ppm: &[f64],
143    mask: &[u8],
144    grid: &crate::Grid,
145    tile: Option<(usize, usize)>,
146    progress: &mut dyn FnMut(usize, usize),
147) -> Result<Vec<f64>, PipelineError> {
148    let spec = crate::models::find_model("bfrnet")
149        .ok_or_else(|| PipelineError::InvalidConfig("bfrnet not in model registry".into()))?;
150    let bytes = crate::models::primary_weight_bytes(spec)
151        .map_err(PipelineError::InvalidConfig)?;
152    match tile {
153        Some((core, halo)) => crate::bgremove::bfrnet_tiled(
154            field_ppm, mask, grid, &bytes, &crate::inversion::TileConfig { core, halo }, progress,
155        ),
156        None => crate::bgremove::bfrnet(field_ppm, mask, grid, &bytes),
157    }
158    .map_err(|e| PipelineError::AlgorithmError(e.to_string()))
159}
160
161#[cfg(not(feature = "onnx"))]
162fn run_bfrnet(
163    _field_ppm: &[f64],
164    _mask: &[u8],
165    _grid: &crate::Grid,
166    _tile: Option<(usize, usize)>,
167    _progress: &mut dyn FnMut(usize, usize),
168) -> Result<Vec<f64>, PipelineError> {
169    Err(PipelineError::InvalidConfig(
170        "BFRnet requires building qsm-core with the 'onnx' feature".into(),
171    ))
172}
173
174#[cfg(test)]
175mod tests {
176    use super::*;
177
178    #[test]
179    fn test_bg_removal_dispatches_vsharp() {
180        let (nx, ny, nz) = (8, 8, 8);
181        let n = nx * ny * nz;
182        let field = vec![0.1; n];
183        let mask = vec![1u8; n];
184        let meta = ScanMetadata {
185            dims: (nx, ny, nz),
186            voxel_size: (1.0, 1.0, 1.0),
187            echo_times: vec![0.005],
188            field_strength: 3.0,
189            b0_direction: (0.0, 0.0, 1.0),
190            slice_geometry: None,
191            affine: None,
192        };
193
194        let result = run_bg_removal(
195            &field, &mask, &meta, &BgRemovalConfig::default(), None,
196            &mut |_, _| {},
197        );
198        assert!(result.is_ok());
199        let r = result.unwrap();
200        assert_eq!(r.local_field_ppm.len(), n);
201        assert_eq!(r.eroded_mask.len(), n);
202    }
203
204    #[test]
205    fn vsharp_rejects_a_bad_working_grid_spacing() {
206        let (nx, ny, nz) = (8, 8, 8);
207        let n = nx * ny * nz;
208        let meta = ScanMetadata {
209            dims: (nx, ny, nz),
210            voxel_size: (1.0, 1.0, 1.0),
211            echo_times: vec![0.005],
212            field_strength: 3.0,
213            b0_direction: (0.0, 0.0, 1.0),
214            slice_geometry: None,
215            affine: None,
216        };
217        for spacing_mm in [0.0, -1.0, f64::NAN] {
218            let mut config = BgRemovalConfig::default();
219            config.vsharp.grid = crate::bgremove::VsharpGrid::Isotropic { spacing_mm };
220            let r = run_bg_removal(&vec![0.1; n], &vec![1u8; n], &meta, &config, None, &mut |_, _| {});
221            assert!(matches!(r, Err(PipelineError::InvalidConfig(_))), "spacing {spacing_mm}");
222        }
223    }
224
225    /// The coverage test in `config.rs` reads the stage sources; this one actually calls the
226    /// stage, so the wiring is proven end to end rather than by inspection.
227    #[test]
228    fn a_gapped_acquisition_is_refused_before_any_work() {
229        let (nx, ny, nz) = (8, 8, 8);
230        let n = nx * ny * nz;
231        let field = vec![0.1; n];
232        let mask = vec![1u8; n];
233        let mut meta = ScanMetadata {
234            dims: (nx, ny, nz),
235            voxel_size: (1.0, 1.0, 3.0),
236            echo_times: vec![0.005],
237            field_strength: 3.0,
238            b0_direction: (0.0, 0.0, 1.0),
239            slice_geometry: None,
240            affine: None,
241        };
242
243        // undeclared: nothing can be concluded, so it runs
244        assert!(run_bg_removal(&field, &mask, &meta, &BgRemovalConfig::default(), None, &mut |_, _| {})
245            .is_ok());
246
247        // declared contiguous: runs
248        meta.slice_geometry =
249            Some(crate::grid::SliceGeometry { thickness: 3.0, axis: 2 });
250        assert!(run_bg_removal(&field, &mask, &meta, &BgRemovalConfig::default(), None, &mut |_, _| {})
251            .is_ok());
252
253        // declared gapped: refused, for every algorithm the enum offers
254        meta.slice_geometry =
255            Some(crate::grid::SliceGeometry { thickness: 2.0, axis: 2 });
256        for algorithm in [
257            BgRemovalAlgorithm::Vsharp, BgRemovalAlgorithm::Pdf, BgRemovalAlgorithm::Lbv,
258            BgRemovalAlgorithm::Ismv, BgRemovalAlgorithm::Sharp, BgRemovalAlgorithm::Resharp,
259            BgRemovalAlgorithm::Harperella, BgRemovalAlgorithm::Iharperella,
260        ] {
261            let config = BgRemovalConfig { algorithm, ..BgRemovalConfig::default() };
262            match run_bg_removal(&field, &mask, &meta, &config, None, &mut |_, _| {}) {
263                Err(PipelineError::InvalidInput(msg)) => {
264                    assert!(msg.contains("1.0000 mm gap"), "{algorithm:?}: {msg}");
265                }
266                Err(other) => panic!("{algorithm:?} gave {other:?}, not InvalidInput"),
267                Ok(_) => panic!("{algorithm:?} accepted a gapped volume"),
268            }
269        }
270    }
271
272    #[test]
273    fn test_bg_removal_dispatches_pdf() {
274        let (nx, ny, nz) = (8, 8, 8);
275        let n = nx * ny * nz;
276        let field = vec![0.1; n];
277        let mask = vec![1u8; n];
278        let meta = ScanMetadata {
279            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
280            echo_times: vec![0.005], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
281            slice_geometry: None,
282            affine: None,
283        };
284        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Pdf, ..Default::default() };
285        let mag = vec![1.0; n];
286        let r = run_bg_removal(&field, &mask, &meta, &config, Some(&mag), &mut |_, _| {}).unwrap();
287        assert_eq!(r.local_field_ppm.len(), n);
288    }
289
290    /// PDF weights the data by the Fit_ppm_complex noise map, so it needs every echo's
291    /// magnitude and says so rather than running unweighted.
292    #[test]
293    fn pdf_requires_the_echo_magnitudes() {
294        let (nx, ny, nz) = (8, 8, 8);
295        let n = nx * ny * nz;
296        let meta = ScanMetadata {
297            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
298            echo_times: vec![0.005, 0.01], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
299            slice_geometry: None,
300            affine: None,
301        };
302        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Pdf, ..Default::default() };
303        let (field, mask) = (vec![0.1; n], vec![1u8; n]);
304        for mag in [None, Some(vec![1.0; n])] {
305            match run_bg_removal(&field, &mask, &meta, &config, mag.as_deref(), &mut |_, _| {}) {
306                Err(PipelineError::InvalidInput(msg)) => {
307                    assert!(msg.contains("PDF needs the magnitude of every echo (2 volumes"), "{msg}");
308                }
309                Err(other) => panic!("{other:?}"),
310                Ok(_) => panic!("PDF ran without the echo magnitudes"),
311            }
312        }
313        assert!(run_bg_removal(&field, &mask, &meta, &config, Some(&vec![1.0; 2 * n]), &mut |_, _| {}).is_ok());
314    }
315
316    /// The dispatcher checks the echo times with `fieldmap::echo_spacing_problem`, so a train
317    /// that is even up to sidecar rounding runs and an uneven one is refused before any work.
318    #[test]
319    fn pdf_checks_echo_spacing_with_the_public_rule() {
320        let (nx, ny, nz) = (8, 8, 8);
321        let n = nx * ny * nz;
322        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Pdf, ..Default::default() };
323        let (field, mask) = (vec![0.1; n], vec![1u8; n]);
324        for (te, ok) in [(vec![4.92e-3, 9.84e-3, 14.77e-3, 19.68e-3], true), (vec![5e-3, 10e-3, 16e-3, 21e-3], false)] {
325            assert_eq!(crate::fieldmap::echo_spacing_problem(&te).is_none(), ok);
326            let meta = ScanMetadata {
327                dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
328                echo_times: te.clone(), field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
329                slice_geometry: None,
330                affine: None,
331            };
332            let r = run_bg_removal(&field, &mask, &meta, &config, Some(&vec![1.0; te.len() * n]), &mut |_, _| {});
333            match r {
334                Ok(_) => assert!(ok, "{te:?} ran"),
335                Err(PipelineError::InvalidInput(msg)) => assert!(!ok && msg.contains("evenly spaced"), "{msg}"),
336                Err(other) => panic!("{other:?}"),
337            }
338        }
339    }
340
341    /// The noise map reaches PDF in the field's unit: the pipeline gives the same local field as
342    /// `pdf` called with `Fit_ppm_complex`'s N_std divided by the ppm-to-radian factor.
343    #[test]
344    fn pdf_gets_the_noise_map_in_ppm() {
345        let (nx, ny, nz) = (12, 12, 12);
346        let n = nx * ny * nz;
347        let mut field = vec![0.0f64; n];
348        let mut mask = vec![0u8; n];
349        let mut mag = vec![0.0f64; 3 * n];
350        for k in 0..nz {
351            for j in 0..ny {
352                for i in 0..nx {
353                    let idx = i + j * nx + k * nx * ny;
354                    let (x, y, z) = (i as f64 - 5.5, j as f64 - 5.5, k as f64 - 5.5);
355                    mask[idx] = u8::from(x * x + y * y + z * z < 16.0);
356                    field[idx] = 0.01 * (0.3 * x - 0.7 * z + 0.02 * x * y * z);
357                    for e in 0..3 {
358                        mag[idx + e * n] = (1.0 + 0.1 * x.abs()) * (1.0 - 0.2 * e as f64);
359                    }
360                }
361            }
362        }
363        let meta = ScanMetadata {
364            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
365            echo_times: vec![0.004, 0.008, 0.012], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
366            slice_geometry: None,
367            affine: None,
368        };
369        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Pdf, ..Default::default() };
370        let got = run_bg_removal(&field, &mask, &meta, &config, Some(&mag), &mut |_, _| {}).unwrap();
371        let rad_per_ppm = 2.0 * std::f64::consts::PI * crate::constants::GAMMA_BAR_HZ_PER_T * 3.0 * 0.004 * 1e-6;
372        let n_std: Vec<f64> = crate::utils::fit_ppm_complex_noise(&mag, 3).iter().map(|v| v / rad_per_ppm).collect();
373        let want = crate::bgremove::pdf(&field, &n_std, &mask, &meta.grid(), (0.0, 0.0, 1.0), &config.pdf, |_, _| {}).unwrap();
374        assert_eq!(got.local_field_ppm, want);
375    }
376
377    /// PDF is the one background remover that builds a dipole kernel, so the
378    /// pipeline must hand it `metadata.b0_direction` rather than assuming `+z`.
379    /// Hardcoding `(0.0, 0.0, 1.0)` makes both calls below identical.
380    #[test]
381    fn test_bg_removal_pdf_uses_metadata_b0_direction() {
382        let (nx, ny, nz) = (12, 12, 12);
383        let n = nx * ny * nz;
384
385        // Interior mask with real background voxels around it, so PDF has
386        // background sources to project onto, plus a spatially varying field.
387        let mut field = vec![0.0f64; n];
388        let mut mask = vec![0u8; n];
389        for k in 0..nz {
390            for j in 0..ny {
391                for i in 0..nx {
392                    let idx = i + j * nx + k * nx * ny;
393                    let (x, y, z) = (i as f64 - 5.5, j as f64 - 5.5, k as f64 - 5.5);
394                    if x * x + y * y + z * z < 9.0 {
395                        mask[idx] = 1;
396                    }
397                    field[idx] = 0.01 * (0.3 * x + 0.5 * y - 0.7 * z + 0.02 * x * y * z);
398                }
399            }
400        }
401
402        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Pdf, ..Default::default() };
403        let run = |bdir: (f64, f64, f64)| {
404            let meta = ScanMetadata {
405                dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
406                echo_times: vec![0.005], field_strength: 3.0, b0_direction: bdir,
407                slice_geometry: None,
408                affine: None,
409            };
410            // high SNR, so PDF's noise-based tolerance lets CG run
411            let mag = vec![1e4; n];
412            run_bg_removal(&field, &mask, &meta, &config, Some(&mag), &mut |_, _| {}).unwrap().local_field_ppm
413        };
414
415        let axial = run((0.0, 0.0, 1.0));
416        // ~30 degrees off +z in the y-z plane.
417        let theta = std::f64::consts::FRAC_PI_6;
418        let oblique = run((0.0, theta.sin(), theta.cos()));
419
420        assert_eq!(axial.len(), n);
421        assert_eq!(oblique.len(), n);
422
423        let max_diff = axial
424            .iter()
425            .zip(&oblique)
426            .map(|(a, b)| (a - b).abs())
427            .fold(0.0f64, f64::max);
428        assert!(
429            max_diff > 1e-9,
430            "oblique B0 gave the same local field as axial (max diff {max_diff:e}); \
431             b0_direction is not reaching the PDF dipole kernel",
432        );
433    }
434
435    #[test]
436    fn test_bg_removal_dispatches_lbv() {
437        let (nx, ny, nz) = (8, 8, 8);
438        let n = nx * ny * nz;
439        let field = vec![0.1; n];
440        let mask = vec![1u8; n];
441        let meta = ScanMetadata {
442            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
443            echo_times: vec![0.005], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
444            slice_geometry: None,
445            affine: None,
446        };
447        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Lbv, ..Default::default() };
448        let r = run_bg_removal(&field, &mask, &meta, &config, None, &mut |_, _| {}).unwrap();
449        assert_eq!(r.local_field_ppm.len(), n);
450    }
451
452    #[test]
453    fn test_bg_removal_dispatches_sharp() {
454        let (nx, ny, nz) = (8, 8, 8);
455        let n = nx * ny * nz;
456        let field = vec![0.1; n];
457        let mask = vec![1u8; n];
458        let meta = ScanMetadata {
459            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
460            echo_times: vec![0.005], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
461            slice_geometry: None,
462            affine: None,
463        };
464        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Sharp, ..Default::default() };
465        let r = run_bg_removal(&field, &mask, &meta, &config, None, &mut |_, _| {}).unwrap();
466        assert_eq!(r.local_field_ppm.len(), n);
467    }
468
469    #[test]
470    fn test_bg_removal_dispatches_ismv() {
471        let (nx, ny, nz) = (8, 8, 8);
472        let n = nx * ny * nz;
473        let field = vec![0.1; n];
474        let mask = vec![1u8; n];
475        let meta = ScanMetadata {
476            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
477            echo_times: vec![0.005], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
478            slice_geometry: None,
479            affine: None,
480        };
481        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Ismv, ..Default::default() };
482        let r = run_bg_removal(&field, &mask, &meta, &config, None, &mut |_, _| {}).unwrap();
483        assert_eq!(r.local_field_ppm.len(), n);
484    }
485
486    #[test]
487    fn test_bg_removal_dispatches_resharp() {
488        let (nx, ny, nz) = (8, 8, 8);
489        let n = nx * ny * nz;
490        let field = vec![0.1; n];
491        let mask = vec![1u8; n];
492        let meta = ScanMetadata {
493            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
494            echo_times: vec![0.005], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
495            slice_geometry: None,
496            affine: None,
497        };
498        let config = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Resharp, ..Default::default() };
499        let r = run_bg_removal(&field, &mask, &meta, &config, None, &mut |_, _| {}).unwrap();
500        assert_eq!(r.local_field_ppm.len(), n);
501    }
502
503    /// iHARPERELLA reads its own params, not HARPERELLA's (#169).
504    #[test]
505    fn iharperella_uses_its_own_params() {
506        let (nx, ny, nz) = (16, 16, 16);
507        let n = nx * ny * nz;
508        let meta = ScanMetadata {
509            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0),
510            echo_times: vec![0.005], field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0),
511            slice_geometry: None,
512            affine: None,
513        };
514        // a blob inside a 10-voxel cube: big enough to leave an interior after the ~3-voxel erosion
515        let field: Vec<f64> = (0..n)
516            .map(|i| {
517                let (x, y, z) = ((i % nx) as f64 - 7.5, ((i / nx) % ny) as f64 - 7.5, (i / (nx * ny)) as f64 - 7.5);
518                (-(x * x + y * y + z * z) / 10.0).exp() + 0.01 * x
519            })
520            .collect();
521        let mask: Vec<u8> = (0..n)
522            .map(|i| {
523                let (x, y, z) = (i % nx, (i / nx) % ny, i / (nx * ny));
524                u8::from((3..13).contains(&x) && (3..13).contains(&y) && (3..13).contains(&z))
525            })
526            .collect();
527        let run = |config: &BgRemovalConfig| {
528            run_bg_removal(&field, &mask, &meta, config, None, &mut |_, _| {}).unwrap().local_field_ppm
529        };
530        let base = BgRemovalConfig { algorithm: BgRemovalAlgorithm::Iharperella, ..Default::default() };
531        let reference = run(&base);
532        assert!(reference.iter().any(|v| v.abs() > 1e-6), "nothing came through");
533        // HARPERELLA's settings must not reach iHARPERELLA
534        let mut other = base.clone();
535        other.harperella = crate::bgremove::HarperellaParams { radius: 2.0, max_iter: 1, tol: 0.5, pad: [3, 3, 3] };
536        assert_eq!(run(&other), reference);
537        // its own settings must
538        let mut own = base.clone();
539        own.iharperella.max_iter = 1;
540        assert_ne!(run(&own), reference);
541    }
542
543    #[test]
544    fn test_bg_removal_validates_dims() {
545        let meta = ScanMetadata {
546            dims: (4, 4, 4),
547            voxel_size: (1.0, 1.0, 1.0),
548            echo_times: vec![0.005],
549            field_strength: 3.0,
550            b0_direction: (0.0, 0.0, 1.0),
551            slice_geometry: None,
552            affine: None,
553        };
554        let field = vec![0.0; 32]; // wrong size, should be 64
555        let mask = vec![1u8; 64];
556
557        let result = run_bg_removal(
558            &field, &mask, &meta, &BgRemovalConfig::default(), None,
559            &mut |_, _| {},
560        );
561        assert!(result.is_err());
562    }
563}
qsm_core/pipeline/bg_removal.rs

qsm_core/pipeline/
bg_removal.rs
