Module vsharp 

V-SHARP background field removal

Variable kernel SHARP uses multiple SMV kernel radii to preserve more brain tissue at edges while still removing background fields.

Each voxel takes the SMV high-pass φ − S_r∗φ of the largest radius r whose ball fits inside the mask. By default that high-passed field is returned as is, which is what the authors’ code, STI Suite’s V_SHARP, does. With VsharpParams::deconvolve on it is then deconvolved by the largest kernel, 1/(1 − Ŝ) with TSVD at threshold (as QSM.jl does). The deconvolution restores low spatial frequencies; on real data those are dominated by residual background and field-map error, so leaving it off gives a more reproducible local field at the cost of attenuating genuinely low-frequency tissue contrast.

Defaults

The defaults are STI Suite 3.0’s V_SHARP (Liu lab, the authors of the V-SHARP papers), which is P-code and was probed black-box (ukb-qsmxt docs/08_vsharp_sti.md): crop to the mask’s bounding box + 12 voxels, linear resampling to a 1 mm grid, balls |x| ≤ r + 1/6 (centre excluded) for r = 1…smvsize mm, the largest fitting ball’s background per voxel, no deconvolution, the background interpolated back; smvsize defaults to 12 mm. Here that is deconvolve: false, an Isotropic 1 mm grid (same 12-voxel crop margin, CROP_MARGIN), max_radius 12 and STI’s radius schedule, 12 → 2 mm in 1 mm steps (VsharpParams::radius_step), and STI’s ball, VsharpBall::StiSuite (|x| ≤ r + 1/6 without its centre voxel, a ball fitting only where all its voxels are in the mask), and STI’s crop: its MaskBoundingBox makes the bounding box even along every axis before the margin (sti_even_box).

Measured (QSM-CI parity, astewartau/QSM.rs#329): on the KKI travelling head (256×256×128, 1 mm, STI Suite V_SHARP in the verified vsharp-sti v2 container, the same total field and mask), the STI ball gives the same eroded mask and NRMSE 2e-6 relative (STI’s float32 output), where the earlier ball (radius r, centre included, VsharpBall::QsmJl) gave 4.4 %. tests/vsharp_sti_parity.rs holds a small case (40×44×36, smvsize 12 and 6, from the same container, regenerated in Actions by tests/data/vsharp_sti_small.py): the same eroded mask and max |Δ| < 5e-8 on a field of up to 0.15, STI’s single-precision filtering. The papers and the code disagree on the deconvolution: Li et al. 2011 (NeuroImage 55:1645, Image Analysis and Results) and Wu et al. 2012 (Background phase removal: “a deconvolution is carried out to restore the local phase”, with the kernel of the largest diameter used) both deconvolve the SMV-filtered phase, without giving the kernel’s truncation; STI Suite’s V_SHARP does not (its impulse response is exactly δ − S₁₂). The code is followed. On 150 UK Biobank subjects × 2 visits (ukb-qsmxt docs/08_vsharp_sti.md, PR #160) deconvolving lowered deep-grey-matter test-retest ICC (0.922 → 0.915) and removed the detectable longitudinal change (SRM putamen + caudate 0.161 → −0.019); on the noise-free CI phantom it is closer to the true local field, since there the low frequencies are tissue. Measured with the pipeline’s ROMEO field and RTS against the previous defaults (deconvolution, native grid, radii 12 → 1 mm), on the voxels both keep: qsm-forward phantom χ XSIM 0.261 → 0.256, NRMSE 84.4 → 85.7 %, r 0.546 → 0.521; in-vivo 3T Philips scan, agreement with the UK Biobank pipeline’s χ (STI Suite V-SHARP + iLSQR) XSIM 0.641 → 0.714, NRMSE 75.2 → 69.4 %, r 0.781 → 0.790. On these 1 mm isotropic inputs the 1 mm grid only adds the crop.

Anisotropic voxels (astewartau/QSM.rs#329, measured 2026-10-30): on UK Biobank’s 1×1×1.5 mm (a 32×32×20 phantom through UKBiobank_QSM_core.m in the verified ukb-qsm v2 container, identical to licensed MATLAB R2023b on it) the eroded mask was STI’s but the local field 4.1 % NRMSE from it. The cause was the crop, not the resampling: STI’s MaskBoundingBox turns the mask’s 11 slices into 12, and that slice’s phase (outside the mask) reaches the 1 mm working grid through the z interpolation. Spied in licensed MATLAB, STI’s resampling is exactly this module’s (interp3 linear, native samples at (1:N)·vs mm, working samples at 1:ceil(N·vs) mm, the mask > 0.5 both ways, the background interpolated back in float32), and its working-grid background is the phase outside the eroded mask and the largest fitting ball’s S_r∗φ inside it. With the even box that phantom agrees to 1.2e-6. A second, smaller cause showed on non-integer resamplings (1.2×0.8×1.5 mm): STI interpolates the mask in single precision, and exact ties at 0.5 fall outside its > 0.5 test where f64 put some just inside (trilinear_mask_f32; 0.96 % and 13 output voxels before). With both, tests/vsharp_sti_aniso_parity.rs (four phantoms from 1×1×1.5 to 0.7×0.7×3 mm, the verified vsharp-sti container) agrees to 7e-7–1.1e-6 NRMSE (STI’s float32) with the same masks; before, 6e-4 to 2 % and up to 22 mask voxels apart.

QSM.jl’s vsharp

The filtering (the per-voxel choice of radius and the deconvolution, and with VsharpBall::QsmJl the balls and the edge test) is a port of QSM.jl’s vsharp (kamesy/QSM.jl fde4b92, _sharp! for a radius vector in src/bgremove/sharp.jl), including its working box: the mask’s bounding box zero-padded to QSM.jl’s fast FFT size for the largest radius, so the deconvolution is periodic on that box (as in [super::sharp]), and its edge test S_r ∗ mask > 1 − √eps(Float64). QSM.jl’s own defaults are reachable with

use qsm_core::bgremove::{VsharpBall, VsharpGrid, VsharpParams};
let vsz = (1.0_f64, 1.0_f64, 1.0_f64);
let (lo, hi) = (vsz.0.min(vsz.1).min(vsz.2), vsz.0.max(vsz.1).max(vsz.2));
let qsmjl = VsharpParams {
    max_radius: 18.0 * lo,   // QSM.jl: r = 18*minimum(vsz):-2*maximum(vsz):2*maximum(vsz)
    min_radius: 2.0 * hi,
    radius_step: 2.0 * hi,
    threshold: 0.05,         // QSM.jl's thr (also this module's default)
    deconvolve: true,        // QSM.jl always deconvolves
    ball: VsharpBall::QsmJl,  // QSM.jl's ball and edge test
    grid: VsharpGrid::Native,
};

The deconvolution threshold default, 0.05, is QSM.jl’s vsharp thr and QSM.m’s vsharp.m thr (kamesy/QSM.m f22dc40, src/bgremove/vsharp.m:15,32): STI Suite’s V_SHARP (the original authors’ code, whose other defaults this module keeps) does not deconvolve, so the threshold has no STI value. QSM.rs defaulted to 0.2 from PR #72 (2026-09-07), a QSM.rs tuning: on the QSM Challenge 2.0 phantom (ROMEO field, V-SHARP, several dipole inversions) 0.2 gave a lower NRMSE on every inversion tested, and was identical to 0.05 on a noise-free field, because the higher threshold leaves out more of the low-frequency band where field-mapping error is amplified. It remains available as an explicit threshold.

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) these settings and QSM.jl’s vsharp in the QSM-CI julia-qsm image (sha256:b3c048ac…) agree to float rounding: NRMSE 2.6e-10 %, max |Δ| 5.8e-11 ppm, identical masks (both outputs saved as float32); with radii 12 → 1 mm in 1 mm steps, NRMSE 2.1e-10 %. Before the working box and edge test were QSM.jl’s (the deconvolution was periodic on the whole field of view, and the edge test was > 1 − √1e-7), the 12 → 1 mm comparison gave NRMSE 15.5 %. tests/sharp_qsmjl_parity.rs checks the port against QSM.jl’s output on small phantoms (max |Δ| < 1e-15 ppm).

The working box depends on the voxel size at the float rounding level: its padding is ⌈r_max / min(vsz)⌉, and NIfTI headers often store 1 mm as 0.99999994 (one float32 ulp below), for which a 12 mm radius pads 13 voxels, not 12. On the QSM-CI KKI travelling head (Philips 3 T, 256×256×128, header z spacing 0.99999994, params.json 1.0) that alone gave NRMSE 1.19 % between QSM.rs (header voxel size) and QSM.jl (params.json’s) with radii 12 → 1 mm (QSMxT/QSM-CI parity run 41029942914); given the same voxel size the outputs are identical in float32. The port keeps QSM.jl’s rule (no rounding of the voxel size), so the two agree whenever they are given the same voxel size; tests/sharp_qsmjl_parity.rs has a case at that spacing (0.75 % from QSM.jl if 1.0 is used instead).

The old edge test 1 − √1e-7 also let a ball count as fitting when it missed a few of its voxels (up to 2 of the 7,153 voxels of a 12 mm ball at 1 mm), so near the rim some voxels took a larger radius than fits. The QSM.jl test fixes that for every setting, including the STI Suite defaults: on the travelling-head scan the default output moved by NRMSE 1.8 % (same output mask), all of it from the edge test (the working box alone changes it by 2e-7 %).

Working grid

The SMV balls live on a voxel grid. On strongly anisotropic voxels (e.g. 0.8 × 0.8 × 3 mm) the native grid makes them crude in the coarse direction: a 1–2 mm ball is a single-slice disk. VsharpParams::grid can instead estimate the background on an isotropic working grid, at the finest input voxel size (VsharpGrid::IsotropicFinest) or at a fixed spacing (VsharpGrid::Isotropic):

1. crop the field and mask to the mask’s bounding box, made even along each axis as STI’s MaskBoundingBox does (sti_even_box), with a zero margin of CROP_MARGIN native voxels per side (this also fixes where the working grid sits relative to the anatomy, so it lands the same way on repeat scans);
2. resample both trilinearly onto the isotropic grid (mask: interpolated value > 0.5), axes whose voxel size already equals the spacing are left as they are;
3. run V-SHARP there (same balls, radii, edge test and deconvolution as on the native grid);
4. interpolate the background (field − local field on the working grid) back to the native grid trilinearly, and return local = native field − background, on the voxels where the back-interpolated working-grid support is > 0.5.

The local field thus keeps native resolution; only the smooth background estimate is made on the working grid. The default is VsharpGrid::Isotropic at 1 mm (STI Suite’s grid; see Defaults). VsharpGrid::IsotropicFinest makes the balls spherical without coarsening any axis, at the cost of a working grid that is larger by the anisotropy ratio. A fixed VsharpGrid::Isotropic spacing gives the same working grid whatever the acquisition resolution, but coarsens axes finer than it, including the mask edge the background is estimated from.

Reference pipeline

STI Suite 3.0’s V_SHARP(φ, mask, 'voxelsize', vs, 'smvsize', R) (Liu lab, the authors of the V-SHARP papers) is P-code; black-box probing (ukb-qsmxt docs/08_vsharp_sti.md, scripts in vsharp_probe/) found: crop to the mask’s bounding box, zero-pad 12 voxels per side, resample linearly to a 1 mm grid, ball |x| ≤ r + 1/6 without its centre for r = 1…R mm, each voxel takes the background S_r∗φ of the largest ball that fits, no deconvolution, background interpolated back and subtracted from the native field. Voxels where only r = 1 fits are dropped, so the radii in effect are 2…R mm and the output mask is the voxels a 2 mm (+1/6) ball fits around, back-interpolated > 0.5; this module’s mask is the union down to VsharpParams::min_radius, which by default is the same 2 mm (see Defaults), and the ball is STI’s (VsharpBall::StiSuite). smvsize defaults to 12 (probed: bitwise equal to an explicit 12), which is also what STI Suite’s GUI passes (its “Radius” field) and what QSMART’s paper used for its V-SHARP comparison (Yaghmaie et al. 2021, §3.4).

In STI Suite’s GUI (probed headless; see crate::inversion::star), V-SHARP’s input is the sum of the echoes’ unwrapped phases from MRPhaseUnwrap on the raw, unmasked phase, its mask is the GUI’s magnitude-threshold mask (not BET), and its eroded mask and output go to QSM_star or QSM_iLSQR with TE = the sum of the echo times.

Reference: Wu, B., Li, W., Guidon, A., Liu, C. (2012). “Whole brain susceptibility mapping using compressed sensing.” Magnetic Resonance in Medicine, 67(1):137-147. https://doi.org/10.1002/mrm.23000

Implementation follows QSM.jl (https://github.com/kamesy/QSM.jl; see QSM.jl’s vsharp); defaults follow STI Suite 3.0’s V_SHARP (see Defaults), except the deconvolution threshold, which only QSM.jl has.

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): sti.unwrap-laplacian, sti.laplacian-solver, sti.echo-combination, sti.echo-weights, sti.bg-vsharp, sti.vsharp-no-deconvolve, sti.vsharp-radius, sti.vsharp-grid, sti.vsharp-min-radius, sti.vsharp-radius-step.

Structs

VsharpParams

V-SHARP algorithm parameters

Enums

VsharpBall

The SMV ball V-SHARP filters with. See the module docs.

VsharpGrid

Grid V-SHARP estimates the background field on. See the module docs.

Constants

CROP_MARGIN

Zero margin (native voxels per side) around the mask’s bounding box when V-SHARP runs on an isotropic working grid. See the module docs.

STI_GRID_MM

Working-grid spacing of STI Suite’s V_SHARP, in mm: the default of VsharpParams::grid.

STI_MIN_RADIUS_MM

Smallest radius in effect in STI Suite’s V_SHARP, in mm (it drops voxels only r = 1 fits): the default of VsharpParams::min_radius.

STI_RADIUS_STEP_MM

Step between STI Suite V_SHARP’s radii, in mm: the default of VsharpParams::radius_step.

Functions

bounding_box 🔒

Inclusive bounding box of the nonzero mask voxels, or None for an empty mask.

sti_even_box 🔒

STI Suite’s MaskBoundingBox(mask, voxelsize) (V_SHARP’s crop): the mask’s bounding box made even along every axis. An odd extent grows by one voxel at its high end, or at its low end when the high end is the array’s last voxel; an odd axis the mask fills entirely loses its last voxel. The voxel size plays no part. Probed black-box in licensed MATLAB R2023b on 17,952 boxes (lengths 8–40, eleven voxel sizes 0.5–3 mm, every position), all reproduced.

trilinear 🔒

Separable trilinear interpolation (MATLAB interp3(..., 'linear')): coords[d] are 1-based fractional source indices along axis d, one output sample per entry, output in the same Fortran order. Samples outside [1, n] on any axis are 0.

trilinear_mask_f32 🔒

trilinear of a 0/1 mask as MATLAB’s interp3(..., 'linear') computes it on a single array, which is how STI Suite’s V_SHARP resamples its mask both ways: each fraction rounded to f32, then (1 − f)·v₀ + f·v₁ in f32, along axis 1 first (interp3’s X, the second array dimension), then axis 0, then axis 2.

vsharp

V-SHARP background field removal

vsharp_isotropic 🔒

V-SHARP with the background estimated on an isotropic grid of spacing h mm anchored at the mask’s bounding box (see the module docs).

vsharp_padded 🔒

V-SHARP on QSM.jl’s padded box.

vsharp_radii 🔒

Build the descending list of SMV kernel radii (mm) from the params.

vsharp_with_radii 🔒

V-SHARP with an explicit list of SMV kernel radii (mm).