In [1]:
%load_ext autoreload
%autoreload 2

In [2]:
import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)

In [3]:
import json
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import GPy
import dipy.reconst.dti as dti
import dipy.reconst.mapmri as mapmri
from diGP.preprocessing_pipelines import get_SPARC_train_and_test, preprocess_SPARC
from diGP.model import GaussianProcessModel
from diGP.evaluation import get_SPARC_metrics

%matplotlib inline
sns.set_style('dark')

In [4]:
with open('../config.json', 'r') as json_file:
    conf = json.load(json_file)
data_paths = conf['SPARC']['data_paths']
q_test_path = conf['SPARC']['q_test_path']

Load data to use for prediction.


In [5]:
source = 'gradient_20'
gtab, data, voxel_size = get_SPARC_train_and_test(data_paths[source], data_paths['goldstandard'], q_test_path)

In [6]:
from dipy.core.gradients import gradient_table
bval_threshold = 2500
cutoff = gtab['train'].bvals < bval_threshold
gtab['cutoff'] = gradient_table(bvals=gtab['train'].bvals[cutoff],
                                bvecs=gtab['train'].bvecs[cutoff])
S0_tenmodel = dti.TensorModel(gtab['cutoff'], return_S0_hat=True)
S0_tenfit = S0_tenmodel.fit(data['train'][:, :, None, cutoff])

b0_DTI = np.squeeze(S0_tenfit.S0_hat)
b0 = np.squeeze(data['train'][:, :, gtab['train'].b0s_mask])

#data['train'] /= b0_DTI[:, :, None]
#data['test'] /= b0_DTI[:, :, None]

GP predictions


In [7]:
gp_model = GaussianProcessModel(gtab['train'], spatial_dims=2,
                                q_magnitude_transform=lambda x:x**0.1, box_cox_lambda=None)

gp_fit = gp_model.fit(data['train'], voxel_size=voxel_size[0:2], retrain=True)

mu = gp_fit.predict(gtab['test'], spatial_shape=data['test'].shape[0:2], voxel_size=voxel_size[0:2], compute_var=False)

pred = {'GP': mu}

Reference predictions with dipy


In [8]:
tenmodel = dti.TensorModel(gtab['train'])
tenfit = tenmodel.fit(data['train'])

pred['DTI'] = tenfit.predict(gtab['test'])

In [9]:
map_model = mapmri.MapmriModel(gtab['train'], positivity_constraint=False, laplacian_weighting='GCV',
                               radial_order=6)
mapfit = map_model.fit(data['train'])

pred['MAP'] = mapfit.predict(gtab['test'])

Evaluation


In [10]:
idx = 50 * np.arange(9)
f, axs = plt.subplots(len(idx), 4)
f.set_figheight(len(idx)*2)
f.set_figwidth(10)
for i, index in enumerate(idx):
    axs[i, 0].imshow(data['test'][:, :, index], vmin=0, vmax=1, cmap="GnBu_r")
    axs[i, 0].axis('off')
    axs[i, 0].set_title('Ground truth,\nb = {}'.format(np.round(gtab['test'].bvals[index])))
    axs[i, 1].imshow(pred['GP'][:, :, index], vmin=0, vmax=1, cmap="GnBu_r")
    axs[i, 1].axis('off')
    axs[i, 1].set_title('GP prediction, \nb = {}'.format(np.round(gtab['test'].bvals[index])))
    axs[i, 2].imshow(pred['DTI'][:, :, index], vmin=0, vmax=1, cmap="GnBu_r")
    axs[i, 2].axis('off')
    axs[i, 2].set_title('DTI prediction, \nb = {}'.format(np.round(gtab['test'].bvals[index])))
    axs[i, 3].imshow(pred['MAP'][:, :, index], vmin=0, vmax=1, cmap="GnBu_r")
    axs[i, 3].axis('off')
    axs[i, 3].set_title('MAP-MRI prediction, \nb = {}'.format(np.round(gtab['test'].bvals[index])))

plt.tight_layout()



In [11]:
for key, value in pred.items():
    print("\n{} model:".format(key))
    get_SPARC_metrics(gtab['test'], data['test'], value, verbose=True)


MAP model:
NMSE low: 0.031040994634508445
NMSE high: 0.09041411861417692
NMSE all: 0.05479024422637584

GP model:
NMSE low: 0.030520581443648297
NMSE high: 0.15438421337450697
NMSE all: 0.08006603421599176

DTI model:
NMSE low: 0.07865178619130589
NMSE high: 0.24003453733688107
NMSE all: 0.14320488664953593