Module pdf 

Projection onto Dipole Fields (PDF) background field removal, a port of PDF.m from the MEDI toolbox (release 2024.11.26)

PDF fits the field inside the mask with the field of dipoles outside it, by weighted least squares,

min_x || W (f − D M_bg x) ||²₂ over the mask, W = 1 / N_std,

and returns what the fit does not explain, (f − D M_bg x) in the mask, as the local field. D is the dipole kernel, M_bg restricts the sources to outside the mask, and N_std is the noise standard deviation of the field.

The port follows PDF.m and its helpers line for line:

• the volume is cropped to the bounding box of the mask and zero padded after it to floor((extent − 1 + n_pad) / 2) · 2 voxels per axis (n_pad = 40), so the background sources live in the box plus the padding, and the result is put back into the full volume;
• W = 1/N_std, set to 0 where N_std is 0, and to 0 outside the mask;
• the normal equations M_bg D W² D M_bg x = M_bg D W² f (dipole_term.m) are solved by conjugate gradients (cgsolve.m, starting from 0, at most n_cg = 30 iterations, returning the last iterate; more iterations are a tuning option, see PdfParams::n_cg), stopping when the relative residual falls below tol · ||D W outside the mask|| / ||M_bg D W² f|| with tol = 0.1;
• the dipole kernel is dipole_kernel.m’s k-space kernel when B0 is along an axis (norm(B0_dir, 1) < 1.01). For an oblique B0 it is the kernel named by space, by default dipole_kernel.m’s 'imagespace' kernel, the FFT of the sampled dipole field (3 (r·B0)² − |r|²) / (4π |r|⁵) (0 at r = 0).

Units and the noise map

Like the other background removals, pdf takes the total field in ppm and returns the local field in ppm. PDF.m works in whatever unit the field is in (MEDI’s pipeline passes radians over the echo spacing), and the weighting is relative, but the absolute scale of N_std still matters: the CG tolerance is tol · ||D W|| / ||D W² f||, which scales with N_std / f. So n_std must be in the same unit as the field, here ppm. MEDI’s pipeline (README.m) takes N_std from Fit_ppm_complex.m (fit_ppm_complex_noise, radians over the echo spacing); divide it by 2π · γ · B0 · ΔTE · 1e-6 for ppm, as the pipeline does.

Deviation from PDF.m

With n_pad = 0, an axis with an odd number of voxels is not padded to an even size, and dipole_kernel.m then builds its grid as -N/2:N/2-1, half a voxel off the sample points. pdf uses the kernel centred on the voxel grid instead. With n_pad > 0 (the default) every axis is even and the kernels are dipole_kernel.m’s.

Parity

Against PDF.m as the verified QSM-CI pdf-cornell container runs it (MEDI 2024.11.26, identical to licensed MATLAB R2023b), at PDF.m’s defaults, local field NRMSE: KKI travelling head (oblique B0, imagespace) 2e-5 %; with the reference’s own ifreq and N_std, QSM-CI’s simulated head (B0 along z, so the k-space kernel) 0.007 % and KKI with space = kspace 0.026 %. At the 30-iteration cap the k-space path is rounding-sensitive (a 1-ulp change of the input moves PDF.m itself by up to 0.03 %), so the CG’s inner products are summed with compensation (see dot); with a plain sum the port was 0.61 % / 0.74 % off (astewartau/QSM.rs#405). Fixture: tests/pdf_cornell_parity.rs.

Reference pipeline

MEDI’s pipeline script (README.m of the 2024.11.26 release, line 124) runs PDF(iFreq, N_std, Mask, matrix_size, voxel_size, B0_dir) with PDF.m’s defaults (tol 0.1, n_CG 30, 'imagespace', n_pad 40) on the field from Fit_ppm_complex unwrapped by unwrapPhase, inside the BET mask of the root-sum-of-squares magnitude (README.m:118, BET defaults), and passes the local field to MEDI_L1. The mSMV authors (same group) run PDF(..., tol = 0, n_CG = 100) after ROMEO unwrapping (github.com/agr78/mSMV, README.md and data/healthy_subjects/proc.m:15). Nothing in PDF.m depends on the field strength.

Reference: Liu, T., Khalidov, I., de Rochefort, L., Spincemaille, P., Liu, J., Tsiouris, A.J., Wang, Y. (2011). “A novel background field removal method for MRI using projection onto dipole fields.” NMR in Biomedicine, 24(9):1129-1136. https://doi.org/10.1002/nbm.1670

de Rochefort, L., Liu, T., Kressler, B., Liu, J., Spincemaille, P., Lebon, V., Wu, J., Wang, Y. (2010). “Quantitative susceptibility map reconstruction from MR phase data using bayesian regularization: validation and application to brain imaging.” Magnetic Resonance in Medicine, 63(1):194-206.

Reference implementation: MEDI toolbox release 2024.11.26 (Cornell MRI Research Lab, functions/PDF.m, dipole_term.m, cgsolve.m, dipole_kernel.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): pdf.mask-bet, pdf.field-complex-fit, and the measured suggestion pdf.cg-iterations (see PdfParams::n_cg).

Structs

PdfParams

PDF parameters (PDF.m defaults)

Enums

DipoleKernelSpace

Which dipole kernel pdf uses when B0 is oblique (PDF.m’s space argument). With B0 along an axis the k-space kernel is used whatever this says, as in PDF.m.

Functions

cgsolve 🔒

cgsolve.m (Romberg; as modified for MEDI) from x0 = 0: conjugate gradients on A x = b, stopping after maxiter iterations or when ||r|| / ||b|| <= tol, with the residual recomputed from x every 50 iterations. Returns the last iterate, as cgsolve.m does.

dot 🔒

Inner product with Neumaier-compensated summation, so the sum is accurate to about one rounding whatever the length (a plain left-to-right sum over the millions of background voxels carries ~1e-13 relative error).

image_space_dipole_kernel 🔒

dipole_kernel.m’s 'imagespace' kernel: the dipole field (3 (r·B0)² − |r|²) / (4π |r|⁵) sampled at the voxel positions r (mm), 0 at r = 0, with the origin at index 0 (fftshift of the -N/2:N/2-1 grid), Fourier transformed.

pdf

PDF background field removal (PDF.m)

zero_inside 🔒