Skip to content

tensor_math

tit.pre.qsi.tensor_math

Tensor and transform algebra for the DTI -> SimNIBS path (numpy only).

Conventions:

  • "world" is scanner RAS in mm (the NIfTI sform), the frame SimNIBS meshes live in. A tensor field held in memory is always a world-frame (..., 3, 3) array.
  • SimNIBS cond2elmdata(correct_FSL=True) reads a stored 6-vector S as T_world = M S M^T with M from :func:simnibs_fsl_matrix. We store S = inv(M) T_world inv(M)^T with that exact M, so the round trip is exact for any output affine, isotropic or not.
  • An ITK/ANTs transform named from-A_to-B maps points of B (the fixed image) to A. ITK stores it in LPS as y = R (x - c) + c + t.

sym 

sym(t6: ndarray) -> ndarray

(..., 6) upper triangle in FSL order -> (..., 3, 3) symmetric tensors.

Source code in tit/pre/qsi/tensor_math.py 
25
26
27
28
29
30
31
def sym(t6: np.ndarray) -> np.ndarray:
    """``(..., 6)`` upper triangle in FSL order -> ``(..., 3, 3)`` symmetric tensors."""
    t6 = np.asarray(t6)
    out = np.empty(t6.shape[:-1] + (3, 3), dtype=np.float64)
    for k, (i, j) in enumerate(FSL_ORDER):
        out[..., i, j] = out[..., j, i] = t6[..., k]
    return out

unsym 

unsym(tensors: ndarray) -> ndarray

(..., 3, 3) -> (..., 6) upper triangle in FSL order.

Source code in tit/pre/qsi/tensor_math.py 
34
35
36
def unsym(tensors: np.ndarray) -> np.ndarray:
    """``(..., 3, 3)`` -> ``(..., 6)`` upper triangle in FSL order."""
    return np.stack([tensors[..., i, j] for i, j in FSL_ORDER], axis=-1)

rotate 

rotate(tensors: ndarray, matrix: ndarray) -> ndarray

matrix @ T @ matrix.T for every tensor in a field.

Source code in tit/pre/qsi/tensor_math.py 
39
40
41
def rotate(tensors: np.ndarray, matrix: np.ndarray) -> np.ndarray:
    """``matrix @ T @ matrix.T`` for every tensor in a field."""
    return np.einsum("ij,...jk,lk->...il", matrix, tensors, matrix)

simnibs_fsl_matrix 

simnibs_fsl_matrix(affine: ndarray) -> ndarray

Exactly SimNIBS 4.6 cond_utils.cond2elmdata with correct_FSL=True.

M = A[:3, :3] / colnorm[:, None] (rows divided by the column norms, as SimNIBS writes it), then M = M @ diag(-1, 1, 1) when det(M) > 0.

Source code in tit/pre/qsi/tensor_math.py 
44
45
46
47
48
49
50
51
52
53
54
def simnibs_fsl_matrix(affine: np.ndarray) -> np.ndarray:
    """Exactly SimNIBS 4.6 ``cond_utils.cond2elmdata`` with ``correct_FSL=True``.

    ``M = A[:3, :3] / colnorm[:, None]`` (rows divided by the column norms, as SimNIBS
    writes it), then ``M = M @ diag(-1, 1, 1)`` when ``det(M) > 0``.
    """
    lin = np.asarray(affine, dtype=np.float64)[:3, :3]
    matrix = lin / np.linalg.norm(lin, axis=0)[:, None]
    if np.linalg.det(matrix) > 0:
        matrix = matrix @ np.diag([-1.0, 1.0, 1.0])
    return matrix

world_to_simnibs 

world_to_simnibs(tensors_world: ndarray, affine: ndarray) -> ndarray

World tensors -> the 6-vectors SimNIBS will read back as those world tensors.

Source code in tit/pre/qsi/tensor_math.py 
57
58
59
def world_to_simnibs(tensors_world: np.ndarray, affine: np.ndarray) -> np.ndarray:
    """World tensors -> the 6-vectors SimNIBS will read back as those world tensors."""
    return unsym(rotate(tensors_world, np.linalg.inv(simnibs_fsl_matrix(affine))))

simnibs_to_world 

simnibs_to_world(t6: ndarray, affine: ndarray) -> ndarray

Stored 6-vectors -> the world tensors SimNIBS builds from them.

Source code in tit/pre/qsi/tensor_math.py 
62
63
64
def simnibs_to_world(t6: np.ndarray, affine: np.ndarray) -> np.ndarray:
    """Stored 6-vectors -> the world tensors SimNIBS builds from them."""
    return rotate(sym(t6), simnibs_fsl_matrix(affine))

itk_euler_to_ras 

itk_euler_to_ras(params, center) -> ndarray

ITK transform parameters -> 4x4 RAS point map.

params is an Euler3DTransform (6 values: angles x, y, z in radians, then the translation; rotation R = Rz @ Rx @ Ry, ITK's default ZXY order) or an AffineTransform (12 values: the matrix row-major, then the translation). center is the ITK fixed parameter. Both are in LPS; the result is in RAS.

Source code in tit/pre/qsi/tensor_math.py 
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def itk_euler_to_ras(params, center) -> np.ndarray:
    """ITK transform parameters -> 4x4 RAS point map.

    *params* is an ``Euler3DTransform`` (6 values: angles x, y, z in radians, then the
    translation; rotation ``R = Rz @ Rx @ Ry``, ITK's default ZXY order) or an
    ``AffineTransform`` (12 values: the matrix row-major, then the translation).
    *center* is the ITK fixed parameter. Both are in LPS; the result is in RAS.
    """
    p = np.asarray(params, dtype=np.float64).ravel()
    c = np.asarray(center, dtype=np.float64).ravel()[:3]
    if p.size == 6:
        ax, ay, az = p[:3]
        rx = np.array(
            [[1, 0, 0], [0, np.cos(ax), -np.sin(ax)], [0, np.sin(ax), np.cos(ax)]]
        )
        ry = np.array(
            [[np.cos(ay), 0, np.sin(ay)], [0, 1, 0], [-np.sin(ay), 0, np.cos(ay)]]
        )
        rz = np.array(
            [[np.cos(az), -np.sin(az), 0], [np.sin(az), np.cos(az), 0], [0, 0, 1]]
        )
        rot, trans = rz @ rx @ ry, p[3:6]
    elif p.size == 12:
        rot, trans = p[:9].reshape(3, 3), p[9:12]
    else:
        raise ValueError(
            f"expected 6 (Euler) or 12 (affine) ITK parameters, got {p.size}"
        )
    lps = np.eye(4)
    lps[:3, :3] = rot
    lps[:3, 3] = c + trans - rot @ c
    return _LPS_FLIP @ lps @ _LPS_FLIP

qsiprep_anat_affines 

qsiprep_anat_affines(affine_raw: ndarray, shape) -> tuple[ndarray, ndarray]

Reproduce QSIPrep's Conform(deoblique_header=True): (A_lps, A_anat).

QSIPrep reorients the input T1w to LPS (A_lps), then replaces an oblique affine with diag(zooms * sign(diag)) keeping the translation column (A_anat). The voxel array is unchanged, so raw world = A_lps @ inv(A_anat) @ anat world. No transform file records this step. Assumes one input T1w and no resize, QSIPrep's behaviour for a single T1w.

Source code in tit/pre/qsi/tensor_math.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
def qsiprep_anat_affines(
    affine_raw: np.ndarray, shape
) -> tuple[np.ndarray, np.ndarray]:
    """Reproduce QSIPrep's ``Conform(deoblique_header=True)``: ``(A_lps, A_anat)``.

    QSIPrep reorients the input T1w to LPS (``A_lps``), then *replaces* an oblique
    affine with ``diag(zooms * sign(diag))`` keeping the translation column
    (``A_anat``). The voxel array is unchanged, so
    ``raw world = A_lps @ inv(A_anat) @ anat world``. No transform file records this
    step. Assumes one input T1w and no resize, QSIPrep's behaviour for a single T1w.
    """
    affine_raw = np.asarray(affine_raw, dtype=np.float64)
    target_sign = (-1, -1, 1)  # L, P, S
    new_to_old = np.zeros((4, 4))
    new_to_old[3, 3] = 1.0
    for old_axis, (world, sign) in enumerate(_axis_codes(affine_raw)):
        if sign == target_sign[world]:
            new_to_old[old_axis, world] = 1.0
        else:
            new_to_old[old_axis, world] = -1.0
            new_to_old[old_axis, 3] = shape[old_axis] - 1
    a_lps = affine_raw @ new_to_old
    zooms = np.linalg.norm(a_lps[:3, :3], axis=0)
    a_anat = a_lps.copy()
    a_anat[:3, :3] = np.diag(zooms * np.sign(np.diag(a_lps[:3, :3])))
    return a_lps, a_anat

acpc_to_t1_world 

acpc_to_t1_world(affine_raw: ndarray, shape, acpc_to_anat_ras) -> ndarray

4x4 map from QSIPrep ACPC world to raw-T1 (= charm m2m T1) world.

G = A_lps @ inv(A_anat) @ inv(X), where acpc_to_anat_ras is X, the RAS form of QSIPrep's from-ACPC_to-anat transform (it maps anat points to ACPC).

Source code in tit/pre/qsi/tensor_math.py 
149
150
151
152
153
154
155
156
def acpc_to_t1_world(affine_raw: np.ndarray, shape, acpc_to_anat_ras) -> np.ndarray:
    """4x4 map from QSIPrep ACPC world to raw-T1 (= charm m2m T1) world.

    ``G = A_lps @ inv(A_anat) @ inv(X)``, where *acpc_to_anat_ras* is ``X``, the RAS
    form of QSIPrep's ``from-ACPC_to-anat`` transform (it maps anat points to ACPC).
    """
    a_lps, a_anat = qsiprep_anat_affines(affine_raw, shape)
    return a_lps @ np.linalg.inv(a_anat) @ np.linalg.inv(acpc_to_anat_ras)

polar_rotation 

polar_rotation(matrix: ndarray) -> ndarray

Rotation factor of a 3x3 (or 4x4) linear map (exact for a rigid map).

Source code in tit/pre/qsi/tensor_math.py 
159
160
161
162
def polar_rotation(matrix: np.ndarray) -> np.ndarray:
    """Rotation factor of a 3x3 (or 4x4) linear map (exact for a rigid map)."""
    u, _, vt = np.linalg.svd(np.asarray(matrix)[:3, :3])
    return u @ vt

eig_desc 

eig_desc(tensors: ndarray) -> tuple[ndarray, ndarray]

Eigenvalues (descending) and matching eigenvectors (columns).

Source code in tit/pre/qsi/tensor_math.py 
165
166
167
168
def eig_desc(tensors: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Eigenvalues (descending) and matching eigenvectors (columns)."""
    vals, vecs = np.linalg.eigh(tensors)
    return vals[..., ::-1], vecs[..., :, ::-1]

fa_md 

fa_md(eigenvalues: ndarray) -> tuple[ndarray, ndarray]

Fractional anisotropy and mean diffusivity from (..., 3) eigenvalues.

Source code in tit/pre/qsi/tensor_math.py 
171
172
173
174
175
176
177
def fa_md(eigenvalues: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Fractional anisotropy and mean diffusivity from ``(..., 3)`` eigenvalues."""
    md = eigenvalues.mean(-1)
    num = np.sqrt(((eigenvalues - md[..., None]) ** 2).sum(-1))
    den = np.sqrt((eigenvalues**2).sum(-1))
    fa = np.where(den > 0, np.sqrt(1.5) * num / np.where(den > 0, den, 1.0), 0.0)
    return fa, md