seq_UI_lfmri/LF_scanner/pypulseq/utils/safe_pns_prediction.py

# 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()