Several bug fixes in Sequence, add_gradients, make_trapezoid, a…

…nd `rotate`:

- Bypass code for `restoreAdditionalShapeSamples` did not include `grad.first` and `grad.last`
- `make_arbitrary_grad` was not available as `pp.make_arbitrary_grad`
- `add_gradients`' detection of arbitrary gradients was incorrect for Python's 0-based indexing
- `add_gradients` handling of non-unique times was broken
- `make_trapezoid` incorrectly checked and reported the minimum duration for gradients specifying area and duration (also fixed in local upstream MATLAB pulseq)
- `rotate` did not have a system parameter
- `rotate` used the incorrect channel for gradient on the second rotation axis

- Small change in `add_gradients` to accept a single gradient (i.e. for cases where a variable number of gradients are added), which also allowed a minor cleanup in `rotate`

pypulseq/Sequence/sequence.py

@@ -1392,13 +1392,13 @@ def waveforms_and_times(
                             # TODO: Implement restoreAdditionalShapeSamples
                             #       https://github.com/pulseq/pulseq/blob/master/matlab/%2Bmr/restoreAdditionalShapeSamples.m
 
-                            out_len[j] += len(grad.tt)
-                            shape_pieces[j].append(np.array(
+                            out_len[j] += len(grad.tt)+2
+                            shape_pieces[j, block_counter] = np.array(
                                 [
-                                    curr_dur + grad.delay + grad.tt,
-                                    grad.waveform,
+                                    curr_dur + grad.delay + np.concatenate(([0], grad.tt, [grad.tt[-1] + self.grad_raster_time/2])),
+                                    np.concatenate(([grad.first], grad.waveform, [grad.last]))
                                 ]
-                            ))
+                            )
                         else:  # Extended trapezoid
                             out_len[j] += len(grad.tt)
                             shape_pieces[j].append(np.array(
@@ -1509,7 +1509,6 @@ def waveforms_and_times(
 
             wave_data.append(np.concatenate(shape_pieces[j], axis=1))
 
-        for j in range(shape_channels):
             rftdiff = np.diff(wave_data[j][0])
             if np.any(rftdiff < eps):
                 raise Warning(

pypulseq/__init__.py

@@ -30,6 +30,7 @@ def round_half_up(n, decimals=0):
 from pypulseq.calc_rf_center import calc_rf_center
 from pypulseq.make_adc import make_adc
 from pypulseq.make_adiabatic_pulse import make_adiabatic_pulse
+from pypulseq.make_arbitrary_grad import make_arbitrary_grad
 from pypulseq.make_arbitrary_rf import make_arbitrary_rf
 from pypulseq.make_block_pulse import make_block_pulse
 from pypulseq.make_sigpy_pulse import *

pypulseq/add_gradients.py

@@ -45,9 +45,11 @@ def add_gradients(
     if max_slew <= 0:
         max_slew = system.max_slew
 
-    if len(grads) < 2:
-        raise ValueError("Cannot add less than two gradients")
-
+    if len(grads) == 0:
+        raise ValueError("No gradients specified")
+    if len(grads) == 1:
+        return grads[0]
+
     # First gradient defines channel
     channel = grads[0].channel
 
@@ -65,7 +67,7 @@ def add_gradients(
         if is_trap[-1]:
             is_arb.append(False)
         else:
-            tt_rast = grads[ii].tt / system.grad_raster_time + 0.5
+            tt_rast = grads[ii].tt / system.grad_raster_time - 0.5
             is_arb.append(np.all(np.abs(tt_rast - np.arange(len(tt_rast)))) < eps)
 
     # Convert to numpy.ndarray for fancy-indexing later on
@@ -103,13 +105,13 @@ def add_gradients(
 
         times = np.sort(np.unique(times))
         dt = times[1:] - times[:-1]
-        ieps = dt < eps
+        ieps = np.flatnonzero(dt < eps)
         if np.any(ieps):
-            dtx = [times[0], *dt]
+            dtx = np.array([times[0], *dt])
             dtx[ieps] = (
                 dtx[ieps] + dtx[ieps + 1]
             )  # Assumes that no more than two too similar values can occur
-            dtx[ieps + 1] = []
+            dtx = np.delete(dtx, ieps + 1)
             times = np.cumsum(dtx)
 
         amplitudes = np.zeros_like(times)

pypulseq/make_trapezoid.py

@@ -4,6 +4,35 @@
 
 from pypulseq.opts import Opts
 
+def calculate_shortest_params_for_area(area, max_slew, max_grad, grad_raster_time):
+    rise_time = (
+        np.ceil(np.sqrt(np.abs(area) / max_slew) / grad_raster_time)
+        * grad_raster_time
+    )
+    if rise_time < grad_raster_time:  # Area was almost 0 maybe
+        rise_time = grad_raster_time
+    amplitude = np.divide(area, rise_time)  # To handle nan
+    t_eff = rise_time
+
+    if np.abs(amplitude) > max_grad:
+        t_eff = (
+            np.ceil(np.abs(area) / max_grad / grad_raster_time)
+            * grad_raster_time
+        )
+        amplitude = area / t_eff
+        rise_time = (
+            np.ceil(np.abs(amplitude) / max_slew / grad_raster_time)
+            * grad_raster_time
+        )
+
+        if rise_time == 0:
+            rise_time = grad_raster_time
+
+    flat_time = t_eff - rise_time
+    fall_time = rise_time
+
+    return amplitude, rise_time, flat_time, fall_time
+
 
 def make_trapezoid(
     channel: str,
@@ -116,12 +145,14 @@ def make_trapezoid(
     elif duration > 0:
         if amplitude == 0:
             if rise_time == 0:
-                dC = 1 / np.abs(2 * max_slew) + 1 / np.abs(2 * max_slew)
-                possible = duration**2 > 4 * np.abs(area) * dC
-                assert possible, (
+                _, rise_time, flat_time, fall_time = calculate_shortest_params_for_area(area, max_slew, max_grad, system.grad_raster_time)
+                min_duration = rise_time + flat_time + fall_time
+                assert duration >= min_duration, (
                     f"Requested area is too large for this gradient. Minimum required duration is "
-                    f"{np.round(np.sqrt(4 * np.abs(area) * dC) * 1e6)} uss"
+                    f"{np.round(min_duration * 1e6)} uss"
                 )
+
+                dC = 1 / np.abs(2 * max_slew) + 1 / np.abs(2 * max_slew)
                 amplitude2 = (
                     duration - np.sqrt(duration**2 - 4 * np.abs(area) * dC)
                 ) / (2 * dC)
@@ -156,33 +187,9 @@ def make_trapezoid(
         if area == None:
             raise ValueError("Must supply area or duration.")
         else:
-            # Find the shortest possible duration. First check if the area can be realized as a triangle.
-            # If not, then it must be a trapezoid.
-            rise_time = (
-                np.ceil(np.sqrt(np.abs(area) / max_slew) / system.grad_raster_time)
-                * system.grad_raster_time
-            )
-            if rise_time < system.grad_raster_time:  # Area was almost 0 maybe
-                rise_time = system.grad_raster_time
-            amplitude2 = np.divide(area, rise_time)  # To handle nan
-            t_eff = rise_time
+            # Find the shortest possible duration.           
+            amplitude2, rise_time, flat_time, fall_time = calculate_shortest_params_for_area(area, max_slew, max_grad, system.grad_raster_time)
 
-            if np.abs(amplitude2) > max_grad:
-                t_eff = (
-                    np.ceil(np.abs(area) / max_grad / system.grad_raster_time)
-                    * system.grad_raster_time
-                )
-                amplitude2 = area / t_eff
-                rise_time = (
-                    np.ceil(np.abs(amplitude2) / max_slew / system.grad_raster_time)
-                    * system.grad_raster_time
-                )
-
-                if rise_time == 0:
-                    rise_time = system.grad_raster_time
-
-            flat_time = t_eff - rise_time
-            fall_time = rise_time
 
     if np.abs(amplitude2) > max_grad:
         raise ValueError("Amplitude violation.")

pypulseq/rotate.py

@@ -5,6 +5,7 @@
 
 from pypulseq.add_gradients import add_gradients
 from pypulseq.scale_grad import scale_grad
+from pypulseq.opts import Opts
 
 
 def __get_grad_abs_mag(grad: SimpleNamespace) -> np.ndarray:
@@ -17,6 +18,7 @@ def rotate(
     *args: SimpleNamespace,
     angle: float,
     axis: str,
+    system=Opts()
 ) -> List[SimpleNamespace]:
     """
     Rotates the corresponding gradient(s) about the given axis by the specified amount. Gradients parallel to the
@@ -85,7 +87,7 @@ def rotate(
         max_mag = np.max((max_mag, __get_grad_abs_mag(g)))
         rotated2.append(scale_grad(grad=g, scale=np.cos(angle)))
         g = scale_grad(grad=g, scale=-np.sin(angle))
-        g.channel = axes_to_rotate[1]
+        g.channel = axes_to_rotate[0]
         rotated1.append(g)
 
     # Eliminate zero-amplitude gradients
@@ -99,17 +101,11 @@ def rotate(
 
     # Add gradients on the corresponding axis together
     g = []
-    if len(rotated1) > 1:
-        g.append(add_gradients(grads=rotated1))
-    else:
-        if len(rotated1) != 0:
-            g.append(rotated1[0])
-
-    if len(rotated2) > 1:
-        g.append(add_gradients(grads=rotated2))
-    else:
-        if len(rotated2) != 0:
-            g.append(rotated2[0])
+    if len(rotated1) != 0:
+        g.append(add_gradients(grads=rotated1, system=system))
+
+    if len(rotated2) != 0:
+        g.append(add_gradients(grads=rotated2, system=system))
 
     # Eliminate zero amplitude gradients
     for i in range(len(g) - 1, -1, -1):