T1 accuracy & inter-vendor agreement¶

T1 plate of the ISMRM/NIST system phantom was scanned with the following T1 mapping protocol:

Click here to reveal the variable flip angle (VFA) protocol.

Parameter (PDw/T1w)	G1_NATIVE	S1_NATIVE	S2_NATIVE	VENUS
Sequence Name	SPGR	FLASH	FLASH	mt_sat v1.0
Flip Angle (°)	6/20	6/20	6/20	6/20
TR (ms)	32/18	32/18	32/18	32/18
TE (ms)	4	4	4	4
FOV (cm)	25.6	25.6	25.6	25.6
Acquisition Matrix	256x256	256x256	256x256	256x256
Receiver BW (kHz)	62.5	62.5	62.5	62.5
RF Phase Increment (°)	115.4	50	50	117

The nominal T1 values (s) of the system phantom (SN = 42) were 1.98, 1.45, 0.98, 0.71, 0.5, 0.35, 0.24, 0.17, 0.13, 0.09 in 10 reference spheres (R1- R10) of the phantom. These values are shown using white cross markers throughout this notebook.

Note

The following code cell loads all the necessary Python packages, defines common variables and downloads the data to /tmp/VENUS directory. If data already exists, the download will be skipped.

import plotly.graph_objs as go
import plotly.express as px
from plotly import tools, subplots
import numpy as np
import pandas as pd
import ipywidgets as widgets
import math
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from ipywidgets import interact, interactive, fixed, interact_manual
init_notebook_mode(connected=True)
from IPython.core.display import display, HTML
import scipy.io as sio
import random
random.seed(123)

from repo2data.repo2data import Repo2Data
import os 
data_req_path = os.path.join("..", "binder", "data_requirement.json")
repo2data = Repo2Data(data_req_path)
data_path = repo2data.install()[0]
derivativesDir = os.path.join(data_path,"qMRFlow")

spheres = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
refMean = [1.989,1.454,0.9841,0.706,0.4967,0.3515,0.24713,0.1753,0.1259,0.089]
sessions = ['rth750rev','rthPRIrev','rthSKYrev',
           'vendor750rev','vendorPRIrev','vendorSKYrev']
convert_tag = {'rth750':'G1<sub>neutral</sub>','rthPRI':'S1<sub>neutral</sub>','rthSKY':'S2<sub>neutral</sub>',
               'vendor750':'G1<sub>native</sub>','vendorPRI':'S1<sub>native</sub>','vendorSKY':'S2<sub>native</sub>'
              }
colors = {'rth750':'rgba(255,8,30,1)','rthPRI':'rgba(255,117,0,1)','rthSKY':'rgba(255,220,0,1)',
         'vendor750':'rgba(13,186,233,1)','vendorPRI':'rgba(2.5,187,180,1)','vendorSKY':'rgba(0,248,205,1)'}
colors_line = {'rth750':'rgba(255,8,30,0.2)','rthPRI':'rgba(255,117,0,0.2)','rthSKY':'rgba(255,220,0,0.2)',
         'vendor750':'rgba(13,186,233,0.2)','vendorPRI':'rgba(2.5,187,180,0.2)','vendorSKY':'rgba(0,248,205,0.2)'}
line_dash = {'rth750':'dash','rthPRI':'dash','rthSKY':'dash',
         'vendor750':'solid','vendorPRI':'solid','vendorSKY':'solid'}
line_sym = {'rth750':'circle','rthPRI':'circle','rthSKY':'circle',
         'vendor750':'square','vendorPRI':'square','vendorSKY':'square'}

---- repo2data starting ----
/Users/agah/opt/anaconda3/envs/jbnew/lib/python3.6/site-packages/repo2data
Config from file :
../binder/data_requirement.json
Destination:
/tmp/VENUS

Info : /tmp/VENUS already downloaded

Peak SNR values¶

Note

Reproduces Figure 3a from the article.

Calculations of signal (average value of the highest signal sphere) and noise (background standard deviation) were performed manually using 3D Slicer.

VENUS PSNR values are on a par with those of vendor-native T1w and PDw images (Fig. 3a).

colors = {'rth750':'rgba(255,8,30,1)','rthPRI':'rgba(255,117,0,1)','rthSKY':'rgba(255,220,0,1)',
         'vendor750':'rgba(0,75,255,1)','vendorPRI':'rgba(24,231,234,1)','vendorSKY':'rgba(24,234,141,1)'}
df = pd.DataFrame()
df['signal'] = [1649,240,264,793,2330,2379,767,110,117,362,1023,1016]
df['noise'] = [1.68,0.34,0.33,0.66,1.64,1.65,1.66,0.32,0.34,0.6,1.72,1.71]
df['session'] = ["rth750","rthPRI","rthSKY","vendor750","vendorPRI","vendorSKY"]*2
df['labels'] = ["G1<sub>VENUS</sub>","S1<sub>VENUS</sub>","S2<sub>VENUS</sub>","G1<sub>NATIVE</sub>","S1<sub>NATIVE</sub>","S2<sub>NATIVE</sub>"]*2
df['acq'] = ["T1w"]*6 + ["PDw"]*6 
df['SNR'] = 10*np.log10(df['signal']/df['noise'])
fig = go.Figure()
base = np.core.defchararray.add(['<br>']*6,df['acq'][0:6])
text = np.core.defchararray.add(list(map(str,np.round(df['SNR'][0:6],2))), base)
fig.add_trace(go.Bar(
        y=df['SNR'][0:6],
        x=df["labels"][0:6],
        showlegend=False,
        name="T1w",
        text=text,
        textposition='auto',
        marker_color= [colors[x] for x in df['session']]))
#aa = list(map(str,np.round(df['SNR'][6:12],2)))
base = np.core.defchararray.add(['<br>']*6,df['acq'][6:12])
text = np.core.defchararray.add(list(map(str,np.round(df['SNR'][6:12],2))), base)
fig.add_trace(go.Bar(
        y=df['SNR'][6:12],
        x=df["labels"][6:12],
        textposition='auto',
        showlegend=False,
        text=text,
        name="PDw",
        marker_color=[colors[x] for x in df['session']]))
fig.update_layout(barmode='group')
fig.update_layout(height=300, width=900)
fig.update_yaxes(title="PSNR")

fig.show()

VENUS vs NATIVE `T1` estimations¶

Note

Reproduces Figure 4b from the article.

T1 values from the vendor-native acquisitions are represented by solid lines and square markers in cold colors, and those from VENUS attain dashed lines and circle markers in hot colors.

Vendor-native measurements, especially G1_NATIVE and S2_NATIVE, overestimate T1. G1_VENUS and S1-2_VENUS remain closer to the reference.

fig = go.Figure()
for jj in range(0,len(sessions),1):
    avg_array  = []
    std_array  =  []
    cur_tag = sessions[jj][0:-3]
    cur_csv = pd.read_csv(derivativesDir + '/sub-phantom/' + 'ses-' + sessions[jj] + '/stat/sub-phantom_ses-' + sessions[jj] + '_stat_summary.csv')
    avg_array = np.array(cur_csv['T1 (mean)'])
    std_array = np.array(cur_csv['T1 (std)'])
    upper = list(np.squeeze(np.array(avg_array)) + np.squeeze(np.array(std_array)))
    lower = list(np.squeeze(np.array(avg_array)) - np.squeeze(np.array(std_array)))
    # Add std lines
    fig.add_trace(go.Scatter(
        x= refMean + refMean[::-1],
        y=upper+lower[::-1],
        fill='toself',
        fillcolor= colors_line[cur_tag],
        line_color='rgba(255,0,0,0)',
        name= convert_tag[cur_tag]
    ))
    # Add markers
    fig.add_trace(go.Scatter(
        x=refMean, y=avg_array,
        line_color=colors[cur_tag],
        name= convert_tag[cur_tag],
        marker = dict(size=20,symbol=line_sym[cur_tag],line=dict(color='black',width=1.5),opacity=0.85),
        line=dict(width=2, dash=line_dash[cur_tag])
    ))

fig.update_layout(height=1000, width=800)
fig.update_layout(legend=dict(x=1, y=1,tracegroupgap=300,font=dict(color="black",size=14)))
axis_template = dict(linecolor = 'white', 
                     showticklabels = True,
                     tickfont=dict(color="black"), 
                     gridcolor = 'lightsteelblue',
                     tickvals = refMean[::-1])
fig.update_xaxes(axis_template,range=[0,2.1], title = "Reference T1 (s)")

fig.update_yaxes(axis_template,tickmode="array",tickvals=refMean,ticktext=["R10","R9","R8","R7","R6","R5","R4","R3","R2","R1"])
fig.update_layout(title='',title_font_color="white",margin=dict(l=0,r=0,t=0,b=0))
fig.update_layout(hovermode = "x")
fig.add_trace(go.Scatter(
    x=refMean, y=refMean,
    line_color='black',
    mode = 'markers',
    name = 'Reference',
    marker = dict(size=30,symbol='cross-open',line = dict(width=2))

))


fig.show()

Percent measurement errors (`∆T1`)¶

Note

Reproduces Figure 4c from the article.

For VENUS, ∆T1 remains low for the physiologically relevant range (0.7 to 2s), whereas deviations reach up to 30.4% for vendor-native measurements.

df = pd.DataFrame()
fig = go.Figure()
for jj in range(0,len(sessions),1):
    cur_csv = pd.read_csv(derivativesDir + '/sub-phantom/' + 'ses-' + sessions[jj] + '/stat/sub-phantom_ses-' + sessions[jj] + '_stat_summary.csv')
    avg_array = np.array(cur_csv['T1 (mean)'])
    std_array = np.array(cur_csv['T1 (std)'])
    cur_tag = sessions[jj][0:-3]
    # Measurement error
    pcterr = ((np.array(avg_array) - np.array(refMean))/np.array(refMean))*100
    fig.add_trace(go.Scatter(
        x=pcterr, y=spheres[::-1],
        line_color=colors[cur_tag],
        name= convert_tag[cur_tag],
        marker = dict(size=20,symbol=line_sym[cur_tag],line=dict(color='black',width=1.5),opacity=0.75),
        line=dict(width=2, dash=line_dash[cur_tag]),
        line_shape='vhv'
    ))
fig.update_layout(height=1000, width=800)
fig.update_layout(legend=dict(x=1, y=1,tracegroupgap=300,font=dict(color="black")),margin=dict(l=0,r=0,t=0,b=0))
axis_template = dict(linecolor = 'black', 
                     showticklabels = True,
                     gridcolor = 'lightsteelblue'
                    )
fig.update_xaxes(axis_template,zerolinecolor="black",zerolinewidth=3)
fig.update_yaxes(axis_template,zerolinecolor="black")
fig.update_layout(hovermode = "y")
fig.update_xaxes(range=[-18,70],title="∆T1 (%)")
fig.update_yaxes(tickmode="array",tickvals=spheres[::-1],ticktext=["R10","R9","R8","R7","R6","R5","R4","R3","R2","R1"],title = "Reference Spheres (R1 - R10)")
fig.update_layout(title_font_color="black")

fig.show()

Averaged `∆T1` comparison¶

Note

Reproduces Figure 4d from the article.

T1 values are averaged over S1-2 (SNATIVE and SVENUS, green square and orange circle) and according to the acquisition type (NATIVE and VENUS, black square and black circle). Inter-vendor percent differences are displayed on hover.

In addition to the prominent improvement in G1 accuracy, S_VENUS is closer to the reference than S_NATIVE for most of the relevant range (∆T1 of 7.6, 3.5, 5.4, 0.7% and 3.2, 0.9, 2, 1.3% for S_NATIVE and S_VENUS, respectively).

In conclusion, VENUS reduces between-vendor differences with an overall accuracy improvement.

You can change the range (lastN) (up to 9) in the following code cell.

sessions = ['vendor750rev','vendorPRIrev','vendorSKYrev',
           'rth750rev','rthPRIrev','rthSKYrev']

# Set to R7-R10
lastN = 4

df = pd.DataFrame()
fig = go.Figure()

rth750 = derivativesDir + '/sub-phantom/' + 'ses-' + sessions[3] + '/stat/sub-phantom_ses-' + sessions[3] + '_stat_summary.csv'
rthPRI = derivativesDir + '/sub-phantom/' + 'ses-' + sessions[4] + '/stat/sub-phantom_ses-' + sessions[4] + '_stat_summary.csv'
rthSKY = derivativesDir + '/sub-phantom/' + 'ses-' + sessions[5] + '/stat/sub-phantom_ses-' + sessions[5] + '_stat_summary.csv'
rth750 = pd.read_csv(rth750)
rthPRI = pd.read_csv(rthPRI)
rthSKY = pd.read_csv(rthSKY)
    
avg_array_ge = rth750['T1 (mean)']
std_array_ge = rth750['T1 (std)']
cur_tag_ge = 'neutral G1'

avg_array_sie = np.mean([rthPRI['T1 (mean)'],rthSKY['T1 (mean)']],axis=0)
std_array_sie = np.mean([rthPRI['T1 (std)'],rthSKY['T1 (std)']],axis=0)
cur_tag = 'neutral S1-2'

avg_array_all = np.mean([rth750['T1 (mean)'],rthPRI['T1 (mean)'],rthSKY['T1 (mean)']],axis=0)
std_array_all = np.mean([rthPRI['T1 (std)'],rthSKY['T1 (std)']],axis=0)

df[cur_tag + '_mean'] = avg_array
df[cur_tag + '_std'] = std_array
upper_sie = list(np.squeeze(np.array(avg_array_sie)) + np.squeeze(np.array(std_array_sie)))
lower_sie = list(np.squeeze(np.array(avg_array_sie)) - np.squeeze(np.array(std_array_sie)))
upper_ge = list(np.squeeze(np.array(avg_array_ge)) + np.squeeze(np.array(std_array_ge)))
lower_ge = list(np.squeeze(np.array(avg_array_ge)) - np.squeeze(np.array(std_array_ge)))


fig.add_trace(go.Scatter(
    x=refMean[:lastN], y=refMean[:lastN],
    line_color='black',
    mode = 'markers',
    name = 'Reference',
    marker = dict(size=30,symbol='cross-open',line = dict(width=2))

))
fig.add_trace(go.Scatter(
    x=refMean[:lastN], y=avg_array_sie[:lastN],
    line_color='orange',
    name= 'neutral S1-2',
    marker = dict(size=20,symbol='circle',line=dict(color='black',width=1.5),opacity=0.9),
    line=dict(width=2, dash='dash'),
    line_shape='spline',
    customdata= np.round(((np.array(avg_array_sie[:lastN]) - np.array(refMean[:lastN]))/np.array(refMean[:lastN]))*100,2),
    hovertemplate='<br>%{y},<b>∆T1:</b>%{customdata}',
    hoverlabel = dict(namelength = -1)
))
fig.add_trace(go.Scatter(
    x=refMean[:lastN], y=avg_array_ge[:lastN],
    line_color='rgba(202,73,68,1)',
    name= 'neutral G1',
    marker = dict(size=20,symbol='circle',line=dict(color='black',width=1.5),opacity=0.9),
    line=dict(width=2, dash='dash'),
    line_shape='spline',
    customdata= np.round(((np.array(avg_array_ge[:lastN]) - np.array(refMean[:lastN]))/np.array(refMean[:lastN]))*100,2),
    hovertemplate='<br>%{y},<b>∆T1:</b>%{customdata}',
    hoverlabel = dict(namelength = -1)
))
fig.add_trace(go.Scatter(
    x=refMean[:lastN], y=avg_array_all[:lastN],
    line_color='rgba(0,0,0,1)',
    name= 'neutral',
    marker = dict(size=5,symbol='circle',line=dict(color='black',width=1.5),opacity=0.9),
    line=dict(width=2, dash='dash'),
    line_shape='spline',
    customdata= np.round(((np.array(avg_array_all[:lastN]) - np.array(refMean[:lastN]))/np.array(refMean[:lastN]))*100,2),
    hovertemplate='<br>%{y},<b>∆T1:</b>%{customdata}',
    hoverlabel = dict(namelength = -1)
))
vendor750 = derivativesDir + '/sub-phantom/' + 'ses-' + sessions[0] + '/stat/sub-phantom_ses-' + sessions[0] + '_stat_summary.csv'
vendorPRI = derivativesDir + '/sub-phantom/' + 'ses-' + sessions[1] + '/stat/sub-phantom_ses-' + sessions[1] + '_stat_summary.csv'
vendorSKY = derivativesDir + '/sub-phantom/' + 'ses-' + sessions[2] + '/stat/sub-phantom_ses-' + sessions[2] + '_stat_summary.csv'
vendor750 = pd.read_csv(vendor750)
vendorPRI = pd.read_csv(vendorPRI)
vendorSKY = pd.read_csv(vendorSKY)
    
avg_array_ge = vendor750['T1 (mean)']
std_array_ge = vendor750['T1 (std)']

avg_array_sie = np.mean([vendorPRI['T1 (mean)'],vendorSKY['T1 (mean)']],axis=0)
std_array_sie = np.mean([vendorPRI['T1 (std)'],vendorSKY['T1 (std)']],axis=0)
cur_tag = 'neutral'

avg_array_all = np.mean([vendor750['T1 (mean)'],vendorPRI['T1 (mean)'],vendorSKY['T1 (mean)']],axis=0)

df[cur_tag + '_mean'] = avg_array
df[cur_tag + '_std'] = std_array
upper_sie = list(np.squeeze(np.array(avg_array_sie)) + np.squeeze(np.array(std_array_sie)))
lower_sie = list(np.squeeze(np.array(avg_array_sie)) - np.squeeze(np.array(std_array_sie)))
upper_ge = list(np.squeeze(np.array(avg_array_ge)) + np.squeeze(np.array(std_array_ge)))
lower_ge = list(np.squeeze(np.array(avg_array_ge)) - np.squeeze(np.array(std_array_ge)))

fig.add_trace(go.Scatter(
    x=refMean[:lastN], y=avg_array_sie[:lastN],
    line_color='green',
    name= 'native S1-2',
    marker = dict(size=20,symbol='square',line=dict(color='black',width=1.5),opacity=0.9),
    line=dict(width=2, dash='solid'),
    line_shape='spline',
    customdata= np.round(((np.array(avg_array_sie[:lastN]) - np.array(refMean[:lastN]))/np.array(refMean[:lastN]))*100,2),
    hovertemplate='<br>%{y}, <b>∆T1:</b>%{customdata}',
    hoverlabel = dict(namelength = -1)

))
fig.add_trace(go.Scatter(
    x=refMean[:lastN], y=avg_array_ge[:lastN],
    line_color='rgba(31,119,180,0.8)',
    name= 'native G1',
    marker = dict(size=20,symbol='square',line=dict(color='black',width=1.5),opacity=0.9),
    line=dict(width=2, dash='solid'),
    line_shape='spline',
    customdata= np.round(((np.array(avg_array_ge[:lastN]) - np.array(refMean[:lastN]))/np.array(refMean[:lastN]))*100,2),
    hovertemplate='<br>%{y}, <b>∆T1:</b>%{customdata}',
    hoverlabel = dict(namelength = -1)
))
fig.add_trace(go.Scatter(
    x=refMean[:lastN], y=avg_array_all[:lastN],
    line_color='rgba(0,0,0,1)',
    name= 'native',
    marker = dict(size=5,symbol='square',line=dict(color='black',width=1.5),opacity=0.9),
    line=dict(width=2, dash='solid'),
    line_shape='spline',
    customdata= np.round(((np.array(avg_array_all[:lastN]) - np.array(refMean[:lastN]))/np.array(refMean[:lastN]))*100,2),
    hovertemplate='<br>%{y},<b>∆T1:</b>%{customdata}',
    hoverlabel = dict(namelength = -1)
))
fig.update_layout(height=1100, width=800)
axis_template = dict(linecolor = 'black', 
                     showticklabels = True,
                     tickfont=dict(color="black"), 
                     gridcolor = 'lightsteelblue',
                     tickvals = refMean[::-1],)
tickvals=refMean
ticktext=["R10","R9","R8","R7","R6","R5","R4","R3","R2","R1"]

fig.update_xaxes(axis_template,title="Reference T1 (s)")
fig.update_yaxes(axis_template,tickmode="array",tickvals=tickvals[:lastN],ticktext=ticktext[:lastN],title="Reference Spheres")

fig.update_layout(hovermode = "x")


fig.show()

Vendor-neutral sequences (VENUS) and multi-center reproducibility of qMRI

T1 accuracy & inter-vendor agreement¶

Peak SNR values¶

VENUS vs NATIVE T1 estimations¶

Percent measurement errors (∆T1)¶

Averaged ∆T1 comparison¶

VENUS vs NATIVE `T1` estimations¶

Percent measurement errors (`∆T1`)¶

Averaged `∆T1` comparison¶