#%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, Bounds
import ipywidgets as widgets
from IPython.display import display, clear_output
check_FAplot = widgets.Checkbox(value=False,description="show FA-plot")
check_TIplot = widgets.Checkbox(value=False,description="show TI-plot")
def reset(T1):
Mz = np.ones(len(T1))
return Mz
def invert(Mz,eff=1.0):
Mz = Mz * -eff
return Mz
def relax(Mz,T1,t):
Mz = 1 - (1-Mz) * np.exp(-t/T1);
return Mz
def read(Mz,FA,N,T1,TR1):
for i in range(N):
sig = np.sin(FA / 180 * np.pi) * Mz;
Mz = Mz * np.cos(FA / 180 * np.pi);
Mz = relax(Mz,T1,TR1)
return Mz,sig
def singleTR(Mz,param):
# calculate helper times
readDur = param['N'] * param['TR1']
TIfill = param['TI'] - readDur/2
TRfill = param['TR0'] - readDur - TIfill
N1 = int(param['N']/2)
N2 = param['N'] - N1
Mz = invert(Mz)
Mz = relax(Mz,param['T1'],TIfill)
Mz,sig = read(Mz,FA=param['FA'],N=N1,T1=param['T1'],TR1=param['TR1'])
Mz = read(Mz,FA=param['FA'],N=N2,T1=param['T1'],TR1=param['TR1'])[0]
Mz = relax(Mz,param['T1'],TRfill)
return Mz,sig
# run enough times to approximate steady state
def runParameters(param):
Mz = reset(param['T1'])
niter=7
for i in range(niter):
Mz,sig = singleTR(Mz,param)
return sig
def plotFAdependence(param):
param2 = param.copy()
FAs = np.linspace(param['FA']-2,param['FA']+2,30)
cont = np.zeros(len(FAs))
sG = np.zeros(len(FAs))
sW = np.zeros(len(FAs))
for i,f in enumerate(FAs):
param2['FA'] = f
sig = runParameters(param2)
cont[i] = sig[1]/sig[0]
sG[i] = sig[0]
sW[i] = sig[1]
fig, ax1 = plt.subplots()
ax1.plot(FAs,cont,color='C0')
ax1.set_title(f'FA variation (FA={param["FA"]:.1f} deg)')
ax1.set_ylabel('WM/GM contrast', color='C0')
ax1.vlines(param['FA'],ymin=min(cont),ymax=max(cont),color='gray')
ax1.set_xlabel('FA [deg]')
ax2 = ax1.twinx()
ax2.plot(FAs,sG,color='C1')
ax2.plot(FAs,sW,color='C2')
plt.grid()
ax2.set_ylabel('GM signal', color='C1')
ax2.legend(['GM','WM'],loc='right')
plt.show()
def plotTIdependence(param):
param2 = param.copy()
TIs = np.linspace(param["TI"]-200,param["TI"]+200,30)
cont = np.zeros(len(TIs))
sG = np.zeros(len(TIs))
sW = np.zeros(len(TIs))
for i,t in enumerate(TIs):
param2['TI'] = t
sig = runParameters(param2)
cont[i] = sig[1]/sig[0]
sG[i] = sig[0]
sW[i] = sig[1]
fig, ax1 = plt.subplots()
ax1.plot(TIs,cont,color='C0')
ax1.set_title(f'TI variation (TI={param["TI"]:.0f} ms)')
ax1.set_ylabel('WM/GM contrast', color='C0')
ax1.vlines(x=param["TI"],ymin=min(cont),ymax=max(cont),color='gray')
ax2 = ax1.twinx()
ax2.plot(TIs,sG,color='C1')
ax2.plot(TIs,sW,color='C2')
plt.grid()
ax1.set_xlabel('TI [ms]')
ax2.set_ylabel('GM signal', color='C1')
ax2.legend(['GM','WM'],loc='right')
plt.show()
def optimize_values(param,target_cont,target_sig):
def cost(x,verbose=False):
param['FA'] = x[0]
param['TI'] = x[1]
sig = runParameters(param)
cont = sig[1]/sig[0]
cost_cont = abs(cont-target_cont)**2
cost_sig = abs(sig[0]-target_sig)**2 * 1e1
cost_sig2 = abs(1/sig[0])* 1e-3
cost_total = cost_cont + cost_sig
if verbose:
print(f'cont: {cont:6.3f}')
print(f'GM-sig: {sig[0]*100:6.3f}')
return cost_total
x0 = [param['FA'],param['TI']]
lb = [1,np.ceil(param['TR1']*param['N']/2)]
ub = [90,param['TR0']-np.ceil(param['TR1']*param['N']/2)]
bounds = Bounds(lb,ub)
x1 = minimize(cost,x0,method='L-BFGS-B',bounds=bounds)
with output:
clear_output()
print(f'Duration: {param["N_part"]*param["TR0"]/1000/60:.2f} min')
# print('\nOld:')
# print(f'FA: {x0[0]:7.2f} deg')
# print(f'TI: {x0[1]:7.2f} ms')
# cost(x0,verbose=True)
print('\nNew:')
print(f'FA: {x1.x[0]:7.2f} deg')
print(f'TI: {x1.x[1]:7.2f} ms')
cost(x1.x,verbose=True)
if check_FAplot.value:
plotFAdependence(param)
if check_TIplot.value:
plotTIdependence(param)
TR_slider = widgets.IntSlider(
value=3360,
min=1000,
max=8000,
step=10,
description='Segment-TR:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d'
)
TR1_slider = widgets.FloatSlider(
value=6.1,
min=1,
max=20,
step=0.1,
description='Readout-TR:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='.1f'
)
N_slider = widgets.IntSlider(
value=132,
min=20,
max=500,
step=1,
description='Seg Length:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d'
)
Npart_slider = widgets.IntSlider(
value=120,
min=10,
max=500,
step=1,
description='N_segments:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='d'
)
Cont_slider = widgets.FloatSlider(
value=1.7,
min=1.0,
max=3.0,
step=0.1,
description='contrast:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='.1f'
)
sig_slider = widgets.FloatSlider(
value=4.0,
min=0.0,
max=4.0,
step=0.1,
description='Signal target:',
disabled=False,
continuous_update=False,
orientation='horizontal',
readout=True,
readout_format='.1f'
)
output = widgets.Output()
def update_values(change):
param = dict()
param['T1'] = np.asarray([2002,1425])
param['N'] = 132
param['TR1'] = 6.06
param['TR0'] = 3360
param['FA'] = 9.0
param['TI'] = 1340
param['N_part'] = Npart_slider.value
param['TR0'] = TR_slider.value
param['TR1'] = TR1_slider.value
param['N'] = N_slider.value
target_cont = Cont_slider.value
target_sig = sig_slider.value/100
optimize_values(param,target_cont,target_sig)
Assuming T1=2002 for GM and T1=1425 for WM
Set sequence parameters (TR,segment length, etc) first. Then select target WM/GM contrast.
The displayed signal is sin(alpha) * Mz * 100, so a 90° pulse after full relaxation would result in signal 100. This is a slightly arbitrary value and only gives an idea about the performance.
sig_slider.observe(update_values, 'value')
TR_slider.observe(update_values, 'value')
N_slider.observe(update_values, 'value')
Cont_slider.observe(update_values, 'value')
Npart_slider.observe(update_values,'value')
TR1_slider.observe(update_values,'value')
check_FAplot.observe(update_values,'value')
check_TIplot.observe(update_values,'value')
display(TR_slider)
display(TR1_slider)
display(N_slider)
display(Npart_slider)
display(Cont_slider)
display(check_FAplot)
display(check_TIplot)
#display(sig_slider)
display(output)
update_values(0)
IntSlider(value=3360, continuous_update=False, description='Segment-TR:', max=8000, min=1000, step=10)
FloatSlider(value=6.1, continuous_update=False, description='Readout-TR:', max=20.0, min=1.0, readout_format='…
IntSlider(value=132, continuous_update=False, description='Seg Length:', max=500, min=20)
IntSlider(value=120, continuous_update=False, description='N_segments:', max=500, min=10)
FloatSlider(value=1.7, continuous_update=False, description='contrast:', max=3.0, min=1.0, readout_format='.1f…
Output()
Old:
FA: 9.00 deg
TI: 1340.00 ms
cont: 1.692
GM-sig: 1.909
param = dict()
param['T1'] = np.asarray([2002,1425]) # GM,WM
param['N'] = 132
param['N_part'] = 120
param['TR1'] = 6.06
param['TR0'] = 3360
param['FA'] = 9.0
param['TI'] = 1340
target_cont = 2.0
target_sig = 0.01
optimize_values(param,target_cont,target_sig)