ready for ismrm

Author: mcencini

examples/scripts/01-simple-simulator.py

@@ -0,0 +1,67 @@
+"""SSFP MR Fingerprinting simulator"""
+
+def mrf(flip, TR, T1, T2, diff=None, device="cpu"):
+    """
+    Simulate an inversion-prepared SSFP sequence with variable flip angles.
+
+    Args:
+        flip (array-like): Flip angle in [deg] of shape (npulses,).
+        TR (float): Repetition time in [ms].
+        T1 (float, array-like): Longitudinal relaxation time in [ms].
+        T2 (float, array-like): Transverse relaxation time in [ms].
+        diff (optional, str, tuple[str]): Arguments to get the signal derivative with respect to.
+            Defaults to None (no differentation).
+        device (optional, str): Computational device. Defaults to "cpu".
+    """
+    # initialize simulator
+    simulator = SSFPMRF(T1=T1, T2=T2, device=device, diff=diff)
+            
+    # run simulator
+    if diff:
+        # actual simulation
+        sig, dsig = simulator(flip=flip, TR=TR)
+        return sig.detach().cpu().numpy(), dsig.detach().cpu().numpy()
+
+    else:
+        # actual simulation
+        sig = simulator(flip=flip, TR=TR)
+        return sig.cpu().numpy()
+
+# Simulator
+from epgtorchx import base
+from epgtorchx import ops
+
+class SSFPMRF(base.BaseSimulator):
+    """Class to simulate inversion-prepared (variable flip angle) SSFP."""
+
+    @staticmethod
+    def sequence(flip, TR, T1, T2, states, signal):
+        
+        # get number of frames and echoes
+        device = flip.device
+        npulses = flip.shape[0]
+
+        # define operators
+        # preparation
+        InvPulse = ops.RFPulse(device, alpha=180.0)
+        Crusher = ops.Spoil()
+        Prep = ops.CompositeOperator(Crusher, InvPulse)
+
+        # readout
+        RF = ops.RFPulse(device) # excitation
+        E = ops.Relaxation(device, TR, T1, T2) # relaxation until TR
+        S = ops.Shift() # gradient spoil
+    
+        # actual sequence loop
+        for n in range(npulses):
+            # apply pulse
+            states = RF(states, flip[n])
+
+            # record signal
+            signal[n] = ops.observe(states)
+
+            # relax, recover and spoil
+            states = E(states)
+            states = S(states)
+
+        return signal

examples/scripts/02-derivatives.py

@@ -104,18 +104,53 @@ def crlb_finitediff_cost(flip, ESP, T1, T2):
 tcost = t1 - t0
 
 # plot derivative
+#%%
+fsz = 20
 plt.figure()
-plt.rcParams.update({'font.size': 22})
+plt.subplot(2,2,1)
+plt.rcParams.update({'font.size': 0.5 * fsz})
+plt.plot(angles, '.')
+plt.xlabel("Echo #", fontsize=fsz)
+plt.xlim([-1, 51])
+plt.ylabel("Flip Angle [deg]", fontsize=fsz)
+
+plt.subplot(2,2,2)
+plt.rcParams.update({'font.size': 0.5 * fsz})
 plt.plot(abs(grad), '-k'), plt.plot(abs(grad0), '*r')
-plt.xlabel("Echo #", fontsize=40)
+plt.xlabel("Echo #", fontsize=fsz)
 plt.xlim([-1, 51])
-plt.ylabel(r"$\frac{\partial signal}{\partial T2}$ [a.u.]", fontsize=40)
+plt.ylabel(r"$\frac{\partial signal}{\partial T2}$ [a.u.]", fontsize=fsz)
 plt.legend(["Finite Diff", "Auto Diff"])
 
-plt.figure()
-plt.rcParams.update({'font.size': 22})
+
+plt.subplot(2,2,3)
+plt.rcParams.update({'font.size': 0.5 * fsz})
 plt.plot(abs(dcost), '-k'), plt.plot(abs(dcost0), '*r')
-plt.xlabel("Echo #", fontsize=40)
+plt.xlabel("Echo #", fontsize=fsz)
 plt.xlim([-1, 51])
-plt.ylabel(r"$\frac{\partial CRLB}{\partial FA}$ [a.u.]", fontsize=40)
-plt.legend(["Finite Diff", "Auto Diff"])
+plt.ylabel(r"$\frac{\partial CRLB}{\partial FA}$ [a.u.]", fontsize=fsz)
+plt.legend(["Finite Diff", "Auto Diff"])
+
+plt.subplot(2,2,4)
+
+# define labels
+# plot results
+labels = ['derivative of signal', 'CRLB objective gradient']
+time_finite = [round(tgrad0, 2), round(tcost0, 2)]
+time_auto = [round(tgrad, 2), round(tcost, 2)]
+
+
+x = np.arange(len(labels))  # the label locations
+width = 0.35  # the width of the bars
+rects1 = plt.bar(x + width/2, time_finite, width, label='Finite Diff')
+rects2 = plt.bar(x - width/2, time_auto, width, label='Auto Diff')
+
+# Add some text for labels, title and custom x-axis tick labels, etc.
+plt.ylabel('Execution Time [s]', fontsize=fsz)
+plt.xticks(x, labels, fontsize=fsz)
+plt.ylim([0, 25])
+plt.legend()
+
+plt.bar_label(rects1, padding=3, fontsize=fsz)
+plt.bar_label(rects2, padding=3, fontsize=fsz)
+

src/epgtorchx/__init__.py

@@ -23,6 +23,8 @@
 from . import regression
 
 from .bloch import *  # noqa
+from .bloch import base
+from .bloch import ops
 from .phantoms import *  # noqa
 
 __all__ = []

src/epgtorchx/bloch/__init__.py

@@ -10,6 +10,7 @@
 from . import model as _model
 from . import ops
 from . import blocks
+from .model import base
 
 from .model import *  # noqa
 

src/epgtorchx/bloch/model/base.py

@@ -160,7 +160,6 @@ def simulate_ssfp(flip, TR, T1, T2, diff=None, device="cpu"):
         B1phase (optional, Union[float, npt.NDArray[float], torch.FloatTensor]): B1 relative phase in [deg]. (0.0 := nominal rf phase). Defaults to None.
 
     """
-
     # main properties
     T1: Union[float, npt.NDArray[float], torch.FloatTensor]  # ms
     T2: Union[float, npt.NDArray[float], torch.FloatTensor]  # ms