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]
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}
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'])
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)