Notebook

Imports¶

In [1]:

from pathlib import Path
import pandas as pd
import json
import nibabel as nib
import numpy as np

from src.database import *
from src.nist import get_reference_NIST_values, get_NIST_ids
from src.tools import calc_error

import matplotlib.pyplot as plt
plt.style.use('custom_matplotlibrc')
plt.rcParams["figure.figsize"] = (10,10)
fig_id = 0

Configurations¶

In [2]:

database_path = Path('databases/3T_NIST_T1maps_general_fitting_model_database.pkl')
output_folder = Path("plots/06_alldatasets_scatter_NIST_general_fitting_model/")

estimate_type = 'median' # median or mean

Define functions¶

In [3]:

def plot_single_scatter(x, y, y_std,
                        title, x_label, y_label,
                        file_prefix, folder_path, fig_id,
                        y_type='linear'):
    if y_type is 'linear':
        plt.errorbar(x,y, y_std, fmt='o', solid_capstyle='projecting')
        ax = plt.gca()
        ax.axline((1, 1), slope=1, linestyle='dashed')
        ax.set_ylim(ymin=0)
        ax.set_xlim(xmin=0)
    if y_type is 'log':
        plt.loglog(x,y,'o')
        ax = plt.gca()
        ax.set_ylim(ymin=20)
        ax.set_xlim(xmin=20)
    if y_type is 'error_t1':
        plt.scatter(x,calc_error(x,y))
        ax = plt.gca()
        ax.axline((1, 0), slope=0, color='k')
        ax.axline((1, -10), slope=0, linestyle='dashed', color='k')
        ax.axline((1, 10), slope=0, linestyle='dashed', color='k')
        ax.set_ylim(ymin=-100, ymax=100)
    
    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)

    fig = plt.gcf()
    

    folder_path.mkdir(parents=True, exist_ok=True)

    if fig_id<10:
        filename = "0" + str(fig_id) + "_" + file_prefix
    else:
        filename = str(fig_id) + "_" + file_prefix

    fig.savefig(folder_path / (str(filename) + '.svg'), facecolor='white')
    fig.savefig(folder_path / (str(filename) + '.png'), facecolor='white')
    fig_id = fig_id + 1
    plt.show()
    return fig_id

<>:5: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:11: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:16: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:5: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:11: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:16: SyntaxWarning: "is" with a literal. Did you mean "=="?
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/2473446924.py:5: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'linear':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/2473446924.py:11: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'log':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/2473446924.py:16: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'error_t1':

Load database¶

In [4]:

df = pd.read_pickle(database_path)

Initialize variables, store mean and standard deviations of T1 in the ROIs¶

In [5]:

dataset_mean = np.zeros((1,14))
dataset_std = np.zeros((1,14))
version = np.array([])
ref_values = np.zeros((1,14))


ii=0
for index, row in df.iterrows():
    if type(df.loc[index]['T1 - NIST sphere 1']) is np.ndarray:

        version = np.append(version,df.loc[index]['phantom serial number'])
        if version[ii] is None:
            version[ii] = 999 # Missing version, only known case is one where we have version > 42 right now.

        if ii==0:
            ref_values = get_reference_NIST_values(version[ii])
        else:
            ref_values = np.vstack((ref_values, get_reference_NIST_values(version[ii])))
        
        tmp_dataset_estimate = np.array([])
        tmp_dataset_std = np.array([])

        for key in get_NIST_ids():
            if estimate_type is 'mean':
                tmp_dataset_estimate = np.append(tmp_dataset_estimate, np.mean(df.loc[index][key]))
            elif estimate_type is 'median':
                tmp_dataset_estimate = np.append(tmp_dataset_estimate, np.median(df.loc[index][key]))
            else:
                Exception('Unsupported dataset estimate type.')

            tmp_dataset_std = np.append(tmp_dataset_std, np.std(df.loc[index][key]))

        if ii==0:
            dataset_estimate = tmp_dataset_estimate  
            dataset_std = tmp_dataset_std
        else:
            dataset_estimate = np.vstack((dataset_estimate, tmp_dataset_estimate))
            dataset_std = np.vstack((dataset_std, tmp_dataset_std))

        ii=ii+1

<>:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:26: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:26: SyntaxWarning: "is" with a literal. Did you mean "=="?
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/21101739.py:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if estimate_type is 'mean':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/21101739.py:26: SyntaxWarning: "is" with a literal. Did you mean "=="?
  elif estimate_type is 'median':

Setup for plots¶

In [6]:

fig_id = 0
dims=ref_values.shape
file_prefix = 'alldatasets'

Plot figures¶

In [7]:

y_type = 'linear'
title = 'Linear axes'


for ii in range(dims[0]):
    
    x=ref_values[ii,:]
    y=dataset_estimate[ii,:]
    y_std=dataset_std[ii,:]
    x_label='Reference T$_1$ (ms)'
    y_label="T$_1$ estimate (ms)"
    folder_path=output_folder

    if y_type is 'linear':
        plt.errorbar(x,y, y_std, fmt='o', solid_capstyle='projecting')
        ax = plt.gca()
        ax.axline((1, 1), slope=1, linestyle='dashed')
        ax.set_ylim(ymin=0, ymax=2500)
        ax.set_xlim(xmin=0, xmax=2500)
    if y_type is 'log':
        plt.loglog(x,y,'o')
        ax = plt.gca()
        ax.set_ylim(ymin=20, ymax=2500)
        ax.set_xlim(xmin=20, xmax=2500)
    if y_type is 'error_t1':
        plt.errorbar(x,calc_error(x,y), fmt='o')
        ax = plt.gca()
        ax.axline((1, 0), slope=0, color='k')
        ax.axline((1, -10), slope=0, linestyle='dashed', color='k')
        ax.axline((1, 10), slope=0, linestyle='dashed', color='k')
        ax.set_ylim(ymin=-100, ymax=100)
        ax.set_xlim(xmin=0, xmax=2500)


    plt.xlabel(x_label)
    plt.ylabel(y_label)

plt.title(title, fontsize=30)
fig = plt.gcf()
    

folder_path.mkdir(parents=True, exist_ok=True)


if fig_id<10:
    filename = "0" + str(fig_id) + "_" + file_prefix
else:
    filename = str(fig_id) + "_" + file_prefix

fig.savefig(folder_path / (str(filename) + '.svg'), facecolor='white')
fig.savefig(folder_path / (str(filename) + '.png'), facecolor='white')
fig_id = fig_id + 1

plt.show()

<>:14: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:20: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:25: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:14: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:20: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:25: SyntaxWarning: "is" with a literal. Did you mean "=="?
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/428255775.py:14: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'linear':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/428255775.py:20: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'log':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/428255775.py:25: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'error_t1':

In [8]:

y_type = 'log'
title = 'Log-log axes'

for ii in range(dims[0]):
    
    x=ref_values[ii,:]
    y=dataset_estimate[ii,:]
    y_std=dataset_std[ii,:]
    x_label='Reference T$_1$ (ms)'
    y_label="T$_1$ estimate (ms)"
    folder_path=output_folder

    if y_type is 'linear':
        plt.errorbar(x,y, y_std, fmt='o', solid_capstyle='projecting')
        ax = plt.gca()
        ax.axline((1, 1), slope=1, linestyle='dashed')
        ax.set_ylim(ymin=0, ymax=2500)
        ax.set_xlim(xmin=0, xmax=2500)
    if y_type is 'log':
        plt.loglog(x,y,'o')
        ax = plt.gca()
        ax.set_ylim(ymin=20, ymax=2500)
        ax.set_xlim(xmin=20, xmax=2500)
    if y_type is 'error_t1':
        plt.errorbar(x,calc_error(x,y), fmt='o')
        ax = plt.gca()
        ax.axline((1, 0), slope=0, color='k')
        ax.axline((1, -10), slope=0, linestyle='dashed', color='k')
        ax.axline((1, 10), slope=0, linestyle='dashed', color='k')
        ax.set_ylim(ymin=-100, ymax=100)
        ax.set_xlim(xmin=0, xmax=2500)


    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)

fig = plt.gcf()
    

folder_path.mkdir(parents=True, exist_ok=True)


if fig_id<10:
    filename = "0" + str(fig_id) + "_" + file_prefix
else:
    filename = str(fig_id) + "_" + file_prefix

fig.savefig(folder_path / (str(filename) + '.svg'), facecolor='white')
fig.savefig(folder_path / (str(filename) + '.png'), facecolor='white')
fig_id = fig_id + 1
plt.show()

<>:13: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:19: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:13: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:19: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/398414603.py:13: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'linear':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/398414603.py:19: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'log':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/398414603.py:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'error_t1':

In [9]:

y_type = 'error_t1'
title = 'Error relative to reference'

for ii in range(dims[0]):
    
    x=ref_values[ii,:]
    y=dataset_estimate[ii,:]
    y_std=dataset_std[ii,:]
    x_label='Reference T$_1$ (ms)'
    y_label="Error (%)"
    folder_path=output_folder

    if y_type is 'linear':
        plt.errorbar(x,y, y_std, fmt='o', solid_capstyle='projecting', capsize=5)
        ax = plt.gca()
        ax.axline((1, 1), slope=1, linestyle='dashed')
        ax.set_ylim(ymin=0, ymax=2500)
        ax.set_xlim(xmin=0, xmax=2500)
    if y_type is 'log':
        plt.loglog(x,y,'o')
        ax = plt.gca()
        ax.set_ylim(ymin=20, ymax=2500)
        ax.set_xlim(xmin=20, xmax=2500)
    if y_type is 'error_t1':
        plt.errorbar(x,calc_error(x,y), fmt='o')
        ax = plt.gca()
        ax.axline((1, 0), slope=0, color='k')
        ax.axline((1, -10), slope=0, linestyle='dashed', color='k')
        ax.axline((1, 10), slope=0, linestyle='dashed', color='k')
        ax.set_ylim(ymin=-100, ymax=100)
        ax.set_xlim(xmin=0, xmax=2500)

    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)
    plt.ylim(top=100)

fig = plt.gcf()
    

folder_path.mkdir(parents=True, exist_ok=True)


if fig_id<10:
    filename = "0" + str(fig_id) + "_" + file_prefix
else:
    filename = str(fig_id) + "_" + file_prefix

fig.savefig(folder_path / (str(filename) + '.svg'), facecolor='white')
fig.savefig(folder_path / (str(filename) + '.png'), facecolor='white')
fig_id = fig_id + 1
plt.show()

<>:13: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:19: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:13: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:19: SyntaxWarning: "is" with a literal. Did you mean "=="?
<>:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/3489514954.py:13: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'linear':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/3489514954.py:19: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'log':
/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_1578/3489514954.py:24: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if y_type is 'error_t1':

In [ ]: