Module sti 

STI: susceptibility tensor imaging

COSMOS combines orientations to make a scalar susceptibility well-posed. STI keeps the orientations but changes the model: susceptibility becomes a symmetric rank-2 tensor, so tissue whose magnetic response depends on direction — white matter, whose myelin sheaths are diamagnetic along the fibre — is described rather than averaged away.

Each orientation measures one projection of the tensor onto its own B0 direction:

    delta_t(k) = (1/3) H_t^T X(k) H_t  -  (H_t . k) (k^T X(k) H_t) / |k|^2

With X symmetric there are six unknowns per voxel — [X11, X12, X13, X22, X23, X33] — so at least six orientations are needed, and in practice ten or more are acquired for conditioning. Setting X = chi I collapses the expression to chi (1/3 - (H.k)^2/|k|^2), the ordinary dipole kernel, which is the sense in which STI contains scalar QSM.

Why this is a per-voxel solve rather than a global one

The reference implementations take two routes that look very different. STI Suite solves a small least-squares problem at each k-point; Bilgic’s toolbox builds a six-channel operator and runs LSQR over the whole volume. They solve the same system: the operator couples the six tensor channels to the N orientations at each k-point independently, with no spatial coupling anywhere, so the global system is block-diagonal in k and separates into N_k independent six-unknown problems. The per-k solve is therefore the exact answer, and the global LSQR an early-stopped approximation of it.

This module takes the direct route: form the 6x6 normal equations at each k-point and solve by Cholesky. Exact, embarrassingly parallel, and with no iteration count to tune. lambda adds a ridge term for the ill-conditioned case.

Units and gauge

Local (tissue) field in ppm in, tensor components in ppm out. As with the scalar dipole kernel there is no DC response, so the k = 0 term is null and every component is recovered zero-mean.

A difference from the scalar dipole kernel

[crate::kernels::dipole::dipole_kernel] keeps the Nyquist planes of an even-length axis; this module drops them (see is_nyquist_bin below). The asymmetry is deliberate. Scalar inversion divides by the kernel pointwise, so whatever convention the kernel picks at an ambiguous bin cancels exactly between forward and inverse. STI instead solves a six-way coupled system at each bin, where the same ambiguity leaves the components inseparable, so the bins have to be excluded rather than carried through. An isotropic tensor therefore reproduces the scalar dipole kernel exactly on an odd grid, and everywhere except those planes on an even one.

Memory

Every orientation’s k-space is needed simultaneously, so this holds N_orientations complex volumes at once — for twelve orientations of a 164x205x205 volume, about 1.3 GB.

References

Liu C. “Susceptibility tensor imaging.” Magnetic Resonance in Medicine. 2010;63(6):1471-7. https://doi.org/10.1002/mrm.22482

Li W, Liu C, Duong TQ, van Zijl PCM, Li X. “Susceptibility tensor imaging (STI) of the brain.” NMR in Biomedicine. 2017;30(4):e3540. https://doi.org/10.1002/nbm.3540

Reference implementations: STI_Parfor.m (STI Suite v3.0, Chunlei Liu) and apply_STI.m (COSMOS_STI_Toolbox, Berkin Bilgic). Ported from the published equations; the code was used only to resolve component ordering and sign conventions, which the two agree on.

Structs

StiParams

STI parameters.

SusceptibilityTensor

A symmetric susceptibility tensor field, in ppm.

TensorMaps

Rotationally invariant maps derived from the tensor’s eigen-decomposition.

Constants

N_TENSOR

Number of independent components of a symmetric 3x3 tensor.

Functions

cholesky6 🔒

Cholesky factorization of a 6x6 symmetric positive-definite matrix, in place. Returns false if the matrix is not positive definite.

cholesky6_solve 🔒

Solve L L^T x = b in place given the lower factor from cholesky6.

is_nyquist_bin 🔒

Whether a bin sits on a Nyquist plane, which only exists on even-length axes.

sti

Reconstruct the susceptibility tensor from multi-orientation local field maps.

sti_forward

Forward STI model: the local field each orientation would measure from tensor.

sti_row 🔒

Build the row of the STI system for one orientation at one k-point.

tensor_maps

Eigen-decompose the tensor and derive MMS, MSA and the principal eigenvector.