This notebook demonstrates the fitting of the RS-SFM model to the data from one voxel.
We start by importing general python tools for analysis and visualization.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
We use the scikit-learn
implementation of the Elastic Net algorithm
import sklearn.linear_model as lm
We use functions from dipy
to calculate model goodness-of-fit, and to represent the gradient table:
from dipy.reconst.cross_validation import coeff_of_determination
from dipy.core import gradients as grad
Finally, we import the local utility functions for this data set, as well as our implementation of the RS-SFM model (which extends the implementation that exists in dipy
)
import utils
from model import Model, BiExponentialIsotropicModel
A utility function reads the data from the files:
data = utils.read_data()
The next cell can be changed to choose data from another voxel (there are 6 voxels from each one of the ROIs).
vox_idx = 0
The fornix ROI can be chosen by changing seen
to seenX
in the following cell
signal = data['seen']['signal'][:, vox_idx]
bvals = data['seen']['bvals']
bvecs = data['seen']['bvecs']
delta = data['seen']['delta']
Delta = data['seen']['Delta']
te = data['seen']['TE']
g = data['seen']['g']
This is the measurement gradient table (b-values and b-vectors)
gtab = grad.gradient_table(bvals,
bvecs,
big_delta=Delta,
small_delta=delta)
The following sets up the model, specifying the measurement scheme (gtab
), as well as the model of the isotropic signal (BiExponentialIsotropicModel
) and the ElasticNet
solver parameters:
model = Model(gtab,
isotropic=BiExponentialIsotropicModel,
solver=lm.ElasticNet(l1_ratio=0.4,
alpha=5e-7,
positive=True,
warm_start=True,
fit_intercept=True,
normalize=True))
We fit the model with the signal provided for this voxel, and with TE and G, used to fit T2 effects, and to weight the measurements
fit = model.fit(signal, te, g)
/Users/arokem/anaconda/lib/python3.5/site-packages/sklearn/linear_model/coordinate_descent.py:466: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations ConvergenceWarning)
We predict back the signal using the derived model parameters
predict = fit.predict(gtab, te)
The first visualization provides an overall sense of the goodness of fit, and calculates the coefficient of determination ($R^2$) of the fit
fig, ax = plt.subplots(1)
plt.scatter(predict, signal)
plt.plot([0, 500], [0, 500], 'k--')
print('R-squared = %s'%coeff_of_determination(signal, predict))
R-squared = 98.7835216409
The next visualization examines the fit in different gradient strength settings. Each plot corresponds to a gradient strength range. Data is plotted so that the signal measured is on the left side of each plot, and the signal predicted is on the right side. The data are plotted as a function of the angle of the direction of the diffusion gradient for each data-point, relative to the sagittal direction ([1, 0, 0]). Different color points in each figure come from measurements with different diffusion weighting strength (b-value).
g_tiers = [0, 85, 185, 265, 350]
g_labels = []
for g_group in range(len(g_tiers)-1):
g_idx = np.where(np.logical_and(g > g_tiers[g_group], g < g_tiers[g_group+1]))
this_predict = predict[g_idx]
this_signal = signal[g_idx]
this_bvec = bvecs[g_idx]
this_bval = bvals[g_idx]
x = []
y1 = []
y2 = []
b = []
for idx, bv in enumerate(this_bvec):
if this_bval[idx]>0:
x.append(np.abs(np.dot(bv, [1, 0, 0])))
b.append(this_bval[idx])
y1.append(this_predict[idx])
y2.append(this_signal[idx])
fig, ax = plt.subplots(1)
pal = np.array(sns.color_palette('spectral', n_colors=len(np.unique(b))))
for ii, bb in enumerate(np.unique(b)):
bb_idx = np.where(b==bb)
ax.scatter(np.array(x)[bb_idx], np.array(y1)[bb_idx], color=pal[np.mod(ii, len(pal))])
ax.scatter(-np.array(x)[bb_idx], np.array(y2)[bb_idx], color=pal[np.mod(ii, len(pal))])
# This one's just so we can get a nice legend:
ax.plot(-1.2, 0, '.', color=pal[np.mod(ii, len(pal))], label='%s'%np.int(bb))
ax.set_xlabel(' Data Prediction\n Angle (relative to saggital)')
ax.set_ylabel('MRI signal (a.u.)')
ax.set_ylim([-10, None])
ax.set_xlim([-1.1, 1.1])
ax.set_xticklabels([0, 90, 45, 0, 45, 90])
ax.set_title('Gradient strength (average) = %3.2f mT/m'%np.mean(g[g_idx]))
plt.legend(loc=1,bbox_to_anchor=(1.2, 1), title='b-value')