account for slice dir inversion

small fixes

please the linter

the linter is stupid

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,11 +58,15 @@ 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.
@@ -78,51 +84,56 @@ def nifti_affine(twixObj):
     """
     # 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 = siemens_quat_to_rot_mat(twixObj.image.slicePos[0][-4:])
+        "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 = 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,
-            ])
+        ]
+    )
+
+    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]
 
-    center = [-fov[0]/2,
-              -fov[1]/2,
-              -(fov[2] - (my[sos] * dthick if sos in my else 0))/2,
-              1]
+    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.