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.
Here, I present various variations of bSSFP: standard bSSFP, fluctuating equilibrium, and alternating steady states. Each variation manipulates sequence parameters such as TR, flip angles, and RF phase to generate unique profiles optimized for specific tasks.
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
[2]
Author: Krishna S. Nayak
Title: Wideband SSFP: Alternating repetition time balanced steady state free precession with increased band spacing
Link: https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.21296
[3]
Author: S S Vasanawala
Title: Fluctuating equilibrium MRI
Link: https://pubmed.ncbi.nlm.nih.gov/10542345/In [59]:
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 *
In [60]:
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()
In [61]:
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¶
In [62]:
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")
Fluctuating equilibrium MRI¶
In [63]:
size = 200 # resolution
Nr = 150
alpha = np.ones(Nr) * np.pi/6
T1 = 170
T2 = 40
M0 = 1
TR = np.ones(Nr) * 2.2
TE = 2.2/2
phi = np.ones(size)
phi[::2] = 0
phi[1::2] = np.pi/2
dphi = np.ones(Nr) * 0
f = 2
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, "FEMR")
Alternating steady state 2-1-1¶
In [64]:
size = 100
TR = [3.1, 0.9, 0.9] * 1000 + [3.1]
Nr = len(TR)
alpha = [np.pi/3] * Nr
T1 = 1000
T2 = 200
M0 = 1
TE = TR[-1]/2
M0 = 1
phi = [0] * Nr
dphi = [0] * Nr
f = 2
b = np.linspace(-f*np.pi, f*np.pi, size)
for dphi in np.linspace(0, np.pi, 2):
dphi = [dphi] * Nr
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, "alternating steady state 2-1-1")
Alternating steady state 2-2-1-1¶
In [65]:
size = 100
TR = [3.45, 1.725, 1.725, 1.725] * 200 + [3.45]
Nr = len(TR)
alpha = [np.pi/3] * Nr
T1 = 1000
T2 = 200
M0 = 1
TE = TR[-1]/2
M0 = 1
phi = [0] * Nr
dphi = [0] * Nr
f = 2
b = np.linspace(-f*np.pi, f*np.pi, size)
for dphi in np.linspace(0, np.pi, 2):
dphi = [dphi] * Nr
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,"")
