Module tgv 

TGV-QSM: Total Generalized Variation for Quantitative Susceptibility Mapping

Single-step QSM: background removal and dipole inversion (and, for a wrapped echo, phase unwrapping) in one primal-dual iteration with TGV regularisation.

References: Langkammer, C., Bredies, K., Poser, B.A., et al. (2015). “Fast quantitative susceptibility mapping using 3D EPI and total generalized variation.” NeuroImage, 111:622-630. https://doi.org/10.1016/j.neuroimage.2015.02.041

Chatnuntawech, I., McDaniel, P., et al. (2017). “Single-step quantitative susceptibility mapping with variational penalties.” NMR in Biomedicine, 30(4):e3570. https://doi.org/10.1002/nbm.3570

Reference implementation: https://github.com/korbinian90/QuantitativeSusceptibilityMappingTGV.jl (qsm_tgv). The iteration, step sizes, dipole stencil, mask erosions, Laplacian mean correction and the defaults (α from get_default_alpha(2), iterations from get_default_iterations, erosions = 3, step_size = 3) are that package’s.

Units

The reference solves

  min_{χ,w,ψ}  α₁‖∇χ − w‖₁ + α₀‖ε(w)‖₁ + ½‖ψ‖²   s.t.  ∇²φ = (⅓∇² − ∂²_z)χ + ∇²ψ

on the phase Laplacian ∇²φ in rad/mm², with α in those units, and divides χ by c = 2π·γ·B0·TE (rad_per_ppm) at the end. QSM.rs does the same: the input (a field in ppm, or one wrapped echo in rad) becomes phase in rad at the given TE and B0 as the first step, and χ is divided by c as the last. α are therefore in the reference’s rad units and mean what they mean there, at the scan’s actual TE and B0. What the caller has decides how ∇²φ is formed (see TgvInput):

• TgvInput::FieldPpm: an unwrapped field in ppm (for example a multi-echo field map), scaled to phase at the echo time TgvParams::te picks. ∇²φ is the ordinary discrete Laplacian: the reference’s get_laplace_phase_romeo with the unwrapping already done.
• TgvInput::WrappedPhase: one wrapped echo in rad, the original single-step TGV-QSM. ∇²φ is the wrap-insensitive sin/cos (DEL) Laplacian, the reference’s default get_laplace_phase_del.

The output is χ in ppm.

Agreement with the reference

On the QSM-CI in-vivo travelling-head scan (Siemens Prisma 2.89 T, 1 mm, 184×256×128, ROMEO total field, brain mask touching the first z plane) against the package in the QSM-CI julia-qsm image (sha256:b3c048ac…, 6c0886a, 4 threads), both at the package’s defaults (1655 iterations) and with both TgvParams::reference_divergence and TgvParams::reference_edge_erosion, over the voxels both return: NRMSE 4.0e-4 % for a wrapped echo (the field as phase at TE1, the package’s get_laplace_phase_del), and 3.7e-4 % for TgvInput::FieldPpm against the package with its own discrete Laplacian kernel on the unwrapped phase. Both are float32 rounding over 1655 iterations. tests/tgv_tgvjl_parity.rs checks both inputs, anisotropic voxels and oblique B0 against the package’s output on small phantoms (relative L2 ≤ 3e-6), with all three reference options on (the comparison covers the package’s whole output box). The defaults differ; see below.

Against the package’s get_laplace_phase_romeo, which ROMEO-unwraps the given phase again before its discrete Laplacian, the in-vivo field differs by NRMSE 12 %: on that scan the re-unwrap adds ±2π at 2 voxels inside the data-term mask, which spreads through the inversion. TgvInput::FieldPpm takes the Laplacian of the field as given.

Where QSM.rs differs from the reference

Two details of the reference code are fixed by default, and one of its outputs is trimmed. Each is an independent option that restores the reference’s behaviour (QSM.rs’s reference rule: unequivocal fixes are the default, and the original’s behaviour stays available for parity):

• Divergence in the w update (TgvParams::reference_divergence). The reference scales each row of q by that row’s own resolution for all three of its differences (update_w_kernel!; the same in the original Python tgv_deconv_update_w). Chambolle–Pock converges to the problem the paper states only when the operator in the primal update is the exact adjoint of the one in the dual update. On anisotropic voxels the reference’s is not the adjoint of its symmetric gradient; a unit test checks the identity for the default and its failure for the reference’s. By default the w update uses the exact adjoint. This is a pure correctness fix: it affects only anisotropic voxels, and on isotropic voxels the two are the same operator.
• Erosion at the array’s edge (TgvParams::reference_edge_erosion). The reference skips neighbours outside the array (erode_mask, src/tgv.jl:218-222; the original Python’s erode_mask treats them as inside), so a mask that reaches an array face is not eroded from it. There the Laplacian and gradient stencils have no valid neighbour (the reference replicates the edge value or drops the term), so the data term at those voxels is not the operator it is elsewhere. By default a neighbour outside the array counts as outside the mask, which removes those voxels. This affects only masks that reach an array face, and costs coverage at those faces.
• Output support (TgvParams::reference_output_box). χ is returned inside the eroded mask (erosions erosions of mask) and is zero elsewhere. The reference returns its whole working box (the eroded mask’s bounding box plus one voxel, reduce_to_mask_box), whose χ outside that mask is the stencils’ spill-over where the data term does not act: a rim of up to two voxels around the eroded mask (no data there; e.g. 44,530 voxels of the in-vivo brain mask). Inside the eroded mask the two are the same. By default QSM.rs trims that rim, since no data constrain it; the option returns it.
  • What it changes in a comparison with the package. Inside an uneroded brain mask the rim is part of the comparison. QSM-CI’s tgv pair (QSMxT/QSM-CI tgv-tgvjl vs tgv-qsmrs, compared inside the dataset’s mask, QSM.rs main 7f8271e, run 41019310531) measured NRMSE 3.84 % (sim, 41,907 rim voxels) and 6.74 % (KKI travelling head, 46,101): all of it the rim. Over the voxels both return, the same outputs differ by 3.8e-4 % (sim) and 4.9e-3 % (KKI): float32 rounding over 1655 iterations (#369).

Measured against the package and, on the QSM-CI phantom, against the true χ. The total field is the input; the defaults are QSM.rs’s, and “both” turns both options on. NRMSE against the package is over the voxels both return. NRMSE and corr against the truth are over the voxels every setting returns.

case	mask reaches	default	divergence only	edge erosion only	both
in vivo, 1 mm (travelling head, brain mask)	z = 0	11.0 %	11.0 %	3.7e-4 %	3.7e-4 %
in vivo, 1 × 1 × 2 mm (same, z pairs averaged)	z = 0	14.6 %	14.4 %	2.3 %	2.0e-4 %
phantom, 1 mm (CI brain mask)	no face	3.8e-4 %	3.8e-4 %	3.8e-4 %	3.8e-4 %
phantom slab, 48 slices at 1 mm	both z faces	99.8 %	99.8 %	4.4e-4 %	4.4e-4 %
phantom slab, 48 slices at 1 × 1 × 2 mm	both z faces	97.9 %	97.9 %	1.7 %	1.1e-4 %

Against the true χ:

case	default	divergence only	edge erosion only	both (= the package)
phantom, 1 mm	59.4 % / 0.807	same	same	same
phantom slab, 1 mm	68.5 % / 0.730	68.5 % / 0.730	1334 % / 0.034	1334 % / 0.034
phantom slab, 1 × 1 × 2 mm	65.6 % / 0.757	65.5 % / 0.757	349 % / 0.082	347 % / 0.082

• Coverage. The edge erosion’s cost is the voxels it removes at the faces: about 1,010 of the in-vivo scan’s 1.5 million (one face), and 70,418 of 617,917 on the 1 mm phantom slab (two faces, 11 %). On a head mask (χ ≠ 0) of the phantom reaching three faces, it removes 308,014 voxels.
• Why it is still the default. When the mask is cut by a face, the reference’s unmasked stencils there spread through the whole slab: χ against the truth goes from NRMSE 69 % to 1334 %.
• Divergence. It changes χ only on anisotropic voxels: by 1.1 % (relative L2) on a 0.9 × 1.0 × 1.3 mm phantom. On the 1 × 1 × 2 mm in-vivo scan with the reference’s erosion, the package’s divergence accounts for the 2.3 % that separates “edge erosion only” from the package. Against the truth it does not help either way: on the 1 × 1 × 2 mm slab the package’s divergence scores 0.1 points lower NRMSE (65.5 % vs 65.6 %, corr equal). The default stays the exact adjoint because that is the operator the TGV problem needs, not for a measured gain.

The phantom cases use QSM-CI’s in-silico phantom (QSM Challenge 2.0 head model, 164 × 205 × 205 at 1 mm, 7 T) with its true total field. They were run in the QSM.rs Actions on the CI’s data (QSM.rs#324). The slabs are 48 slices centred on the brain, cut from the field, the mask and the truth; the 1 × 1 × 2 mm slab averages pairs of 96 slices.

Optimizations:

• Bounding box reduction: only process the region containing the mask
• Pre-allocated buffers: all temporary arrays allocated once outside the loop

Reference pipeline

QuantitativeSusceptibilityMappingTGV.jl, commit 6c0886a:

• Input: one echo’s phase in rad, its TE and the field strength, plus a brain mask (README “Example to run TGV”; qsm_tgv(phase, mask, res; TE, fieldstrength, ...), src/tgv.jl L10). No unwrapping or background removal beforehand: the Laplacian of the phase does both (default get_laplace_phase_del; get_laplace_phase_romeo is the option for an unwrapped phase). χ is scaled by 1/(2π·TE·B0·42.5781) at the end (L224-228).
• Mask: erosions = 3 erosions of the given mask with a diamond (6-neighbour) element, then one more for the data-term mask mask0 (L14-20, L218-222). The README’s “Masking for QSM” points to phase-quality masking (ROMEO voxel quality, robustmask(romeovoxelquality(...))) combined with a brain mask.
• Resolution and orientation: no resampling; anisotropic voxels and oblique B0_dir are handled directly (README intro; the 27-point oblique stencil, L160). The default iteration count scales with the voxel volume (get_default_iterations, L108-115), and the README says the defaults “should lead to good results independently of the acquired resolution”. dedimensionalize = false by default.
• Defaults: regularization = 2 (α₀ = 0.002, α₁ = 0.003), step_size = 3, correct_laplacian = true. The settings that reproduce the original Python version are step_size = 1, iterations = 1000, alpha = [0.003, 0.001], erosions = 5, laplacian = get_laplace_phase3, correct_laplacian = false (README).
• Langkammer et al. 2015 (abstract): 3D EPI at 3T and 7T, 1 mm isotropic at 3T.

Structs

BoundingBox 🔒

Bounding box for mask region

TgvParams

TGV parameters. alpha0/alpha1 are in the reference’s rad units (see the module docs).

TgvWorkspace 🔒

Pre-allocated workspace for TGV iteration

Enums

TgvInput

What TGV-QSM reconstructs from. Both carry the echo time and field strength that turn the input into phase in rad, the units the reference (and α) work in.

Constants

TGV_GAMMA_MHZ_PER_T

γ/2π in MHz/T, as used by TGV-QSM (Julia and Python reference implementations): the reference implementation’s value, not crate::constants. Keep it.

Functions

apply_laplacian_inplace 🔒

Apply the Laplacian (Neumann boundary). With mask, the output is zero outside it; without, it is computed at every voxel.

apply_stencil 🔒

Apply dipole stencil to a 3D volume Uses Neumann BC at boundaries (matching Julia’s wave_local)

compute_dipole_stencil

Compute dipole stencil (27-point spatial kernel)

compute_field_laplacian

Laplacian of an unwrapped field (any units), zero outside mask.

compute_laplacian 🔒

Compute discrete Laplacian of a 3D array

compute_oblique_stencil

Compute SVD-fitted oblique 27-point dipole stencil via DST-based Poisson solves.

compute_operator_norm_sqr 🔒

Compute the squared spectral norm of the TGV operator matrix via power iteration

compute_phase_laplacian

Compute Laplacian of wrapped phase (rad) using the DEL method.

dst3d 🔒

3D separable DST-I, unnormalised: X[k] = Σ x[j]·sin(π·(j+1)·(k+1)/(N+1)) along each axis.

dst_sin_table 🔒

Build sin table for DST-I of given length. sin_table[n][k] = sin(π*(n+1)*(k+1)/(N+1))

erode_mask 🔒

One erosion with the 6-neighbour (strel_diamond) element. A neighbour outside the array counts as outside the mask, so a mask that reaches the array’s edge is eroded from it (the default: the Laplacian and gradient stencils have no neighbour there). With skip_outside, neighbours outside the array are skipped instead, as the reference’s erode_mask does (src/tgv.jl:218-222, 6c0886a, checkbounds), so such a mask is not eroded from the edge (TgvParams::reference_edge_erosion).

extract_subvolume 🔒

Extract sub-volume from full volume

frobenius_norm 🔒

Frobenius norm of symmetric 3x3 tensor (6 components)

get_default_alpha

Default (alpha0, alpha1) for a regularization strength (the reference’s get_default_alpha): max([0.001, 0.001] + [0.001, 0.002]·(regularization − 1), 0). 2 is the default, giving α₀ = 0.002, α₁ = 0.003.

get_default_iterations

Get default number of iterations based on voxel size and step size.

grad_norm_sq 🔒

Compute gradient norm squared

insert_subvolume 🔒

Insert sub-volume back into full volume

norm3 🔒

L2 norm of 3-component vector

project_linf3 🔒

L-infinity projection for 3-component vector

project_linf6 🔒

L-infinity projection for 6-component symmetric tensor

rad_per_ppm

Phase (rad) per ppm of field at echo time te_s (s) and field strength b0_t (T): 2π·γ·B0·TE with γ/2π = TGV_GAMMA_MHZ_PER_T.

solve_symmetric_pseudoinverse 🔒

Solve a symmetric positive semi-definite system via eigendecomposition with thresholding.

symdiv_rows_f32 🔒

The divergence of the symmetric tensor q in the w update, as the reference computes it (update_w_kernel!, src/tgv_helper.jl:88-105, 6c0886a; the same in the original Python tgv_deconv_update_w, qsm_tgv_cython_helper.pyx): row d of the tensor ((qxx, qxy, qxz) for w_x, …) takes its backward differences along x, y and z all scaled by 1/res_d, the row’s own resolution, not by the resolution of the axis each difference is along. On isotropic voxels that is the ordinary divergence; on anisotropic voxels it is not the adjoint of the symmetric gradient used in the q update. Used only with TgvParams::reference_divergence; the default is the adjoint ([crate::utils::gradient::symdiv_inplace_f32]). See the module docs.

tgv_qsm

Main TGV-QSM reconstruction