Module tfi 

TFI (preconditioned Total Field Inversion)

Single-step QSM: inverts the total field, before background removal, over the whole field of view. Background sources outside the mask are part of the solution, made cheap to reach by a preconditioner P, and the susceptibility is χ = P·y:

min_y ||wG ∇(P y)||₁ + 2λ ||m (D P y − f)||²₂

What follows which source

Cornell’s own TFI script (TFI_L1.m) is not public. The 2024.11.26 MEDI toolbox has the TFI solver inside MEDI_L1.m: when RDF.mat holds a preconditioner P (parse_QSM_input.m), MEDI_L1 solves for y with χ = P·y. This port follows that code line for line:

• 2024 MEDI_L1.m and its helpers (fgrad.m, bdiv.m, gradient_mask.m, dataterm_mask.m, cgsolve.m, dipole_kernel.m): the 'linear' data term (and the toolbox’s default exponential one, as an option), the data weight m from N_std, the edge mask wG and then wG.*Mask (no regularisation outside the mask), the lagged-diffusivity weight 1/sqrt(|wG ∇(P y)|² + 1e-6·P²) (epsilon scaled by P²), the operator PH(reg(P dx)) + 2λ·PH(Dconv(|w|² Dconv(P dx))), the right-hand side, the Gauss-Newton stop ||dy||/||y||, MERIT, the output P y in ppm masked by Mask, and the defaults of parse_QSM_input.m (λ = 1000, MERIT on, 10 Gauss-Newton and 100 CG iterations, tolerances 0.1 and 0.01, percentage 0.9).
• Liu et al. 2017, for what the toolbox does not settle for TFI: P = 1 inside the mask and P_B = 30 outside (eq. 3; 30 is the paper’s in-vivo value); the input is the total field over the whole field of view, so there is no SMV step (an SMV high-pass would remove the background TFI is meant to fit); and the data term is the linear one of eq. 5, which is MEDI_L1’s 'linear' option rather than its exponential default. (With the exponential term, MEDI_L1 itself fails on in-vivo total fields in radians: the field wraps inside exp(i·). Set linear: false to run the toolbox default.)

Not ported: the R2*-based preconditioner (eq. 4), zero padding and the CSF zero reference.

Parity

Against TFI as the verified QSM-CI tfi-cornell container runs it (MEDI 2024.11.26 MEDI_L1.m with a preconditioner, identical to licensed MATLAB R2023b), astewartau/QSM.rs#407. Two departures are deliberate, each with the original’s behaviour as an option:

• Dipole kernel on odd grids (TfiKernel). dipole_kernel.m builds its k grid as -N/2:N/2-1, half a frequency step off the DFT grid on an odd axis, so it is not the dipole kernel there. On QSM-CI’s simulated head (164×205×205) that alone made the port’s χ 0.77× the original’s (42 % NRMSE); with TfiKernel::MediOddGrid 2.5 %. The DFT-grid kernel fits the true χ better (demeaned, in the mask: NRMSE 62.1 % against 71.4 %, r 0.785 against 0.718), so it stays the default. On even grids the two are the same.
• Precision (TfiPrecision). MEDI_L1.m runs in single precision (tempn = single(N_std) makes the data weight single, and MATLAB’s mixed arithmetic carries that into everything). On the fixture’s small odd phantom (tests/tfi_cornell_parity.rs) the port with both options is 2.9e-6 relative after one Gauss-Newton step and 4.7e-4 with all defaults; in double precision 4.7e-3. The container itself differs by 3.0e-4 there on another CPU.

Beyond that the original is not reproducible at full size: its 100-iteration CG solves at cg_tol 0.01 in single precision amplify rounding. In the same container, 1e-7 relative noise on the simulated head’s total field moves its χ by 1.6 % NRMSE, and the container on another CPU (QSM-CI’s runner against this VM) differs from itself by 2.8 % on the KKI travelling head. On the reference’s own inputs the port with both options is 3.6 % (sim) and 1.26 % (KKI) from it; with the DFT-grid kernel in double precision, 42 % and 1.25 %.

Reference pipeline

From Liu et al. 2017 (Methods, in-vivo experiments); not stated above:

• Brain mask (the M where P = 1): BET. Field map: nonlinear (complex) field estimation followed by graph-cut phase unwrapping. No background removal or high-pass before TFI.
• In vivo the paper uses the R2*-based preconditioner of eq. 4, P = 1 in M ∩ {R2* < 30 s⁻¹} and P_B elsewhere, with P_B = 30 (10 for the phantom). Only eq. 3 is ported here.
• Data: 3 T, 6 echoes, TE₁ 3.5 ms, ΔTE 3.6 ms, 1×1×3 mm³, monopolar readout.
• The paper’s solver settings differ from MEDI_L1’s defaults, which this port uses: λ = 1e-3 on the L1 term (not the data term), Gauss-Newton stop ||dy||/||y|| < 0.01, up to 300 CG iterations. CSF is the reference for reported values.

Units

Like the other inversions, tfi takes the field in ppm and returns ppm. MEDI_L1 works on the field in radians accumulated over the echo spacing (delta_TE) at CF = γ·B0; the field is converted with 2π·CF·delta_TE·1e-6 first and back at the end, as MEDI_L1.m does. The result depends on this scale, so b0 and echo_times are scan parameters that the pipeline fills from metadata, and TfiParams::te picks the echo time the conversion uses (by default the reference’s choice; see there).

Magnitude and data weight

As for medi, the echo magnitudes (EchoMagnitudes) are required: the edge mask uses their root sum of squares (iMag = sqrt(sum(abs(iField).^2,4))) and the data weight m = Mask./N_std (mean 1 in the mask) uses Fit_ppm_complex.m’s noise map (fit_ppm_complex_noise), as MEDI’s pipeline makes them, so TFI always weights the data as the reference does.

Reference: Liu, Z., Kee, Y., Zhou, D., Wang, Y., Spincemaille, P. (2017). “Preconditioned total field inversion (TFI) method for quantitative susceptibility mapping.” Magnetic Resonance in Medicine, 78(1):303-315. https://doi.org/10.1002/mrm.26331

Reference implementation: MEDI toolbox 2024.11.26 (Weill Cornell), functions/MEDI_L1.m.

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): tfi.mask-bet, tfi.te-spacing.

Structs

TfiParams

TFI parameters (2024 MEDI_L1.m / parse_QSM_input.m defaults; P_B and the linear data term from Liu et al. 2017)

Enums

TfiKernel

Which dipole kernel tfi uses: they differ only on grids with an odd number of voxels along some axis.

TfiPrecision

Arithmetic of tfi’s Gauss-Newton iterations.

Functions

tfi

Preconditioned total field inversion (2024 MEDI_L1.m with a preconditioner)

tfi_l1 🔒

MEDI_L1.m with the preconditioner, on an iMag and an N_std

tfi_single 🔒

tfi_l1 in MATLAB’s single precision (TfiPrecision::Single): MEDI_L1.m line for line with MATLAB’s types. tempn = single(N_std) makes the data weight m, and with it the data b0 = m.*RDF, single; from there every product with a single operand is single, so the operator, the right-hand side, cgsolve.m’s vectors and scalars, dx and the iterate x (x + dx) are single. Double operands (P, the edge mask, ε, the kernel, and after MERIT the data weight and data, which dataterm_mask remakes from double(N_std)) are rounded to single where they meet a single one. FFTs are single precision.