Added PNS calculations

This is a direct translation of the MATLAB Pulseq calcPNS function to python, renamed to calculate_pns. This includes a partial translation of the MATLAB safe_pns_prediction repository (https://github.com/filip-szczepankiewicz/safe_pns_prediction), put into utils/safe_pns_prediction.py. To test the functionality, this adds the pnsTest sequence "example" from MATLAB Pulseq.

pypulseq/Sequence/calc_pns.py

@@ -0,0 +1,99 @@
+from types import SimpleNamespace
+from typing import Tuple
+import matplotlib.pyplot as plt
+import pypulseq as pp
+import numpy as np
+
+from pypulseq import Sequence
+from pypulseq.utils.safe_pns_prediction import safe_gwf_to_pns, safe_plot
+
+from pypulseq.utils.siemens.readasc import readasc
+from pypulseq.utils.siemens.asc_to_hw import asc_to_hw
+
+
+def calc_pns(
+        obj : Sequence, hardware : SimpleNamespace, do_plots: bool = True
+        ) -> Tuple[bool, np.array, np.ndarray, np.array]:
+    """
+    Calculate PNS using safe model implementation by Szczepankiewicz and Witzel
+    See http://github.com/filip-szczepankiewicz/safe_pns_prediction
+    
+    Returns pns levels due to respective axes (normalized to 1 and not to 100#)
+    
+    Parameters
+    ----------
+    hardware : SimpleNamespace
+        Hardware specifications. See safe_example_hw() from
+        the safe_pns_prediction package. Alternatively a text file
+        in the .asc format (Siemens) can be passed, e.g. for Prisma
+        it is MP_GPA_K2309_2250V_951A_AS82.asc (we leave it as an
+        exercise to the interested user to find were these files
+        can be acquired from)
+    do_plots : bool, optional
+        Plot the results from the PNS calculations. The default is True.
+
+    Returns
+    -------
+    ok : bool
+        Boolean flag indicating whether peak PNS is within acceptable limits
+    pns_norm : numpy.array [N]
+        PNS norm over all gradient channels, normalized to 1
+    pns_components : numpy.array [Nx3]
+        PNS levels per gradient channel
+    t_pns : np.array [N]
+        Time axis for the pns_norm and pns_components arrays
+    """
+
+    # acquire the entire gradient wave form
+    gw = obj.waveforms_and_times()[0]
+    if do_plots:
+        plt.figure()
+        plt.plot(gw[0][0], gw[0][1], gw[1][0], gw[1][1], gw[2][0], gw[2][1]) # plot the entire gradient shape
+        plt.title('gradient wave form, in T/m')
+
+    # find beginning and end times and resample GWs to a regular sampling raster
+    tf = []
+    tl = []
+    for i in range(3):
+        if gw[i].shape[1] > 0:
+            tf.append(gw[i][0,0])
+            tl.append(gw[i][0,-1])
+
+    nt_min = np.floor(min(tf) / obj.grad_raster_time + pp.eps) 
+    nt_max = np.ceil(max(tl) / obj.grad_raster_time - pp.eps)
+
+    # shift raster positions to the centers of the raster periods
+    nt_min = nt_min + 0.5
+    nt_max = nt_max - 0.5
+    if nt_min < 0.5:
+        nt_min = 0.5
+
+    t_axis = (np.arange(0,np.floor(nt_max-nt_min) + 1) + nt_min) * obj.grad_raster_time
+
+    gwr = np.zeros((t_axis.shape[0],3))
+    for i in range(3):
+        if gw[i].shape[1] > 0:
+            gwr[:,i] = np.interp(t_axis, gw[i][0], gw[i][1])
+
+    if type(hardware) == str:
+        # this loads the parameters from the provided text file
+        asc, _ = readasc(hardware)
+        hardware = asc_to_hw(asc)
+
+    # use the Szczepankiewicz' and Witzel's implementation
+    [pns_comp,res] = safe_gwf_to_pns(gwr/obj.system.gamma, np.nan*np.ones(t_axis.shape[0]), obj.grad_raster_time, hardware) # the RF vector is unused in the code inside but it is zeropaded and exported ... 
+
+    # use the exported RF vector to detect and undo zero-padding
+    pns_comp = 0.01 * pns_comp[~np.isfinite(res.rf[1:]),:]
+
+    # calc pns_norm and the final ok/not_ok
+    pns_norm = np.sqrt((pns_comp**2).sum(axis=1))
+    ok = all(pns_norm<1)
+
+    # ready
+    if do_plots:
+        # plot results
+        plt.figure()
+        safe_plot(pns_comp*100, obj.grad_raster_time)
+
+    return ok, pns_norm, pns_comp, t_axis

pypulseq/Sequence/sequence.py

@@ -17,6 +17,7 @@
 from pypulseq.Sequence.ext_test_report import ext_test_report
 from pypulseq.Sequence.read_seq import read
 from pypulseq.Sequence.write_seq import write as write_seq
+from pypulseq.Sequence.calc_pns import calc_pns
 from pypulseq.calc_rf_center import calc_rf_center
 from pypulseq.check_timing import check_timing as ext_check_timing
 from pypulseq.decompress_shape import decompress_shape
@@ -570,6 +571,40 @@ def fnint(arr_poly):
         k_traj_adc = k_traj[:, i_adc]
 
         return k_traj_adc, k_traj, t_excitation, t_refocusing, t_adc
+
+    def calculate_pns(
+            self, hardware: SimpleNamespace, do_plots: bool = True
+            ) -> Tuple[bool, np.array, np.ndarray, np.array]:
+        """
+        Calculate PNS using safe model implementation by Szczepankiewicz and Witzel
+        See http://github.com/filip-szczepankiewicz/safe_pns_prediction
+        
+        Returns pns levels due to respective axes (normalized to 1 and not to 100#)
+        
+        Parameters
+        ----------
+        hardware : SimpleNamespace
+            Hardware specifications. See safe_example_hw() from
+            the safe_pns_prediction package. Alternatively a text file
+            in the .asc format (Siemens) can be passed, e.g. for Prisma
+            it is MP_GPA_K2309_2250V_951A_AS82.asc (we leave it as an
+            exercise to the interested user to find were these files
+            can be acquired from)
+        do_plots : bool, optional
+            Plot the results from the PNS calculations. The default is True.
+
+        Returns
+        -------
+        ok : bool
+            Boolean flag indicating whether peak PNS is within acceptable limits
+        pns_norm : numpy.array [N]
+            PNS norm over all gradient channels, normalized to 1
+        pns_components : numpy.array [Nx3]
+            PNS levels per gradient channel
+        t_pns : np.array [N]
+            Time axis for the pns_norm and pns_components arrays
+        """
+        return calc_pns(self, hardware, do_plots=do_plots)
 
     def check_timing(self) -> Tuple[bool, List[str]]:
         """

pypulseq/seq_examples/scripts/pns_test.py

@@ -0,0 +1,111 @@
+import pypulseq as pp
+from pypulseq.utils.safe_pns_prediction import safe_example_hw
+import numpy as np
+
+from copy import copy
+
+# Set system limits
+sys = pp.Opts(max_grad=32, grad_unit='mT/m',
+              max_slew=130, slew_unit='T/m/s',
+              rf_ringdown_time=20e-6,
+              rf_dead_time=100e-6,
+              adc_dead_time=20e-6,
+              B0=2.89)
+
+seq = pp.Sequence()      # Create a new sequence object
+
+## prepare test objects
+# pns is induced by the ramps, so we use long gradients to isolate the
+# effects of the ramps
+gpt = 10e-3
+delay = 30e-3
+rt_min = sys.grad_raster_time
+rt_test = np.floor(sys.max_grad/sys.max_slew/sys.grad_raster_time)*sys.grad_raster_time
+ga_min = sys.max_slew*rt_min
+ga_test = sys.max_slew*rt_test
+
+gx_min = pp.make_trapezoid(channel='x',system=sys,amplitude=ga_min,rise_time=rt_min,fall_time=2*rt_min,flat_time=gpt)
+gy_min = copy(gx_min)
+gy_min.channel = 'y'
+gz_min = copy(gx_min)
+gz_min.channel = 'z'
+
+gx_test = pp.make_trapezoid(channel='x',system=sys,amplitude=ga_test,rise_time=rt_test,fall_time=2*rt_test,flat_time=gpt)
+gy_test = copy(gx_test)
+gy_test.channel = 'y'
+gz_test = copy(gx_test)
+gz_test.channel = 'z'
+
+g_min = [gx_min,gy_min,gz_min]
+g_test = [gx_test,gy_test,gz_test]
+
+# dummy FID sequence
+# Create non-selective pulse 
+rf = pp.make_block_pulse(np.pi/2,duration=0.1e-3, system=sys)
+# Define delays and ADC events
+adc = pp.make_adc(512,duration=6.4e-3, system=sys)
+
+
+## Define sequence blocks
+seq.add_block(pp.make_delay(delay))
+for a in range(3):
+    seq.add_block(g_min[a])
+    seq.add_block(pp.make_delay(delay))
+    seq.add_block(g_test[a])
+    seq.add_block(pp.make_delay(delay))
+    for b in range(a+1, 3):
+        seq.add_block(g_min[a],g_min[b])
+        seq.add_block(pp.make_delay(delay))
+        seq.add_block(g_test[a],g_test[b])
+        seq.add_block(pp.make_delay(delay))
+
+seq.add_block(*g_min)
+seq.add_block(pp.make_delay(delay))
+seq.add_block(*g_test)
+seq.add_block(pp.make_delay(delay))
+seq.add_block(*g_min)
+seq.add_block(rf)
+seq.add_block(adc)
+
+## check whether the timing of the sequence is correct
+ok, error_report = seq.check_timing()
+
+if (ok):
+    print('Timing check passed successfully')
+else:
+    print('Timing check failed! Error listing follows:')
+    print(error_report)
+    print('\n')
+
+
+## do some visualizations
+seq.plot()             # Plot all sequence waveforms
+
+## 'install' to the IDEA simulator
+# seq.write('idea/external.seq')
+
+## PNS calc
+pns_ok, pns_n, pns_c, tpns = seq.calculate_pns(safe_example_hw(), do_plots=True) # Safe example HW
+
+# pns_ok, pns_n, pns_c, tpns = seq.calculate_pns('idea/asc/MP_GPA_K2309_2250V_951A_AS82.asc', do_plots=True) # prisma
+# pns_ok, pns_n, pns_c, tpns = seq.calculate_pns('idea/asc/MP_GPA_K2309_2250V_951A_GC98SQ.asc', do_plots=True) # aera-xq
+
+# ## load simulation results 
+
+# #[sll,~,~,vscale]=dsv_read('idea/dsv/prisma_pulseq_SLL.dsv')
+# #[sll,~,~,vscale]=dsv_read('idea/dsv/aera_pulseq_SLL.dsv')
+# [sll,~,~,vscale]=dsv_read('idea/dsv/terra_pulseq_SLL.dsv')
+# sll=cumsum(sll/vscale)
+
+# ## plot
+# figureplot(sll(104:end)) # why 104? good question
+# hold on
+# plot(tpns*1e5-0.5,pns_n)
+# title('comparing internal and IDEA PNS predictions')
+
+# ## manual time alignment to calculate relative differences
+# ssl_s=104+tpns(1)*1e5-1.5
+# ssl_e=ssl_s+length(pns_n)-1
+# #figureplot(sll(ssl_s:ssl_e))hold on plot(pns_n)
+# figureplot((sll(ssl_s:ssl_e)-pns_n)./pns_n*100)
+# title('relative difference in #')