Module ismv 

Iterative Spherical Mean Value (iSMV) background field removal

A harmonic (background) field equals its own spherical mean. iSMV holds the field fixed on a thin band at the edge of the mask and, inside the mask eroded by the sphere, repeatedly replaces the field by its spherical mean (SMV), F ← SMV(F)·M + (mask − M)·field. That converges to the harmonic field matching the edge band, i.e. the background; the local field is (field − F)·M. No deconvolution follows, so the sphere can be small and the mask is eroded less than by SHARP.

Which implementation (decided 2026-10-29, astewartau/QSM.rs#322)

The default, IsmvVariant::Wen2014, is the authors’ paper and code together:

• Code: the iSMV authors’ lab (Cornell, Yi Wang’s group) ships iSMV.m with its kernel tdcsmv.m in github.com/agr78/mSMV (42a1f18, code/dependencies/MEDI_functions/). This module ports its loop, kernel and erosion line by line (see iSMV.m, step by step). Nothing in that repository (or in MEDI 2024.11.26, which ships no iSMV) calls it, so the code fixes no radius, and its hard-coded tolerance, 0.001, under-converges (below).
• Paper: Wen, Zhou, Liu, Spincemaille, Wang 2014 (Methods, PMC4026360) stop “when the relative change between successive iterates falls below 5 × 10⁻⁵” and compare radii of 1 and 6 voxels, recommending neither (“much less dependent on the choice of radius”). So the defaults are the paper’s: tol 5 × 10⁻⁵ on iSMV.m’s relative change, and a radius of 6 voxels. The paper’s other radius, 1 voxel, degenerates with tdcsmv.m: its in-plane loop runs over -(n−1):(n−1) with n = ceil(r/vs₁) = 1, so the “sphere” is a 3-voxel column along the third axis.
• 6 voxels in mm is 6 × the smallest voxel dimension (IsmvParams::radius_mm): the paper’s in-vivo voxels were 0.94 × 0.94 × 2 mm, and tdcsmv.m counts the sphere’s in-plane extent in units of voxel_size(1), the in-plane (smallest) voxel size. On isotropic voxels the two readings agree. iSMV.m itself takes the radius in mm.
• No iteration cap by default (IsmvParams::max_iter None), as in iSMV.m (the paper reports none). The iteration is a contraction on the eroded mask (the SMV of every voxel of M reaches the fixed edge band or the zero outside the mask), so it always terminates for tol > 0; a cap is available.

QSM.jl’s ismv, which this module followed from 2026-10-28 (the designated implementation of Kames’s SMV methods when the authors have no code), is kept as IsmvVariant::QsmJl for the QSM.jl parity pair (IsmvParams::qsmjl); see QSM.jl’s ismv.

Why the paper’s tolerance (measured, 2026-10-28)

The stopping tolerance matters more than the radius. At tol 1e-3 (iSMV.m as shipped, or QSM.jl’s default) the iteration stops after 9–15 iterations and leaves background in the field, which shrinks deep-grey-matter contrast about threefold in vivo. CI phantom (qsm-forward, 1 mm; local field scored over the brain mask as CI scores it, and over the 941,844 voxels every setting keeps; RTS χ over those voxels) and the QSM-CI travelling heads (10 acquisitions of 4 subjects, ROMEO field, SynthSeg regions; |slope − 1| of each scan’s regional χ against the subject’s mean across scanners):

setting	phantom field, CI scoring: NRMSE % / r	common voxels: NRMSE % / r	RTS χ: NRMSE % / r	in vivo median |slope−1|	pallidum χ (ppb)
QSM.jl: 2 × max voxel, 1e-3 of the first residual, 500	60.9 / 0.800	53.3 / 0.853	68.4 / 0.740	0.111	41
iSMV.m as shipped, 2 voxels, 0.001	63.0 / 0.778	53.7 / 0.847	68.9 / 0.731	0.106	44
iSMV.m as shipped, 6 voxels, 0.001	57.9 / 0.817	36.2 / 0.934	59.9 / 0.806	0.055	97
iSMV.m, 6 voxels, 5e-5 (Wen 2014; the default)	53.0 / 0.848	25.2 / 0.968	53.6 / 0.848	0.036	129
QSM.m ismv.m: 3 mm, √eps of the first residual, 500	48.1 / 0.878	27.4 / 0.963	51.2 / 0.862	0.038	134
QSM.rs before 2026-10-28: 5 mm, 1e-6, 50 (no reference)	51.4 / 0.858	26.2 / 0.966	54.2 / 0.845	0.024	124

The default took 94 iterations on the phantom. Its cost is coverage: a 6-voxel sphere keeps fewer edge voxels (in vivo, on average 1.06 M brain voxels, against 1.25 M at 3 mm and 1.33 M at 2 voxels).

iSMV.m, step by step

0. Odd sizes (QSM.rs’s fix). Every odd axis is zero-padded by one slice at its end and the result cropped back, as SEPIA does before its FFT-based background removal (zeropad_odd_dimension(…, 'pre'), kschan0214/sepia db939cc; see super::sharp). tdcsmv.m places the kernel at f(end/2+1+k, …), which is not an integer index on an odd axis, so iSMV.m errors there; padding is the unequivocal fix, and on even sizes it is a no-op.
1. Kernel a = tdcsmv(matrix_size, voxel_size, radius): a sphere of radius mm with partial-volume weights from sub-voxel samples (tdcsmv), normalised to sum 1, fftn(fftshift(SMV)).
2. Erosion M = ifftn(fftn(Mask).*a) > .999 (MATLAB’s > compares real parts), periodic on the (padded) array. BV = (Mask − M).*iFreq.
3. Loop from F1 = iFreq (the whole input array, also outside the mask, as iSMV.m): F2 = ifftn(fftn(F1).*a).*M + BV; d = norm(F1(:)−F2(:))/norm(F1(:)); F1 = F2; while d ≥ tol (iSMV.m: while (d>=.001); here IsmvParams::tol, default 5e-5). The iterates are kept complex, as MATLAB’s are (their imaginary part is rounding, ~1e-17). A NaN d (an all-zero input) stops the loop, as in MATLAB.
4. Output L = (iFreq − F2).*M (its real part; iSMV.m returns the complex array) and M. The local field is zero outside the eroded mask M, not only outside the input mask.

tdcsmv.m has two more places where it cannot run as written; both changes leave every input it handles unchanged:

• It uses voxel_size(1) for the second axis too (sub-voxel grid, offsets and extent), so on voxels with vs₁ ≠ vs₂ its sphere has the wrong size along that axis. Here the second axis uses voxel_size(2); with vs₁ = vs₂ (every case tdcsmv.m is right for) nothing changes.
• It picks the sub-voxel factor w as the smallest integer with vs₁·w integral in single precision (while (int16(g*w)~=single(g*w))), unless vs₁ = vs₃ (then 10), then doubles it. For a voxel size with no such w (a header value like 0.93750006) that loop never ends in MATLAB (int16 saturates at 32767). Here the search stops at TDCSMV_MAX_SUBVOXEL and falls back to w = 10, tdcsmv.m’s own isotropic choice.

A sphere that does not fit the array (MATLAB would grow f and then fail on the size mismatch) wraps around it periodically.

QSM.jl’s ismv (IsmvVariant::QsmJl)

A port of QSM.jl’s ismv (kamesy/QSM.jl fde4b92, src/bgremove/ismv.jl), with its defaults from IsmvParams::qsmjl: radius r = 2·max(vsz) mm, stop when the residual ‖fᵢ − fᵢ₋₁‖ falls to tol = 1e-3 of the first one, at most maxit = 500 iterations. Its ball is the voxel centres within r (no partial volume), its erosion S ∗ mask > 1 − √eps(Float64), and it returns mask·(field − f) over the input mask with the eroded mask.

QSM.jl crops to the mask’s bounding box and zero-pads to a fast FFT size before filtering; this port filters on the whole array, periodically. iSMV only applies the SMV kernel forwards, at voxels whose ball lies inside the mask, so the two are the same unless the mask reaches two opposite faces of the array (then a ball can wrap around here and count as fitting). QSM.jl’s edge test is S ∗ mask > 1 − √eps(Float64), this port’s > 1 − 1e-10: the same for any ball of fewer than 6·10⁷ voxels.

On the QSM-CI in-vivo travelling-head scan (Siemens Prisma 2.89 T, 1 mm, 184×256×128, ROMEO total field and brain mask), QSM.jl’s ismv in the QSM-CI julia-qsm image (sha256:b3c048ac…, QSM.jl fde4b92) and this port agree to float rounding (both outputs saved as float32), with identical masks: at QSM.jl’s defaults NRMSE 1.4e-9 %, max |Δ| 2.3e-10 ppm; at radius 5 mm, tol 1e-6, 50 iterations NRMSE 3.1e-9 %, max |Δ| 4.7e-10 ppm.

Other implementations, for the record

implementation	radius	stopping rule	cap
Wen et al. 2014 + lab iSMV.m (the default, IsmvVariant::Wen2014)	6 voxels (paper; also 1)	relative change < 5e-5 (paper); iSMV.m: < 0.001	none
QSM.jl ismv (fde4b92; IsmvVariant::QsmJl)	2·max(vsz) mm	residual < 1e-3 of the first	500
QSM.m ismv.m (kamesy/QSM.m f22dc40, src/bgremove/ismv.m:12-17,37-39)	3 mm	residual < √eps of the first	500

SEPIA does not offer iSMV.

Pipeline rules

Choosing IsmvVariant::QsmJl keeps the tolerance and iteration cap that are set, which by default are the paper’s (5e-5, no cap). QSM.jl’s own are published as Reference rules in crate::pipeline::rules::RULES, so a host can offer them with the variant: ismv-qsmjl.tolerance (1e-3) and ismv-qsmjl.max-iterations (500), QSM.jl fde4b92 src/bgremove/ismv.jl:39-40 (tol = 1e-3, maxit = 500; IsmvParams::qsmjl). Its radius needs no rule: radius None is already the variant’s own (2 · max(voxel size)).

Tests

tests/ismv_wen2014_parity.rs checks the default against the lab’s iSMV.m/tdcsmv.m run in GNU Octave (tests/data/ismv_wen2014_small.m): at the paper’s defaults on 1 mm and anisotropic voxels, iSMV.m as shipped (tol 0.001), and an odd-sized array padded as above. tests/ismv_qsmjl_parity.rs checks IsmvVariant::QsmJl against QSM.jl’s output.

Reference: Wen, Y., Zhou, D., Liu, T., Spincemaille, P., Wang, Y. (2014). “An iterative spherical mean value method for background field removal in MRI.” Magnetic Resonance in Medicine, 72(4):1065-1071. https://doi.org/10.1002/mrm.24998

Reference implementation: https://github.com/agr78/mSMV (iSMV.m, tdcsmv.m, 42a1f18); QSM.jl’s: https://github.com/kamesy/QSM.jl

Structs

IsmvParams

iSMV parameters. The defaults are Wen et al. 2014’s with the lab’s iSMV.m (see the module docs); IsmvParams::qsmjl gives QSM.jl’s.

Enums

IsmvVariant

Which implementation of iSMV to reproduce. See the module docs.

Constants

LAB_ERODE_THRESHOLD 🔒

iSMV.m’s erosion threshold: M = ifftn(fftn(Mask).*a) > .999.

TDCSMV_MAX_SUBVOXEL

Largest sub-voxel factor tdcsmv.m’s search for an integral voxel_size(1)·w tries before falling back to 10 (MATLAB’s loop would not end). See the module docs.

Functions

expected_iterations 🔒

Iterations an uncapped geometric convergence from d_prev to d needs to reach tol, for the progress callback (at least one more than iter).

ismv

iSMV background field removal (see the module docs)

ismv_qsmjl 🔒

QSM.jl’s ismv with an explicit radius (mm); see the module docs.

ismv_wen2014 🔒

The lab’s iSMV.m with Wen et al. 2014’s stopping rule (steps 0–4 of the module docs).

tdcsmv

tdcsmv.m (agr78/mSMV 42a1f18): the SMV sphere of radius mm on a dims array, centred at the FFT origin, normalised to sum 1 (image space; tdcsmv.m returns its fftn).