Merge pull request #266 from j-obriot/affine

Nifti affine matrix

Author: j-obriot

src/mrinufft/io/siemens.py

@@ -3,7 +3,7 @@
 from __future__ import annotations
 
 import numpy as np
-from .utils import siemens_quat_to_rot_mat
+from .utils import siemens_quat_to_rot_mat, nifti_affine
 
 
 def read_siemens_rawdat(
@@ -76,6 +76,7 @@ def read_siemens_rawdat(
         "n_slices": int(twixObj.image.NSli),
         "n_average": int(twixObj.image.NAve),
         "orientation": siemens_quat_to_rot_mat(twixObj.image.slicePos[0][-4:]),
+        "affine": nifti_affine(twixObj),
         "acs": None,
     }
     if "refscan" in twixObj.keys():

src/mrinufft/io/utils.py

@@ -38,14 +38,16 @@ def add_phase_to_kspace_with_shifts(kspace_data, kspace_loc, normalized_shifts):
     return kspace_data * phase
 
 
-def siemens_quat_to_rot_mat(quat):
+def siemens_quat_to_rot_mat(quat, return_det=False):
     """
     Calculate the rotation matrix from Siemens Twix quaternion.
 
     Parameters
     ----------
     quat : np.ndarray
         The quaternion from the Siemens Twix file.
+    return_det : bool
+        Whether to return the determinent of the rotation before norm
 
     Returns
     -------
@@ -56,12 +58,82 @@ def siemens_quat_to_rot_mat(quat):
     R = np.zeros((4, 4))
     R[:3, :3] = Rotation.from_quat([quat[1], quat[2], quat[3], quat[0]]).as_matrix()
     R[:, (0, 1)] = R[:, (1, 0)]
-    if np.linalg.det(R[:3, :3]) < 0:
+    det = np.linalg.det(R[:3, :3])
+    if det < 0:
         R[2] = -R[2]
     R[-1, -1] = 1
+    if return_det:
+        return R, det
     return R
 
 
+def nifti_affine(twixObj):
+    """
+    Calculate the affine transformation matrix from Siemens Twix object.
+
+    Parameters
+    ----------
+    twixObj : twixObj
+        The twix object returned by mapVBVD.
+
+    Returns
+    -------
+    np.ndarray
+        The affine transformation matrix which is a 4x4 matrix.
+        This can be passed as input to `affine` parameter in `nibabel`.
+    """
+    # required keys
+    keys = {
+        "dthick": ("sSliceArray", "asSlice", "0", "dThickness"),
+        "dread": ("sSliceArray", "asSlice", "0", "dReadoutFOV"),
+        "dphase": ("sSliceArray", "asSlice", "0", "dPhaseFOV"),
+        "lbase": ("sKSpace", "lBaseResolution"),
+        "lphase": ("sKSpace", "lPhaseEncodingLines"),
+        "ucdim": ("sKSpace", "ucDimension"),
+    }
+    sos = ("sKSpace", "dSliceOversamplingForDialog")
+    rot, det = siemens_quat_to_rot_mat(twixObj.image.slicePos[0][-4:], True)
+    my = twixObj.hdr.MeasYaps
+
+    for k in keys.keys():
+        if keys[k] not in my:
+            return rot
+
+    dthick = my[keys["dthick"]]
+    fov = np.array(
+        [
+            my[keys["dread"]],
+            my[keys["dphase"]],
+            dthick * (1 + my[sos] if sos in my else 1),
+        ]
+    )
+
+    lpart = ("sKSpace", "lPartitions")
+    res = np.array(
+        [
+            my[keys["lbase"]],
+            my[keys["lphase"]],
+            my[lpart] if my[keys["ucdim"]] == 4 and lpart in my else 1,
+        ]
+    )
+
+    scale = np.diag([*(fov / res), 1])
+
+    offset = twixObj.image.slicePos[0][:3]
+
+    fovz = fov[2] - (my[sos] * dthick if sos in my else 0)
+    center = [-fov[0] / 2, -fov[1] / 2, -fovz / 2, 1]
+
+    t = (rot @ center)[:3] - offset
+    if det < 0:
+        t[2] = (rot @ center)[2] * 2 - t[2]
+
+    full_mat = rot @ scale
+    full_mat[:3, 3] = t
+
+    return full_mat
+
+
 def remove_extra_kspace_samples(kspace_data, num_samples_per_shot):
     """Remove extra samples from k-space data.