Module heidi 

HEIDI — Homogeneity Enabled Incremental Dipole Inversion

HEIDI splits k-space at the dipole cone and treats the two halves differently:

• Well-conditioned coefficients (|D(k)| > threshold) are taken straight from an already-computed χ map — normally the minimally regularised LSQR solution of [super::lsqr_qsm] — and held as exact constraints.
• Ill-conditioned coefficients inside the cone, where direct division amplifies noise without bound, are recovered by minimising a weighted total variation subject to those constraints.

minimise  Σ_v ‖ (M_x ∂_x χ, M_y ∂_y χ, M_z ∂_z χ) ‖   subject to   (F χ)|cone = (F χ_init)|cone

The “homogeneity” of the name is the weighting M. Each direction gets its own 0/1 mask, built by thresholding the gradient of the field map: where the field is locally flat the susceptibility is expected to be flat too, and TV is applied; across an anatomical edge the mask opens and the gradient goes unpenalised. A Laplacian test additionally opens the mask wherever the field curves sharply. Voxels that are an edge in every direction get the weight HeidiParams::gradient_mask_floor, 0 in the shipped code (see that field).

The constrained minimisation is solved with NESTA (Nesterov’s smoothing with continuation): the TV term is smoothed by a parameter μ that is annealed down over several outer steps, each of which runs a fixed number of accelerated gradient iterations. Because the data constraint is an orthogonal projection (A·Aᵀ = I for a restricted FFT), the projection step is exact and costs one FFT round trip.

Ported from the authors’ code, qsm_heidi’s inversion_heidi.m and computesusceptibility_heidi.m (https://gitlab.com/R01NS114227/qsm_heidi), run after its inversion_lsqr.m ([super::lsqr_qsm]):

• the seed is the LSQR solution on LSQR’s grid (the mask’s bounding box with a 40-voxel anti-aliasing frame, odd extents), unmasked, and HEIDI runs on that same grid;
• the cone is |D| > 0.1 (the driver’s PostProcCone.threshold) with D the spatial-domain unit dipole response (super::lsqr_qsm::spatial_dipole_kernel);
• NESTA runs with the driver’s tolerance, machine epsilon, so every continuation step uses all of its inner iterations (8 × 500);
• the homogeneity masks are built from the phase denoised by ITK 3.20’s gradient anisotropic diffusion (the executable diffusionimagefilter.m runs, 5 iterations, time step 0.05, conductance 1), in the phase’s units. That executable works in float and sums the mean squared gradient that sets its conductance scale in float, which underestimates it (by about 1 % on a whole brain); QSM.rs computes it in double by default, and [AnisotropicDiffusionParams::reference_itk] reproduces the executable bit for bit. The two give the same masks but for the voxels at the gradient thresholds (about 0.07 % on an in-vivo brain, QSM-CI parity case s005001-cbic-prisma-bridge-run1).

Units: the field is in ppm and χ comes back in ppm. The gradient and Laplacian thresholds below are therefore in ppm/mm and ppm/voxel², independent of echo time and field strength — upstream divides the phase by TE[ms]·B₀[T] before thresholding, which is the same quantity. The cone fill is not scale-free: NESTA’s smoothing ends at an absolute μ (machine epsilon) and upstream’s χ is in units of the phase (ppm · γ·TE·B₀·1e-6, γ in rad/s/T), so the seed is converted to those units before NESTA and back after, with HeidiParams::b0 and HeidiParams::echo_times (scan parameters the pipeline fills from metadata) at the echo time HeidiParams::te picks.

Reference: Schweser, F., Sommer, K., Deistung, A., Reichenbach, J.R. (2012). “Quantitative susceptibility mapping for investigating subtle susceptibility variations in the human brain.” NeuroImage, 62(3):2083-2100. https://doi.org/10.1016/j.neuroimage.2012.05.067

NESTA: Becker, S., Bobin, J., Candès, E.J. (2011). “NESTA: A fast and accurate first-order method for sparse recovery.” SIAM Journal on Imaging Sciences, 4(1):1-39. https://doi.org/10.1137/090756855

Any seed: HEIDI as streaking reduction

The cone fill does not depend on where the seed came from. The authors’ computesusceptibility_heidi(suscInterim, …) (qsm_heidi 4592631) takes the initial χ map from outside and only refills its cone (the changelog calls the step “PostProcCone”), and SEPIA’s HEIDI add-on (Salman, Schweser and Chan, 2025) offers it after any dipole inversion as “Streaking reduction by HEIDI”: Wrapper_QSM_HEIDI4all(localField, chiRaw, …), called from wrapper/QSMMacro.m on the method’s result. heidi_cone_fill is that step:

• Grid: as Wrapper_QSM_HEIDI4all, the field, mask and seed are cropped to the mask’s bounding box with the 40-voxel anti-aliasing frame (the grid heidi runs on), and χ is masked at the end (chi .* mask). The seed is the method’s (masked) result as it is.
• Units: the qsm_heidi pipeline’s, phase at the chosen echo time, for every seed, so heidi and heidi_cone_fill fill the cone the same way. Wrapper_QSM_HEIDI4all instead passes the field and the seed in Hz/2π, a factor (2π)²·TE (≈ 0.8 at TE = 20 ms) from those; the factor moves the gradient thresholds and NESTA’s absolute final μ only.
• Seed suitability. The paper takes the well-conditioned region from a minimally regularised inversion and only the cone from the TV prior. A seed whose region outside the cone is close to data-consistent fits that: LSQR (the reference, which heidi runs), and TKD, TSVD, small-λ Tikhonov, NDI and RTS (RTS already splits k-space the same way). A regularised seed (MEDI, TV and the FANSI family, NLTV, iLSQR, STAR, AMP-PE, the deep-learning methods) keeps its regularisation outside the cone, as only the cone is replaced. The pipeline publishes this as the heidi.seed-* rules (crate::pipeline::rules::RULES).
• Truncated seeds. TKD and TSVD zero or clip |D| < threshold; with that threshold above the cone threshold, truncated coefficients would be kept as if they were data. Upstream stops with an error there (“DipoleFilter parameter must not be larger that Cone threshold”, computesusceptibility_heidi.m L922-923), and so does the pipeline (crate::pipeline::heidi_seed_problem).
• Not offered: methods without a local field (the total-field and phase methods: the homogeneity weights come from the local field), χ-separation and multi-orientation reconstructions (upstream refuses more than one field: “You cannot use HEIDI with COSMOS”).

Reference pipeline

Same inputs as super::lsqr_qsm (qsm_heidi commit 4592631): a SHARP-corrected local field (makefile L14-20, background_correction_sharp.mat), the magnitude, the LSQR result, and the mask mask_tissue .* mask_phase .* phsInput.mask (inversion_heidi.m). The echo time for the phase units is the last one (phase2susceptibilityfactor.m L70-71), as the pipeline sets.

Pipeline rules

The facts above that a configuration can depart from are published as rules in crate::pipeline::rules::RULES (each with its strength, basis, reason and source): lsqr.bg-sharp, heidi.te-last.

Structs

HeidiParams

Parameters for heidi.

HomogeneityMasks

The three per-direction homogeneity weights, as used by the TV term.

TvWorkspace 🔒

Scratch buffers for the NESTA loop, so the inner iterations allocate nothing.

Constants

DEFAULT_GRADIENT_THRESHOLD_PPM

Upstream DEFAULT_GRADIENTMASK_THRESHOLD, converted to ppm/mm.

DEFAULT_LAPLACIAN_THRESHOLD_PPM

Upstream DEFAULT_GRADIENTMASK_LAPLACIANTRHESHOLD, converted to ppm/voxel².

GAMMA 🔒

Gyromagnetic ratio of the proton, in rad/s/T: the reference implementation’s value (getconstants.m), not crate::constants. Keep it.

UPSTREAM_PPM_SCALE 🔒

Upstream expresses its mask thresholds per unit of phase / (TE[ms] · B₀[T]). One ppm of field is γ · 1e-3 · 1e-6 of that, so dividing an upstream threshold by this constant converts it to ppm.

Functions

cone_fill_iterations

Progress steps of the cone fill: heidi_cone_fill, and heidi’s second stage.

cone_mask 🔒

|D(k)| > threshold: the k-space region the dipole model conditions well, with D the spatial-domain unit dipole response, as upstream (generateunitdipoleresponse).

discrete_gradient 🔒

Allocating wrapper around discrete_gradient_into.

discrete_gradient_into 🔒

Backward difference along axis, scaled by voxel size, zero on the leading face.

heidi

HEIDI dipole inversion: LSQR ([super::lsqr_qsm]) followed by HEIDI’s cone fill, as qsm_heidi’s inversion_lsqr.m then inversion_heidi.m.

heidi_cone_fill

HEIDI’s cone fill on any seed χ map: streaking reduction after a dipole inversion, as the authors’ computesusceptibility_heidi(suscInterim, …) and SEPIA’s Wrapper_QSM_HEIDI4all.

heidi_on_grid 🔒

HEIDI’s cone fill on a grid with odd extents (upstream’s), from the seed chi_init (unmasked). Returns χ unmasked.

homogeneity_masks

Build the three per-direction homogeneity weights from the field map.

nesta_cone_fill 🔒

NESTA with continuation: fill the dipole cone under a weighted-TV prior while holding the well-conditioned coefficients of chi_init fixed.

threshold_gradient 🔒

1.0 where |∂_axis field| < threshold, else 0.0.

tv_objective 🔒

Smoothed weighted-TV objective and its gradient.