Module msmv 

mSMV (maximum spherical mean value) filtering, the field step of MEDI-mSMV

mSMV removes residual background field near the brain boundary, the main source of shadow artifacts, without eroding the mask (Roberts et al. 2024). It is not a background field removal step on its own: the authors run it inside MEDI, as MEDI_L1(..., 'msmv', radius), on the local field from a conventional background removal (PDF, V-SHARP or LBV in the paper), and the inversion then fits the SMV-filtered field with the SMV-filtered dipole kernel (δ − S_r1) d (the paper’s Eq. 12, “MEDI-mSMV”). Its output is a high-passed field, (δ − S_r1) applied to the local field plus a boundary correction, so inverting it with the plain dipole kernel d, as a background removal output would be, is the wrong forward model. In QSM.rs it is therefore an option of medi: MediParams::msmv. This module holds the filter itself, a port of msmv.m and kernel_lim.m, for MEDI and for parity checks.

The authors’ method, as their code runs it

Reference code: agr78/mSMV (commit 42a1f18), prerelease/MEDI_toolbox-mSMV100125.zip, the newest MEDI toolbox with mSMV: functions/msmv.m, kernel_lim.m, MaskErode.m, SMV.m, sphere_kernel.m, MEDI_L1.m and parse_QSM_input.m. MEDI_L1(..., 'msmv', r) (the 'msmv' keyword takes the radius, like 'smv', and turns SMV on):

1. reads R2s from the input file and errors if it is missing;
2. does not erode the mask (smv_shrink_mask is forced off) and calls RDF = msmv(RDF, Mask, R2s, voxel_size, r, B0_mag, prefilter, maxk, vessel_radius);
3. inverts with D = (1 − SphereK).*D and the SMV-propagated noise tempn = sqrt(SMV(tempn.^2) + tempn.^2), exactly as MEDI-SMV does, over the full mask.

msmv.m itself, in MEDI’s units (the field in radians over the echo spacing delta_TE):

• SphereK = single(sphere_kernel(..., radius)), the partial-volume sphere of MEDI; Mne = SMV(Mask, SphereK) > 0.999, Me = Mask − Mne (the boundary band of width ~radius).
• Prefilter (default on): RDF_s = Mask.*(RDF − SMV(RDF, SphereK)) (the paper’s Eq. 4). With prefilter = 0 the input is used as it is (not masked); the README says this is “typically used with SHARP, RESHARP, VSHARP”, i.e. background removals that already applied an SMV. The paper ran it with the prefilter on after PDF, V-SHARP and LBV.
• Threshold t = kernel_lim(RDF, ...) on the unfiltered input: starting from a rendered sphere of half the smallest voxel, the radius grows in 1e-3 mm steps and t is the first max|Mask.*(RDF − SMV(RDF, K))| that reaches tmin = 0.01·B0_mag/3 rad. At half a voxel the rendered sphere is the centre voxel alone (the residual is 0); the residual then grows in steps as sub-voxel samples of the neighbours enter the sphere, so t is either just above tmin (smooth fields) or the first nonzero step, which strong local fields put well above it (0.138 rad against tmin = 0.023 rad on the 7 T qsm-forward phantom after PDF): the code’s form of the paper’s max(t_min, lim_{r'→0} ...).
• Vessel mask: Mask_ev = Mask − MaskErode(Mask, radius + 1), Mv = imbinarize(fibermetric(Mask_ev.*R2s, 1:vessel_radius, 'ObjectPolarity', 'bright'), 0) with vessel_radius = max(15, 8·max(voxel_size)) by default. fibermetric takes the thicknesses in voxels; the threshold 0 keeps every voxel with any bright-tubular response at any of the 15+ scales (the Frangi sign condition λ2, λ3 < 0), not only veins: 70 % of the band on the qsm-forward phantom (in MATLAB as here), so the correction acts on the rest.
• Mb = imbinarize(|Me.*RDF_s|, t) & ~Mv; while sum(Mb)/sum(Mask) > 1e-6: recompute Mb from the current field, RDF_s = Mask.*(RDF_s − SMV(Mb.*RDF_s, r2)) with the minimum radius r2 = min(voxel_size)/2 + 0.05 mm, stop after maxk − 1 passes (default maxk 5, so at most 4 passes). The loop test uses the Mb of the previous pass.

Where the paper and the code differ (the code is followed)

• Threshold. The paper (Eq. 7) defines t = max(t_min, lim_{r'→0} max_M |b_L0 − S_r' b_L0|) on the prefiltered field with t_min = 0.3 Hz at 3 T, “scaled accordingly for other field strengths”. kernel_lim.m uses the unfiltered input field and tmin = 0.01·B0/3 radians over delta_TE, which is 0.3 Hz only for delta_TE ≈ 5.3 ms (0.49 Hz at the paper’s in-vivo 3.25 ms). The field strength comes from B0_mag, which MEDI_L1 reads from the input file and otherwise leaves at msmv.m’s default of 3 T. QSM.rs passes the scan’s field strength, as the paper’s scaling asks; at 3 T this is the code’s value.
• Iterations. The paper’s i_max = 5; the code’s maxk = 5 gives at most 4 passes.
• Vessel mask. The paper says “a Frangi filter … followed by thresholding” without the scales or threshold; the code’s are above.

Ported here

msmv follows msmv.m line for line, with MEDI’s rendered sphere_kernel for every SMV (the mSMV release’s version, which takes the real part of the kernel’s FFT), kernel_lim.m as written, MaskErode.m, and MATLAB’s fibermetric (R2023b, default StructureSensitivity 0.01 for double input, ported in fibermetric). R2* is a required input, as it is for MEDI_L1‘s mSMV path; MEDI’s pipeline (and the authors’ SEPIA patch) computes it with ARLO from the echo magnitudes, which is what medi does (arlo_unclipped: arlo.m exactly, negative values kept, since the vessel filter sees them).

Differences from MATLAB, none of which changes the method:

• Arithmetic is in double. msmv.m stores SphereK as single, which makes SMV(RDF, SphereK) and, with the prefilter on, the whole filtered field single precision; here the kernel is rounded to single as MATLAB stores it but the field stays double. Voxels within single rounding of the threshold t can be classified differently.
• SMVs with small kernels (kernel_lim’s growing sphere and the minimum-radius r2 kernel) are applied as a direct circular convolution with the same kernel instead of by FFT (same operator, equal to rounding; kernel_lim runs dozens to hundreds of them).
• fibermetric is evaluated only on the boundary band Me, the only place Mv is used (Mb is zero outside Me), so MsmvOutput::vessel_mask is false outside it.
• kernel_lim.m loops forever if the residual never reaches tmin (a field below tmin/2 everywhere); here the threshold is then tmin, which flags no voxel (see kernel_lim).

Parity with the authors’ MATLAB

Against msmv.m and MEDI_L1(..., 'msmv', 5) (R2023b, the prerelease above):

• Small phantom (tests/msmv_medi_parity.rs, 32×36×28, anisotropic voxels): t identical, Me and the vessel mask identical (5615 of 6272 band voxels), the filtered field to 5e-8 rad (relative 1e-7; single rounding) with the prefilter and to 1e-16 without; one MEDI-mSMV Gauss-Newton step to 3e-7 ppm (relative 6e-7).
• qsm-forward phantom after MATLAB’s PDF.m (164×206×206 once padded, 7 T): t 0.138 rad equal to 1e-13 relative, Me (321,821 voxels) and Mv (225,330) identical, the field to 1.5e-7 relative with the prefilter and 3e-16 without. MEDI-mSMV with MEDI’s defaults agrees to 4e-3 relative L2 (1e-2 ppm max) and one step solved to cg_tol 1e-4 to 2.4e-4 (4e-4 ppm max), as close as MEDI-SMV gets in the same setup (1.9e-4): MATLAB’s MEDI runs in single precision, and at the default cg_tol 0.01 CG can stop one iteration apart, which moves χ by ~1e-3 (see the parity test’s docs).

Measured

QSM.rs MEDI with its defaults (λ 1000, MERIT, MEDI+0), after QSM.rs PDF, V-SHARP and LBV, comparing MEDI without SMV (the paper’s Eq. 14), MEDI-SMV (Eq. 13, mask eroded by 5 mm) and MEDI-mSMV (Eq. 12). “GM variance” is the paper’s shadow score, the χ variance in grey matter (×1e-4 ppm²; lower means less shadow, down to the true variation).

qsm-forward phantom (7 T, TE 4/12/20/28 ms, 1 mm; ground-truth GM variance 1.61), XSIM and NRMSE over the brain mask (MEDI-SMV is zero on the third of the mask it erodes):

field	BFR	MEDI XSIM / NRMSE / GM var	MEDI-SMV	MEDI-mSMV
noise-free	PDF	0.597 / 65 % / 3.64	0.563 / 63 % / 2.37	0.561 / 65 % / 2.92
noise-free	V-SHARP	0.548 / 63 % / 2.88	0.510 / 66 % / 2.46	0.547 / 63 % / 2.42
noise-free	LBV	0.612 / 59 % / 3.32	0.529 / 65 % / 2.40	0.606 / 59 % / 2.46
SNR 50	PDF	0.305 / 103 % / 9.20	0.352 / 76 % / 3.55	0.305 / 90 % / 5.58
SNR 50	V-SHARP	0.321 / 80 % / 3.87	0.335 / 77 % / 3.12	0.333 / 78 % / 3.07
SNR 50	LBV	0.296 / 98 % / 9.84	0.336 / 77 % / 3.48	0.358 / 78 % / 3.52

In vivo, four 3 T multi-echo scans (EMTP travelling heads in QSM-CI: Philips Achieva and Ingenia, Siemens Skyra and Prisma), GM variance in the SynthSeg cortex, mean of the four (MEDI-SMV covers only 55–67 % of the cortex; MEDI and MEDI-mSMV 93–100 %):

BFR	MEDI	MEDI-SMV	MEDI-mSMV
PDF	15.1	6.0	7.3
V-SHARP	7.0	5.5	4.5
LBV	17.0	5.8	6.2

MEDI-mSMV lowered the GM variance against plain MEDI for every background removal and every scan (PDF by half, LBV to a third, V-SHARP by a third, the paper’s direction and size), to about MEDI-SMV’s level without its erosion, and lowered the χ spread in the 5 mm edge band. On the phantom it did not raise XSIM in general: noise-free it was level with MEDI after V-SHARP and LBV and lower after PDF (0.561 vs 0.597); at SNR 50 it was level after PDF and higher after V-SHARP (0.333 vs 0.321) and LBV (0.358 vs 0.296), with lower NRMSE. It helps most where the local field carries residual background at the edge (PDF, LBV, noisy data) and little where V-SHARP’s erosion already removed it.

Prefilter off (msmv.m‘s prefilter = 0) breaks the forward model in MEDI_L1: the field is then the background removal’s local field, unfiltered, but MEDI_L1 still fits it with (δ − S) d. On the phantom that gave NRMSE 460 % (PDF), 280 % (V-SHARP) and 740 % (LBV), and the authors’ MATLAB gives the same (481 % for PDF + MEDI_L1(..., 'msmv', 5) with prefilter = 0). It is only consistent for a field that is already high-passed by the same SMV; V-SHARP’s deconvolved output is not. The option is kept because the authors’ code has it; keep the default.

Use with other inversions

Only MEDI is offered. The authors made and tested MEDI-mSMV only (with PDF, V-SHARP and LBV before it). An inversion whose forward model is an explicit dipole convolution could in principle fit the mSMV field with (δ − S_r1) d instead of d, but none was made or evaluated by the authors, and the result would be a new method, so QSM.rs does not offer one.

Reference: Roberts, A.G., Romano, D.J., Şişman, M., Dimov, A.V., Nguyen, T.D., Kovanlikaya, I., Gauthier, S.A., Wang, Y., Spincemaille, P. (2024). “Maximum spherical mean value filtering for whole-brain QSM.” Magnetic Resonance in Medicine, 91(4):1586-1597. https://doi.org/10.1002/mrm.29963

Reference implementation: https://github.com/agr78/mSMV (prerelease/MEDI_toolbox-mSMV100125.zip)

Structs

MsmvOutput

msmv’s output with the intermediate results, for parity checks and evaluation

MsmvParams

mSMV options (msmv.m’s prefilter, maxk and vessel_radius, with its defaults). The prefilter radius is MEDI’s SMV radius and the field strength the scan’s (see MediParams::msmv).

Enums

SparseKernel 🔒

SMV(x, matrix_size, voxel_size, r) for a small r: MEDI’s rendered sphere applied as a circular convolution with real(fftn(...))’s kernel, i.e. the sphere symmetrised ((S(o) + S(−o))/2), which is the sphere itself on even grids. Large spheres go through the FFT.

Functions

kernel_lim

kernel_lim.m: the shadow threshold t in radians (see the module docs).

msmv

msmv.m: mSMV-filter a local field (see the module docs).