Data Fitting

At first glance, one could be tempted to fit VFA data using a non-linear least squares fitting algorithm such as Levenberg-Marquardt with Eq. 1, which typically only has two free fitting variables (T1 and M0). Although this is a valid way of estimating T1 from VFA data, it is rarely done in practice because a simple refactoring of Equation 1 allows T1 values to be estimated with a linear least square fitting algorithm, which substantially reduces the processing time. Without any approximations, Equation 1 can be rearranged into the form y = mx+b (Gupta 1977):

As the third term does not change between measurements (it is constant for each θn), it can be grouped into the constant for a simpler representation:

With this rearranged form of Equation 1, T1 can be simply estimated from the slope of a linear regression calculated from Sn/sin(θn) and Sn/tan(θn) values:

If data were acquired using only two flip angles – a very common VFA acquisition protocol – then the slope can be calculated using the elementary slope equation. Figure 5 displays both Equation 1 and 4 plotted for a noisy dataset.

Figure 5. Mean and standard deviation of the VFA signal plotted using the nonlinear form (Equation 1 – blue) and linear form (Equation 4 – red). Monte Carlo simulation details: SNR = 25, N = 1000. VFA simulation details: TR = 25 ms, T1 = 900 ms. 

%use octave

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

params.EXC_FA = [1:4,5:5:90];

%% Calculate signals
%
% To see all the options available, run `help vfa_t1.analytical_solution`

params.TR = 0.025;
params.EXC_FA = [2:9,10:5:90];

% White matter
x.M0 = 1;
x.T1 = 0.900; % in milliseconds

Model = vfa_t1; 

Opt.SNR = 25;
Opt.TR = params.TR;
Opt.T1 = x.T1;

clear Model.Prot.VFAData.Mat(:,1) 
Model.Prot.VFAData.Mat = zeros(length(params.EXC_FA),2);
Model.Prot.VFAData.Mat(:,1) = params.EXC_FA';
Model.Prot.VFAData.Mat(:,2) = Opt.TR;

for jj = 1:1000
    [FitResult{jj}, noisyData{jj}] = Model.Sim_Single_Voxel_Curve(x,Opt,0); 
    fittedT1(jj) = FitResult{jj}.T1;
    noisyData_array(jj,:) = noisyData{jj}.VFAData;
    noisyData_array_div_sin(jj,:) = noisyData_array(jj,:) ./ sind(Model.Prot.VFAData.Mat(:,1))';
    noisyData_array_div_tan(jj,:) = noisyData_array(jj,:) ./ tand(Model.Prot.VFAData.Mat(:,1))';
end
        
for kk=1:length(params.EXC_FA)
    data_mean(kk) = mean(noisyData_array(:,kk));
    data_std(kk) = std(noisyData_array(:,kk));
    
    data_mean_div_sin(kk) = mean(noisyData_array_div_sin(:,kk));
    data_std_div_sin(kk) = std(noisyData_array_div_sin(:,kk));
    
    data_mean_div_tan(kk) = mean(noisyData_array_div_tan(:,kk));
    data_std_div_tan(kk) = std(noisyData_array_div_tan(:,kk));
end


%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

params_highres.EXC_FA = [2:1:90];

%% Calculate signals
%
% To see all the options available, run `help vfa_t1.analytical_solution`

params_highres.TR = params.TR * 1000; % in milliseconds
    
% White matter
params_highres.T1 = x.T1*1000; % in milliseconds

signal_WM = vfa_t1.analytical_solution(params_highres);
signal_WM_div_sin = signal_WM ./ sind(params_highres.EXC_FA);
signal_WM_div_tan = signal_WM ./ tand(params_highres.EXC_FA);

(View simulation code) 

%use sos
%get params --from Octave
%get data_mean --from Octave
%get data_mean_div_sin --from Octave
%get data_mean_div_tan --from Octave
%get data_std --from Octave
%get data_std_div_sin --from Octave
%get data_std_div_tan --from Octave
%get params_highres --from Octave
%get signal_WM --from Octave
%get signal_WM_div_sin --from Octave
%get signal_WM_div_tan --from Octave

import matplotlib.pyplot as plt
import plotly.plotly as py
import plotly.graph_objs as go
import numpy as np
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
config={'showLink': False, 'displayModeBar': False}

init_notebook_mode(connected=True)

from IPython.core.display import display, HTML

data1 = dict(
        visible = True,
        x = params_highres["EXC_FA"],
        y = signal_WM,
        name = 'Analytical Solutions',
        text = params["EXC_FA"],
        mode = 'lines', 
        line = dict(
            color = ('rgb(0, 0, 0)'),
            dash = 'dot'),
        hoverinfo='none')

data2 = dict(
        visible = True,
        x = signal_WM_div_tan,
        y = signal_WM_div_sin,
        name = 'Analytical Solutions',
        text = params_highres["EXC_FA"],
        mode = 'lines',
        xaxis='x2',
        yaxis='y2',
        line = dict(
            color = ('rgb(0, 0, 0)'),
            dash = 'dot'
            ),
        hoverinfo='none',
        showlegend=False)

data3 = dict(
        visible = True,
        x = params["EXC_FA"],
        y = data_mean,
        name = 'Nonlinear Form - Noisy',
        text = ["Flip angle: " + str(x) + "°" for x in params["EXC_FA"]],
        mode = 'markers',
        hoverinfo = 'y+text',
        line = dict(
            color = ('rgb(22, 96, 167)'),
            ),
        error_y=dict(
            type='data',
            array=data_std,
            visible=True,
            color = ('rgb(142, 192, 240)')
        ))

data4 = dict(
        visible = True,
        x = data_mean_div_tan,
        y = data_mean_div_sin,
        name = 'Linear Form - Noisy',
        text = ["Flip angle: " + str(x) + "°" for x in params["EXC_FA"]],
        mode = 'markers',
        xaxis='x2',
        yaxis='y2',
        hoverinfo = 'x+y+text',
        line = dict(
            color = ('rgb(205, 12, 24)'),
            ),
        error_x=dict(
            type='data',
            array=data_std_div_tan,
            visible=True,
            color = ('rgb(248, 135, 142)')
        ),
        error_y=dict(
            type='data',
            array=data_std_div_sin,
            visible=True,
            color = ('rgb(248, 135, 142)')
        ))

data = [data1, data2, data3, data4]

layout = go.Layout(
    width=580,
    height=450,
    margin=go.layout.Margin(
        l=80,
        r=80,
        b=60,
        t=60,
    ),
    annotations=[
        dict(
            x=0.5004254919715793,
            y=-0.14,
            showarrow=False,
            text='Excitation Flip Angle (<i>θ<sub>n</sub></i>)',
            font=dict(
                family='Times New Roman',
                size=22,
                color=('rgb(21, 91, 158)')
            ),
            xref='paper',
            yref='paper'
        ),
        dict(
            x=-0.17,
            y=0.5,
            showarrow=False,
            text='Signal (<i>S<sub>n</sub></i>)',
            font=dict(
                family='Times New Roman',
                size=22,
                color=('rgb(21, 91, 158)')
            ),
            textangle=-90,
            xref='paper',
            yref='paper'
        ),
        dict(
            x=0.5004254919715793,
            y=1.15,
            showarrow=False,
            text='<i>S<sub>n</sub></i> / tan(<i>θ<sub>n</sub></i>)',
            font=dict(
                family='Times New Roman',
                size=22,
                color=('rgb(169, 10, 20)') 
            ),
            xref='paper',
            yref='paper'
        ),
        dict(
            x=1.16,
            y=0.5,
            showarrow=False,
            text='<i>S<sub>n</sub></i> / sin(<i>θ<sub>n</sub></i>)',
            font=dict(
                family='Times New Roman',
                size=22,
                color=('rgb(169, 10, 20)') 
            ),
            xref='paper',
            yref='paper',
            textangle=-90,
        ),
    ],
    xaxis=dict(
        autorange=False,
        range=[params['EXC_FA'][0], params['EXC_FA'][-1]],
        showgrid=False,
        linecolor='black',
        linewidth=2
    ),
    yaxis=dict(
        autorange=True,
        showgrid=False,
        linecolor='black',
        linewidth=2
    ),
    xaxis2=dict(
        autorange=False,
        range=[0, 1],
        showgrid=False,
        mirror=True,
        overlaying= 'x',
        anchor= 'y2',
        side= 'top',
        linecolor='black',
        linewidth=2
    ),
    yaxis2=dict(
        autorange=False,
        range=[0, 1],
        showgrid=False,
        overlaying= 'y',
        anchor= 'x',
        side= 'right',
        linecolor='black',
        linewidth=2
    ),
    legend=dict(
        x=0.32,
        y=0.98,
        traceorder='normal',
        font=dict(
            family='Times New Roman',
            size=12,
            color='#000'
        ),
        bordercolor='#000000',
        borderwidth=2
    ), 
)

fig = dict(data=data, layout=layout)

plot(fig, filename = 'vfa_fig_5.html', config = config)
display(HTML('vfa_fig_5.html'))

(View plot code)