From c251d7799352f39dc671f9bc249f0d00e181bd1b Mon Sep 17 00:00:00 2001 From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com> Date: Mon, 11 Sep 2023 13:50:52 +0200 Subject: [PATCH] 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 | 99 ++++++ pypulseq/Sequence/sequence.py | 35 ++ pypulseq/seq_examples/scripts/pns_test.py | 111 ++++++ pypulseq/utils/safe_pns_prediction.py | 411 ++++++++++++++++++++++ pypulseq/utils/siemens/asc_to_hw.py | 72 ++++ pypulseq/utils/siemens/readasc.py | 97 +++++ 6 files changed, 825 insertions(+) create mode 100644 pypulseq/Sequence/calc_pns.py create mode 100644 pypulseq/seq_examples/scripts/pns_test.py create mode 100644 pypulseq/utils/safe_pns_prediction.py create mode 100644 pypulseq/utils/siemens/asc_to_hw.py create mode 100644 pypulseq/utils/siemens/readasc.py diff --git a/pypulseq/Sequence/calc_pns.py b/pypulseq/Sequence/calc_pns.py new file mode 100644 index 00000000..6ab848bc --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py index 23f64d3f..a509f637 100755 --- a/pypulseq/Sequence/sequence.py +++ b/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]]: """ diff --git a/pypulseq/seq_examples/scripts/pns_test.py b/pypulseq/seq_examples/scripts/pns_test.py new file mode 100644 index 00000000..bd07509f --- /dev/null +++ b/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 #') \ No newline at end of file diff --git a/pypulseq/utils/safe_pns_prediction.py b/pypulseq/utils/safe_pns_prediction.py new file mode 100644 index 00000000..0883b5e6 --- /dev/null +++ b/pypulseq/utils/safe_pns_prediction.py @@ -0,0 +1,411 @@ +# This code is a direct Python translation of the relevant functions in +# https://github.com/filip-szczepankiewicz/safe_pns_prediction/ to perform +# PNS calculations with pypulseq +# +# A small modification was made to safe_plot to plot long sequences better + + +# BSD 3-Clause License + +# Copyright (c) 2018, Filip Szczepankiewicz and Thomas Witzel +# All rights reserved. + +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: + +# 1. Redistributions of source code must retain the above copyright notice, this +# list of conditions and the following disclaimer. + +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. + +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. + +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +from types import SimpleNamespace + +import numpy as np +import matplotlib.pyplot as plt + + +def safe_example_hw(): + # function hw = safe_example_hw() + # + # SAFE model parameters for EXAMPLE scanner hardware (not a real scanner). + # See comments for units. + + hw = SimpleNamespace() + hw.name = 'MP_GPA_EXAMPLE' + hw.checksum = '1234567890' + hw.dependency = '' + + hw.x = SimpleNamespace() + hw.x.tau1 = 0.20 # ms + hw.x.tau2 = 0.03 # ms + hw.x.tau3 = 3.00 # ms + hw.x.a1 = 0.40 + hw.x.a2 = 0.10 + hw.x.a3 = 0.50 + hw.x.stim_limit = 30.0 # T/m/s + hw.x.stim_thresh = 24.0 # T/m/s + hw.x.g_scale = 0.35 # 1 + + hw.y = SimpleNamespace() + hw.y.tau1 = 1.50 # ms + hw.y.tau2 = 2.50 # ms + hw.y.tau3 = 0.15 # ms + hw.y.a1 = 0.55 + hw.y.a2 = 0.15 + hw.y.a3 = 0.30 + hw.y.stim_limit = 15.0 # T/m/s + hw.y.stim_thresh = 12.0 # T/m/s + hw.y.g_scale = 0.31 # 1 + + hw.z = SimpleNamespace() + hw.z.tau1 = 2.00 # ms + hw.z.tau2 = 0.12 # ms + hw.z.tau3 = 1.00 # ms + hw.z.a1 = 0.42 + hw.z.a2 = 0.40 + hw.z.a3 = 0.18 + hw.z.stim_limit = 25.0 # T/m/s + hw.z.stim_thresh = 20.0 # T/m/s + hw.z.g_scale = 0.25 # 1 + return hw + + +def safe_example_gwf(): + # function function [gwf, rf, dt] = safe_example_gwf() + # Waveform with some frequency matching by Filip Szczepankiewicz. + # + # Waveform was optimized in the NOW framework by Jens Sjölund et al. + # https://github.com/jsjol/NOW + # + # Optimization was Maxwell-compensated to remove effects of concomitant + # gradients. + # https://arxiv.org/ftp/arxiv/papers/1903/1903.03357.pdf + + ## STE + dt = 1e-3 # ms + + # T/m + gwf = 0.08 * np.array([ + [0, 0, 0], + [-0.2005, 0.9334, 0.3029], + [-0.2050, 0.9324, 0.3031], + [-0.2146, 0.9302, 0.3032], + [-0.2313, 0.9263, 0.3030], + [-0.2589, 0.9193, 0.3019], + [-0.3059, 0.9060, 0.2980], + [-0.3892, 0.8767, 0.2883], + [-0.3850, 0.7147, 0.3234], + [-0.3687, 0.5255, 0.3653], + [-0.3509, 0.3241, 0.4070], + [-0.3323, 0.1166, 0.4457], + [-0.3136, -0.0906, 0.4783], + [-0.2956, -0.2913, 0.5019], + [-0.2790, -0.4793, 0.5139], + [-0.2642, -0.6491, 0.5118], + [-0.2518, -0.7957, 0.4939], + [-0.2350, -0.8722, 0.4329], + [-0.2187, -0.9111, 0.3541], + [-0.2063, -0.9409, 0.2747], + [-0.1977, -0.9627, 0.1933], + [-0.1938, -0.9768, 0.1080], + [-0.1967, -0.9820, 0.0159], + [-0.2114, -0.9751, -0.0883], + [-0.2292, -0.9219, -0.2150], + [-0.2299, -0.8091, -0.3561], + [-0.2290, -0.6748, -0.5011], + [-0.2253, -0.5239, -0.6460], + [-0.2178, -0.3620, -0.7868], + [-0.2056, -0.1948, -0.9194], + [-0.1391, -0.0473, -0.9908], + [-0.0476, 0.0607, -0.9987], + [ 0.0215, 0.1452, -0.9909], + [ 0.0725, 0.2136, -0.9759], + [ 0.1114, 0.2709, -0.9579], + [ 0.1426, 0.3204, -0.9383], + [ 0.1690, 0.3641, -0.9177], + [ 0, 0, 0], + [ 0, 0, 0], + [ 0, 0, 0], + [ 0, 0, 0], + [ 0, 0, 0], + [ 0, 0, 0], + [ 0, 0, 0], + [-0.3734, -0.1768, 0.9125], + [-0.3825, -0.2310, 0.8965], + [-0.3919, -0.2895, 0.8752], + [-0.4015, -0.3543, 0.8465], + [-0.4108, -0.4290, 0.8065], + [-0.4182, -0.5202, 0.7469], + [-0.4178, -0.6423, 0.6451], + [-0.3855, -0.8173, 0.4321], + [-0.3110, -0.9418, 0.1401], + [-0.2526, -0.9669, -0.0674], + [-0.2100, -0.9541, -0.2213], + [-0.1766, -0.9227, -0.3474], + [-0.1491, -0.8788, -0.4570], + [-0.1258, -0.8239, -0.5555], + [-0.1056, -0.7583, -0.6459], + [-0.0882, -0.6809, -0.7293], + [-0.0734, -0.5900, -0.8061], + [-0.0615, -0.4830, -0.8753], + [-0.0533, -0.3556, -0.9349], + [-0.0506, -0.2005, -0.9801], + [-0.0575, -0.0019, -1.0000], + [-0.0909, 0.2976, -0.9521], + [-0.3027, 0.9509, -0.0860], + [-0.2737, 0.9610, -0.0692], + [-0.2524, 0.9675, -0.0596], + [-0.2364, 0.9719, -0.0533], + [-0.2245, 0.9749, -0.0490], + [-0.2158, 0.9770, -0.0459], + [-0.2097, 0.9785, -0.0439], + [-0.2058, 0.9794, -0.0426], + [-0.2039, 0.9798, -0.0420], + [ 0, 0, 0] + ]) + + rf = np.ones(gwf.shape[0]) + rf[40:] = -1 + + return gwf, rf, dt + + +def safe_hw_check(hw): + # function safe_hw_check(hw) + # + # Make sure that all is well with the hardware configuration. + + if abs(hw.x.a1 + hw.x.a2 + hw.x.a3 - 1) > 0.001 or \ + abs(hw.y.a1 + hw.y.a2 + hw.y.a3 - 1) > 0.001 or \ + abs(hw.z.a1 + hw.z.a2 + hw.z.a3 - 1) > 0.001: + raise ValueError('Hardware specification a1+a2+a3 must be equal to 1!') + + axl = ['x', 'y', 'z'] + fnl = ['stim_limit', 'stim_thresh', 'tau1', 'tau2', 'tau3', 'a1', 'a2', 'a3', 'g_scale'] + + for axn in axl: + if not hasattr(hw, axn): + raise ValueError(f"'{axn}' missing in hardware specification") + + hw_ax = getattr(hw, axn) + for par in fnl: + if not hasattr(hw_ax, par): + raise ValueError(f"'{axn}.{par}' missing in hardware specification") + + +def safe_longest_time_const(hw): + # function ltau = safe_longest_time_const(hw) + # Get the longest time constant. Can be used to estimate the size of zero + # padding. + + return max([hw.x.tau1, hw.x.tau2, hw.x.tau3, + hw.y.tau1, hw.y.tau2, hw.y.tau3, + hw.z.tau1, hw.z.tau2, hw.z.tau3]) + + +def safe_pns_model(dgdt, dt, hw): + # function stim = safe_pns_model(dgdt, dt, hw) + # + # dgdt (nx3) is in T/m/s + # dt (1x1) is in s + # All time coefficients (a1 and tau1 etc.) are in ms. + # + # This PNS model is based on the SAFE-abstract + # SAFE-Model - A New Method for Predicting Peripheral Nerve Stimulations in MRI + # by Franz X. Herbank and Matthias Gebhardt. Abstract No 2007. + # Proc. Intl. Soc. Mag. Res. Med. 8, 2000, Denver, Colorado, USA + # https://cds.ismrm.org/ismrm-2000/PDF7/2007.PDF + # + # The main SAFE-model was coded by Thomas Witzel @ Martinos Center, + # MGH, HMS, Boston, MA, USA. + # + # The code was adapted/expanded/corrected by Filip Szczepankiewicz @ LMI + # BWH, HMS, Boston, MA, USA, and Lund University, Sweden. + + stim1 = hw.a1 * abs( safe_tau_lowpass(dgdt , hw.tau1, dt * 1000) ) + stim2 = hw.a2 * safe_tau_lowpass(abs(dgdt), hw.tau2, dt * 1000) + stim3 = hw.a3 * abs( safe_tau_lowpass(dgdt , hw.tau3, dt * 1000) ) + + stim = (stim1 + stim2 + stim3) / hw.stim_limit * hw.g_scale * 100 + + return stim + + # Not sure where something goes awry, probably in the lowpass filter, but + # compared to the Siemens simulator we are exactly a factor of pi off, so + # I'm dividing the final result by pi. + # Note also that the final result is essentially some kind of arbitrary + # unit. - TW + + # UPDATE 210720 - The pi factor was not quite correct. Instead, the correct + # factor was determined by the gradient scale factor (hw.g_scale, defined + # in the .asc file). Thanks to Maxim Zaitsev for supporting this buggfix and + # validating that the updated code is accurate. - FSz + + +def safe_tau_lowpass(dgdt, tau, dt, eps=1e-16): + # function fw = safe_tau_lowpass(dgdt, tau, dt) + # + # Apply a RC lowpass filter with time constant tau = RC to data with sampling + # interval dt. NOTE tau and dt need to be in the same unit (i.e. s or ms) + # The SAFE model abstract by Hebrank et.al. just says "Lowpass with time-constant tau", + # so I decided to make the most simple filter possible here. + # The RC lowpass is also appealing because its something Siemens could have + # easily implemented on their hardware stimulation monitors, so I'm probably + # pretty close. - TW + # + # UPDATE 230206 - There was a factor alpha missing on the first sample it + # has now been corrected. Thanks to Oliver Schad for finding this error. + # - FSz + + alpha = dt / (tau + dt) + + # Calculate number of elements in filter to reach desired accuracy (eps) + n = min(round(np.log(eps) / np.log(1-alpha)), dgdt.shape[0]) + filt = (1-alpha)**np.arange(n) + + # Implements lowpass filter using convolution to get rid of for loop in original code + return alpha * np.convolve(dgdt, filt)[:dgdt.shape[0]] + + +def safe_gwf_to_pns(gwf, rf, dt, hw, do_padding=True): + # function [pns, res] = safe_gwf_to_pns(gwf, rf, dt, hw, doPadding) + # + # gwf (nx3) in T/m + # dt (1x1) in s + # hw (struct) is structure that describes the hardware configuration and PNS + # response. Example: hw = safe_example_hw(). + # doPadding adds zeropadding based on the decay time. + # + # This PNS model is based on the SAFE-abstract + # SAFE-Model - A New Method for Predicting Peripheral Nerve Stimulations in MRI + # by Franz X. Herbank and Matthias Gebhardt. Abstract No 2007. + # Proc. Intl. Soc. Mag. Res. Med. 8, 2000, Denver, Colorado, USA + # https://cds.ismrm.org/ismrm-2000/PDF7/2007.PDF + # + # The main SAFE-model was coded by Thomas Witzel @ Martinos Center, + # MGH, HMS, Boston, MA, USA. + # + # The code was adapted/expanded by Filip Szczepankiewicz @ LMI + # BWH, HMS, Boston, MA, USA. + + if do_padding: + zpt = safe_longest_time_const(hw) * 4 / 1000 # s + pad1 = round(zpt/4/dt) + pad2 = round(zpt/1/dt) + + gwf = np.pad(gwf, ((pad1, pad2), (0,0))) + rf = np.pad(rf, (pad1, pad2)) + + safe_hw_check(hw) + + dgdt = np.diff(gwf, axis=0) / dt + pns = np.zeros(dgdt.shape) + + pns[:,0] = safe_pns_model(dgdt[:,0], dt, hw.x) + pns[:,1] = safe_pns_model(dgdt[:,1], dt, hw.y) + pns[:,2] = safe_pns_model(dgdt[:,2], dt, hw.z) + + # Export relevant paramters + res = SimpleNamespace() + res.pns = pns + res.gwf = gwf + res.rf = rf + res.dgdt = dgdt + res.dt = dt + res.hw = hw + + return pns, res + +def safe_plot(pns, dt=None, envelope=True, envelope_points=500): + # function h = safe_plot(pns, dt) + # pns is relative PNS waveform (nx3) + # dt is time step size in seconds. + + pnsnorm = np.sqrt((pns**2).sum(axis=1)) + + # FZ: Added option to plot the moving maximum of pns and pnsnorm to keep + # plots for long sequences intelligible + if envelope and pns.shape[0] > envelope_points: + N = int(np.ceil(pns.shape[0] / envelope_points)) + if dt != None: + dt *= N + + if pns.shape[0] % N != 0: + pns = np.concatenate((pns, np.zeros((N - pns.shape[0] % N, pns.shape[1])))) + pnsnorm = np.concatenate((pnsnorm, np.zeros((N - pnsnorm.shape[0] % N)))) + + pns = pns.reshape(pns.shape[0]//N, N, pns.shape[1]) + pns = pns.max(axis=1) + pnsnorm = pnsnorm.reshape(pnsnorm.shape[0]//N, N) + pnsnorm = pnsnorm.max(axis=1) + + if dt == None: + ttot = 1 # au + xlabstr = 'Time [a.u.]' + else: + ttot = pns.shape[0] * dt * 1000 # ms + xlabstr = 'Time [ms]' + + + t = np.linspace(0, ttot, pns.shape[0]) + + plt.plot(t, pns[:,0], 'r-', + t, pns[:,1], 'g-', + t, pns[:,2], 'b-', + t, pnsnorm , 'k-') + + plt.ylim([0, 120]) + plt.xlim([min(t), max(t)]) + + plt.title(f'Predicted PNS ({max(pnsnorm):0.0f}%)') + + plt.xlabel(xlabstr) + plt.ylabel('Relative stimulation [%]') + + plt.plot([0, max(t)], [max(pnsnorm), max(pnsnorm)], 'k:') + + plt.legend([f'X ({max(pns[:,0]):0.0f}%)', + f'Y ({max(pns[:,1]):0.0f}%)', + f'Z ({max(pns[:,2]):0.0f}%)', + f'nrm ({max(pnsnorm):0.0f}%)'], loc='best') + + +def safe_example(): + # Load an exampe gradient waveform + [gwf, rf, dt] = safe_example_gwf() + + # Load reponse parameters for example hardware + hw = safe_example_hw() + + # Check if hardware parameters are consistent + safe_hw_check(hw) + + # Check if this hw is part of the library (validate hw) + # safe_hw_verify(hw) + + # Predict PNS levels + pns, res = safe_gwf_to_pns(gwf, rf, dt, hw, 1) + + # Plot some results + safe_plot(pns, dt) + + +if __name__ == '__main__': + safe_example() diff --git a/pypulseq/utils/siemens/asc_to_hw.py b/pypulseq/utils/siemens/asc_to_hw.py new file mode 100644 index 00000000..e3bef621 --- /dev/null +++ b/pypulseq/utils/siemens/asc_to_hw.py @@ -0,0 +1,72 @@ +from types import SimpleNamespace +import numpy as np + + +def asc_to_hw(asc : dict) -> SimpleNamespace: + """ + Convert ASC dictionary from readasc to SAFE hardware description. + + Parameters + ---------- + asc : dict + ASC dictionary, see readasc + + Returns + ------- + SimpleNamespace + SAFE hardware description + """ + hw = SimpleNamespace() + + if 'asCOMP' in asc and 'tName' in asc['asCOMP']: + hw.name = asc['asCOMP']['tName'] + else: + hw.name = 'unknown' + + if 'GradPatSup' in asc: + asc_pns = asc['GradPatSup']['Phys']['PNS'] + else: + asc_pns = asc + + #hw.look_ahead = 1.0 # MZ: this is not a real hardware parameter but a coefficient, with which the final result is multiplied + hw.x = SimpleNamespace() + hw.x.tau1 = asc_pns['flGSWDTauX'][0] # ms + hw.x.tau2 = asc_pns['flGSWDTauX'][1] # ms + hw.x.tau3 = asc_pns['flGSWDTauX'][2] # ms + hw.x.a1 = asc_pns['flGSWDAX'][0] + hw.x.a2 = asc_pns['flGSWDAX'][1] + hw.x.a3 = asc_pns['flGSWDAX'][2] + hw.x.stim_limit = asc_pns['flGSWDStimulationLimitX'] # T/m/s + hw.x.stim_thresh = asc_pns['flGSWDStimulationThresholdX'] # T/m/s + + hw.y = SimpleNamespace() + hw.y.tau1 = asc_pns['flGSWDTauY'][0] # ms + hw.y.tau2 = asc_pns['flGSWDTauY'][1] # ms + hw.y.tau3 = asc_pns['flGSWDTauY'][2] # ms + hw.y.a1 = asc_pns['flGSWDAY'][0] + hw.y.a2 = asc_pns['flGSWDAY'][1] + hw.y.a3 = asc_pns['flGSWDAY'][2] + hw.y.stim_limit = asc_pns['flGSWDStimulationLimitY'] # T/m/s + hw.y.stim_thresh = asc_pns['flGSWDStimulationThresholdY'] # T/m/s + + hw.z = SimpleNamespace() + hw.z.tau1 = asc_pns['flGSWDTauZ'][0] # ms + hw.z.tau2 = asc_pns['flGSWDTauZ'][1] # ms + hw.z.tau3 = asc_pns['flGSWDTauZ'][2] # ms + hw.z.a1 = asc_pns['flGSWDAZ'][0] + hw.z.a2 = asc_pns['flGSWDAZ'][1] + hw.z.a3 = asc_pns['flGSWDAZ'][2] + hw.z.stim_limit = asc_pns['flGSWDStimulationLimitZ'] # T/m/s + hw.z.stim_thresh = asc_pns['flGSWDStimulationThresholdZ'] # T/m/s + + if 'asGPAParameters' in asc: + hw.x.g_scale = asc['asGPAParameters'][0]['sGCParameters']['flGScaleFactorX'] + hw.y.g_scale = asc['asGPAParameters'][0]['sGCParameters']['flGScaleFactorY'] + hw.z.g_scale = asc['asGPAParameters'][0]['sGCParameters']['flGScaleFactorZ'] + else: + print('Warning: Gradient scale factors not in ASC file: assuming 1/pi') + hw.x.g_scale = 1/np.pi + hw.y.g_scale = 1/np.pi + hw.z.g_scale = 1/np.pi + + return hw diff --git a/pypulseq/utils/siemens/readasc.py b/pypulseq/utils/siemens/readasc.py new file mode 100644 index 00000000..1f5602de --- /dev/null +++ b/pypulseq/utils/siemens/readasc.py @@ -0,0 +1,97 @@ +import re +from typing import Tuple + +def readasc(filename : str) -> Tuple[dict, dict]: + """ + Reads Siemens ASC ascii-formatted textfile and returns a dictionary + structure. + E.g. a[0].b[2][3].c = "string" + parses into: + asc['a'][0]['b'][2][3]['c'] = "string" + + Parameters + ---------- + filename : str + Filename of the ASC file. + + Returns + ------- + asc : dict + Dictionary of ASC part of file. + extra : dict + Dictionary of other fields after "ASCCONV END" + """ + + asc, extra = {}, {} + + # Read asc file and convert it into a dictionary structure + with open(filename, 'r') as fp: + end_of_asc = False + + for next_line in fp: + next_line = next_line.strip() + + if next_line == '### ASCCONV END ###': # find end of mrProt in the asc file + end_of_asc = True + + if next_line == '' or next_line[0] == '#': + continue + + # regex wizardry: Matches lines like 'a[0].b[2][3].c = "string" # comment' + # Note this assumes correct formatting, e.g. does not check whether + # brackets match. + match = re.match('^\s*([a-zA-Z0-9\[\]\._]+)\s*\=\s*((\"[^\"]*\"|\\\'[^\\\']\\\')|(\d+)|([0-9\.e\-]+))\s*((#|\/\/)(.*))?$', next_line) + + if match: + field_name = match[1] + + # Keep track of where to put the value: base[assign_to] = value + if end_of_asc: + base = extra + else: + base = asc + + assign_to = None + + # Iterate over every segment of the field name + parts = field_name.split('.') + for p in parts: + # Update base so final assignement is like: base[assign_to][p] = value + if assign_to != None and assign_to not in base: + base[assign_to] = {} + if assign_to != None: + base = base[assign_to] + + # Iterate over brackets + start = p.find('[') + if start != -1: + name = p[:start] + assign_to = name + + while start != -1: + stop = p.find(']', start) + index = int(p[start+1:stop]) + + # Update base so final assignement is like: base[assign_to][p][index] = value + if assign_to not in base: + base[assign_to] = {} + base = base[assign_to] + assign_to = index + + start = p.find('[', stop) + else: + assign_to = p + + # Depending on which regex section matched we can infer the value type + if match[3]: + base[assign_to] = match[3][1:-1] + elif match[4]: + base[assign_to] = int(match[4]) + elif match[5]: + base[assign_to] = float(match[5]) + else: + raise RuntimeError('This should not be reached') + elif next_line.find('=') != -1: + raise RuntimeError(f'Bug: ASC line with an assignment was not parsed correctly: {next_line}') + + return asc, extra