Module medi 

MEDI (Morphology Enabled Dipole Inversion), a port of MEDI_L1.m from the MEDI toolbox (release 2024.11.26)

Solves

min_χ ||wG ∇χ||₁ + 2λ ||m (exp(i Dχ) − exp(i f))||²₂

by Gauss-Newton, with the L1 term handled by lagged diffusivity and each Gauss-Newton step solved by conjugate gradients. f is the local field in radians, m the data weight (dataterm_mask.m) and wG the binary edge mask from the magnitude (gradient_mask.m): 0 at the strongest magnitude edges, 1 elsewhere, per gradient direction.

The port follows MEDI_L1.m and its helpers line for line: fgrad.m/bdiv.m (forward differences, zero at the last voxel, and their adjoint), the per-direction weight 1/sqrt(|wG ∇χ|² + 1e-6), cgsolve.m, sphere_kernel.m (the partial-volume sphere) for the SMV step, and the MERIT reweighting. Defaults are those of MEDI_L1.m and parse_QSM_input.m: λ = 1000, 10 Gauss-Newton iterations, 100 CG iterations, tolerances 0.1 and 0.01, percentage 0.9, SNR data weighting, SMV (radius 5 mm, shrinking the mask), MERIT and the CSF zero reference (MEDI+0, lambda_CSF 100) on. Not ported: the preconditioner P, the linear formulation and the Lreg term (all off unless requested).

MEDI+0: CSF zero reference

MEDI’s pipeline (README.m) computes R2* from the echo magnitudes with arlo.m (arlo), the ventricular CSF mask from it with extract_CSF.m (extract_csf), and passes the mask to MEDI_L1, which then adds reg_CSF.m’s term lambda_CSF/2 · ||M_CSF (χ − mean_CSF(χ))||² to the cost and subtracts the CSF mean from the result. medi does the same from the echo magnitudes and echo times it requires (MediParams::csf, default on, as parse_QSM_input sets lambda_CSF = 100). Differences from MATLAB, which only matter where MATLAB fails or gives no answer:

• Fewer than three echoes: arlo.m gives no R2* for one echo (so MEDI runs without the CSF term) and 0 everywhere for two (which would make the whole central brain “CSF”). medi runs without the CSF term in both cases.
• A CSF mask with no voxels makes MEDI_L1 return NaN (the mean over no voxels); medi runs without the CSF term instead.
• extract_CSF.m errors when the central sphere holds fewer than three low-R2* regions; extract_csf uses those it finds.
• MEDI_L1 subtracts the CSF mean from every voxel, so the result is −mean outside the mask; medi keeps the result zero outside the mask, like the other inversions. Inside the mask the values are the same.

Units

Like the other inversions, medi takes the local field in ppm and returns ppm. MEDI works on the field in radians accumulated over the echo spacing (delta_TE, TE(2) − TE(1) for multi-echo data and TE for one echo, as MEDI’s DICOM readers set it) at the centre frequency CF = γ·B0. The field is converted with 2π·CF·delta_TE·1e-6 first and the result converted back at the end, as MEDI_L1.m does. The result depends on this scale (the data term is nonlinear and λ is fixed), so b0 and echo_times are scan parameters that the pipeline fills from metadata, and MediParams::te picks the echo time (by default delta_TE).

Magnitude and data weight

MEDI_L1.m weights the data by 1/N_std, the noise standard deviation of the field (normalised to mean 1 in the mask), and builds the edge mask from iMag. MEDI’s pipeline (README.m) takes both from the complex echoes: iMag = sqrt(sum(abs(iField).^2,4)) and N_std from Fit_ppm_complex.m, which depends only on the echo magnitudes (fit_ppm_complex_noise). medi requires the echo magnitudes (EchoMagnitudes) and their echo times (MediParams::echo_times) and does the same, so it always weights the data as MEDI_L1 does in MEDI’s pipeline; one echo is fitted as Fit_ppm_complex.m fits one echo. Fit_ppm_complex.m assumes evenly spaced echoes.

Zero padding

zeropad pads the volume with zeros at the end of each axis before the inversion and crops the result back. That is what MEDI_L1’s 'zeropad' option is meant to do, but in MEDI 2024.11.26 the option fails: parse_QSM_input pads the volumes but leaves matrix_size and the masks (Mask, Mask_CSF) at the unpadded size, so MEDI_L1 stops on mismatched arrays. The QSM-CI medi-cornell container pads the inputs itself before calling MEDI_L1 instead (QSMxT/QSM-CI#309).

Input: the local field

medi takes the local field from a background removal, with or without SMV or mSMV. The MEDI toolbox’s pipeline (README.m) runs RDF = PDF(iFreq, N_std, Mask, ...) and then MEDI_L1('lambda',1000,'lambda_CSF',100,'merit','smv',5) on RDF. Its naming convention calls RDF the “Relative Difference Field, aka local field”, and MEDI_L1 reads iFreq (the total field) without using it. The mSMV paper (Roberts et al. 2024, Eqs. 4 and 13) defines MEDI-SMV on the local field after PDF, b_L0 = b̂_L − S·M·b̂_L on the eroded mask, and SEPIA’s MEDI wrapper passes MEDI its local field. MEDI’s SMV step only removes the residual background that a background removal leaves near the edge. The pipeline therefore routes MEDI through a background removal (InversionConfig::input), PDF in the authors’ pipeline (reference_bg_removal).

Until #287, QSM.rs gave MEDI-SMV the total field and the brain mask instead. The two routes were compared with MEDI’s defaults (λ 1000, MERIT, MEDI+0) and QSM.rs’s PDF. Data: the qsm-forward 7 T head phantom (TE 4/12/20/28 ms), noise-free from its true total field and at SNR 50 through QSM.rs’s field mapping, scored over the brain mask; and four 3 T EMTP travelling-head scans in QSM-CI (Philips Achieva, Siemens Skyra, Siemens Prisma, Philips Ingenia), scored by the χ variance in the SynthSeg cortex (the mSMV paper’s shadow score):

input → MEDI	phantom, noise-free: XSIM / NRMSE / corr	phantom, SNR 50	in vivo cortex χ variance (×1e-4 ppm²), 4 scans
total field → MEDI-SMV (before #287)	0.563 / 62.6 % / 0.780	0.352 / 76.2 % / 0.657	4.76 / 8.99 / 4.12 / 5.71
PDF → MEDI-SMV (authors)	0.563 / 62.6 % / 0.780	0.352 / 76.2 % / 0.658	4.86 / 9.09 / 4.14 / 5.76
total field → MEDI-SMV, no shrinking	0.233 / 140.6 % / 0.458	0.079 / 315.0 % / 0.207	850 / 383 / 875 / 983
PDF → MEDI-SMV, no shrinking	0.558 / 65.1 % / 0.760	0.292 / 94.0 % / 0.511	6.9 / 12.7 / 9.1 / 11.4

With the default mask shrinking the two routes give nearly the same χ: the SMV high-pass removes any harmonic field inside the eroded mask, and the total field differs from PDF’s local field by a field that is harmonic there. The maps differ by 0.5 ppb RMS on the phantom (relative L2 1.3 %) and 1.5–4.9 ppb in vivo (4–14 %). Without shrinking, the edge band keeps the background field on the total-field route, and the result is unusable. The authors’ route costs a PDF run, 10–80 s on these volumes.

SMV options

What MEDI_L1.m does with each combination of 'smv'/'nosmv', smv_shrink_mask and 'msmv' (SphereK is sphere_kernel.m’s kernel, SMV(x) = ifftn(SphereK·fftn(x))):

options	mask	field RDF	kernel D	noise tempn
'nosmv'	unchanged	unchanged	D	N_std
'smv', shrink on (default)	SMV(Mask) > 0.999	Mask·(RDF − SMV(RDF))	(1 − SphereK)·D	sqrt(SMV(tempn²) + tempn²)
'smv', shrink off, release 2024.11.26	unchanged	Mask·(RDF − SMV(RDF))	(1 − SphereK)·D	sqrt(SMV(tempn²) + tempn²)
'smv', shrink off, mSMV prerelease	unchanged	unchanged	(1 − SphereK)·D	sqrt(SMV(tempn²) + tempn²)
'msmv' (prerelease only; forces shrink off)	unchanged	msmv(RDF, Mask, R2s, ...)	(1 − SphereK)·D	sqrt(SMV(tempn²) + tempn²)

medi follows the 2024.11.26 release for the first three rows (MediParams::smv, MediParams::smv_shrink_mask) and the prerelease for 'msmv' (MediParams::msmv), which only the prerelease has. The prerelease (MEDI_toolbox-mSMV100125.zip) moved the field high-pass inside if opts.smv_shrink_mask && ~opts.msmv, so with SMV on and shrinking off it fits the unfiltered field with the SMV-filtered kernel. Its comments and the mSMV paper do not describe this case: the paper’s MEDI-SMV (Eq. 13) is the shrinking one. The two releases agree in every other row and, apart from the mSMV code, in the rest of MEDI_L1.m. In both, 'msmv'/'smv' set smv = 1. In the prerelease, 'nosmv' given after 'msmv' turns SMV off and with it mSMV, a combination medi does not have: MediParams::msmv implies SMV.

Checked against both toolboxes in MATLAB R2023b (tests/medi_smv_parity.rs, a 32×36×28 phantom, no CSF term). The two MEDI_L1.m versions give identical χ with 'nosmv' and with 'smv' and shrinking. With shrinking off they differ by 101 % relative L2 after one step and 176 % with all defaults. medi against the release, χ over the mask:

options	one Gauss-Newton step: max abs diff (ppm), rel. L2	all defaults (10 steps)
'nosmv'	2.5e-7, 1.1e-6	1.4e-3, 1.9e-3
'smv', shrink on	1.1e-8, 2.1e-6	2.9e-8, 5.8e-7
'smv', shrink off	1.7e-7, 5.7e-7	2.0e-4, 2.1e-4

After one step the differences are single (MATLAB) against double (here) rounding. Over ten steps they also include CG stopping one iteration apart at cg_tol 0.01.

tests/medi_cornell_parity.rs checks the same against the verified QSM-CI medi-cornell container (MEDI 2024.11.26, QSMxT/QSM-CI#309; regenerated in Actions by tests/data/medi_cornell_small.py), measured 2026-10-30: one step 2.7e-7 / 1.3e-8 / 1.5e-7 ppm max abs diff for the three rows above (the container writes χ as float32); and MEDI_L1’s defaults end to end (SMV with shrinking, MERIT, 10 steps, MEDI+0 with the CSF mask from arlo and extract_csf on the phantom with three CSF-like blobs) 2.9e-8 ppm, relative L2 6.0e-7, every CG solve stopping where MATLAB’s did. On the way, fit_ppm_complex_noise (N_std), arlo (R2*) and extract_csf are bit-identical to the container’s.

MEDI-mSMV

MediParams::msmv is MEDI_L1(..., 'msmv', radius) from the mSMV authors’ MEDI toolbox (agr78/mSMV, prerelease/MEDI_toolbox-mSMV100125.zip; Roberts et al. 2024, Eq. 12): SMV is on with radius MediParams::smv_radius, but instead of eroding the mask and high-passing the field, MEDI_L1 passes the local field and the full mask to msmv.m (msmv), then inverts the result with the SMV-filtered dipole kernel (1 − SphereK)·D and the SMV-propagated noise, as MEDI-SMV does. The mask stays whole, so the result reaches the brain edge. msmv.m needs R2* for its vessel mask, and MEDI_L1 errors without it; medi computes it from the echo magnitudes with arlo.m (arlo_unclipped), as MEDI’s pipeline and the authors’ SEPIA patch do, so MEDI-mSMV needs at least three echoes (and panics with fewer; the pipeline returns an error). msmv.m‘s field strength B0_mag is MediParams::b0. MEDI-mSMV takes the local field from a background removal (the paper: PDF, V-SHARP and LBV), not the total field. See msmv for what the filter does and how it compares with the authors’ MATLAB.

Deliberate deviation from MEDI_L1.m: with SMV on, sphere_kernel.m builds its grid as -N/2:N/2-1, which along an axis with an odd number of voxels centres the sphere half a voxel off the grid. The SMV operator is then not symmetric, and MEDI_L1 diverges (in MATLAB as here: NRMSE above 10000% on a 164×205×205 phantom). So when SMV is on, medi pads every axis that is still odd after zeropad by one more voxel and crops the result back. On grids that are even after zeropad, and with SMV off, nothing changes and the result is MEDI_L1’s.

Reference pipeline

What MEDI’s own pipeline script (README.m of the 2024.11.26 release; the release has no MEDI_full_process.m) runs around MEDI_L1. Inputs, units and the CSF term are above.

• Field map: Fit_ppm_complex(iField) on the complex echoes (evenly spaced; README.m:87), which also gives N_std; Fit_ppm_complex_bipolar for bipolar readouts (README.m:98). Unwrapping: unwrapPhase(iMag, iFreq_raw, matrix_size), MEDI’s region-growing unwrapper (README.m:93; graph-cut and Laplacian unwrapping are commented alternatives).
• Brain mask: FSL BET via the bundled bet2 MEX, BET(iMag, matrix_size, voxel_size) on the root-sum-of-squares magnitude (README.m:118), with BET’s defaults (fractional threshold 0.5, gradient threshold 0; functions/bet2/bet2.cpp:136). genMask.m (magnitude threshold) is a commented alternative.
• Background removal: PDF, PDF(iFreq, N_std, Mask, matrix_size, voxel_size, B0_dir) with PDF.m’s defaults (README.m:124). LBV with tol 0.005 is the commented alternative (README.m:129). MEDI_L1 then inverts this local field (RDF), so with SMV on, the SMV high-pass is applied on top of the PDF local field, not to the total field.
• Inversion: MEDI_L1('lambda',1000,'lambda_CSF',100,'merit','smv',5) (README.m:144), which equals the 2024 defaults (parse_QSM_input.m:10-20: MERIT and SMV on, radius 5 mm). The 2020 release’s parse_QSM_input.m defaulted to MERIT and SMV off.
• ppm conversion: delta_TE is TE(2) − TE(1), or TE for one echo (Read_DICOM.m:166-168), and CF is read from the DICOM header; QSM.rs uses CF = γ·B0.
• Field strength: nothing in MEDI_L1, PDF, extract_CSF or arlo depends on it beyond CF; λ = 1000 and extract_CSF’s 5 Hz R2* threshold are the same at 3 T and 7 T.

Reference: Liu, T., Liu, J., de Rochefort, L., Spincemaille, P., Khalidov, I., Ledoux, J.R., Wang, Y. (2011). “Morphology enabled dipole inversion (MEDI) from a single-angle acquisition: comparison with COSMOS in human brain imaging.” Magnetic Resonance in Medicine, 66(3):777-783. https://doi.org/10.1002/mrm.22816

Liu, J., Liu, T., de Rochefort, L., Ledoux, J., Khalidov, I., Chen, W., Tsiouris, A.J., Wisnieff, C., Spincemaille, P., Prince, M.R., Wang, Y. (2012). “Morphology enabled dipole inversion for quantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map.” NeuroImage, 59(3):2560-2568.

Reference implementation: MEDI toolbox release 2024.11.26 (Cornell MRI Research Lab, functions/MEDI_L1.m)

The double-precision helpers below are shared with tfi, which is MEDI_L1.m run with a preconditioner, and with the χ-separation ports.

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): medi.mask-bet, medi.field-complex-fit, medi.te-spacing, medi-msmv.three-echoes, medi-msmv.bg-removal, medi-msmv.prefilter, medi.bg-pdf.

Structs

CgBuffers 🔒

cgsolve.m’s vectors, reused across Gauss-Newton iterations

EchoMagnitudes

The magnitude of each echo, n_echoes volumes stacked with the echo slowest. MEDI’s edge mask uses their root sum of squares and its data weight the Fit_ppm_complex.m noise.

MediParams

MEDI parameters (MEDI_L1.m / parse_QSM_input.m defaults)

Functions

bounding_box 🔒

The smallest box holding every true voxel (the whole grid if there is none)

cgsolve 🔒

cgsolve.m from a zero start; the last iterate is left in cg.x. progress receives the iteration count. The dot products are summed in index order, as before.

dataterm_mask 🔒

dataterm_mask.m into w: 0 gives 1 everywhere; 1 gives Mask./N_std (0 where that is not finite) normalised to mean 1 over the mask.

dconv 🔒

Dconv: out = real(ifftn(D.*fftn(x)))

dconv_weighted 🔒

Dconv(w2.*Dconv(p)), left in the real part of buf. w2 is zero outside mbox, so the inner Dconv is computed only on mbox and the outer one reads input only there; the values are those of the two plain Dconvs (see [Fft3dWorkspace::fft3d_with]).

extract_csf

extract_CSF.m (MEDI toolbox 2024.11.26): the ventricular CSF mask for MEDI+0.

fgrad 🔒

fgrad.m: forward differences, zero at the last voxel along each axis

gradient_mask 🔒

gradient_mask.m (binary mode): 1 where |∇(iMag.*Mask)| is at or below a threshold, per direction, over the whole volume. The threshold starts at 1% of max(iMag) and is scaled by 1.05 or 0.95 until the number of gradient components above it, over the number of voxels in the mask, crosses percentage.

label6 🔒

6-connected components (bwconncomp(·, 6)), numbered in the order of their first voxel in memory order (x fastest, as MATLAB’s column-major order). Returns each voxel’s component and the component sizes.

load_real 🔒

Loads a real array into the transform’s rows as complex(x, 0)

medi

MEDI dipole inversion (MEDI_L1.m)

medi_l1 🔒

MEDI_L1.m on an iMag, an N_std and, for MEDI+0, a Mask_CSF

mul_kernel 🔒

Multiplies the transform’s rows by k as they are loaded (K.*buf, fused into the first pass of the next transform)

par_fill 🔒

out[i] = f(i) for every i (over threads with the parallel feature)

par_update 🔒

f(i, &mut out[i]) for every i (over threads with the parallel feature)

reg 🔒

reg(dx) = div(wG.*(Vr.*(wG.*grad(real(dx))))) into out, with wvw = wG.*Vr.*wG: bdiv of the weighted fgrad, a z plane at a time so the gradients never leave the cache. Each value is computed with the same operations as fgrad followed by bdiv.

smv 🔒

SMV.m with a precomputed kernel: ifftn(fftn(x).*K), left in buf

smv_noise 🔒

The SMV noise propagation in MEDI_L1.m, in place: tempn = sqrt(SMV(tempn.^2, SphereK) + tempn.^2)

sphere_kernel 🔒

sphere_kernel.m as MEDI_L1.m and extract_CSF.m use it here: fftn(fftshift(Sphere)) in double, complex (the MEDI toolbox’s sphere_kernel.m before the mSMV release added real(...); the two differ only on odd grids, which medi pads away)

sphere_kernel_real 🔒

sphere_kernel(matrix_size, voxel_size, radius) with the mSMV release’s sphere_kernel.m (real(fftn(fftshift(Sphere)))): the real part of the FFT of sphere_kernel_sparse, in double, as SMV.m builds it from a radius (SMV(x, matrix_size, voxel_size, r)).

sphere_kernel_single 🔒

sphere_kernel_real rounded to single precision: single(sphere_kernel(...)), as msmv.m (and MEDI_L1.m) store the prefilter kernel SphereK.

sphere_kernel_sparse 🔒

sphere_kernel.m’s sphere in image space: a sphere of radius mm with partial-volume edge voxels (each split into 20³ sub-voxels), on MATLAB’s -N/2:N/2-1 grid, fftshifted (circshift by floor(N/2)) and normalised to sum 1. Returned sparsely as (flat index, value) for the nonzero voxels, in column-major order. Only voxels within one voxel of the sphere are visited, so small radii (mSMV’s kernel_lim.m and its minimum-radius kernel) are cheap; the values are those of the full-grid loop.