Module lsqr_qsm 

Minimally regularised LSQR dipole inversion (Schweser 2010)

Solves the dipole inversion as one large least-squares system and leans on LSQR’s minimum-norm property instead of an explicit regularisation term — hence “minimally regularised”.

Two extra sets of unknowns make that work:

• a residual field r, added to the data equation with weight w, which absorbs field contributions the dipole model cannot explain (unmodelled background, noise, non-local sources). Because LSQR converges to the minimum-norm solution, ‖r‖ is implicitly penalised, and w sets how cheaply the fit may spend field on r rather than on χ.
• a global offset c, one scalar, absorbing any constant field shift left over from background removal.

With χ over the whole (cropped and padded) grid, r and the data over the mask only:

minimise ‖ W · ( M[ F⁻¹(D · F χ) + c ] + w·r − b ) ‖₂   over  [χ; r; c]

W is the optional magnitude image (used as is, as upstream does) and M restricts to the mask.

Ported from the authors’ code, inversion_lsqr.m and computesusceptibility_lsqr.m of qsm_heidi (https://gitlab.com/R01NS114227/qsm_heidi), with its settings:

• the grid: field, mask and magnitude are cropped to the mask’s bounding box, each odd extent grown by one voxel, zero-padded by 20 voxels per side and the last voxel dropped (zealouscrop, prepareforconvolution with its 40-voxel anti-aliasing frame), so every extent is odd and the crop starts 20 voxels in;
• the dipole kernel: the unit dipole response (3cos²θ − 1)/(4πr³) sampled in the spatial domain on that grid and Fourier transformed (generateunitdipoleresponse with isFourierDomainFormula = false), see spatial_dipole_kernel;
• w = 0.2/9.4 · B0, the global offset, tolerance 1e-5 and at most 400 iterations of MATLAB’s lsqr (its stopping rules, via super::ilsqr::lsqr_matlab).

Units

The problem is linear and homogeneous in the field: scaling the field scales χ, r and c alike and leaves LSQR’s iterates and stopping tests unchanged. Upstream solves for the phase in radians and converts χ to ppm at the end; here the field goes in ppm and χ comes out in ppm, which is the same solution (on a 3T scan the two agree to rounding).

This map is what HEIDI consumes: [super::heidi] keeps its well-conditioned k-space coefficients and re-derives the rest under a weighted-TV prior.

Reference: Schweser, F., Deistung, A., Lehr, B.W., Reichenbach, J.R. (2010). “Differentiation between diamagnetic and paramagnetic cerebral lesions based on magnetic susceptibility mapping.” Medical Physics, 37(10):5165-5178. https://doi.org/10.1118/1.3481505

Reference pipeline

qsm_heidi (https://gitlab.com/R01NS114227/qsm_heidi, commit 4592631): the makefile (L14-20) and def.sh feed inversion_lsqr.m a SHARP-corrected local field (background_correction_sharp.mat, phase in rad), the magnitude, and masks made by an upstream step not in the repository (create_masks.mat). The mask the inversion uses is the product of a tissue mask, a phase mask and the background-removal mask (inversion_lsqr.m, suscCompMask = mask_tissue .* mask_phase .* phsInput.mask). The SHARP settings are not in the repository. The phase-to-ppm factor uses the last echo time (phase2susceptibilityfactor.m L70-71), which does not change the LSQR solution (see Units).

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.

Structs

LsqrQsmParams

Parameters for lsqr_qsm.

Constants

ANTIALIASING_FRAME

Width of upstream’s anti-aliasing frame (DEFAULT_ANTIALIASINGFRAMEWIDTH), in voxels, split evenly over both sides of each axis.

FRAME_OFFSET 🔒

Where the crop starts on the padded grid: half the anti-aliasing frame.

RESIDUAL_WEIGHTING_PER_TESLA

Residual weighting per tesla: upstream’s 0.2 / 9.4 * B0 (“0.2 optimal in QUASAR paper”).

Functions

crop_grid 🔒

The grid upstream solves on: the mask’s bounding box (zealouscrop), each extent made even by one trailing voxel and then grown by the anti-aliasing frame less one voxel (prepareforconvolution: makeeven, 20 voxels per side, makeodd).

crop_pad 🔒

Copy the crop box of x (on dims) into a zero grid of the padded size, FRAME_OFFSET voxels in along each axis.

lsqr_qsm

Minimally regularised LSQR dipole inversion.

lsqr_qsm_core 🔒

lsqr_qsm_padded 🔒

The solve on upstream’s grid. Returns the cropped-and-padded grid, its crop box, and χ on it (unmasked), or None for an empty mask.

spatial_dipole_kernel

The dipole kernel in k-space (FFT order, real) as qsm_heidi’s generateunitdipoleresponse(dims, true, Options) builds it with isFourierDomainFormula = false and a voxel size given: the unit dipole response (3(b·r̂)² − 1)/(4π|r|³) (zero at r = 0) sampled on an isotropic grid of the smallest voxel size, 2⌊n·a/2⌋ + 1 voxels along an axis of n voxels and aspect a = vs/min(vs), Fourier transformed, cropped to the central n frequencies, real part, DC set to zero. For isotropic voxels the fine grid is the grid itself.

uncrop 🔒

Inverse of crop_pad: the crop box back in place on dims, zero elsewhere.

Type Aliases

CropBox 🔒

Per axis: (first voxel of the crop, its extent, the padded extent). See crop_grid.