SLR isodelay

docs/mri_rf.rst

@@ -89,6 +89,7 @@ SLR Pulse Design Functions
     :nosignatures:
 
     sigpy.mri.rf.slr.dzrf
+    sigpy.mri.rf.slr.get_isodelay
     sigpy.mri.rf.slr.root_flip
     sigpy.mri.rf.slr.dz_gslider_rf
     sigpy.mri.rf.slr.dz_gslider_b

sigpy/mri/rf/slr.py

@@ -12,7 +12,8 @@
 
 __all__ = ['dzrf', 'dzls', 'msinc', 'dzmp', 'fmp', 'dzlp',
            'b2rf', 'b2a', 'mag2mp', 'ab2rf', 'dz_gslider_b', 'dz_gslider_rf',
-           'root_flip', 'dz_recursive_rf', 'dz_hadamard_b', 'calc_ripples']
+           'root_flip', 'dz_recursive_rf', 'dz_hadamard_b', 'calc_ripples',
+           'get_isodelay']
 
 """ Functions for SLR pulse design
     SLR algorithm simplifies the solution of the Bloch equations
@@ -513,6 +514,27 @@ def ab2rf(a, b):
     return rf
 
 
+def get_isodelay(rf, dt):
+    r"""Function to compute isodelay of a selective RF pulse, for calulation of
+    refocusing gradient area. Approximates isodelay as the time interval
+    between peak RF energy and end of pulse, neglecting the (small) nonlinear
+    dependence on flip angle.
+
+    Args:
+        rf (array): input RF pulse.
+        dt (float): hardware dwell time (s).
+
+    Returns:
+        RF pulse isodelay. Always positive. (s).
+
+    References:
+        Bernstein, M.A. King, K.F. and Zhou, X.J. (2004). Handbook of MRI pulse
+        sequences. Amsterdam: Academic Press.
+    """
+
+    return (np.size(rf) - np.argmax(abs(rf)))*dt
+
+
 def root_flip(b, d1, flip, tb, verbose=False):
     r"""Exhaustive root-flip pattern search for min-peak b1
 

tests/mri/rf/test_ptx.py

@@ -2,7 +2,7 @@
 import numpy as np
 import sigpy as sp
 import numpy.testing as npt
-import scipy.ndimage.filters as filt
+from scipy.ndimage import gaussian_filter
 
 from sigpy.mri import rf, linop, sim
 
@@ -23,7 +23,7 @@ def problem_2d(dim):
         circle = x * x + y * y <= int(img_shape[0] / 6) ** 2
         target = np.zeros(img_shape)
         target[circle] = 1
-        target = filt.gaussian_filter(target, 1)
+        target = gaussian_filter(target, 1)
         target = target.astype(np.complex64)
 
         sens = sim.birdcage_maps(sens_shape)
@@ -43,7 +43,7 @@ def problem_3d(dim, Nz):
         circle = x * x + y * y + z * z <= int(img_shape[0] / 5) ** 2
         target = np.zeros(img_shape)
         target[circle] = 1
-        target = filt.gaussian_filter(target, 1)
+        target = gaussian_filter(target, 1)
         target = target.astype(np.complex64)
         sens = sp.mri.sim.birdcage_maps(sens_shape)
 

tests/mri/rf/test_slr.py

@@ -17,7 +17,7 @@ def test_st(self):
         N = 128
         tb = 16
         filts = ['ls', 'ms', 'pm', 'min', 'max']
-        for idx, filt in enumerate(filts):
+        for filt in filts:
             pulse = sp.mri.rf.dzrf(N, tb, ptype='st', ftype=filt,
                                    d1=0.01, d2=0.01)
 
@@ -36,8 +36,10 @@ def test_inv(self):
         ptype = 'ex'
         filts = ['min', 'max']  # filts produce inconsistent inversions
 
-        for idx, filt in enumerate(filts):
+        isodelays = []
+        for filt in filts:
             pulse = rf.slr.dzrf(N, tb, ptype, filt, d1, d2)
+            isodelays.append(rf.slr.get_isodelay(pulse, dt=4e-6))
 
             [_, b] = rf.sim.abrm(pulse, np.arange(-2 * tb, 2 * tb, 0.01))
             mz = 1 - 2 * np.abs(b) ** 2
@@ -47,6 +49,8 @@ def test_inv(self):
                             mz[int(len(mz) / 2 + len(mz)/3)]])
 
             npt.assert_almost_equal(pts, np.array([1, -0.2, 1]), decimal=1)
+        # test that isodelay is shorter for minphase inv than maxphase inv
+        npt.assert_array_less(isodelays[0], isodelays[1])
 
     def test_root_flipped(self):
         tb = 12