Module tgv 

TGV-QSM: Total Generalized Variation for Quantitative Susceptibility Mapping

Single-step QSM reconstruction from wrapped phase data using TGV regularization.

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

The algorithm solves: min_χ ||∇²(phase) - D*χ||₂² + α₁||∇χ - w||₁ + α₀||ε(w)||₁

where:

• χ is the susceptibility map
• w is an auxiliary vector field (velocity)
• ε(w) is the symmetric gradient of w
• D is the dipole kernel
• α₀, α₁ are TGV regularization parameters

Optimizations:

• Bounding box reduction: only process the region containing the mask
• Pre-allocated buffers: all temporary arrays allocated once outside the loop
• Early termination: convergence check every 100 iterations

Structs

BoundingBox 🔒

Bounding box for mask region

TgvParams

TGV parameters

TgvWorkspace 🔒

Pre-allocated workspace for TGV iteration

Functions

apply_laplacian_inplace 🔒

Apply Laplacian with mask

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_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 using the DEL method

compute_relative_change 🔒

Compute relative change for convergence check

dst1 🔒

1D DST-I: X[k] = Σ x[n] * sin(π*(n+1)*(k+1)/(N+1))

dst3d 🔒

3D separable DST-I (Type I Discrete Sine Transform).

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

Erode mask by one voxel (6-connected)

extract_subvolume 🔒

Extract sub-volume from full volume

frobenius_norm 🔒

Frobenius norm of symmetric 3x3 tensor (6 components)

get_default_alpha

Get default alpha values based on regularization level (1-4)

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

solve_symmetric_pseudoinverse 🔒

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

tgv_qsm

Main TGV-QSM reconstruction

tgv_qsm_with_progress

TGV-QSM with progress callback (optimized version)