SSFP simualtions

Author: Zimu Huo

Date: 02.2023

In bSSFP, the entire spin evolution can be modeled mathematically through bloch equations, which can be broken down into three components: excitation E, relaxation R, and precession P. These three components are matrix operators representing the sequence events. Immediately after the RF pulse R, the spin undergoes T1 and T2 relaxation E and off-resonance precession P. The bSSFP sequence achieves the steady state by repetitively executing these operations, while ensuring that the gradients are refocused within each TR.

---

References

[1] 
Author: Dr Neal K Bangerter
Title: Contrast enhancement and artifact reduction in steady state magnetic resonance imaging
Link: https://www.proquest.com/openview/41a8dcfb0f16a1289210b3bd4f9ea82b/1.pdf?cbl=18750&diss=y&pq-origsite=gscholar

import sys
sys.path.insert(1, '../')
from matplotlib import pyplot as plt
import numpy as np
import sys
sys.path.insert(1, '../')
sys.path.insert(1, '../../')
from simulator import *

def plot_bssfp(x, data, title):
    magnitude, phase = np.absolute(data), np.angle(data)
    plt.figure()
    plt.subplot(211)
    plt.ylabel('Magitude')
    plt.title(title)
    plt.grid(True)
    plt.plot(x, magnitude)
    plt.subplot(212)
    plt.xlabel('Off-Resonance (Hz)')
    plt.ylabel('Phase')
    plt.grid(True)
    plt.plot(x, phase)
    plt.show()

def rotMat(alpha, phi):
    # RF rotation matrix 
    # Author Zimu Huo
    rotation = np.array([
    [np.cos(alpha)*np.sin(phi)**2 + np.cos(phi)**2,          
         (1-np.cos(alpha))*np.cos(phi)*np.sin(phi),         
     -np.sin(alpha)*np.sin(phi)],
        
    [    (1-np.cos(alpha))*np.cos(phi)*np.sin(phi),        
       np.cos(alpha)*np.cos(phi)**2+np.sin(phi)**2,         
                        np.sin(alpha)*np.cos(phi)],
        
    [                    np.sin(alpha)*np.sin(phi),                      
                        -np.sin(alpha)*np.cos(phi),                    
                                     np.cos(alpha)]
    ])
    return rotation
    
def iterative_SSFP(M0, alpha, phi, dphi, beta, TR, TE, T1, T2):
    # SSFP simulator
    # Author Zimu Huo
    signal = np.zeros([1,3])
    signal = np.asarray([0, 0, M0])
    Nr = len(TR)
    for i in range(Nr-1):
        signal = np.matmul(rotMat(alpha[i], phi[i]+np.sum(dphi[0:i])),signal) 
        signal[0] = signal[0]*np.exp(-TR[i]/T2)
        signal[1] = signal[1]*np.exp(-TR[i]/T2)
        signal[2] = M0+(signal[2]-M0)*np.exp(-TR[i]/T1)
        P = np.array([
            [ np.cos(beta[i]),  np.sin(beta[i]),   0],
            [-np.sin(beta[i]),  np.cos(beta[i]),   0],
            [            0,             0,   1]
            ])
        signal = np.matmul(P,signal)
        
    signal = np.matmul(rotMat(alpha[-1], phi[-1]+np.sum(dphi[:-1])),signal) 
    
    signal[0] = signal[0]*np.exp(-TE/T2)
    signal[1] = signal[1]*np.exp(-TE/T2)
    signal[2] = M0+(signal[2]-M0)*np.exp(-TE/T1) 
    P = np.array([
            [ np.cos(beta[-1]),  np.sin(beta[-1]),   0],
            [-np.sin(beta[-1]),  np.cos(beta[-1]),   0],
            [            0,             0,      1]
            ])
    signal = np.matmul(P,signal)
    return signal
def bssfp_bloch(b, size,  M0, alpha, phi, dphi, TR, TE, T1, T2):
    # SSFP simulator
    # Author Zimu Huo
    samples = np.zeros([size,3])
    Nr = len(TR)
    for index, beta in enumerate((b)):
        # some safty checks 
        betas = np.ones(Nr) * beta * TR[:Nr] / np.max(TR[:Nr])
        betas[-1] = beta * TE / np.max(TR[:Nr])
        ntime = int(Nr)
        mx = np.zeros(Nr)
        my = np.zeros(Nr)
        mz = np.zeros(Nr)
        alpha = np.asarray(alpha).astype("float64")[:Nr]
        beta = np.asarray(betas).astype("float64")[:Nr]
        TR = np.asarray(TR).astype("float64")[:Nr]
        phi = np.asarray(phi).astype("float64")[:Nr]
        dphi = np.asarray(dphi).astype("float64")[:Nr]
        
        samples[index,:] = iterative_SSFP(M0 = M0 , alpha = alpha, phi = phi, dphi = dphi, beta = betas, TR= TR, TE= TE, T1 = T1, T2 = T2)
    data = np.zeros([size,1], dtype = complex)
    data.real = samples[:,1].reshape(-1,1)
    data.imag = samples[:,0].reshape(-1,1)
    axis = b / (2*np.pi*TR[0]/1000)
    return axis, data

# bSSFP

size = 200
Nr = 2000
alpha = np.ones(Nr) * np.pi/6
T1 = 1000
T2 = 200
M0 = 1
TR = np.ones(Nr) * 10
TE = 10/2
M0 = 2
phi = np.ones(Nr) * 0
dphi = np.ones(Nr) * 0
f = 4
b = np.linspace(-f*np.pi, f*np.pi, size)
x, data = bssfp_bloch(b, size, M0 = M0 , alpha = alpha, phi = phi, dphi = dphi, TR= TR, TE= TE, T1 = T1, T2 = T2)
plot_bssfp(x, data, "bSSFP")

size = 200
Nr = 2000

T1 = 1000
T2 = 200
M0 = 1
TR = np.ones(Nr) * 10
TE = 10/2
M0 = 2
phi = np.ones(Nr) * 0
dphi = np.ones(Nr) * 0
f = 2
b = np.linspace(-f*np.pi, f*np.pi, size)

signal = np.zeros([size, 4], dtype = complex)
for idx, flip_angle in enumerate (np.asarray([np.pi/ 9, np.pi / 6, np.pi / 4, np.pi / 2])):
    alpha = np.ones(Nr) * flip_angle
    x, data = bssfp_bloch(b, size, M0 = M0 , alpha = alpha, phi = phi, dphi = dphi, TR= TR, TE= TE, T1 = T1, T2 = T2)
    signal[:,idx] = data.flatten()

plt.figure()
plt.ylabel('Magitude')
plt.title("bSSFP profile vs flip angle")
plt.xlabel('Off-Resonance (Hz)')
plt.grid(True)
plt.plot(x, np.abs(signal))
plt.legend([20, 30, 45, 60])
plt.show()

## Different TR

size = 200
Nr = 2000

T1 = 1000
T2 = 200

M0 = 1
TE = 10/2
M0 = 2
phi = np.ones(Nr) * 0
dphi = np.ones(Nr) * 0
f = 2
b = np.linspace(-f*np.pi, f*np.pi, size)

signal = np.zeros([size, 4], dtype = complex)
xs = np.zeros([size, 4], dtype = complex)
for idx, TR in enumerate (np.asarray([2, 4, 6, 8])):
    TE = TR/2
    TR = np.ones(Nr) * TR
    x, data = bssfp_bloch(b, size, M0 = M0 , alpha = alpha, phi = phi, dphi = dphi, TR= TR, TE= TE, T1 = T1, T2 = T2)
    signal[:,idx] = data.flatten()
    xs[:,idx] = x.flatten()

plt.figure()
plt.ylabel('Magitude')
plt.title("bSSFP profile vs TR")
plt.grid(True)
plt.plot(xs, np.abs(signal))
plt.xlabel('Off-Resonance (Hz)')
plt.legend([2, 4, 6, 8])
plt.show()

num_phasecyle = 2
size = 200
Nr = size
alpha = np.ones(Nr) * np.pi/6
T1 = 1000
T2 = 100
M0 = 10
TR = np.ones(Nr) * 10
TE = 10/2
M0 = 2
phi = np.ones(Nr) * 0
dphi = np.ones(Nr) * 0
f = 2
b = np.linspace(-np.pi*f, np.pi*f, Nr)
signal = []
for dphi in np.linspace(0, 2*np.pi, num_phasecyle , endpoint=False):
    dphi = np.ones(Nr) * dphi
    x, data = bssfp_bloch(b, size, M0 = M0 , alpha = alpha, phi = phi, dphi = dphi, TR= TR, TE= TE, T1 = T1, T2 = T2)
    signal.append(data)

plt.figure()
plt.ylabel('Magitude')
plt.title("phase cycled bSSFP profile")
plt.xlabel('Off-Resonance (Hz)')
plt.grid(True)
plt.plot(x, np.abs(np.squeeze(signal)).T)
plt.show()