Merge branch 'master' into add_prephasors

Author: chaithyagr

src/mrinufft/io/nsp.py

@@ -368,6 +368,7 @@ def read_trajectory(
     raster_time: float = DEFAULT_RASTER_TIME,
     read_shots: bool = False,
     normalize_factor: float = KMAX,
+    pre_skip: int = 0,
 ):
     """Get k-space locations from gradient file.
 
@@ -389,9 +390,13 @@ def read_trajectory(
     read_shots : bool, optional
         Whether in read shots configuration which accepts an extra
         point at end, by default False
-    normalize : float, optional
+    normalize_factor : float, optional
         Whether to normalize the k-space locations, by default 0.5
         When None, normalization is not done.
+    pre_skip: int, optional
+        Number of samples to skip from the start of each shot,
+        by default 0. This is useful when we want to avoid artifacts
+        from ADC switching in UTE sequences.
 
     Returns
     -------
@@ -522,6 +527,15 @@ def read_trajectory(
         if normalize_factor is not None:
             Kmax = img_size / 2 / fov
             kspace_loc = kspace_loc / Kmax * normalize_factor
+        if pre_skip > 0:
+            if pre_skip >= num_samples_per_shot:
+                raise ValueError(
+                    "skip_first_Nsamples should be less than num_adc_samples"
+                )
+            oversample_factor = num_adc_samples / num_samples_per_shot
+            skip_samples = pre_skip * int(oversample_factor)
+            kspace_loc = kspace_loc[:, skip_samples:]
+            params["num_adc_samples"] = num_adc_samples - skip_samples
         return kspace_loc, params
 
 
@@ -532,6 +546,7 @@ def read_arbgrad_rawdat(
     squeeze: bool = True,
     slice_num: int | None = None,
     contrast_num: int | None = None,
+    pre_skip: int = 0,
     data_type: str = "ARBGRAD_VE11C",
 ):  # pragma: no cover
     """Read raw data from a Siemens MRI file.
@@ -550,6 +565,10 @@ def read_arbgrad_rawdat(
         The slice to read, by default None. This applies for 2D data.
     contrast_num: int, optional
         The contrast to read, by default None.
+    pre_skip : int, optional
+        Number of samples to skip from the start of each shot,
+        by default 0. This is useful when we want to avoid artifacts
+        from ADC switching in UTE sequences.
     data_type : str, optional
         The type of data to read, by default 'ARBGRAD_VE11C'.
 
@@ -598,4 +617,12 @@ def read_arbgrad_rawdat(
                 "Phoenix", ("sFastImaging", "lTurboFactor")
             )[0]
             hdr["type"] = "ARBGRAD_MP2RAGE"
+    if pre_skip > 0:
+        samples_to_skip = int(hdr["oversampling_factor"] * pre_skip)
+        if samples_to_skip >= hdr["n_adc_samples"]:
+            raise ValueError(
+                "Samples to skip should be less than n_samples in the data"
+            )
+        data = data[:, :, samples_to_skip:]
+        hdr["n_adc_samples"] -= samples_to_skip
     return data, hdr

src/mrinufft/operators/interfaces/cufinufft.py

@@ -64,7 +64,7 @@ def __init__(
         self._kz = cp.array(samples[:, 2], copy=False) if self.ndim == 3 else None
         for i in [1, 2]:
             self._make_plan(i, **kwargs)
-            self._set_pts(i, samples)
+            self._set_pts(i)
 
     @property
     def dtype(self):
@@ -84,13 +84,15 @@ def _make_plan(self, typ, **kwargs):
             **kwargs,
         )
 
-    def _set_pts(self, typ, samples):
+    def _set_kxyz(self, samples):
+        self._kx.set(samples[:, 0])
+        self._ky.set(samples[:, 1])
+        if self.ndim == 3:
+            self._kz.set(samples[:, 2])
+
+    def _set_pts(self, typ):
         plan = self.grad_plan if typ == "grad" else self.plans[typ]
-        plan.setpts(
-            cp.array(samples[:, 0], copy=False),
-            cp.array(samples[:, 1], copy=False),
-            cp.array(samples[:, 2], copy=False) if self.ndim == 3 else None,
-        )
+        plan.setpts(self._kx, self._ky, self._kz)
 
     def _destroy_plan(self, typ):
         if self.plans[typ] is not None:
@@ -272,10 +274,11 @@ def samples(self, new_samples):
                 np.float32, copy=False
             )
         )
+        self.raw_op._set_kxyz(self._samples)
         for typ in [1, 2, "grad"]:
             if typ == "grad" and not self._grad_wrt_traj:
                 continue
-            self.raw_op._set_pts(typ, samples=self._samples)
+            self.raw_op._set_pts(typ)
         self.compute_density(self._density_method)
 
     @FourierOperatorBase.density.setter
@@ -808,7 +811,7 @@ def _make_plan_grad(self, **kwargs):
             isign=1,
             **kwargs,
         )
-        self.raw_op._set_pts(typ="grad", samples=self.samples)
+        self.raw_op._set_pts(typ="grad")
 
     def get_lipschitz_cst(self, max_iter=10, **kwargs):
         """Return the Lipschitz constant of the operator.

src/mrinufft/operators/off_resonance.py

@@ -198,6 +198,9 @@ class MRIFourierCorrected(FourierOperatorBase):
         Must have same shape as ``b0_map``.
         The default is ``None`` (purely imaginary field).
         Also supports Cupy arrays and Torch tensors.
+    negate: bool, optional, default=False
+        If True, negate the field map. Useful for matching the convention of
+        your field map generation.
     backend: str, optional
         The backend to use for computations. Either 'cpu', 'gpu' or 'torch'.
         The default is `cpu`.
@@ -221,6 +224,7 @@ def __init__(
         r2star_map=None,
         B=None,
         tl=None,
+        negate=False,
         backend="cpu",
     ):
         if backend == "gpu" and not CUPY_AVAILABLE:
@@ -235,7 +239,9 @@ def __init__(
             raise ValueError("Unsupported backend.")
 
         self._fourier_op = fourier_op
-
+        if not isinstance(negate, bool):
+            raise ValueError("negate must be a boolean value.")
+        self.isign = -1 if negate else 1
         self.n_coils = fourier_op.n_coils
         self.shape = fourier_op.shape
         self.smaps = fourier_op.smaps
@@ -275,7 +281,7 @@ def __init__(
             self.C = None
             self.field_map = field_map
         else:
-            self.C = _get_spatial_coefficients(field_map, self.tl)
+            self.C = _get_spatial_coefficients(field_map, self.tl, isign=self.isign)
             self.field_map = None
 
     def op(self, data, *args):
@@ -377,11 +383,10 @@ def _get_complex_fieldmap(b0_map, r2star_map=None):
     return field_map
 
 
-def _get_spatial_coefficients(field_map, tl):
+def _get_spatial_coefficients(field_map, tl, isign=-1):
     xp = get_array_module(field_map)
-
     # get spatial coeffs
-    C = xp.exp(-tl * field_map[..., None])
+    C = xp.exp(isign * tl * field_map[..., None])
     C = C[None, ...].swapaxes(0, -1)[
         ..., 0
     ]  # (..., n_time_segments) -> (n_time_segments, ...)

tests/test_offres_exp_approx.py

@@ -12,7 +12,7 @@
 from mrinufft._array_compat import CUPY_AVAILABLE
 from mrinufft._utils import get_array_module
 from mrinufft.operators.off_resonance import MRIFourierCorrected
-
+from mrinufft import get_operator
 
 from helpers import to_interface, assert_allclose
 from helpers.factories import _param_array_interface
@@ -106,3 +106,34 @@ def test_zmap_coeff(zmap, mask, array_interface):
     )
     actual = calculate_approx_offresonance_term(B, C)
     assert_allclose(actual, expected, atol=1e-3, rtol=1e-3, interface=array_interface)
+
+
+def test_b0_map_upsampling_warns_and_matches_shape():
+    """Test that MRIFourierCorrected upscales the b0_map and warns if shape mismatch exists."""
+
+    shape_target = (16, 16, 16)
+    b0_shape = (8, 8, 8)
+
+    b0_map = np.ones(b0_shape, dtype=np.float32)
+    kspace = np.zeros((10, 3), dtype=np.float32)
+    smaps = np.ones((1, *shape_target), dtype=np.complex64)
+    readout_time = np.ones(10, dtype=np.float32)
+
+    nufft = get_operator("finufft")(
+        samples=kspace,
+        shape=shape_target,
+        n_coils=1,
+        smaps=smaps,
+        density=False,
+    )
+
+    with pytest.warns(UserWarning):
+        op = MRIFourierCorrected(
+            nufft,
+            b0_map=b0_map,
+            readout_time=readout_time,
+        )
+
+        # check that no exception is raised and internal shape matches
+        assert op.B.shape[1] == len(readout_time)
+        assert op.shape == shape_target