diff --git a/Dockerfile b/Dockerfile index d89ef26..2de89a3 100644 --- a/Dockerfile +++ b/Dockerfile @@ -20,7 +20,7 @@ RUN apt-get update && \ wget \ gfortran \ cmake \ - bsdtar \ + bsdtar \ rsync \ imagemagick \ gnuplot-x11 \ @@ -47,6 +47,7 @@ RUN cd $HOME/work;\ pip install pyarrow>=0.4.0; \ pip install octave_kernel \ sos==0.17.7 \ + sos-r \ sos-notebook==0.17.2 \ sos-python==0.9.12.1 \ sos-bash==0.12.3 \ @@ -59,12 +60,12 @@ RUN cd $HOME/work;\ scipy \ plotly==3.10.0 \ flask; \ - python -m sos_notebook.install;\ - git clone --single-branch -b master https://github.com/qMRLab/t1_notebooks.git; \ + python -m sos_notebook.install; \ + git clone --single-branch -b blog_afi https://github.com/qMRLab/t1_notebooks.git; \ cd t1_notebooks;\ git clone https://github.com/neuropoly/qMRLab.git; \ cd qMRLab; \ - git checkout d15a553f9d93457c3ed59861380852c54458c2b4; \ + git checkout mb/afi_b1_filter_patch; \ cd ..; \ chmod -R 777 $HOME/work/t1_notebooks; \ octave --eval "cd qMRLab; \ diff --git a/afi_blog/ActualFlipAngleImaging.html b/afi_blog/ActualFlipAngleImaging.html new file mode 100644 index 0000000..be545d4 --- /dev/null +++ b/afi_blog/ActualFlipAngleImaging.html @@ -0,0 +1,15212 @@ + + + + + + +ActualFlipAngleImaging + + +
+ +
+ + + +
+ + +
+ +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ Display content:
+
+ + +
+ + +
+ + +
+
+
+ + + + + +
+
+ + +
+ +
+
+
In [1]:
+ +
+
+
% **Blog post code introduction**
+% 
+% Congrats on activating the "All cells" option in this interactive blog post =D
+%
+% Below, several new HTML blocks have appears prior to the figures, displaying the Octave/MATLAB code that was used to generate the figures in this blog post.
+%
+% If you want to reproduce the data on your own local computer, you simply need to have qMRLab installed in your Octave/MATLAB path and run the "startup.m" file, as is shown below.
+%
+% If you want to get under the hood and modify the code right now, you can do so in the Jupyter Notebook of this blog post hosted on MyBinder. The link to it is in the introduction above.
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+
+
In [2]:
+
+ +
+
+ + +
+
+ +
+
+ + +
+
+
+
+

Actual flip-angle imaging B1 mapping

+

+

+
+
+
+
+
+
+
+

+Transmit radiofrequency field maps (B1+, or B1 for short) are used in diverse applications in MRI including: the study of electrical properties in tissues in vivo (Sled & Pike 1998; Katscher et al. 2009), specific absorption rate (SAR) calculations (Ibrahim et al. 2001), the calibration of quantitative T1 (Deoni 2007; Boudreau et al. 2017) and T2 (Sled and Pike 2000) maps, better parameter estimation from magnetization transfer measurements (Ropele et al. 2005; Boudreau et al. 2018), B1 shimming to improve image quality at whole-body ultra high fields (van den Bergen et al. 2007), or quality control of RF coils (Yarnykh 2007). Several B1 mapping techniques have been developed, and they can be broadly divided as magnitude-based and phase-based methods. The double angle method (DAM) is a saturation-recovery magnitude-based method that takes the ratio of the signal intensity of two magnitude images measured with different excitation flip angles (Insko & Bolinger 1993; Stollberger and Wach 1996). The Bloch-Siegert shift technique is a rapid phase-based method that encodes the B1 information into phase signal (Sacolick et al. 2010). The actual flip-angle imaging (AFI) is a magnitude-based B1 mapping method that consists of a 3D acquisition that benefits from good anatomical coverage. In addition, this technique allows the acquisitions of whole-body (~7 min) and brain (~3 min) B1 maps leading to a feasible implementation in clinics (Yarnykh 2004; Yarnykh 2007). On the other hand, the AFI pulse sequence has certain constraints that need to be considered for this B1 mapping method to be widely deployed. Some of the limitations include the use of spoiler gradients that can give rise to prohibitive SAR values (Sacolick et al. 2010), and the pulse sequence modifications on the MRI machine to implement the AFI method. +

+ +

+In this blog post, we will focus on presenting details about the AFI B1 mapping method. We will cover signal modeling, data fitting, the benefits and the pitfalls of the technique. The figures are generated using the qMRLab module for this method. +

+
+
+
+
+
+
+
+
+

Signal Modelling

+
+
+
+
+
+
+
+
+

+The pulse sequence of the AFI method (Figure 1) is composed of two identical RF pulses and two different delays (TR1 < TR2). After each RF pulse, the signal intensity is acquired followed by a spoiler to destroy the residual transverse magnetization next to the following RF pulse. This method implements a pulsed steady-state signal with a gradient-echo acquisition, thus preventing the use of long repetition times (Yarnykh 2007). It has been demonstrated that if the delays TR1 and TR2 are sufficiently short (e.g. TR1/TR2 = 20 ms/100 ms), and the transverse magnetization is completely spoiled, the ratio of signal intensities (r = S2/S1) depends on the flip angle of applied pulses and is highly insensitive to T1 (Yarnykh 2007). +

+
+
+
+
+
+
+
+
+
+

+ +Figure 1. Simplified pulse sequence diagram of an actual flip-angle imaging (AFI) pulse sequence with a gradient echo readout. TR1: repetition time 1, TR2: repetition time 2, θ: excitation flip angle for the measurement, IMG: image acquisition (k-space readout), SPOIL: spoiler gradient. + +

+

+
+
+
+
+
+
+
+
+
+

+The magnetization of an AFI experiment can be modeled under steady-state conditions by the implementation of a fast repetition of the sequence (TR1 < TR2 < T1). The solution of the Bloch equations for the AFI method is given by Equations 1 and 2 that represent the longitudinal magnetization before the application of the RF pulses: +

+ +

+

+

+ +

+

+

+ +

+Mz1,2 is the longitudinal magnetization of both pulses, M0 is the magnetization at thermal equilibrium, TR1 is the delay time after the first pulse, TR2 is the delay time after the second identical pulse (Figure 1), and θ is the excitation flip angle. The steady-state longitudinal magnetization Mz curves for different T1 values for a range of θn and TR values are shown in Figure 2. +

+ +
+
+
+
+
+
+
+
+
+

+ +Figure 2. Longitudinal magnetization before the first radiofrequency pulse (Equation 1, solid lines) and before the second identical pulse (Equation 2, dashed lines) for three different T1 values. + +

+
+
+
+
+ +
+ +
+
+
In [3]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Adds qMRLab to the path of the environment
+
+cd ../qMRLab
+startup
+
+ +
+
+
+
+ +
+
+ +
+ +
+ +
+ + +
+
warning: function /home/jovyan/work/t1_notebooks/qMRLab/External/NODDI_toolbox_v1.0/ParforProgMonv2/example.m shadows a core library function
+warning: called from
+    startup at line 2 column 1
+-----------------------------------------------
+Release v2.4.2
+Hawdy developer!  The version specified in version.txt is ahead of the latest published release v2.4.1.
+Please do not forget pushing a new commit tag upon merge or publish the new release.
+-----------------------------------------------
+loading struct
+loading io
+loading statistics
+loading optim
+loading image
+
+
+
+ +
+
+
+ +
+
+ + + +
+ +
+
+
In [4]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to generate the data required for Figure 2 of the blog post
+
+clear all
+
+%% Setup parameters
+% All times are in milliseconds
+% All flip angles are in degrees
+
+TR1_range = 5:5:50;
+n = 5;
+TR2_range = n*TR1_range;
+
+params.EXC_FA = 1:90;
+
+%% Calculate signals
+%
+% To see all the options available, run `help b1_afi.analytical_solution`
+
+for ii = 1:length(TR1_range)
+    params.TR1 = TR1_range(ii);
+    params.TR2 = TR2_range(ii);
+    
+    % White matter
+    params.T1 = 900; % in milliseconds
+    [Mz1, Mz2] = b1_afi.analytical_solution(params);
+    signal1_WM(ii,:) = Mz1;
+    signal2_WM(ii,:) = Mz2;
+
+    % Grey matter
+    params.T1 = 1500;  % in milliseconds
+    [Mz1, Mz2] = b1_afi.analytical_solution(params);
+    signal1_GM(ii,:) = Mz1;
+    signal2_GM(ii,:) = Mz2;
+
+    % CSF
+    params.T1 = 4000;  % in milliseconds
+    [Mz1, Mz2] = b1_afi.analytical_solution(params);
+    signal1_CSF(ii,:) = Mz1;
+    signal2_CSF(ii,:) = Mz2;
+end
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+
+
In [5]:
+
+ +
+
+ + + +
+ +
+
+
In [6]:
+ +
+
+
# PYTHON CODE
+
+init_notebook_mode(connected=True)
+
+data1_1 = [dict(
+        visible = False,
+        line=dict(color='royalblue'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_WM[ii]))),
+        name = 'S<sub>1</sub>: T<sub>1</sub> = 0.9 s (White Matter)',
+        text = 'S<sub>1</sub>: T<sub>1</sub> = 0.9 s (White Matter)',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data1_1[4]['visible'] = True
+
+data1_2 = [dict(
+        visible = False,
+        line=dict(color='royalblue', dash='dash'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal2_WM[ii]))),
+        name = 'S<sub>2</sub>: T<sub>1</sub> = 0.9 s (White Matter)',
+        text = 'S<sub>2</sub>: T<sub>1</sub> = 0.9 s (White Matter)',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data1_2[4]['visible'] = True
+
+data2_1 = [dict(
+        visible = False,
+        line=dict(color='firebrick'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_GM[ii]))),
+        name = 'S<sub>1</sub>: T<sub>1</sub> = 1.5 s (Grey Matter)',
+        text = 'S<sub>1</sub>: T<sub>1</sub> = 1.5 s (Grey Matter)',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data2_1[4]['visible'] = True
+
+data2_2 = [dict(
+        visible = False,
+        line=dict(color='firebrick', dash='dash'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal2_GM[ii]))),
+        name = 'S<sub>2</sub>: T<sub>1</sub> = 1.5 s (Grey Matter)',
+        text = 'S<sub>2</sub>: T<sub>1</sub> = 1.5 s (Grey Matter)',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data2_2[4]['visible'] = True
+
+data3_1 = [dict(
+        visible = False,
+        line=dict(color='orange'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_CSF[ii]))),
+        name = 'S<sub>1</sub>: T<sub>1</sub> = 4.0 s (Cerebrospinal Fluid)',
+        text = 'S<sub>1</sub>: T<sub>1</sub> = 4.0 s (Cerebrospinal Fluid)',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data3_1[4]['visible'] = True
+
+data3_2 = [dict(
+        visible = False,
+        line=dict(color='orange', dash='dash'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal2_CSF[ii]))),
+        name = 'S<sub>2</sub>: T<sub>1</sub> = 4.0 s (Cerebrospinal Fluid)',
+        text = 'S<sub>2</sub>: T<sub>1</sub> = 4.0 s (Cerebrospinal Fluid)',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data3_2[4]['visible'] = True
+
+data = data1_1 + data1_2 + data2_1 +data2_2 + data3_1 + data3_2
+
+steps = []
+for i in range(len(TR1_range)):
+    step = dict(
+        method = 'restyle',  
+        args = ['visible', [False] * len(data1_1)],
+        label = str(TR1_range[i])
+        )
+    step['args'][1][i] = True # Toggle i'th trace to "visible"
+    steps.append(step)
+
+sliders=[
+    dict(
+        x = 0,
+        y = -0.09,
+        active = 3,
+        currentvalue = {"prefix": "TR1 value (ms): <b>"},
+        pad = {"t": 50, "b": 10},
+        steps = steps)]
+
+layout = go.Layout(
+    width=650,
+    height=520,
+    margin=go.layout.Margin(
+        l=80,
+        r=40,
+        b=60,
+        t=10,
+    ),
+    annotations=[
+        dict(
+            x=0.5004254919715793,
+            y=-0.18,
+            showarrow=False,
+            text='Excitation Flip Angle (°)',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=-0.15,
+            y=0.5,
+            showarrow=False,
+            text='Long. Magnetization (M<sub>z</sub>)',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            textangle=-90,
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=-0.007,
+            y=-0.25,
+            showarrow=False,
+            text='TR<sub>2</sub> = 5TR<sub>1</sub>',
+            font=dict(
+                family='Times New Roman',
+                size=12
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis=dict(
+        autorange=False,
+        range=[0, params['EXC_FA'][-1]],
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    yaxis=dict(
+        autorange=True,
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    legend=dict(
+        x=0.5,
+        y=0.9,
+        traceorder='normal',
+        font=dict(
+            family='Times New Roman',
+            size=12,
+            color='#000'
+        ),
+        bordercolor='#000000',
+        borderwidth=2
+    ), 
+    sliders=sliders
+)
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-line', config = config)
+
+ +
+
+
+
+ +
+
+ + + + +
+ +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+
+ +

+The analytical solution of the Bloch equations in a steady-state experiment (Equation 1 and Equation 2) makes several assumptions leading to practical challenges. First, it is assumed that the longitudinal magnetization has reached a steady state after a sufficiently large number of repetition times (TR), and that the transverse magnetization is perfectly spoiled prior to each pulse. To explore these properties, a numerical approach known as Bloch simulations is used to estimate the signal from an MRI experiment given a set of sequence parameters. Here, the Bloch simulations allow us to estimate the magnetization using a different number of sequence repetitions, and look at a special case when the steady-state is not achieved (due to a small number of sequence repetitions). As can be seen in Figure 3, the number of repetitions required to reach a steady-state depends on T1 and the flip angle. +

+
+
+
+
+
+
+
+
+
+

+ +Figure 3. Signal 1 (blue) and Signal 2 (red) curves simulated using Bloch simulations (solid lines) for a number of repetitions ranging from 1 to 150, plotted against the ideal case (Equations 1 and 2 – dashed lines). Simulation details: TR1 = 20 ms, TR2 = 100 ms, T1 = 900 ms, 100 spins. Ideal spoiling was used for this set of Bloch simulations (transverse magnetization was set to 0 at the end of each TR1,2). + +

+
+
+
+
+ +
+ +
+
+
In [7]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to generate the data required for Figure 3 of the blog post
+
+clear all
+
+%% Setup parameters
+% All times are in milliseconds
+% All flip angles are in degrees
+
+% White matter
+params.T1 = 900; % in milliseconds
+params.T2 = 10000;
+params.TR1 = 20;
+params.TR2 = 100;
+params.TE = 5;
+params.EXC_FA = 1:90;
+Nex_range = 1:1:150;
+
+%% Calculate signals
+% This results were precomputed. Uncomment the lines below to run the simulation.
+signal1_bloch = load('../afi_blog/blochsim/signal1_blochsim.mat');
+signal2_bloch = load('../afi_blog/blochsim/signal2_blochsim.mat');
+signal1_analytical = load('../afi_blog/blochsim/signal1_analytical.mat');
+signal2_analytical = load('../afi_blog/blochsim/signal2_analytical.mat');
+signal1_blochsim = signal1_bloch.signal1_blochsim;
+signal2_blochsim = signal2_bloch.signal2_blochsim;
+signal1_analytical = signal1_analytical.signal1_analytical;
+signal2_analytical = signal2_analytical.signal2_analytical;
+
+%{
+% Bloch simulations
+% To see all the options available, run `help b1_afi.analytical_solution`
+
+for ii = 1:length(Nex_range)
+    params.Nex = Nex_range(ii);
+    
+    [Mz1, Mz2] = b1_afi.analytical_solution(params);
+    signal1_analytical(ii,:) = Mz1;
+    signal2_analytical(ii,:) = Mz2;
+
+    [~, Msig1, Msig2] = b1_afi.bloch_sim(params);
+    signal1_blochsim(ii,:) = abs(complex(Msig1));
+    signal2_blochsim(ii,:) = abs(complex(Msig2));
+end
+%}
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+
+
In [8]:
+
+ +
+
+ + + +
+ +
+
+
In [9]:
+ +
+
+
# PYTHON CODE
+
+init_notebook_mode(connected=True)
+
+data1_1 = [dict(
+        visible = False,
+        line=dict(color='royalblue', dash='dash'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_analytical[ii]))),
+        name = 'S<sub>1</sub>: Analytical Solution',
+        text = 'S<sub>1</sub>: Analytical Solution',
+        hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]
+
+data1_1[7]['visible'] = True
+
+data1_2 = [dict(
+        visible = False,
+        line=dict(color='firebrick', dash='dash'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal2_analytical[ii]))),
+        name = 'S<sub>2</sub>: Analytical Solution',
+        text = 'S<sub>2</sub>: Analytical Solution',
+        hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]
+
+data1_2[7]['visible'] = True
+
+data2_1 = [dict(
+        visible = False,
+        line=dict(color='royalblue'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_blochsim[ii]))),
+        name = 'S<sub>1</sub>: Bloch Simulation',
+        text = 'S<sub>1</sub>: Bloch Simulation',
+        hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]
+
+data2_1[7]['visible'] = True
+
+data2_2 = [dict(
+        visible = False,
+        line=dict(color='firebrick'),
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal2_blochsim[ii]))),
+        name = 'S<sub>2</sub>: Bloch Simulation',
+        text = 'S<sub>2</sub>: Bloch Simulation',
+        hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]
+
+data2_2[7]['visible'] = True
+
+data = data1_1 + data2_1 + data1_2 + data2_2
+
+steps = []
+for i in range(len(Nex_range)):
+    step = dict(
+        method = 'restyle',  
+        args = ['visible', [False] * len(data1_1)],
+        label = str(Nex_range[i])
+        )
+    step['args'][1][i] = True # Toggle i'th trace to "visible"
+    steps.append(step)
+
+sliders = [dict(
+    x = 0,
+    y = -0.02,
+    active = 5,
+    currentvalue = {"prefix": "n<sup>th</sup> TR1-TR2: <b>"},
+    pad = {"t": 50, "b": 10},
+    steps = steps)]
+
+layout = go.Layout(
+    width=580,
+    height=450,
+    margin=go.layout.Margin(
+        l=80,
+        r=40,
+        b=60,
+        t=10,
+    ),
+    annotations=[
+        dict(
+            x=0.5004254919715793,
+            y=-0.18,
+            showarrow=False,
+            text='Excitation Flip Angle (°)',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=-0.15,
+            y=0.5,
+            showarrow=False,
+            text='Signal',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            textangle=-90,
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis=dict(
+        autorange=False,
+        range=[0, params['EXC_FA'][-1]],
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    yaxis=dict(
+        autorange=True,
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    legend=dict(
+        x=0.5,
+        y=0.9,
+        traceorder='normal',
+        font=dict(
+            family='Times New Roman',
+            size=12,
+            color='#000'
+        ),
+        bordercolor='#000000',
+        borderwidth=2
+    ), 
+    sliders=sliders
+)
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-line', config = config)
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+
+

+In practice, gradient and RF spoiling are important parameters to consider in an AFI experiment. A combination of both (Zur et al. 1991; Bernstein et al. 2004) is typically recommended, and Figure 4 shows how this better approximates the ideal spoiling case. +

+
+
+
+
+
+
+
+
+
+

+ +Figure 4. Signal 1 curves estimated using Bloch simulations for three categories of signal spoiling: (1) ideal spoiling (blue), gradient & RF Spoiling (red), and no spoiling (orange). Simulation details: TR1 = 20 ms, TR2 = 100 ms, T1 = 900 ms, T2 = 100 ms, TE = 5 ms, 100 spins. For the ideal spoiling case, the transverse magnetization is set to zero at the end of each TR. For the gradient & RF spoiling case, each spin is rotated by different increments of phase (2𝜋 / # of spins) to simulate complete dephasing from gradient spoiling, and the RF phase of the excitation pulse is ɸn = ɸn-1 + nɸ0 = ½ ɸ0(n2 + n + 2) (Bernstein et al. 2004) with ɸ0 = 39° (Zur et al. 1991) after each TR. + +

+
+
+
+
+ +
+ +
+
+
In [10]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to generate the data required for Figure 4 of the blog post
+
+clear all
+
+%% Setup parameters
+% All times are in milliseconds
+% All flip angles are in degrees
+
+% White matter
+params.T1 = 900; % in milliseconds
+params.T2 = 100;
+params.TR1 = 20;
+params.TR2 = 100;
+params.TE = 5;
+params.EXC_FA = 1:90;
+Nex_range = [1:9, 10:10:100];
+
+%% Calculate signals
+% This results were precomputed. Uncomment the lines below to run the simulation.
+signal1_ideal = load('../afi_blog/blochsim/signal1_ideal_spoil.mat');
+signal2_ideal = load('../afi_blog/blochsim/signal2_ideal_spoil.mat');
+signal1_optimal = load('../afi_blog/blochsim/signal1_optimal_crush_and_rf_spoil.mat');
+signal2_optimal = load('../afi_blog/blochsim/signal2_optimal_crush_and_rf_spoil.mat');
+signal1_no_spoil = load('../afi_blog/blochsim/signal1_no_gradient_and_rf_spoil.mat');
+signal2_no_spoil = load('../afi_blog/blochsim/signal2_no_gradient_and_rf_spoil.mat');
+signal1_ideal_spoil = signal1_ideal.signal1_ideal_spoil;
+signal2_ideal_spoil = signal2_ideal.signal2_ideal_spoil;
+signal1_optimal_crush_and_rf_spoil = signal1_optimal.signal1_optimal_crush_and_rf_spoil;
+signal2_optimal_crush_and_rf_spoil = signal2_optimal.signal2_optimal_crush_and_rf_spoil;
+signal1_no_gradient_and_rf_spoil = signal1_no_spoil.signal1_no_gradient_and_rf_spoil;
+signal2_no_gradient_and_rf_spoil = signal2_no_spoil.signal2_no_gradient_and_rf_spoil;
+
+%{
+% To see all the options available, run `help b1_afi.analytical_solution`
+
+for ii = 1:length(Nex_range)
+    params.Nex = Nex_range(ii);
+    
+    params.crushFlag = 1;
+    
+    [~, Msig1, Msig2] = b1_afi.bloch_sim(params);
+    signal1_ideal_spoil(ii,:) = abs(Msig1);
+    signal2_ideal_spoil(ii,:) = abs(Msig2);
+    
+    params.inc = 39;
+    params.partialDephasing = 1;
+    params.partialDephasingFlag = 1;
+    params.crushFlag = 0;
+    
+    [~, Msig1, Msig2] = b1_afi.bloch_sim(params);
+    signal1_optimal_crush_and_rf_spoil(ii,:) = abs(Msig1);
+    signal2_optimal_crush_and_rf_spoil(ii,:) = abs(Msig2);
+    
+    params.inc = 0;
+    params.partialDephasing = 0;
+
+    [~, Msig1, Msig2] = b1_afi.bloch_sim(params);
+    signal1_no_gradient_and_rf_spoil(ii,:) = abs(Msig1);
+    signal2_no_gradient_and_rf_spoil(ii,:) = abs(Msig2);
+end
+%}
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+
+
In [11]:
+
+ +
+
+ + + +
+ +
+
+
In [12]:
+ +
+
+
# PYTHON CODE
+
+init_notebook_mode(connected=True)
+
+data1_1 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_ideal_spoil[ii]))),
+        name = 'Ideal Spoiling',
+        text = 'Ideal Spoiling',
+        hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]
+
+data1_1[10]['visible'] = True
+
+data2_1 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_optimal_crush_and_rf_spoil[ii]))),
+        name = 'Gradient & RF Spoiling',
+        text = 'Gradient & RF Spoiling',
+        hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]
+
+data2_1[10]['visible'] = True
+
+data3_1 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(signal1_no_gradient_and_rf_spoil[ii]))),
+        name = 'No Spoiling',
+        text = 'No Spoiling',
+        hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]
+
+data3_1[10]['visible'] = True
+
+data = data1_1 + data2_1+ data3_1
+
+steps = []
+for i in range(len(Nex_range)):
+    step = dict(
+        method = 'restyle',  
+        args = ['visible', [False] * len(data1_1)],
+        label = str(Nex_range[i])
+        )
+    step['args'][1][i] = True # Toggle i'th trace to "visible"
+    steps.append(step)
+
+sliders = [dict(
+    x = 0,
+    y = -0.02,
+    active = 10,
+    currentvalue = {"prefix": "n<sup>th</sup> TR1-TR2: <b>"},
+    pad = {"t": 50, "b": 10},
+    steps = steps
+)]
+
+layout = go.Layout(
+    width=580,
+    height=450,
+    margin=go.layout.Margin(
+        l=80,
+        r=40,
+        b=60,
+        t=10,
+    ),
+    annotations=[
+        dict(
+            x=0.5004254919715793,
+            y=-0.18,
+            showarrow=False,
+            text='Excitation Flip Angle (°)',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=-0.15,
+            y=0.5,
+            showarrow=False,
+            text='Signal',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            textangle=-90,
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis=dict(
+        autorange=False,
+        range=[0, params['EXC_FA'][-1]],
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    yaxis=dict(
+        autorange=True,
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    legend=dict(
+        x=0.5,
+        y=0.9,
+        traceorder='normal',
+        font=dict(
+            family='Times New Roman',
+            size=12,
+            color='#000'
+        ),
+        bordercolor='#000000',
+        borderwidth=2
+    ), 
+    sliders=sliders
+)
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-line', config = config)
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+

Data Fitting

+
+
+
+
+
+
+
+
+

+The ratio of Equations 1 and 2, gives rise to Equation 3 that depends on the parameters T1, TR1, TR2 and the excitation flip angle (θ). +

+ +

+

+

+ +

+Equation 3 can be simplified if the Taylor series expansion of the exponential function is used, followed by a first-order approximation to its terms. For this expansion, short TR1 and TR2 (TR1 < T1 and TR2 < T1) are assumed to approximate the signal intensities ratio (Equation 4) where n = TR2/TR1. +

+ +

+

+

+ +

+Finally, a measure of the actual flip-angle (θ) can be achieved by solving Equation 4 to obtain Equation 5, which only depends on the signal intensities ratio (r = S2/S1) and the parameters TR1 and TR2. +

+ +

+

+

+ +

+The actual flip-angle is estimated using an approximation (Equation 4) of a complete analytical solution (Equation 3), and the nature of this approximation makes it worthwhile to assess the accuracy of the signal intensities ratio between both equations. Next, a set of simulations are displayed to analyze how the choice of r is affected by T1, TR1 and TR2. First, the effect of the relaxation time T1 is simulated in Figure 5 for both the approximation and the complete analytical solution. +

+ +
+
+
+
+
+
+
+
+
+

+ +Figure 5. Effect of the relaxation time T1 on the ratio r. Signal intensities ratio is plotted as a function of the flip angle for the complete analytical solution (Equation 3 - blue) and the first-order approximation (Equation 4 - orange). AFI simulation details: TR1 = 20 ms, TR2 = 100 ms and variable T1. + +

+
+
+
+ +
+ +
+
+
In [13]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to generate the data required for Figure 5 of the blog post
+
+clear all
+
+%% Setup parameters
+% All times are in milliseconds
+% All flip angles are in degrees
+
+%% Calculate signals
+%
+% To see all the options available, run `help b1_afi.analytical_solution`
+
+T1_range = [10:10:200, 200:400:2000];
+
+params.TR1 = 20;
+params.TR2 = 100;
+params.EXC_FA = 1:90;
+n = params.TR2/params.TR1;
+
+for jj=1:length(T1_range)
+    params.T1 = T1_range(jj);
+    
+    % Range of T1
+    [Mz1, Mz2] = b1_afi.analytical_solution(params);
+    signal1_analytical = Mz1;
+    signal2_analytical = Mz2;
+
+    r_analytical(jj,:) = signal2_analytical./signal1_analytical;
+    r_approximation(jj,:) = (1 + n.*cosd(params.EXC_FA))./(n + cosd(params.EXC_FA));
+end
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+
+
In [14]:
+
+ +
+
+ + + +
+ +
+
+
In [15]:
+ +
+
+
# PYTHON CODE
+
+init_notebook_mode(connected=True)
+
+data1 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(r_analytical[ii]))),
+        name = 'Analytical',
+        text = 'Analytical',
+        hoverinfo = 'x+y+text') for ii in range(len(T1_range))]
+
+data1[10]['visible'] = True
+
+data2 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(r_approximation[ii]))),
+        name = 'Approximation',
+        text = 'Approximation',
+        hoverinfo = 'x+y+text') for ii in range(len(T1_range))]
+
+data2[10]['visible'] = True
+
+data = data1 + data2
+
+steps = []
+for i in range(len(T1_range)):
+    step = dict(
+        method = 'restyle',  
+        args = ['visible', [False] * len(data1)],
+        label = str(T1_range[i])
+        )
+    step['args'][1][i] = True # Toggle i'th trace to "visible"
+    steps.append(step)
+
+sliders = [dict(
+    x = 0,
+    y = -0.02,
+    active = 10,
+    currentvalue = {"prefix": "T1 value (ms): <b>"},
+    pad = {"t": 50, "b": 10},
+    steps = steps)]
+
+layout = go.Layout(
+    width=580,
+    height=450,
+    margin=go.layout.Margin(
+        l=80,
+        r=40,
+        b=60,
+        t=10,
+    ),
+    annotations=[
+        dict(
+            x=0.5004254919715793,
+            y=-0.18,
+            showarrow=False,
+            text='Excitation Flip Angle (°)',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=-0.15,
+            y=0.5,
+            showarrow=False,
+            text='r = S<sub>2</sub>/S<sub>1</sub>',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            textangle=-90,
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis=dict(
+        autorange=False,
+        range=[0, params['EXC_FA'][-1]],
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    yaxis=dict(
+        autorange=True,
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    legend=dict(
+        x=0.5,
+        y=0.9,
+        traceorder='normal',
+        font=dict(
+            family='Times New Roman',
+            size=12,
+            color='#000'
+        ),
+        bordercolor='#000000',
+        borderwidth=2
+    ), 
+    sliders=sliders
+)
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-line', config = config)
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+
+ +

+The signal ratio r is highly insensitive to the relaxation time T1, except for the low T1 values and large flip angles (>70°). This shows that the Taylor expansion is a good approximation to the signal ratio r because it is possible to get rid of the inverse quadratic T1 dependance by taking the first-order terms of the expansion. +

+ +

+The effect of the TR1 parameter on the signal ratio is shown in Figure 6. To assess the influence of the repetition time, we fix n=5 and vary the parameter TR1 in accordance to the relation n = TR2/TR1. As TR1 increases (> 50 ms), the approximated ratio r slightly deviates from the analytical approach. Although the deviation is slight only at high flip angles, a good signal ratio approximation can be achieved for a wide range of flip angles and repetition times. +

+ +
+
+
+
+
+
+
+
+
+

+ +Figure 6. Effect of the repetition time TR1 on the ratio r. Signal intensities ratio is plotted as a function of the flip angle for the complete analytical solution (Equation 3 - blue) and the first-order approximation (Equation 4 - orange). AFI simulation details: Variable TR1 ranging from 10 to 60 ms, fixed ratio n = 5 and T1 = 900 ms. + +

+
+
+
+ +
+ +
+
+
In [16]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to generate the data required for Figure 6 of the blog post
+
+clear all
+
+%% Setup parameters
+% All times are in milliseconds
+% All flip angles are in degrees
+
+%% Calculate signals
+%
+% To see all the options available, run `help b1_afi.analytical_solution`
+
+TR1_range = 10:5:60;
+
+params.T1 = 900;
+params.EXC_FA = 1:90;
+n = 5;
+
+for jj=1:length(TR1_range)
+    params.TR1 = TR1_range(jj);
+    params.TR2 = n*params.TR1;
+    
+    % Fixed: T1 = 900 ms, n = 5
+    [Mz1, Mz2] = b1_afi.analytical_solution(params);
+    signal1_analytical = Mz1;
+    signal2_analytical = Mz2;
+
+    r_analytical(jj,:) = signal2_analytical./signal1_analytical;
+    r_approximation(jj,:) = (1 + n.*cosd(params.EXC_FA))./(n + cosd(params.EXC_FA));
+end
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+
+
In [17]:
+
+ +
+
+ + + +
+ +
+
+
In [18]:
+ +
+
+
# PYTHON CODE
+
+init_notebook_mode(connected=True)
+
+data1 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(r_analytical[ii]))),
+        name = 'Analytical',
+        text = 'Analytical',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data1[2]['visible'] = True
+
+data2 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(r_approximation[ii]))),
+        name = 'Approximation',
+        text = 'Approximation',
+        hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]
+
+data2[2]['visible'] = True
+
+data = data1 + data2
+
+steps = []
+for i in range(len(TR1_range)):
+    step = dict(
+        method = 'restyle',  
+        args = ['visible', [False] * len(data1)],
+        label = str(TR1_range[i])
+        )
+    step['args'][1][i] = True # Toggle i'th trace to "visible"
+    steps.append(step)
+
+sliders = [dict(
+    x = 0,
+    y = -0.02,
+    active = 2,
+    currentvalue = {"prefix": "TR1 value (ms): <b>"},
+    pad = {"t": 50, "b": 10},
+    steps = steps)]
+
+layout = go.Layout(
+    width=580,
+    height=450,
+    margin=go.layout.Margin(
+        l=80,
+        r=40,
+        b=60,
+        t=10,
+    ),
+    annotations=[
+        dict(
+            x=0.5004254919715793,
+            y=-0.18,
+            showarrow=False,
+            text='Excitation Flip Angle (°)',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=-0.15,
+            y=0.5,
+            showarrow=False,
+            text='r = S<sub>2</sub>/S<sub>1</sub>',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            textangle=-90,
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis=dict(
+        autorange=False,
+        range=[0, params['EXC_FA'][-1]],
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    yaxis=dict(
+        autorange=True,
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    legend=dict(
+        x=0.5,
+        y=0.9,
+        traceorder='normal',
+        font=dict(
+            family='Times New Roman',
+            size=12,
+            color='#000'
+        ),
+        bordercolor='#000000',
+        borderwidth=2
+    ), 
+    sliders=sliders
+)
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-line', config = config)
+
+ +
+
+
+
+ +
+
+ + + + +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+
+ +

+Finally, the effect of the parameter n on the signal ratio r (Figure 7) does not seem to significantly affect the signal ratio between the approximated equation and the analytical approach. However, the parameter n has a major impact on the sensitivity of the AFI method to variations in the flip angle. Figure 7 shows that the increase of the parameter n (= TR2/TR1) allows for improvement of the dynamic range of flip angles measurements. These simulations have shown that an optimal implementation of the AFI method requires a careful selection of sequence parameters. +

+ +
+
+
+
+
+
+
+
+
+

+ +Figure 7. Effect of n (TR2 to TR1 ratio) on the ratio r. The signal intensities ratio is plotted as a function of the flip angle for the complete analytical solution (Equation 3 - blue) and the first-order approximation (Equation 4 - orange). AFI simulation details: Variable n ranging from 2 to 6, fixed TR1 = 20 ms and T1 = 900 ms. + +

+
+
+
+ +
+ +
+
+
In [19]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to generate the data required for Figure 7 of the blog post
+
+clear all
+
+%% Setup parameters
+% All times are in milliseconds
+% All flip angles are in degrees
+
+%% Calculate signals
+%
+% To see all the options available, run `help b1_afi.analytical_solution`
+
+n_range = 2:1:6;
+
+params.T1 = 900;
+params.TR1 = 20;
+params.EXC_FA = 1:90;
+
+for jj=1:length(n_range)
+    n = n_range(jj);
+    params.TR2 = n*params.TR1;
+    
+    % Fixed: T1 = 900 ms, TR1 = 20
+    [Mz1, Mz2] = b1_afi.analytical_solution(params);
+    signal1_analytical = Mz1;
+    signal2_analytical = Mz2;
+
+    r_analytical(jj,:) = signal2_analytical./signal1_analytical;
+    r_approximation(jj,:) = (1 + n.*cosd(params.EXC_FA))./(n + cosd(params.EXC_FA));
+end
+
+ +
+
+
+
+ +
+
+ + + +
+ +
+
+
In [20]:
+
+ +
+
+ + + +
+ +
+
+
In [21]:
+ +
+
+
# PYTHON CODE
+
+init_notebook_mode(connected=True)
+
+data1 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(r_analytical[ii]))),
+        name = 'Analytical',
+        text = 'Analytical',
+        hoverinfo = 'x+y+text') for ii in range(len(n_range))]
+
+data1[2]['visible'] = True
+
+data2 = [dict(
+        visible = False,
+        mode = 'lines',
+        x = params["EXC_FA"],
+        y = abs(np.squeeze(np.asarray(r_approximation[ii]))),
+        name = 'Approximation',
+        text = 'Approximation',
+        hoverinfo = 'x+y+text') for ii in range(len(n_range))]
+
+data2[2]['visible'] = True
+
+data = data1 + data2
+
+steps = []
+for i in range(len(n_range)):
+    step = dict(
+        method = 'restyle',  
+        args = ['visible', [False] * len(data1)],
+        label = str(n_range[i])
+        )
+    step['args'][1][i] = True # Toggle i'th trace to "visible"
+    steps.append(step)
+
+sliders = [dict(
+    x = 0,
+    y = -0.02,
+    active = 2,
+    currentvalue = {"prefix": "n value: <b>"},
+    pad = {"t": 50, "b": 10},
+    steps = steps)]
+
+layout = go.Layout(
+    width=580,
+    height=450,
+    margin=go.layout.Margin(
+        l=80,
+        r=40,
+        b=60,
+        t=10,
+    ),
+    annotations=[
+        dict(
+            x=0.5004254919715793,
+            y=-0.18,
+            showarrow=False,
+            text='Excitation Flip Angle (°)',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=-0.15,
+            y=0.5,
+            showarrow=False,
+            text='r = S<sub>2</sub>/S<sub>1</sub>',
+            font=dict(
+                family='Times New Roman',
+                size=22
+            ),
+            textangle=-90,
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis=dict(
+        autorange=False,
+        range=[0, params['EXC_FA'][-1]],
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    yaxis=dict(
+        autorange=False,
+        range=[0, 1],
+        showgrid=False,
+        linecolor='black',
+        linewidth=2
+    ),
+    legend=dict(
+        x=0.5,
+        y=0.9,
+        traceorder='normal',
+        font=dict(
+            family='Times New Roman',
+            size=12,
+            color='#000'
+        ),
+        bordercolor='#000000',
+        borderwidth=2
+    ), 
+    sliders=sliders
+)
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-line', config = config)
+
+ +
+
+
+
+ +
+
+ + + + +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+
+

+Figure 8 displays an example AFI dataset and its corresponding field B1 map in a healthy human brain. Although not clearly visible, both AFI images present a small Gibbs ringing artifact that is propagated and amplified due to the AFI calculation consisting of the division of both images (Boudreau et al. 2017). The ringing artifact is clearly seen in the unfiltered/raw B1 field map shown in Figure 8 (right). +

+
+
+
+
+
+
+
+
+
+

+ +Figure 8. Example actual flip-angle imaging dataset (left) and a resulting raw B1 map of a healthy adult brain (right). The relevant VFA protocol parameters used were: TR1 = 20 ms, TR2 = 100 ms and θnominal = 60°. The B1 map (right) was fitted using the approximate r ratio (Equation 5). + +

+
+
+
+ +
+ +
+
+
In [22]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Download Actual Flip-Angle Imaging brain MRI data for Figure 8 of the blog post
+
+cmd = ['curl -L -o b1_afi.zip https://osf.io/csjgx/download?version=6'];
+[STATUS,MESSAGE] = unix(cmd);
+unzip('b1_afi.zip');
+
+ +
+
+
+
+ +
+
+ +
+ +
+ +
+ + +
+
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
+                                 Dload  Upload   Total   Spent    Left  Speed
+100   459  100   459    0     0    739      0 --:--:-- --:--:-- --:--:--   737
+100 5281k  100 5281k    0     0  1620k      0  0:00:03  0:00:03 --:--:-- 4028k
+
+
+
+ +
+
+
+ +
+
+ + + +
+ +
+
+
In [23]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to generate the data required for Figure 8 of the blog post
+
+clear all
+
+% Format qMRLab b1_afi model parameters, and load them into the Model object
+Model = b1_afi; 
+Model.Prot.Sequence.Mat = [60, 20, 100];
+
+% Format data structure so that they may be fit by the model
+data = struct();
+data.AFIData1 = load_nii_data('b1_afi/AFIData1.nii.gz');
+data.AFIData2 = load_nii_data('b1_afi/AFIData2.nii.gz');
+
+FitResults = FitData(data,Model,0); % The '0' flag is so that no wait bar is shown.
+
+ +
+
+
+
+ +
+
+ +
+ +
+ +
+ + +
+
ans = 1
+ans = 1
+warning: load: '/home/jovyan/work/t1_notebooks/qMRLab/Deploy/Documentation/iqmr_gethelp.mat' found by searching load path
+warning: called from
+    qMRinfo at line 15 column 8
+    get_MRIinputs_optional at line 130 column 20
+    sanityCheck at line 73 column 27
+    FitData at line 35 column 1
+=============== qMRLab::Fit ======================
+Operation has been started: b1_afi
+Elapsed time is 0.0901079 seconds.
+Operation has been completed: b1_afi
+==================================================
+
+
+
+ +
+
+
+ +
+
+ + + +
+ +
+
+
In [24]:
+ +
+
+
%% MATLAB/OCTAVE CODE
+% Code used to re-orient the images to make pretty figures, and to assign variables with the axis lengths.
+
+B1map_raw = imrotate(squeeze(FitResults.B1map_raw(:,:,14)),-90);
+
+xAxis = [0:size(B1map_raw,2)-1];
+yAxis = [0:size(B1map_raw,1)-1];
+
+% Raw MRI AFI data
+AFIData1 = imrotate(squeeze(data.AFIData1(:,:,14)),-90);
+AFIData2 = imrotate(squeeze(data.AFIData2(:,:,14)),-90);
+
+% Load mask
+mask = load_nii_data('../afi_blog/fitResults/FitResults_median_filter/mask.nii.gz');
+mask = imrotate(squeeze(mask(:,:,14)),-90);
+
+ +
+
+
+
+ +
+
+ +
+ +
+ +
+ + +
+
ans = 1
+
+
+
+ +
+
+
+ +
+
+ + + +
+ +
+
+
In [25]:
+
+ +
+
+ + + +
+ +
+
+
In [26]:
+ +
+
+
from plotly import tools
+
+# Masking B1 map
+B1map_raw = np.asarray(B1map_raw)
+mask = np.asarray(mask)
+B1map_raw_masked = B1map_raw*mask
+B1map_raw_masked[np.isnan(B1map_raw_masked)] = 0
+
+
+trace1 = go.Heatmap(x = xAxis,
+                   y = yAxis,
+                   z=AFIData1,
+                   colorscale='Greys',
+                   showscale = False,
+                   visible=False,
+                   name = 'Signal1')
+trace2 = go.Heatmap(x = xAxis,
+                   y = yAxis,
+                   z=AFIData2,
+                   colorscale='Greys',
+                   showscale = False,
+                   visible=True,
+                   name = 'Signal2')
+trace3 = go.Heatmap(x = xAxis,
+                   y = yAxis,
+                   z=B1map_raw_masked,
+                   zmin=0.7,
+                   zmax=1.3,
+                   colorscale='RdBu',
+                   xaxis='x2',
+                   yaxis='y2',
+                   visible=True,
+                   name = 'B1 values')
+
+data=[trace1, trace2, trace3]
+
+
+updatemenus = list([
+    dict(active=1,
+         x = 0.09,
+         xanchor = 'left',
+         y = -0.15,
+         yanchor = 'bottom',
+         direction = 'up',
+         font=dict(
+                family='Times New Roman',
+                size=16
+            ),
+         buttons=list([   
+            dict(label = 'Signal 1',
+                 method = 'update',
+                 args = [{'visible': [True, False, True]},
+                         ]),
+            dict(label = 'Signal 2',
+                 method = 'update',
+                 args = [{'visible': [False, True, True]},
+                           ]),
+        ])
+    )
+])
+
+layout = dict(
+    width=560,
+    height=345,
+    margin = dict(
+                t=40,
+                r=50,
+                b=10,
+                l=50),
+    annotations=[
+        dict(
+            x=0.07,
+            y=1.15,
+            showarrow=False,
+            text='Input Data',
+            font=dict(
+                family='Times New Roman',
+                size=26
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=0.60,
+            y=1.15,
+            showarrow=False,
+            text='Raw B<sub>1</sub> map',
+            font=dict(
+                family='Times New Roman',
+                size=26
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=1.12,
+            y=1.15,
+            showarrow=False,
+            text='B<sub>1</sub>',
+            font=dict(
+                family='Times New Roman',
+                size=26
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0, 0.58]),
+    yaxis = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0, 1]),
+    xaxis2 = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0.40, 0.98]),
+    yaxis2 = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0, 1], anchor='x2'),
+    showlegend = False,
+    autosize = False,
+    updatemenus=updatemenus
+)
+
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-heatmap', config = config)
+
+ +
+
+
+
+ +
+
+ + +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+
+

+The ringing artifact shown in Figure 8 can be attenuated by implementing a smoothing process. Figure 9 shows the raw (left) and the filtered (right) B1 map where a median filter was used to smooth the field map. +

+
+
+
+
+
+
+
+
+
+

+ +Figure 9. Raw (left) and filtered (right) B1 map. A median filter of size 7x7x7 pixels was used to attenuate the Gibbs ringing artifact. + +

+
+
+
+ +
+ +
+
+
In [27]:
+ +
+
+
B1map_filtered = load_nii_data('../afi_blog/fitResults/FitResults_median_filter/B1map_filtered.nii.gz');
+B1map_filtered = imrotate(squeeze(B1map_filtered(:,:,14)),-90);
+
+ +
+
+
+
+ +
+
+ +
+ +
+ +
+ + +
+
ans = 1
+
+
+
+ +
+
+
+ +
+
+ + + +
+ +
+
+
In [28]:
+
+ +
+
+ + + +
+ +
+
+
In [29]:
+ +
+
+
from plotly import tools
+
+# Masking B1 map
+B1map_filtered = np.asarray(B1map_filtered)
+mask = np.asarray(mask)
+B1map_filtered_masked = B1map_filtered*mask
+B1map_filtered_masked[np.isnan(B1map_filtered_masked)] = 0
+
+trace1 = go.Heatmap(x = xAxis,
+                   y = yAxis,
+                   z=B1map_raw_masked,
+                   zmin=0.7,
+                   zmax=1.3,
+                   colorscale='RdBu',
+                   showscale = False,
+                   visible=True,
+                   name = 'B1 values')
+trace2 = go.Heatmap(x = xAxis,
+                   y = yAxis,
+                   z=B1map_filtered_masked,
+                   zmin=0.7,
+                   zmax=1.3,
+                   colorscale='RdBu',
+                   xaxis='x2',
+                   yaxis='y2',
+                   visible=True,
+                   name = 'B1 values')
+
+data=[trace1, trace2]
+
+layout = dict(
+    width=560,
+    height=310,
+    margin = dict(
+                t=40,
+                r=50,
+                b=10,
+                l=50),
+    annotations=[
+        dict(
+            x=0.04,
+            y=1.15,
+            showarrow=False,
+            text='Raw B<sub>1</sub> map',
+            font=dict(
+                family='Times New Roman',
+                size=26
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=0.60,
+            y=1.15,
+            showarrow=False,
+            text='Filtered B<sub>1</sub> map',
+            font=dict(
+                family='Times New Roman',
+                size=26
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+        dict(
+            x=1.12,
+            y=1.15,
+            showarrow=False,
+            text='B<sub>1</sub>',
+            font=dict(
+                family='Times New Roman',
+                size=26
+            ),
+            xref='paper',
+            yref='paper'
+        ),
+    ],
+    xaxis = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0, 0.58]),
+    yaxis = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0, 1]),
+    xaxis2 = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0.40, 0.98]),
+    yaxis2 = dict(range = [0,127], autorange = False,
+             showgrid = False, zeroline = False, showticklabels = False,
+             ticks = '', domain=[0, 1], anchor='x2'),
+    showlegend = False,
+    autosize = False,
+)
+
+
+fig = dict(data=data, layout=layout)
+
+iplot(fig, filename = 'basic-heatmap', config = config)
+
+ +
+
+
+
+ +
+
+ + +
+ +
+ + + +
+
+ + +
+ +
+
+ +
+ +
+
+ +
+
+ + +
+
+
+
+

Benefits and Pitfalls

+
+
+
+
+
+
+
+
+

+B1 mapping is of interest for diverse MRI applications, and several mapping techniques have been developed. The DAM method consists of acquiring two scans at two different flip angles. To avoid the dependence of the signal on T1, long repetition times are required to allow the recovery of the longitudinal magnetization between pulses (Yarnykh 2007; Insko & Bolinger 1993). The AFI method overcomes this practical limitation by repeating the pulse sequence at a fast rate to achieve a pulsed state of magnetization and shorter time delays between pulses. In addition, due to scan-time constraints, B1 mapping methods are often implemented in 2D (Chavez & Stanisz 2011). However, the accuracy of the measurements of 2D B1 mapping techniques is compromised by the slice profile effects due to the problem of nonuniform excitation across slices (Yarnykh 2007; Chavez & Stanisz 2011). The AFI method on the other hand, adresses this issue using a fast 3D implementation leading to scans with an excellent anatomical coverage in clinically feasible times, with an increase in signal-to-noise ratio compared to 2D multislice acquisitions. +

+ +

+The performance of the AFI method is based on the following assumptions. First, the two images acquired at different times should be registered to avoid motion effects. It is also assumed that the signal is insensitive to the main magnetic field non-uniformities and chemical shift effects that are canceled out when taking the signal ratio r (Yarnykh 2007). Despite some clear advantages over other B1 mapping techniques, the application of spoiler gradients to mitigate the T1 dependence can be a limitation due to significant levels of RF power depositions (Sacolick et al. 2010). Furthermore, it is necessary to adapt the AFI pulse sequence to different scanner manufacturers, and programming experience is required to bring this technique into the clinic. +

+
+
+
+
+
+
+
+
+

Works Cited

+
+
+
+
+
+
+
+
+

+Bernstein, M., King, K. & Zhou, X., 2004. Handbook of MRI Pulse Sequences, Elsevier. +

+ +

+Boudreau, M., Tardif, C.L., Stikov, N., Sled, J.G., Lee, W. & Pike, G.B., 2017. B1 mapping for bias-correction in quantitative T1 imaging of the brain at 3T using standard pulse sequences. J. Magn. Reson. Imaging, 46, pp.1673-1682. +

+ +

+Boudreau, M., Stikov, N. & Pike, G.B, 2018. B1-sensitivity analysis of quantitative magnetization transfer imaging. Magn. Reson. Med., 79, pp.276-285. +

+ +

+Chavez, S. & Stanisz, G.J., 2012. A novel method for simultaneous 3D B1 and T1 mapping: the method of slopes (MoS). NMR Biomed., 25, pp.1043-1055. +

+ +

+Deoni, S.C., 2007. High-resolution T1 mapping of the brain at 3T with driven equilibrium single pulse observation of T1 with high-speed incorporation of RF field inhomogeneities (DESPOT1-HIFI). J. Magn. Reson. Imaging, 26, pp.1106-1111. +

+ +

+Ibrahim, T.S., Lee, R., Baertlein, B.A. & Robitaille, P.M., 2001. B1 field homogeneity and SAR calculations for the birdcage coil. Phys Med Biol., 46(2), pp.609-619. +

+ +

+Insko, E. & Bolinger, L., 1993. Mapping of the radio frequency field. J. Magn. Reson. Ser. A, 103(1), pp.82–85. +

+ +

+Katscher, U., Voigt, T., Findeklee, C., Vernickel, P., Nehrke, K. & Dössel, O., 2009. Determination of electric conductivity and local SAR via B1 mapping. IEEE Trans Med Imaging., 28(9), pp.1365-1374. +

+ +

+Ropele, S., Filippi, M., Valsasina, P., Korteweg, T., Barkhof, F., Tofts, P.S., Samson, R., Miller, D.H. & Fazekas, F., 2005. Assessment and correction of B1-induced errors in magnetization transfer ratio measurements. Magn. Reson. Med., 53, pp.134-140. +

+ +

+Sacolick, L.I., Wiesinger, F., Hancu, I. & Vogel, M.W., 2010. B1 mapping by Bloch-Siegert shift. Magn. Reson. Med., 63, pp.1315-1322. +

+ +

+Sled, J.G. & Pike, G.B., 1998. Standing-wave and RF penetration artifacts caused by elliptic geometry: an electrodynamic analysis of MRI. IEEE Trans. Med. Imaging., 17(4), pp.653-662. +

+ +

+Sled, J.G. & Pike, G.B., 2000. Correction for B1 and B0 variations in quantitative T2 measurements using MRI. Magn. Reson. Med., 43(4), pp.589–593. +

+ +

+Stollberger, R. & Wach, P., 1996. Imaging of the active B1 field in vivo. Magn. Reson. Med., 35(2), pp.246–251. +

+ +

+van den Bergen, B., Van den Berg, C.A., Bartels, L.W. & Lagendijk, J.J., 2007. 7 T body MRI: B1 shimming with simultaneous SAR reduction. Phys Med Biol., 52(17), pp.5429-5441. +

+ +

+Yarnykh, V.L. & Yuan, C. Actual flip angle imaging in the pulsed steady state. In: Proceedings of the 12th Annual Meeting of ISMRM, Kyoto, Japan, 2004 (Abstract 194). +

+ +

+Yarnykh, V.L., 2007. Actual Flip-Angle Imaging in the pulsed steady state: A method for rapid three-dimensional mapping of the transmitted radiofrequency field. Magn. Reson. Med., 57(1), pp.192-200. +

+ +

+Yarnykh, V.L., 2010. Optimal radiofrequency and gradient spoiling for improved accuracy of T1 and B1 measurements using fast steady-state techniques. Magn. Reson. Med., 63(6), pp.1610–1626. +

+ +

+Zur, Y., Wood, M.L. & Neuringer, L.J., 1991. Spoiling of transverse magnetization in steady-state sequences. Magn. Reson. Med., 21(2), pp.251–263. +

+ +
+
+
+
+ +
+ +
+
+
In [30]:
+ +
+
+
# PYTHON CODE
+
+display(HTML(
+    '<style type="text/css">'
+    '.output_subarea {'
+        'display: block;'
+        'margin-left: auto;'
+        'margin-right: auto;'
+    '}'
+    '.blog_body {'
+        'line-height: 2;'
+        'font-family: timesnewroman;'
+        'font-size: 18px;'
+        'margin-left: 0px;'
+        'margin-right: 0px;'
+    '}'
+    '.biblio_body {'
+        'line-height: 1.5;'
+        'font-family: timesnewroman;'
+        'font-size: 18px;'
+        'margin-left: 0px;'
+        'margin-right: 0px;'
+    '}'
+    '.note_body {'
+        'line-height: 1.25;'
+        'font-family: timesnewroman;'
+        'font-size: 18px;'
+        'margin-left: 0px;'
+        'margin-right: 0px;'
+        'color: #696969'
+    '}'
+    '.figure_caption {'
+        'line-height: 1.5;'
+        'font-family: timesnewroman;'
+        'font-size: 16px;'
+        'margin-left: 0px;'
+        'margin-right: 0px'
+    '</style>'
+))
+
+ +
+
+
+
+ +
+
+ +
+ +
+ +
+ + + +
+ +
+ +
+ +
+
+
+ +
+
+ + +
+ + + + + + + + + + + + + + + + + diff --git a/afi_blog/ActualFlipAngleImaging.ipynb b/afi_blog/ActualFlipAngleImaging.ipynb new file mode 100644 index 0000000..3437852 --- /dev/null +++ b/afi_blog/ActualFlipAngleImaging.ipynb @@ -0,0 +1,2293 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "kernel": "SoS", + "tags": [ + "scratch" + ] + }, + "source": [ + "# Welcome to a qMRLab interactive blog post Jupyter Notebook!\n", + "\n", + "If this is your first time running a Juptyer Notebook, there's a lot of tutorials available online to help. [Here's one](https://www.dataquest.io/blog/jupyter-notebook-tutorial/) for your convenience.\n", + "\n", + "## Introduction\n", + "\n", + "This notebook contains everything needed to reproduce the Actual Flip-Angle B1 blog post on the [qMRLab website](). In fact, this notebook generated the HTML for the blog post too! This notebook is currently running on a MyBinder server that only you can access, but if you want to be kept up-to-date on any changes that the developpers make to this notebook, you should go to it's [GitHub repository](https://github.com/qMRLab/t1_notebooks) and follow it by clicking the \"Watch\" button in the top right (you may need to create a GitHub account, if you don't have one already).\n", + "\n", + "## Tips\n", + "\n", + "Here's a few things you can do in this notebook\n", + "\n", + "### Code\n", + "* Run the entire processing by clicking above on the \"Kernel\" tab, then \"Restart & Run All\". It will be complete when none of the cells have an asterix \"\\*\" in the square brackets.\n", + "* To change the code, you need to click once on code cells. To re-run that cell, click the \"Run\" button above when the cell is selected.\n", + " * **Note:** Cells can depend on previous cells, or even on previous runs of the cell itself, so it's best to run all the previous cells beforehand.\n", + "* This binder runs on SoS, which allows the mixing of Octave (i.e. an open-source MATLAB) and Python cells. Take a look a the drop down menu on the top right of the cells to know which one you are running.\n", + "* To transfer data from cells of one language to another, you need to create a new cell in the incoming language and run `%get (param name) --from (outgoing language)`. See cells below for several examples within this notebook.\n", + "\n", + "### HTML\n", + "* To reproduce the HTML of the blog post, run the entire processing pipeline (see point one in the previous section), then save the notebook (save icon, top left). Now, click on the drop down menu on the left pannel, and select `%sossave --to html --force` . After a few seconds, it should output \"Workflow saved to ActualFlipAngleImaging.html\" – click on the HTML name, and you're done!\n", + "* Cells with tags called \"scratch\" are not displayed in the generated HTML.\n", + "* Cells with the tag \"report_output\" display the output (e.g. figures) in the generated HTML.\n", + "* Currently in an un-run notebook, the HTML is not formatted like the website. To do so, run the Python module import cell (`# Module imports`) and then very last cell (`display(HTML(...)`).\n", + "\n", + "**If you have any other questions or comments, please raise them in a [GitHub issue](https://github.com/qMRLab/t1_notebooks/issues).**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "SoS", + "tags": [ + "scratch" + ] + }, + "source": [ + "# Note\n", + "\n", + "The following cell is meant to be displayed for instructional purposes in the blog post HTML when \"All cells\" gets displayed (i.e. the Octave code)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "% **Blog post code introduction**\n", + "% \n", + "% Congrats on activating the \"All cells\" option in this interactive blog post =D\n", + "%\n", + "% Below, several new HTML blocks have appears prior to the figures, displaying the Octave/MATLAB code that was used to generate the figures in this blog post.\n", + "%\n", + "% If you want to reproduce the data on your own local computer, you simply need to have qMRLab installed in your Octave/MATLAB path and run the \"startup.m\" file, as is shown below.\n", + "%\n", + "% If you want to get under the hood and modify the code right now, you can do so in the Jupyter Notebook of this blog post hosted on MyBinder. The link to it is in the introduction above." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "# Module imports\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import plotly.plotly as py\n", + "import plotly.graph_objs as go\n", + "import numpy as np\n", + "from plotly import __version__\n", + "from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n", + "config={'showLink': False, 'displayModeBar': False}\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "from IPython.core.display import display, HTML" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "SoS" + }, + "source": [ + "

Actual flip-angle imaging B1 mapping

\n", + "

" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "SoS" + }, + "source": [ + "

\n", + "

\n", + "Transmit radiofrequency field maps (B1+, or B1 for short) are used in diverse applications in MRI including: the study of electrical properties in tissues in vivo (Sled & Pike 1998; Katscher et al. 2009), specific absorption rate (SAR) calculations (Ibrahim et al. 2001), the calibration of quantitative T1 (Deoni 2007; Boudreau et al. 2017) and T2 (Sled and Pike 2000) maps, better parameter estimation from magnetization transfer measurements (Ropele et al. 2005; Boudreau et al. 2018), B1 shimming to improve image quality at whole-body ultra high fields (van den Bergen et al. 2007), or quality control of RF coils (Yarnykh 2007). Several B1 mapping techniques have been developed, and they can be broadly divided as magnitude-based and phase-based methods. The double angle method (DAM) is a saturation-recovery magnitude-based method that takes the ratio of the signal intensity of two magnitude images measured with different excitation flip angles (Insko & Bolinger 1993; Stollberger and Wach 1996). The Bloch-Siegert shift technique is a rapid phase-based method that encodes the B1 information into phase signal (Sacolick et al. 2010). The actual flip-angle imaging (AFI) is a magnitude-based B1 mapping method that consists of a 3D acquisition that benefits from good anatomical coverage. In addition, this technique allows the acquisitions of whole-body (~7 min) and brain (~3 min) B1 maps leading to a feasible implementation in clinics (Yarnykh 2004; Yarnykh 2007). On the other hand, the AFI pulse sequence has certain constraints that need to be considered for this B1 mapping method to be widely deployed. Some of the limitations include the use of spoiler gradients that can give rise to prohibitive SAR values (Sacolick et al. 2010), and the pulse sequence modifications on the MRI machine to implement the AFI method.\n", + "

\n", + "\n", + "

\n", + "In this blog post, we will focus on presenting details about the AFI B1 mapping method. We will cover signal modeling, data fitting, the benefits and the pitfalls of the technique. The figures are generated using the qMRLab module for this method. \n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "SoS" + }, + "source": [ + "

Signal Modelling

" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "SoS" + }, + "source": [ + "
\n", + "

\n", + "The pulse sequence of the AFI method (Figure 1) is composed of two identical RF pulses and two different delays (TR1 < TR2). After each RF pulse, the signal intensity is acquired followed by a spoiler to destroy the residual transverse magnetization next to the following RF pulse. This method implements a pulsed steady-state signal with a gradient-echo acquisition, thus preventing the use of long repetition times (Yarnykh 2007). It has been demonstrated that if the delays TR1 and TR2 are sufficiently short (e.g. TR1/TR2 = 20 ms/100 ms), and the transverse magnetization is completely spoiled, the ratio of signal intensities (r = S2/S1) depends on the flip angle of applied pulses and is highly insensitive to T1 (Yarnykh 2007). \n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "SoS" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 1. Simplified pulse sequence diagram of an actual flip-angle imaging (AFI) pulse sequence with a gradient echo readout. TR1: repetition time 1, TR2: repetition time 2, θ: excitation flip angle for the measurement, IMG: image acquisition (k-space readout), SPOIL: spoiler gradient.\n", + "\n", + "

\n", + "
\n", + "\n", + "

\n", + "

" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "The magnetization of an AFI experiment can be modeled under steady-state conditions by the implementation of a fast repetition of the sequence (TR1 < TR2 < T1). The solution of the Bloch equations for the AFI method is given by Equations 1 and 2 that represent the longitudinal magnetization before the application of the RF pulses:\n", + "

\n", + "\n", + "

\n", + "

\n", + "

\n", + " \n", + "

\n", + "

\n", + "

\n", + "\n", + "

\n", + "Mz1,2 is the longitudinal magnetization of both pulses, M0 is the magnetization at thermal equilibrium, TR1 is the delay time after the first pulse, TR2 is the delay time after the second identical pulse (Figure 1), and θ is the excitation flip angle. The steady-state longitudinal magnetization Mz curves for different T1 values for a range of θn and TR values are shown in Figure 2.\n", + "

\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 2. Longitudinal magnetization before the first radiofrequency pulse (Equation 1, solid lines) and before the second identical pulse (Equation 2, dashed lines) for three different T1 values.\n", + "\n", + "

\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Adds qMRLab to the path of the environment\n", + "\n", + "cd ../qMRLab\n", + "startup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to generate the data required for Figure 2 of the blog post\n", + "\n", + "clear all\n", + "\n", + "%% Setup parameters\n", + "% All times are in milliseconds\n", + "% All flip angles are in degrees\n", + "\n", + "TR1_range = 5:5:50;\n", + "n = 5;\n", + "TR2_range = n*TR1_range;\n", + "\n", + "params.EXC_FA = 1:90;\n", + "\n", + "%% Calculate signals\n", + "%\n", + "% To see all the options available, run `help b1_afi.analytical_solution`\n", + "\n", + "for ii = 1:length(TR1_range)\n", + " params.TR1 = TR1_range(ii);\n", + " params.TR2 = TR2_range(ii);\n", + " \n", + " % White matter\n", + " params.T1 = 900; % in milliseconds\n", + " [Mz1, Mz2] = b1_afi.analytical_solution(params);\n", + " signal1_WM(ii,:) = Mz1;\n", + " signal2_WM(ii,:) = Mz2;\n", + "\n", + " % Grey matter\n", + " params.T1 = 1500; % in milliseconds\n", + " [Mz1, Mz2] = b1_afi.analytical_solution(params);\n", + " signal1_GM(ii,:) = Mz1;\n", + " signal2_GM(ii,:) = Mz2;\n", + "\n", + " % CSF\n", + " params.T1 = 4000; % in milliseconds\n", + " [Mz1, Mz2] = b1_afi.analytical_solution(params);\n", + " signal1_CSF(ii,:) = Mz1;\n", + " signal2_CSF(ii,:) = Mz2;\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get params --from Octave\n", + "%get TR1_range --from Octave\n", + "%get TR2_range --from Octave\n", + "%get signal1_WM --from Octave\n", + "%get signal2_WM --from Octave\n", + "%get signal1_GM --from Octave\n", + "%get signal2_GM --from Octave\n", + "%get signal1_CSF --from Octave\n", + "%get signal2_CSF --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "data1_1 = [dict(\n", + " visible = False,\n", + " line=dict(color='royalblue'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_WM[ii]))),\n", + " name = 'S1: T1 = 0.9 s (White Matter)',\n", + " text = 'S1: T1 = 0.9 s (White Matter)',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data1_1[4]['visible'] = True\n", + "\n", + "data1_2 = [dict(\n", + " visible = False,\n", + " line=dict(color='royalblue', dash='dash'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal2_WM[ii]))),\n", + " name = 'S2: T1 = 0.9 s (White Matter)',\n", + " text = 'S2: T1 = 0.9 s (White Matter)',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data1_2[4]['visible'] = True\n", + "\n", + "data2_1 = [dict(\n", + " visible = False,\n", + " line=dict(color='firebrick'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_GM[ii]))),\n", + " name = 'S1: T1 = 1.5 s (Grey Matter)',\n", + " text = 'S1: T1 = 1.5 s (Grey Matter)',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data2_1[4]['visible'] = True\n", + "\n", + "data2_2 = [dict(\n", + " visible = False,\n", + " line=dict(color='firebrick', dash='dash'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal2_GM[ii]))),\n", + " name = 'S2: T1 = 1.5 s (Grey Matter)',\n", + " text = 'S2: T1 = 1.5 s (Grey Matter)',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data2_2[4]['visible'] = True\n", + "\n", + "data3_1 = [dict(\n", + " visible = False,\n", + " line=dict(color='orange'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_CSF[ii]))),\n", + " name = 'S1: T1 = 4.0 s (Cerebrospinal Fluid)',\n", + " text = 'S1: T1 = 4.0 s (Cerebrospinal Fluid)',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data3_1[4]['visible'] = True\n", + "\n", + "data3_2 = [dict(\n", + " visible = False,\n", + " line=dict(color='orange', dash='dash'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal2_CSF[ii]))),\n", + " name = 'S2: T1 = 4.0 s (Cerebrospinal Fluid)',\n", + " text = 'S2: T1 = 4.0 s (Cerebrospinal Fluid)',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data3_2[4]['visible'] = True\n", + "\n", + "data = data1_1 + data1_2 + data2_1 +data2_2 + data3_1 + data3_2\n", + "\n", + "steps = []\n", + "for i in range(len(TR1_range)):\n", + " step = dict(\n", + " method = 'restyle', \n", + " args = ['visible', [False] * len(data1_1)],\n", + " label = str(TR1_range[i])\n", + " )\n", + " step['args'][1][i] = True # Toggle i'th trace to \"visible\"\n", + " steps.append(step)\n", + "\n", + "sliders=[\n", + " dict(\n", + " x = 0,\n", + " y = -0.09,\n", + " active = 3,\n", + " currentvalue = {\"prefix\": \"TR1 value (ms): \"},\n", + " pad = {\"t\": 50, \"b\": 10},\n", + " steps = steps)]\n", + "\n", + "layout = go.Layout(\n", + " width=650,\n", + " height=520,\n", + " margin=go.layout.Margin(\n", + " l=80,\n", + " r=40,\n", + " b=60,\n", + " t=10,\n", + " ),\n", + " annotations=[\n", + " dict(\n", + " x=0.5004254919715793,\n", + " y=-0.18,\n", + " showarrow=False,\n", + " text='Excitation Flip Angle (°)',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=-0.15,\n", + " y=0.5,\n", + " showarrow=False,\n", + " text='Long. Magnetization (Mz)',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " textangle=-90,\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=-0.007,\n", + " y=-0.25,\n", + " showarrow=False,\n", + " text='TR2 = 5TR1',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis=dict(\n", + " autorange=False,\n", + " range=[0, params['EXC_FA'][-1]],\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " yaxis=dict(\n", + " autorange=True,\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " legend=dict(\n", + " x=0.5,\n", + " y=0.9,\n", + " traceorder='normal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12,\n", + " color='#000'\n", + " ),\n", + " bordercolor='#000000',\n", + " borderwidth=2\n", + " ), \n", + " sliders=sliders\n", + ")\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-line', config = config)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "\n", + "

\n", + "The analytical solution of the Bloch equations in a steady-state experiment (Equation 1 and Equation 2) makes several assumptions leading to practical challenges. First, it is assumed that the longitudinal magnetization has reached a steady state after a sufficiently large number of repetition times (TR), and that the transverse magnetization is perfectly spoiled prior to each pulse. To explore these properties, a numerical approach known as Bloch simulations is used to estimate the signal from an MRI experiment given a set of sequence parameters. Here, the Bloch simulations allow us to estimate the magnetization using a different number of sequence repetitions, and look at a special case when the steady-state is not achieved (due to a small number of sequence repetitions). As can be seen in Figure 3, the number of repetitions required to reach a steady-state depends on T1 and the flip angle.\n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 3. Signal 1 (blue) and Signal 2 (red) curves simulated using Bloch simulations (solid lines) for a number of repetitions ranging from 1 to 150, plotted against the ideal case (Equations 1 and 2 – dashed lines). Simulation details: TR1 = 20 ms, TR2 = 100 ms, T1 = 900 ms, 100 spins. Ideal spoiling was used for this set of Bloch simulations (transverse magnetization was set to 0 at the end of each TR1,2).\n", + "\n", + "

\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to generate the data required for Figure 3 of the blog post\n", + "\n", + "clear all\n", + "\n", + "%% Setup parameters\n", + "% All times are in milliseconds\n", + "% All flip angles are in degrees\n", + "\n", + "% White matter\n", + "params.T1 = 900; % in milliseconds\n", + "params.T2 = 10000;\n", + "params.TR1 = 20;\n", + "params.TR2 = 100;\n", + "params.TE = 5;\n", + "params.EXC_FA = 1:90;\n", + "Nex_range = 1:1:150;\n", + "\n", + "%% Calculate signals\n", + "% This results were precomputed. Uncomment the lines below to run the simulation.\n", + "signal1_bloch = load('../afi_blog/blochsim/signal1_blochsim.mat');\n", + "signal2_bloch = load('../afi_blog/blochsim/signal2_blochsim.mat');\n", + "signal1_analytical = load('../afi_blog/blochsim/signal1_analytical.mat');\n", + "signal2_analytical = load('../afi_blog/blochsim/signal2_analytical.mat');\n", + "signal1_blochsim = signal1_bloch.signal1_blochsim;\n", + "signal2_blochsim = signal2_bloch.signal2_blochsim;\n", + "signal1_analytical = signal1_analytical.signal1_analytical;\n", + "signal2_analytical = signal2_analytical.signal2_analytical;\n", + "\n", + "%{\n", + "% Bloch simulations\n", + "% To see all the options available, run `help b1_afi.analytical_solution`\n", + "\n", + "for ii = 1:length(Nex_range)\n", + " params.Nex = Nex_range(ii);\n", + " \n", + " [Mz1, Mz2] = b1_afi.analytical_solution(params);\n", + " signal1_analytical(ii,:) = Mz1;\n", + " signal2_analytical(ii,:) = Mz2;\n", + "\n", + " [~, Msig1, Msig2] = b1_afi.bloch_sim(params);\n", + " signal1_blochsim(ii,:) = abs(complex(Msig1));\n", + " signal2_blochsim(ii,:) = abs(complex(Msig2));\n", + "end\n", + "%}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get params --from Octave\n", + "%get Nex_range --from Octave\n", + "%get signal1_analytical --from Octave\n", + "%get signal2_analytical --from Octave\n", + "%get signal1_blochsim --from Octave\n", + "%get signal2_blochsim --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "data1_1 = [dict(\n", + " visible = False,\n", + " line=dict(color='royalblue', dash='dash'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_analytical[ii]))),\n", + " name = 'S1: Analytical Solution',\n", + " text = 'S1: Analytical Solution',\n", + " hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]\n", + "\n", + "data1_1[7]['visible'] = True\n", + "\n", + "data1_2 = [dict(\n", + " visible = False,\n", + " line=dict(color='firebrick', dash='dash'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal2_analytical[ii]))),\n", + " name = 'S2: Analytical Solution',\n", + " text = 'S2: Analytical Solution',\n", + " hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]\n", + "\n", + "data1_2[7]['visible'] = True\n", + "\n", + "data2_1 = [dict(\n", + " visible = False,\n", + " line=dict(color='royalblue'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_blochsim[ii]))),\n", + " name = 'S1: Bloch Simulation',\n", + " text = 'S1: Bloch Simulation',\n", + " hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]\n", + "\n", + "data2_1[7]['visible'] = True\n", + "\n", + "data2_2 = [dict(\n", + " visible = False,\n", + " line=dict(color='firebrick'),\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal2_blochsim[ii]))),\n", + " name = 'S2: Bloch Simulation',\n", + " text = 'S2: Bloch Simulation',\n", + " hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]\n", + "\n", + "data2_2[7]['visible'] = True\n", + "\n", + "data = data1_1 + data2_1 + data1_2 + data2_2\n", + "\n", + "steps = []\n", + "for i in range(len(Nex_range)):\n", + " step = dict(\n", + " method = 'restyle', \n", + " args = ['visible', [False] * len(data1_1)],\n", + " label = str(Nex_range[i])\n", + " )\n", + " step['args'][1][i] = True # Toggle i'th trace to \"visible\"\n", + " steps.append(step)\n", + "\n", + "sliders = [dict(\n", + " x = 0,\n", + " y = -0.02,\n", + " active = 5,\n", + " currentvalue = {\"prefix\": \"nth TR1-TR2: \"},\n", + " pad = {\"t\": 50, \"b\": 10},\n", + " steps = steps)]\n", + "\n", + "layout = go.Layout(\n", + " width=580,\n", + " height=450,\n", + " margin=go.layout.Margin(\n", + " l=80,\n", + " r=40,\n", + " b=60,\n", + " t=10,\n", + " ),\n", + " annotations=[\n", + " dict(\n", + " x=0.5004254919715793,\n", + " y=-0.18,\n", + " showarrow=False,\n", + " text='Excitation Flip Angle (°)',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=-0.15,\n", + " y=0.5,\n", + " showarrow=False,\n", + " text='Signal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " textangle=-90,\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis=dict(\n", + " autorange=False,\n", + " range=[0, params['EXC_FA'][-1]],\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " yaxis=dict(\n", + " autorange=True,\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " legend=dict(\n", + " x=0.5,\n", + " y=0.9,\n", + " traceorder='normal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12,\n", + " color='#000'\n", + " ),\n", + " bordercolor='#000000',\n", + " borderwidth=2\n", + " ), \n", + " sliders=sliders\n", + ")\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-line', config = config)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "In practice, gradient and RF spoiling are important parameters to consider in an AFI experiment. A combination of both (Zur et al. 1991; Bernstein et al. 2004) is typically recommended, and Figure 4 shows how this better approximates the ideal spoiling case. \n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 4. Signal 1 curves estimated using Bloch simulations for three categories of signal spoiling: (1) ideal spoiling (blue), gradient & RF Spoiling (red), and no spoiling (orange). Simulation details: TR1 = 20 ms, TR2 = 100 ms, T1 = 900 ms, T2 = 100 ms, TE = 5 ms, 100 spins. For the ideal spoiling case, the transverse magnetization is set to zero at the end of each TR. For the gradient & RF spoiling case, each spin is rotated by different increments of phase (2𝜋 / # of spins) to simulate complete dephasing from gradient spoiling, and the RF phase of the excitation pulse is ɸn = ɸn-1 + nɸ0 = ½ ɸ0(n2 + n + 2) (Bernstein et al. 2004) with ɸ0 = 39° (Zur et al. 1991) after each TR.\n", + "\n", + "

\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to generate the data required for Figure 4 of the blog post\n", + "\n", + "clear all\n", + "\n", + "%% Setup parameters\n", + "% All times are in milliseconds\n", + "% All flip angles are in degrees\n", + "\n", + "% White matter\n", + "params.T1 = 900; % in milliseconds\n", + "params.T2 = 100;\n", + "params.TR1 = 20;\n", + "params.TR2 = 100;\n", + "params.TE = 5;\n", + "params.EXC_FA = 1:90;\n", + "Nex_range = [1:9, 10:10:100];\n", + "\n", + "%% Calculate signals\n", + "% This results were precomputed. Uncomment the lines below to run the simulation.\n", + "signal1_ideal = load('../afi_blog/blochsim/signal1_ideal_spoil.mat');\n", + "signal2_ideal = load('../afi_blog/blochsim/signal2_ideal_spoil.mat');\n", + "signal1_optimal = load('../afi_blog/blochsim/signal1_optimal_crush_and_rf_spoil.mat');\n", + "signal2_optimal = load('../afi_blog/blochsim/signal2_optimal_crush_and_rf_spoil.mat');\n", + "signal1_no_spoil = load('../afi_blog/blochsim/signal1_no_gradient_and_rf_spoil.mat');\n", + "signal2_no_spoil = load('../afi_blog/blochsim/signal2_no_gradient_and_rf_spoil.mat');\n", + "signal1_ideal_spoil = signal1_ideal.signal1_ideal_spoil;\n", + "signal2_ideal_spoil = signal2_ideal.signal2_ideal_spoil;\n", + "signal1_optimal_crush_and_rf_spoil = signal1_optimal.signal1_optimal_crush_and_rf_spoil;\n", + "signal2_optimal_crush_and_rf_spoil = signal2_optimal.signal2_optimal_crush_and_rf_spoil;\n", + "signal1_no_gradient_and_rf_spoil = signal1_no_spoil.signal1_no_gradient_and_rf_spoil;\n", + "signal2_no_gradient_and_rf_spoil = signal2_no_spoil.signal2_no_gradient_and_rf_spoil;\n", + "\n", + "%{\n", + "% To see all the options available, run `help b1_afi.analytical_solution`\n", + "\n", + "for ii = 1:length(Nex_range)\n", + " params.Nex = Nex_range(ii);\n", + " \n", + " params.crushFlag = 1;\n", + " \n", + " [~, Msig1, Msig2] = b1_afi.bloch_sim(params);\n", + " signal1_ideal_spoil(ii,:) = abs(Msig1);\n", + " signal2_ideal_spoil(ii,:) = abs(Msig2);\n", + " \n", + " params.inc = 39;\n", + " params.partialDephasing = 1;\n", + " params.partialDephasingFlag = 1;\n", + " params.crushFlag = 0;\n", + " \n", + " [~, Msig1, Msig2] = b1_afi.bloch_sim(params);\n", + " signal1_optimal_crush_and_rf_spoil(ii,:) = abs(Msig1);\n", + " signal2_optimal_crush_and_rf_spoil(ii,:) = abs(Msig2);\n", + " \n", + " params.inc = 0;\n", + " params.partialDephasing = 0;\n", + "\n", + " [~, Msig1, Msig2] = b1_afi.bloch_sim(params);\n", + " signal1_no_gradient_and_rf_spoil(ii,:) = abs(Msig1);\n", + " signal2_no_gradient_and_rf_spoil(ii,:) = abs(Msig2);\n", + "end\n", + "%}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get params --from Octave\n", + "%get Nex_range --from Octave\n", + "%get signal1_ideal_spoil --from Octave\n", + "%get signal1_optimal_crush_and_rf_spoil --from Octave\n", + "%get signal1_no_gradient_and_rf_spoil --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "data1_1 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_ideal_spoil[ii]))),\n", + " name = 'Ideal Spoiling',\n", + " text = 'Ideal Spoiling',\n", + " hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]\n", + "\n", + "data1_1[10]['visible'] = True\n", + "\n", + "data2_1 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_optimal_crush_and_rf_spoil[ii]))),\n", + " name = 'Gradient & RF Spoiling',\n", + " text = 'Gradient & RF Spoiling',\n", + " hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]\n", + "\n", + "data2_1[10]['visible'] = True\n", + "\n", + "data3_1 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(signal1_no_gradient_and_rf_spoil[ii]))),\n", + " name = 'No Spoiling',\n", + " text = 'No Spoiling',\n", + " hoverinfo = 'x+y+text') for ii in range(len(Nex_range))]\n", + "\n", + "data3_1[10]['visible'] = True\n", + "\n", + "data = data1_1 + data2_1+ data3_1\n", + "\n", + "steps = []\n", + "for i in range(len(Nex_range)):\n", + " step = dict(\n", + " method = 'restyle', \n", + " args = ['visible', [False] * len(data1_1)],\n", + " label = str(Nex_range[i])\n", + " )\n", + " step['args'][1][i] = True # Toggle i'th trace to \"visible\"\n", + " steps.append(step)\n", + "\n", + "sliders = [dict(\n", + " x = 0,\n", + " y = -0.02,\n", + " active = 10,\n", + " currentvalue = {\"prefix\": \"nth TR1-TR2: \"},\n", + " pad = {\"t\": 50, \"b\": 10},\n", + " steps = steps\n", + ")]\n", + "\n", + "layout = go.Layout(\n", + " width=580,\n", + " height=450,\n", + " margin=go.layout.Margin(\n", + " l=80,\n", + " r=40,\n", + " b=60,\n", + " t=10,\n", + " ),\n", + " annotations=[\n", + " dict(\n", + " x=0.5004254919715793,\n", + " y=-0.18,\n", + " showarrow=False,\n", + " text='Excitation Flip Angle (°)',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=-0.15,\n", + " y=0.5,\n", + " showarrow=False,\n", + " text='Signal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " textangle=-90,\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis=dict(\n", + " autorange=False,\n", + " range=[0, params['EXC_FA'][-1]],\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " yaxis=dict(\n", + " autorange=True,\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " legend=dict(\n", + " x=0.5,\n", + " y=0.9,\n", + " traceorder='normal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12,\n", + " color='#000'\n", + " ),\n", + " bordercolor='#000000',\n", + " borderwidth=2\n", + " ), \n", + " sliders=sliders\n", + ")\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-line', config = config)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "

Data Fitting

" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "The ratio of Equations 1 and 2, gives rise to Equation 3 that depends on the parameters T1, TR1, TR2 and the excitation flip angle (θ). \n", + "

\n", + "\n", + "

\n", + "

\n", + "

\n", + "\n", + "

\n", + "Equation 3 can be simplified if the Taylor series expansion of the exponential function is used, followed by a first-order approximation to its terms. For this expansion, short TR1 and TR2 (TR1 < T1 and TR2 < T1) are assumed to approximate the signal intensities ratio (Equation 4) where n = TR2/TR1.\n", + "

\n", + "\n", + "

\n", + "

\n", + "

\n", + "\n", + "

\n", + "Finally, a measure of the actual flip-angle (θ) can be achieved by solving Equation 4 to obtain Equation 5, which only depends on the signal intensities ratio (r = S2/S1) and the parameters TR1 and TR2. \n", + "

\n", + "\n", + "

\n", + "

\n", + "

\n", + "\n", + "

\n", + "The actual flip-angle is estimated using an approximation (Equation 4) of a complete analytical solution (Equation 3), and the nature of this approximation makes it worthwhile to assess the accuracy of the signal intensities ratio between both equations. Next, a set of simulations are displayed to analyze how the choice of r is affected by T1, TR1 and TR2. First, the effect of the relaxation time T1 is simulated in Figure 5 for both the approximation and the complete analytical solution.\n", + "

\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 5. Effect of the relaxation time T1 on the ratio r. Signal intensities ratio is plotted as a function of the flip angle for the complete analytical solution (Equation 3 - blue) and the first-order approximation (Equation 4 - orange). AFI simulation details: TR1 = 20 ms, TR2 = 100 ms and variable T1.\n", + "\n", + "

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to generate the data required for Figure 5 of the blog post\n", + "\n", + "clear all\n", + "\n", + "%% Setup parameters\n", + "% All times are in milliseconds\n", + "% All flip angles are in degrees\n", + "\n", + "%% Calculate signals\n", + "%\n", + "% To see all the options available, run `help b1_afi.analytical_solution`\n", + "\n", + "T1_range = [10:10:200, 200:400:2000];\n", + "\n", + "params.TR1 = 20;\n", + "params.TR2 = 100;\n", + "params.EXC_FA = 1:90;\n", + "n = params.TR2/params.TR1;\n", + "\n", + "for jj=1:length(T1_range)\n", + " params.T1 = T1_range(jj);\n", + " \n", + " % Range of T1\n", + " [Mz1, Mz2] = b1_afi.analytical_solution(params);\n", + " signal1_analytical = Mz1;\n", + " signal2_analytical = Mz2;\n", + "\n", + " r_analytical(jj,:) = signal2_analytical./signal1_analytical;\n", + " r_approximation(jj,:) = (1 + n.*cosd(params.EXC_FA))./(n + cosd(params.EXC_FA));\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get params --from Octave\n", + "%get T1_range --from Octave\n", + "%get r_analytical --from Octave\n", + "%get r_approximation --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "data1 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(r_analytical[ii]))),\n", + " name = 'Analytical',\n", + " text = 'Analytical',\n", + " hoverinfo = 'x+y+text') for ii in range(len(T1_range))]\n", + "\n", + "data1[10]['visible'] = True\n", + "\n", + "data2 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(r_approximation[ii]))),\n", + " name = 'Approximation',\n", + " text = 'Approximation',\n", + " hoverinfo = 'x+y+text') for ii in range(len(T1_range))]\n", + "\n", + "data2[10]['visible'] = True\n", + "\n", + "data = data1 + data2\n", + "\n", + "steps = []\n", + "for i in range(len(T1_range)):\n", + " step = dict(\n", + " method = 'restyle', \n", + " args = ['visible', [False] * len(data1)],\n", + " label = str(T1_range[i])\n", + " )\n", + " step['args'][1][i] = True # Toggle i'th trace to \"visible\"\n", + " steps.append(step)\n", + "\n", + "sliders = [dict(\n", + " x = 0,\n", + " y = -0.02,\n", + " active = 10,\n", + " currentvalue = {\"prefix\": \"T1 value (ms): \"},\n", + " pad = {\"t\": 50, \"b\": 10},\n", + " steps = steps)]\n", + "\n", + "layout = go.Layout(\n", + " width=580,\n", + " height=450,\n", + " margin=go.layout.Margin(\n", + " l=80,\n", + " r=40,\n", + " b=60,\n", + " t=10,\n", + " ),\n", + " annotations=[\n", + " dict(\n", + " x=0.5004254919715793,\n", + " y=-0.18,\n", + " showarrow=False,\n", + " text='Excitation Flip Angle (°)',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=-0.15,\n", + " y=0.5,\n", + " showarrow=False,\n", + " text='r = S2/S1',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " textangle=-90,\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis=dict(\n", + " autorange=False,\n", + " range=[0, params['EXC_FA'][-1]],\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " yaxis=dict(\n", + " autorange=True,\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " legend=dict(\n", + " x=0.5,\n", + " y=0.9,\n", + " traceorder='normal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12,\n", + " color='#000'\n", + " ),\n", + " bordercolor='#000000',\n", + " borderwidth=2\n", + " ), \n", + " sliders=sliders\n", + ")\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-line', config = config)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "\n", + "

\n", + "The signal ratio r is highly insensitive to the relaxation time T1, except for the low T1 values and large flip angles (>70°). This shows that the Taylor expansion is a good approximation to the signal ratio r because it is possible to get rid of the inverse quadratic T1 dependance by taking the first-order terms of the expansion.\n", + "

\n", + "\n", + "

\n", + "The effect of the TR1 parameter on the signal ratio is shown in Figure 6. To assess the influence of the repetition time, we fix n=5 and vary the parameter TR1 in accordance to the relation n = TR2/TR1. As TR1 increases (> 50 ms), the approximated ratio r slightly deviates from the analytical approach. Although the deviation is slight only at high flip angles, a good signal ratio approximation can be achieved for a wide range of flip angles and repetition times.\n", + "

\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 6. Effect of the repetition time TR1 on the ratio r. Signal intensities ratio is plotted as a function of the flip angle for the complete analytical solution (Equation 3 - blue) and the first-order approximation (Equation 4 - orange). AFI simulation details: Variable TR1 ranging from 10 to 60 ms, fixed ratio n = 5 and T1 = 900 ms.\n", + "\n", + "

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to generate the data required for Figure 6 of the blog post\n", + "\n", + "clear all\n", + "\n", + "%% Setup parameters\n", + "% All times are in milliseconds\n", + "% All flip angles are in degrees\n", + "\n", + "%% Calculate signals\n", + "%\n", + "% To see all the options available, run `help b1_afi.analytical_solution`\n", + "\n", + "TR1_range = 10:5:60;\n", + "\n", + "params.T1 = 900;\n", + "params.EXC_FA = 1:90;\n", + "n = 5;\n", + "\n", + "for jj=1:length(TR1_range)\n", + " params.TR1 = TR1_range(jj);\n", + " params.TR2 = n*params.TR1;\n", + " \n", + " % Fixed: T1 = 900 ms, n = 5\n", + " [Mz1, Mz2] = b1_afi.analytical_solution(params);\n", + " signal1_analytical = Mz1;\n", + " signal2_analytical = Mz2;\n", + "\n", + " r_analytical(jj,:) = signal2_analytical./signal1_analytical;\n", + " r_approximation(jj,:) = (1 + n.*cosd(params.EXC_FA))./(n + cosd(params.EXC_FA));\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get params --from Octave\n", + "%get TR1_range --from Octave\n", + "%get r_analytical --from Octave\n", + "%get r_approximation --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "data1 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(r_analytical[ii]))),\n", + " name = 'Analytical',\n", + " text = 'Analytical',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data1[2]['visible'] = True\n", + "\n", + "data2 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(r_approximation[ii]))),\n", + " name = 'Approximation',\n", + " text = 'Approximation',\n", + " hoverinfo = 'x+y+text') for ii in range(len(TR1_range))]\n", + "\n", + "data2[2]['visible'] = True\n", + "\n", + "data = data1 + data2\n", + "\n", + "steps = []\n", + "for i in range(len(TR1_range)):\n", + " step = dict(\n", + " method = 'restyle', \n", + " args = ['visible', [False] * len(data1)],\n", + " label = str(TR1_range[i])\n", + " )\n", + " step['args'][1][i] = True # Toggle i'th trace to \"visible\"\n", + " steps.append(step)\n", + "\n", + "sliders = [dict(\n", + " x = 0,\n", + " y = -0.02,\n", + " active = 2,\n", + " currentvalue = {\"prefix\": \"TR1 value (ms): \"},\n", + " pad = {\"t\": 50, \"b\": 10},\n", + " steps = steps)]\n", + "\n", + "layout = go.Layout(\n", + " width=580,\n", + " height=450,\n", + " margin=go.layout.Margin(\n", + " l=80,\n", + " r=40,\n", + " b=60,\n", + " t=10,\n", + " ),\n", + " annotations=[\n", + " dict(\n", + " x=0.5004254919715793,\n", + " y=-0.18,\n", + " showarrow=False,\n", + " text='Excitation Flip Angle (°)',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=-0.15,\n", + " y=0.5,\n", + " showarrow=False,\n", + " text='r = S2/S1',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " textangle=-90,\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis=dict(\n", + " autorange=False,\n", + " range=[0, params['EXC_FA'][-1]],\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " yaxis=dict(\n", + " autorange=True,\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " legend=dict(\n", + " x=0.5,\n", + " y=0.9,\n", + " traceorder='normal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12,\n", + " color='#000'\n", + " ),\n", + " bordercolor='#000000',\n", + " borderwidth=2\n", + " ), \n", + " sliders=sliders\n", + ")\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-line', config = config)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "\n", + "

\n", + "Finally, the effect of the parameter n on the signal ratio r (Figure 7) does not seem to significantly affect the signal ratio between the approximated equation and the analytical approach. However, the parameter n has a major impact on the sensitivity of the AFI method to variations in the flip angle. Figure 7 shows that the increase of the parameter n (= TR2/TR1) allows for improvement of the dynamic range of flip angles measurements. These simulations have shown that an optimal implementation of the AFI method requires a careful selection of sequence parameters.\n", + "

\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 7. Effect of n (TR2 to TR1 ratio) on the ratio r. The signal intensities ratio is plotted as a function of the flip angle for the complete analytical solution (Equation 3 - blue) and the first-order approximation (Equation 4 - orange). AFI simulation details: Variable n ranging from 2 to 6, fixed TR1 = 20 ms and T1 = 900 ms.\n", + "\n", + "

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to generate the data required for Figure 7 of the blog post\n", + "\n", + "clear all\n", + "\n", + "%% Setup parameters\n", + "% All times are in milliseconds\n", + "% All flip angles are in degrees\n", + "\n", + "%% Calculate signals\n", + "%\n", + "% To see all the options available, run `help b1_afi.analytical_solution`\n", + "\n", + "n_range = 2:1:6;\n", + "\n", + "params.T1 = 900;\n", + "params.TR1 = 20;\n", + "params.EXC_FA = 1:90;\n", + "\n", + "for jj=1:length(n_range)\n", + " n = n_range(jj);\n", + " params.TR2 = n*params.TR1;\n", + " \n", + " % Fixed: T1 = 900 ms, TR1 = 20\n", + " [Mz1, Mz2] = b1_afi.analytical_solution(params);\n", + " signal1_analytical = Mz1;\n", + " signal2_analytical = Mz2;\n", + "\n", + " r_analytical(jj,:) = signal2_analytical./signal1_analytical;\n", + " r_approximation(jj,:) = (1 + n.*cosd(params.EXC_FA))./(n + cosd(params.EXC_FA));\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get params --from Octave\n", + "%get n_range --from Octave\n", + "%get r_analytical --from Octave\n", + "%get r_approximation --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "data1 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(r_analytical[ii]))),\n", + " name = 'Analytical',\n", + " text = 'Analytical',\n", + " hoverinfo = 'x+y+text') for ii in range(len(n_range))]\n", + "\n", + "data1[2]['visible'] = True\n", + "\n", + "data2 = [dict(\n", + " visible = False,\n", + " mode = 'lines',\n", + " x = params[\"EXC_FA\"],\n", + " y = abs(np.squeeze(np.asarray(r_approximation[ii]))),\n", + " name = 'Approximation',\n", + " text = 'Approximation',\n", + " hoverinfo = 'x+y+text') for ii in range(len(n_range))]\n", + "\n", + "data2[2]['visible'] = True\n", + "\n", + "data = data1 + data2\n", + "\n", + "steps = []\n", + "for i in range(len(n_range)):\n", + " step = dict(\n", + " method = 'restyle', \n", + " args = ['visible', [False] * len(data1)],\n", + " label = str(n_range[i])\n", + " )\n", + " step['args'][1][i] = True # Toggle i'th trace to \"visible\"\n", + " steps.append(step)\n", + "\n", + "sliders = [dict(\n", + " x = 0,\n", + " y = -0.02,\n", + " active = 2,\n", + " currentvalue = {\"prefix\": \"n value: \"},\n", + " pad = {\"t\": 50, \"b\": 10},\n", + " steps = steps)]\n", + "\n", + "layout = go.Layout(\n", + " width=580,\n", + " height=450,\n", + " margin=go.layout.Margin(\n", + " l=80,\n", + " r=40,\n", + " b=60,\n", + " t=10,\n", + " ),\n", + " annotations=[\n", + " dict(\n", + " x=0.5004254919715793,\n", + " y=-0.18,\n", + " showarrow=False,\n", + " text='Excitation Flip Angle (°)',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=-0.15,\n", + " y=0.5,\n", + " showarrow=False,\n", + " text='r = S2/S1',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=22\n", + " ),\n", + " textangle=-90,\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis=dict(\n", + " autorange=False,\n", + " range=[0, params['EXC_FA'][-1]],\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " yaxis=dict(\n", + " autorange=False,\n", + " range=[0, 1],\n", + " showgrid=False,\n", + " linecolor='black',\n", + " linewidth=2\n", + " ),\n", + " legend=dict(\n", + " x=0.5,\n", + " y=0.9,\n", + " traceorder='normal',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12,\n", + " color='#000'\n", + " ),\n", + " bordercolor='#000000',\n", + " borderwidth=2\n", + " ), \n", + " sliders=sliders\n", + ")\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-line', config = config)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "Figure 8 displays an example AFI dataset and its corresponding field B1 map in a healthy human brain. Although not clearly visible, both AFI images present a small Gibbs ringing artifact that is propagated and amplified due to the AFI calculation consisting of the division of both images (Boudreau et al. 2017). The ringing artifact is clearly seen in the unfiltered/raw B1 field map shown in Figure 8 (right).\n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "\n", + "Figure 8. Example actual flip-angle imaging dataset (left) and a resulting raw B1 map of a healthy adult brain (right). The relevant VFA protocol parameters used were: TR1 = 20 ms, TR2 = 100 ms and θnominal = 60°. The B1 map (right) was fitted using the approximate r ratio (Equation 5).\n", + "\n", + "

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Download Actual Flip-Angle Imaging brain MRI data for Figure 8 of the blog post\n", + "\n", + "cmd = ['curl -L -o b1_afi.zip https://osf.io/csjgx/download?version=6'];\n", + "[STATUS,MESSAGE] = unix(cmd);\n", + "unzip('b1_afi.zip');\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to generate the data required for Figure 8 of the blog post\n", + "\n", + "clear all\n", + "\n", + "% Format qMRLab b1_afi model parameters, and load them into the Model object\n", + "Model = b1_afi; \n", + "Model.Prot.Sequence.Mat = [60, 20, 100];\n", + "\n", + "% Format data structure so that they may be fit by the model\n", + "data = struct();\n", + "data.AFIData1 = load_nii_data('b1_afi/AFIData1.nii.gz');\n", + "data.AFIData2 = load_nii_data('b1_afi/AFIData2.nii.gz');\n", + "\n", + "FitResults = FitData(data,Model,0); % The '0' flag is so that no wait bar is shown." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "%% MATLAB/OCTAVE CODE\n", + "% Code used to re-orient the images to make pretty figures, and to assign variables with the axis lengths.\n", + "\n", + "B1map_raw = imrotate(squeeze(FitResults.B1map_raw(:,:,14)),-90);\n", + "\n", + "xAxis = [0:size(B1map_raw,2)-1];\n", + "yAxis = [0:size(B1map_raw,1)-1];\n", + "\n", + "% Raw MRI AFI data\n", + "AFIData1 = imrotate(squeeze(data.AFIData1(:,:,14)),-90);\n", + "AFIData2 = imrotate(squeeze(data.AFIData2(:,:,14)),-90);\n", + "\n", + "% Load mask\n", + "mask = load_nii_data('../afi_blog/fitResults/FitResults_median_filter/mask.nii.gz');\n", + "mask = imrotate(squeeze(mask(:,:,14)),-90);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get B1map_raw --from Octave\n", + "%get AFIData1 --from Octave\n", + "%get AFIData2 --from Octave\n", + "%get mask --from Octave\n", + "%get xAxis --from Octave\n", + "%get yAxis --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "from plotly import tools\n", + "\n", + "# Masking B1 map\n", + "B1map_raw = np.asarray(B1map_raw)\n", + "mask = np.asarray(mask)\n", + "B1map_raw_masked = B1map_raw*mask\n", + "B1map_raw_masked[np.isnan(B1map_raw_masked)] = 0\n", + "\n", + "\n", + "trace1 = go.Heatmap(x = xAxis,\n", + " y = yAxis,\n", + " z=AFIData1,\n", + " colorscale='Greys',\n", + " showscale = False,\n", + " visible=False,\n", + " name = 'Signal1')\n", + "trace2 = go.Heatmap(x = xAxis,\n", + " y = yAxis,\n", + " z=AFIData2,\n", + " colorscale='Greys',\n", + " showscale = False,\n", + " visible=True,\n", + " name = 'Signal2')\n", + "trace3 = go.Heatmap(x = xAxis,\n", + " y = yAxis,\n", + " z=B1map_raw_masked,\n", + " zmin=0.7,\n", + " zmax=1.3,\n", + " colorscale='RdBu',\n", + " xaxis='x2',\n", + " yaxis='y2',\n", + " visible=True,\n", + " name = 'B1 values')\n", + "\n", + "data=[trace1, trace2, trace3]\n", + "\n", + "\n", + "updatemenus = list([\n", + " dict(active=1,\n", + " x = 0.09,\n", + " xanchor = 'left',\n", + " y = -0.15,\n", + " yanchor = 'bottom',\n", + " direction = 'up',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=16\n", + " ),\n", + " buttons=list([ \n", + " dict(label = 'Signal 1',\n", + " method = 'update',\n", + " args = [{'visible': [True, False, True]},\n", + " ]),\n", + " dict(label = 'Signal 2',\n", + " method = 'update',\n", + " args = [{'visible': [False, True, True]},\n", + " ]),\n", + " ])\n", + " )\n", + "])\n", + "\n", + "layout = dict(\n", + " width=560,\n", + " height=345,\n", + " margin = dict(\n", + " t=40,\n", + " r=50,\n", + " b=10,\n", + " l=50),\n", + " annotations=[\n", + " dict(\n", + " x=0.07,\n", + " y=1.15,\n", + " showarrow=False,\n", + " text='Input Data',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=26\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=0.60,\n", + " y=1.15,\n", + " showarrow=False,\n", + " text='Raw B1 map',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=26\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=1.12,\n", + " y=1.15,\n", + " showarrow=False,\n", + " text='B1',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=26\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0, 0.58]),\n", + " yaxis = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0, 1]),\n", + " xaxis2 = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0.40, 0.98]),\n", + " yaxis2 = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0, 1], anchor='x2'),\n", + " showlegend = False,\n", + " autosize = False,\n", + " updatemenus=updatemenus\n", + ")\n", + "\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-heatmap', config = config)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "The ringing artifact shown in Figure 8 can be attenuated by implementing a smoothing process. Figure 9 shows the raw (left) and the filtered (right) B1 map where a median filter was used to smooth the field map.\n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + " \n", + "Figure 9. Raw (left) and filtered (right) B1 map. A median filter of size 7x7x7 pixels was used to attenuate the Gibbs ringing artifact. \n", + "\n", + "

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Octave" + }, + "outputs": [], + "source": [ + "B1map_filtered = load_nii_data('../afi_blog/fitResults/FitResults_median_filter/B1map_filtered.nii.gz');\n", + "B1map_filtered = imrotate(squeeze(B1map_filtered(:,:,14)),-90);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "scratch" + ] + }, + "outputs": [], + "source": [ + "%get B1map_filtered --from Octave" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3", + "tags": [ + "report_output" + ] + }, + "outputs": [], + "source": [ + "from plotly import tools\n", + "\n", + "# Masking B1 map\n", + "B1map_filtered = np.asarray(B1map_filtered)\n", + "mask = np.asarray(mask)\n", + "B1map_filtered_masked = B1map_filtered*mask\n", + "B1map_filtered_masked[np.isnan(B1map_filtered_masked)] = 0\n", + "\n", + "trace1 = go.Heatmap(x = xAxis,\n", + " y = yAxis,\n", + " z=B1map_raw_masked,\n", + " zmin=0.7,\n", + " zmax=1.3,\n", + " colorscale='RdBu',\n", + " showscale = False,\n", + " visible=True,\n", + " name = 'B1 values')\n", + "trace2 = go.Heatmap(x = xAxis,\n", + " y = yAxis,\n", + " z=B1map_filtered_masked,\n", + " zmin=0.7,\n", + " zmax=1.3,\n", + " colorscale='RdBu',\n", + " xaxis='x2',\n", + " yaxis='y2',\n", + " visible=True,\n", + " name = 'B1 values')\n", + "\n", + "data=[trace1, trace2]\n", + "\n", + "layout = dict(\n", + " width=560,\n", + " height=310,\n", + " margin = dict(\n", + " t=40,\n", + " r=50,\n", + " b=10,\n", + " l=50),\n", + " annotations=[\n", + " dict(\n", + " x=0.04,\n", + " y=1.15,\n", + " showarrow=False,\n", + " text='Raw B1 map',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=26\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=0.60,\n", + " y=1.15,\n", + " showarrow=False,\n", + " text='Filtered B1 map',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=26\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " dict(\n", + " x=1.12,\n", + " y=1.15,\n", + " showarrow=False,\n", + " text='B1',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=26\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + " ],\n", + " xaxis = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0, 0.58]),\n", + " yaxis = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0, 1]),\n", + " xaxis2 = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0.40, 0.98]),\n", + " yaxis2 = dict(range = [0,127], autorange = False,\n", + " showgrid = False, zeroline = False, showticklabels = False,\n", + " ticks = '', domain=[0, 1], anchor='x2'),\n", + " showlegend = False,\n", + " autosize = False,\n", + ")\n", + "\n", + "\n", + "fig = dict(data=data, layout=layout)\n", + "\n", + "iplot(fig, filename = 'basic-heatmap', config = config)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "

Benefits and Pitfalls

" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "B1 mapping is of interest for diverse MRI applications, and several mapping techniques have been developed. The DAM method consists of acquiring two scans at two different flip angles. To avoid the dependence of the signal on T1, long repetition times are required to allow the recovery of the longitudinal magnetization between pulses (Yarnykh 2007; Insko & Bolinger 1993). The AFI method overcomes this practical limitation by repeating the pulse sequence at a fast rate to achieve a pulsed state of magnetization and shorter time delays between pulses. In addition, due to scan-time constraints, B1 mapping methods are often implemented in 2D (Chavez & Stanisz 2011). However, the accuracy of the measurements of 2D B1 mapping techniques is compromised by the slice profile effects due to the problem of nonuniform excitation across slices (Yarnykh 2007; Chavez & Stanisz 2011). The AFI method on the other hand, adresses this issue using a fast 3D implementation leading to scans with an excellent anatomical coverage in clinically feasible times, with an increase in signal-to-noise ratio compared to 2D multislice acquisitions.\n", + "

\n", + "\n", + "

\n", + "The performance of the AFI method is based on the following assumptions. First, the two images acquired at different times should be registered to avoid motion effects. It is also assumed that the signal is insensitive to the main magnetic field non-uniformities and chemical shift effects that are canceled out when taking the signal ratio r (Yarnykh 2007). Despite some clear advantages over other B1 mapping techniques, the application of spoiler gradients to mitigate the T1 dependence can be a limitation due to significant levels of RF power depositions (Sacolick et al. 2010). Furthermore, it is necessary to adapt the AFI pulse sequence to different scanner manufacturers, and programming experience is required to bring this technique into the clinic.\n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "

Works Cited

" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "kernel": "Python3" + }, + "source": [ + "
\n", + "

\n", + "Bernstein, M., King, K. & Zhou, X., 2004. Handbook of MRI Pulse Sequences, Elsevier.\n", + "

\n", + " \n", + "

\n", + "Boudreau, M., Tardif, C.L., Stikov, N., Sled, J.G., Lee, W. & Pike, G.B., 2017. B1 mapping for bias-correction in quantitative T1 imaging of the brain at 3T using standard pulse sequences. J. Magn. Reson. Imaging, 46, pp.1673-1682.\n", + "

\n", + " \n", + "

\n", + "Boudreau, M., Stikov, N. & Pike, G.B, 2018. B1-sensitivity analysis of quantitative magnetization transfer imaging. Magn. Reson. Med., 79, pp.276-285.\n", + "

\n", + " \n", + "

\n", + "Chavez, S. & Stanisz, G.J., 2012. A novel method for simultaneous 3D B1 and T1 mapping: the method of slopes (MoS). NMR Biomed., 25, pp.1043-1055.\n", + "

\n", + " \n", + "

\n", + "Deoni, S.C., 2007. High-resolution T1 mapping of the brain at 3T with driven equilibrium single pulse observation of T1 with high-speed incorporation of RF field inhomogeneities (DESPOT1-HIFI). J. Magn. Reson. Imaging, 26, pp.1106-1111.\n", + "

\n", + " \n", + "

\n", + "Ibrahim, T.S., Lee, R., Baertlein, B.A. & Robitaille, P.M., 2001. B1 field homogeneity and SAR calculations for the birdcage coil. Phys Med Biol., 46(2), pp.609-619.\n", + "

\n", + "\n", + "

\n", + "Insko, E. & Bolinger, L., 1993. Mapping of the radio frequency field. J. Magn. Reson. Ser. A, 103(1), pp.82–85.\n", + "

\n", + " \n", + "

\n", + "Katscher, U., Voigt, T., Findeklee, C., Vernickel, P., Nehrke, K. & Dössel, O., 2009. Determination of electric conductivity and local SAR via B1 mapping. IEEE Trans Med Imaging., 28(9), pp.1365-1374.\n", + "

\n", + " \n", + "

\n", + "Ropele, S., Filippi, M., Valsasina, P., Korteweg, T., Barkhof, F., Tofts, P.S., Samson, R., Miller, D.H. & Fazekas, F., 2005. Assessment and correction of B1-induced errors in magnetization transfer ratio measurements. Magn. Reson. Med., 53, pp.134-140. \n", + "

\n", + " \n", + "

\n", + "Sacolick, L.I., Wiesinger, F., Hancu, I. & Vogel, M.W., 2010. B1 mapping by Bloch-Siegert shift. Magn. Reson. Med., 63, pp.1315-1322.\n", + "

\n", + " \n", + "

\n", + "Sled, J.G. & Pike, G.B., 1998. Standing-wave and RF penetration artifacts caused by elliptic geometry: an electrodynamic analysis of MRI. IEEE Trans. Med. Imaging., 17(4), pp.653-662.\n", + "

\n", + "\n", + "

\n", + "Sled, J.G. & Pike, G.B., 2000. Correction for B1 and B0 variations in quantitative T2 measurements using MRI. Magn. Reson. Med., 43(4), pp.589–593.\n", + "

\n", + "\n", + "

\n", + "Stollberger, R. & Wach, P., 1996. Imaging of the active B1 field in vivo. Magn. Reson. Med., 35(2), pp.246–251.\n", + "

\n", + "\n", + "

\n", + "van den Bergen, B., Van den Berg, C.A., Bartels, L.W. & Lagendijk, J.J., 2007. 7 T body MRI: B1 shimming with simultaneous SAR reduction. Phys Med Biol., 52(17), pp.5429-5441.\n", + "

\n", + "\n", + "

\n", + "Yarnykh, V.L. & Yuan, C. Actual flip angle imaging in the pulsed steady state. In: Proceedings of the 12th Annual Meeting of ISMRM, Kyoto, Japan, 2004 (Abstract 194).\n", + "

\n", + "\n", + "

\n", + "Yarnykh, V.L., 2007. Actual Flip-Angle Imaging in the pulsed steady state: A method for rapid three-dimensional mapping of the transmitted radiofrequency field. Magn. Reson. Med., 57(1), pp.192-200.\n", + "

\n", + " \n", + "

\n", + "Yarnykh, V.L., 2010. Optimal radiofrequency and gradient spoiling for improved accuracy of T1 and B1 measurements using fast steady-state techniques. Magn. Reson. Med., 63(6), pp.1610–1626.\n", + "

\n", + " \n", + "

\n", + "Zur, Y., Wood, M.L. & Neuringer, L.J., 1991. Spoiling of transverse magnetization in steady-state sequences. Magn. Reson. Med., 21(2), pp.251–263.\n", + "

\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3" + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "\n", + "display(HTML(\n", + " ''\n", + "))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "Python3" + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "celltoolbar": "Tags", + "kernelspec": { + "display_name": "SoS", + "language": "sos", + "name": "sos" + }, + "language_info": { + "codemirror_mode": "sos", + "file_extension": ".sos", + "mimetype": "text/x-sos", + "name": "sos", + "nbconvert_exporter": "sos_notebook.converter.SoS_Exporter", + "pygments_lexer": "sos" + }, + "sos": { + "kernels": [ + [ + "Octave", + "octave", + "Octave", + "#dff8fb" + ], + [ + "Python3", + "python3", + "Python3", + "#FFD91A" + ] + ], + "panel": { + "displayed": true, + "height": 0, + "style": "side" + }, + "version": "0.17.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/afi_blog/afi_pulsesequence.png b/afi_blog/afi_pulsesequence.png new file mode 100644 index 0000000..a702fa0 Binary files /dev/null and b/afi_blog/afi_pulsesequence.png differ diff --git a/afi_blog/afi_pulsesequence.svg b/afi_blog/afi_pulsesequence.svg new file mode 100644 index 0000000..1face31 --- /dev/null +++ b/afi_blog/afi_pulsesequence.svg @@ -0,0 +1,676 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + TR1 + + TR2 + + TR + + + IMG + + + + + + + + + + + + + + + + + + + + + + + + + + SPOIL + + α + + + + + + + + + + + + IMG + + + + + + + + + + + + + + + + + + + + + + + + + SPOIL + + α + + + + + + + + + + + + + + diff --git a/afi_blog/blochsim/signal1_analytical.mat b/afi_blog/blochsim/signal1_analytical.mat new file mode 100644 index 0000000..b026104 Binary files /dev/null and b/afi_blog/blochsim/signal1_analytical.mat differ diff --git a/afi_blog/blochsim/signal1_blochsim.mat b/afi_blog/blochsim/signal1_blochsim.mat new file mode 100644 index 0000000..f9464f7 Binary files /dev/null and b/afi_blog/blochsim/signal1_blochsim.mat differ diff --git a/afi_blog/blochsim/signal1_ideal_spoil.mat b/afi_blog/blochsim/signal1_ideal_spoil.mat new file mode 100644 index 0000000..1713d74 Binary files /dev/null and b/afi_blog/blochsim/signal1_ideal_spoil.mat differ diff --git a/afi_blog/blochsim/signal1_no_gradient_and_rf_spoil.mat b/afi_blog/blochsim/signal1_no_gradient_and_rf_spoil.mat new file mode 100644 index 0000000..3a7feb3 Binary files /dev/null and b/afi_blog/blochsim/signal1_no_gradient_and_rf_spoil.mat differ diff --git a/afi_blog/blochsim/signal1_optimal_crush_and_rf_spoil.mat b/afi_blog/blochsim/signal1_optimal_crush_and_rf_spoil.mat new file mode 100644 index 0000000..46cc369 Binary files /dev/null and b/afi_blog/blochsim/signal1_optimal_crush_and_rf_spoil.mat differ diff --git a/afi_blog/blochsim/signal2_analytical.mat b/afi_blog/blochsim/signal2_analytical.mat new file mode 100644 index 0000000..5cab0c0 Binary files /dev/null and b/afi_blog/blochsim/signal2_analytical.mat differ diff --git a/afi_blog/blochsim/signal2_blochsim.mat b/afi_blog/blochsim/signal2_blochsim.mat new file mode 100644 index 0000000..7676c6f Binary files /dev/null and b/afi_blog/blochsim/signal2_blochsim.mat differ diff --git a/afi_blog/blochsim/signal2_ideal_spoil.mat b/afi_blog/blochsim/signal2_ideal_spoil.mat new file mode 100644 index 0000000..e6b9a2b Binary files /dev/null and b/afi_blog/blochsim/signal2_ideal_spoil.mat differ diff --git a/afi_blog/blochsim/signal2_no_gradient_and_rf_spoil.mat b/afi_blog/blochsim/signal2_no_gradient_and_rf_spoil.mat new file mode 100644 index 0000000..a53d7f0 Binary files /dev/null and b/afi_blog/blochsim/signal2_no_gradient_and_rf_spoil.mat differ diff --git a/afi_blog/blochsim/signal2_optimal_crush_and_rf_spoil.mat b/afi_blog/blochsim/signal2_optimal_crush_and_rf_spoil.mat new file mode 100644 index 0000000..2d9cb2a Binary files /dev/null and b/afi_blog/blochsim/signal2_optimal_crush_and_rf_spoil.mat differ diff --git a/afi_blog/equation1.png b/afi_blog/equation1.png new file mode 100644 index 0000000..819afca Binary files /dev/null and b/afi_blog/equation1.png differ diff --git a/afi_blog/equation2.png b/afi_blog/equation2.png new file mode 100644 index 0000000..53aea7a Binary files /dev/null and b/afi_blog/equation2.png differ diff --git a/afi_blog/equation3.png b/afi_blog/equation3.png new file mode 100644 index 0000000..0aa2cd5 Binary files /dev/null and b/afi_blog/equation3.png differ diff --git a/afi_blog/equation4.png b/afi_blog/equation4.png new file mode 100644 index 0000000..0882b96 Binary files /dev/null and b/afi_blog/equation4.png differ diff --git a/afi_blog/equation5.png b/afi_blog/equation5.png new file mode 100644 index 0000000..84e1143 Binary files /dev/null and b/afi_blog/equation5.png differ diff --git a/afi_blog/fitResults/FitResults_median_filter/B1map_filtered.nii.gz b/afi_blog/fitResults/FitResults_median_filter/B1map_filtered.nii.gz new file mode 100644 index 0000000..0ec3acb Binary files /dev/null and b/afi_blog/fitResults/FitResults_median_filter/B1map_filtered.nii.gz differ diff --git a/afi_blog/fitResults/FitResults_median_filter/B1map_raw.nii.gz b/afi_blog/fitResults/FitResults_median_filter/B1map_raw.nii.gz new file mode 100644 index 0000000..b5e2fd1 Binary files /dev/null and b/afi_blog/fitResults/FitResults_median_filter/B1map_raw.nii.gz differ diff --git a/afi_blog/fitResults/FitResults_median_filter/FitResults.mat b/afi_blog/fitResults/FitResults_median_filter/FitResults.mat new file mode 100644 index 0000000..0fdedc0 Binary files /dev/null and b/afi_blog/fitResults/FitResults_median_filter/FitResults.mat differ diff --git a/afi_blog/fitResults/FitResults_median_filter/Spurious.nii.gz b/afi_blog/fitResults/FitResults_median_filter/Spurious.nii.gz new file mode 100644 index 0000000..3028040 Binary files /dev/null and b/afi_blog/fitResults/FitResults_median_filter/Spurious.nii.gz differ diff --git a/afi_blog/fitResults/FitResults_median_filter/mask.nii.gz b/afi_blog/fitResults/FitResults_median_filter/mask.nii.gz new file mode 100644 index 0000000..dbbc572 Binary files /dev/null and b/afi_blog/fitResults/FitResults_median_filter/mask.nii.gz differ