Merge branch 'master' into display_dens

src/mrinufft/io/siemens.py

@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import numpy as np
+from .utils import siemens_quat_to_rot_mat
 
 
 def read_siemens_rawdat(
@@ -73,6 +74,7 @@ def read_siemens_rawdat(
         "n_adc_samples": int(twixObj.image.NCol),
         "n_slices": int(twixObj.image.NSli),
         "n_average": int(twixObj.image.NAve),
+        "orientation": siemens_quat_to_rot_mat(twixObj.image.slicePos[0][-4:]),
     }
     if slice_num is not None and hdr["n_slices"] < slice_num:
         raise ValueError("The slice number is out of bounds.")

src/mrinufft/io/utils.py

@@ -1,6 +1,7 @@
 """Module containing utility functions for IO in MRI NUFFT."""
 
 import numpy as np
+from scipy.spatial.transform import Rotation
 
 
 def add_phase_to_kspace_with_shifts(kspace_data, kspace_loc, normalized_shifts):
@@ -35,3 +36,27 @@ def add_phase_to_kspace_with_shifts(kspace_data, kspace_loc, normalized_shifts):
     phi = np.sum(kspace_loc * normalized_shifts, axis=-1)
     phase = np.exp(-2 * np.pi * 1j * phi)
     return kspace_data * phase
+
+
+def siemens_quat_to_rot_mat(quat):
+    """
+    Calculate the rotation matrix from Siemens Twix quaternion.
+
+    Parameters
+    ----------
+    quat : np.ndarray
+        The quaternion from the Siemens Twix file.
+
+    Returns
+    -------
+    np.ndarray
+        The affine rotation matrix which is a 4x4 matrix.
+        This can be passed as input to `affine` parameter in `nibabel`.
+    """
+    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:
+        R[2] = -R[2]
+    R[-1, -1] = 1
+    return R

src/mrinufft/operators/interfaces/gpunufft.py

@@ -590,12 +590,12 @@ def pipe(
             The oversampling factor the volume shape
         normalize: bool
             Whether to normalize the density compensation.
-            We normalize such that the energy of PSF = 1
         """
         if GPUNUFFT_AVAILABLE is False:
             raise ValueError(
                 "gpuNUFFT is not available, cannot " "estimate the density compensation"
             )
+        original_shape = volume_shape
         volume_shape = (np.array(volume_shape) * osf).astype(int)
         grid_op = MRIGpuNUFFT(
             samples=kspace_loc,
@@ -607,11 +607,10 @@ def pipe(
             max_iter=num_iterations
         )
         if normalize:
-            spike = np.zeros(volume_shape)
-            mid_loc = tuple(v // 2 for v in volume_shape)
-            spike[mid_loc] = 1
-            psf = grid_op.adj_op(grid_op.op(spike))
-            density_comp /= np.linalg.norm(psf)
+            test_op = MRIGpuNUFFT(samples=kspace_loc, shape=original_shape, **kwargs)
+            test_im = np.ones(original_shape, dtype=np.complex64)
+            test_im_recon = test_op.adj_op(density_comp * test_op.op(test_im))
+            density_comp /= np.mean(np.abs(test_im_recon))
         return density_comp.squeeze()
 
     def get_lipschitz_cst(self, max_iter=10, tolerance=1e-5, **kwargs):