Estimating diffusion time dependent q-space indices using qt-dMRI

from dipy.data.fetcher import (fetch_qtdMRI_test_retest_2subjects,
                               read_qtdMRI_test_retest_2subjects)
from dipy.reconst import qtdmri, dti
import matplotlib.pyplot as plt
import numpy as np

fetch_qtdMRI_test_retest_2subjects()
data, cc_masks, gtabs = read_qtdMRI_test_retest_2subjects()

plt.figure()
qtdmri.visualise_gradient_table_G_Delta_rainbow(gtabs[0])
plt.savefig('qt-dMRI_acquisition_scheme.png')

subplot_titles = ["Subject1 Test", "Subject1 Retest",
                  "Subject2 Test", "Subject2 Retest"]
fig = plt.figure()
plt.subplots(nrows=2, ncols=2)
for i, (data_, mask_, gtab_) in enumerate(zip(data, cc_masks, gtabs)):
    # take the middle slice
    data_middle_slice = data_[:, :, 2]
    mask_middle_slice = mask_[:, :, 2]

    # estimate fractional anisotropy (FA) for this slice
    tenmod = dti.TensorModel(gtab_)
    tenfit = tenmod.fit(data_middle_slice, data_middle_slice[..., 0] > 0)
    fa = tenfit.fa

    # set mask color to green with 0.5 opacity as overlay
    mask_template = np.zeros(np.r_[mask_middle_slice.shape, 4])
    mask_template[mask_middle_slice == 1] = np.r_[0., 1., 0., .5]

    # produce the FA images with corpus callosum masks.
    plt.subplot(2, 2, 1 + i)
    plt.title(subplot_titles[i], fontsize=15)
    plt.imshow(fa, cmap='Greys_r', origin='lower', interpolation='nearest')
    plt.imshow(mask_template, origin='lower', interpolation='nearest')
    plt.axis('off')
plt.tight_layout()
plt.savefig('qt-dMRI_datasets_fa_with_ccmasks.png')

tau_min = gtabs[0].tau.min()
tau_max = gtabs[0].tau.max()
taus = np.linspace(tau_min, tau_max, 5)

qtdmri_fits = []
msds = []
rtops = []
rtaps = []
rtpps = []
for i, (data_, mask_, gtab_) in enumerate(zip(data, cc_masks, gtabs)):
    # select the corpus callosum voxel for every dataset
    cc_voxels = data_[mask_ == 1]
    # initialize the qt-dMRI model.
    # recommended basis orders are radial_order=6 and time_order=2.
    # The combined Laplacian and l1-regularization using Generalized
    # Cross-Validation (GCV) and Cross-Validation (CV) settings is most robust,
    # but can be used separately and with weightings preset to any positive
    # value to optimize for speed.
    qtdmri_mod = qtdmri.QtdmriModel(
        gtab_, radial_order=6, time_order=2,
        laplacian_regularization=True, laplacian_weighting='GCV',
        l1_regularization=True, l1_weighting='CV'
    )
    # fit the model.
    # Here we take every 5th voxel for speed, but of course all voxels can be
    # fit for a more robust result later on.
    qtdmri_fit = qtdmri_mod.fit(cc_voxels[::5])
    qtdmri_fits.append(qtdmri_fit)
    # We estimate MSD, RTOP, RTAP and RTPP for the chosen diffusion times.
    msds.append(np.array(list(map(qtdmri_fit.msd, taus))))
    rtops.append(np.array(list(map(qtdmri_fit.rtop, taus))))
    rtaps.append(np.array(list(map(qtdmri_fit.rtap, taus))))
    rtpps.append(np.array(list(map(qtdmri_fit.rtpp, taus))))

def plot_mean_with_std(ax, time, ind1, plotcolor, ls='-', std_mult=1,
                       label=''):
    means = np.mean(ind1, axis=1)
    stds = np.std(ind1, axis=1)
    ax.plot(time, means, c=plotcolor, lw=3, label=label, ls=ls)
    ax.fill_between(time,
                    means + std_mult * stds,
                    means - std_mult * stds,
                    alpha=0.15, color=plotcolor)
    ax.plot(time, means + std_mult * stds, alpha=0.25, color=plotcolor)
    ax.plot(time, means - std_mult * stds, alpha=0.25, color=plotcolor)

# we first generate the data to produce the background index isolines.
Delta_ = np.linspace(0.005, 0.02, 100)
MSD_ = np.linspace(4e-5, 10e-5, 100)
Delta_grid, MSD_grid = np.meshgrid(Delta_, MSD_)
D_grid = MSD_grid / (6 * Delta_grid)
D_levels = np.r_[1, 5, 7, 10, 14, 23, 30] * 1e-4

fig = plt.figure(figsize=(10, 3))
# start with the plot of subject 1.
ax = plt.subplot(1, 2, 1)
# first plot the background
plt.contourf(Delta_ * 1e3, 1e5 * MSD_, D_grid, levels=D_levels, cmap='Greys',
             alpha=.5)

# plot the test-retest mean MSD and standard deviation of subject 1.
plot_mean_with_std(ax, taus * 1e3, 1e5 * msds[0], 'r', 'dashdot',
                   label='MSD Test')
plot_mean_with_std(ax, taus * 1e3, 1e5 * msds[1], 'g', 'dashdot',
                   label='MSD Retest')
ax.legend(fontsize=13)
# plot some text markers to clarify the background diffusivity lines.
ax.text(.0091 * 1e3, 6.33, 'D=14e-4', fontsize=12, rotation=35)
ax.text(.0091 * 1e3, 4.55, 'D=10e-4', fontsize=12, rotation=25)
ax.set_ylim(4, 9.5)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_title(r'Test-Retest MSD(:math:`\tau`) Subject 1', fontsize=15)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_ylabel('MSD (:math:`10^{-5}mm^2`)', fontsize=17)

# then do the same thing for subject 2.
ax = plt.subplot(1, 2, 2)
plt.contourf(Delta_ * 1e3, 1e5 * MSD_, D_grid, levels=D_levels, cmap='Greys',
             alpha=.5)
cb = plt.colorbar()
cb.set_label('Free Diffusivity (:math:`mm^2/s`)', fontsize=18)

plot_mean_with_std(ax, taus * 1e3, 1e5 * msds[2], 'r', 'dashdot')
plot_mean_with_std(ax, taus * 1e3, 1e5 * msds[3], 'g', 'dashdot')
ax.set_ylim(4, 9.5)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_title(r'Test-Retest MSD(:math:`\tau`) Subject 2', fontsize=15)
plt.savefig('qt_indices_msd.png')

# Again, first we define the data for the background illustration.
Delta_ = np.linspace(0.005, 0.02, 100)
RTXP_ = np.linspace(1, 200, 100)
Delta_grid, RTXP_grid = np.meshgrid(Delta_, RTXP_)
D_grid = 1 / (4 * RTXP_grid ** 2 * np.pi * Delta_grid)
D_levels = np.r_[1, 2, 3, 4, 6, 9, 15, 30] * 1e-4
D_colors = np.tile(np.linspace(.8, 0, 7), (3, 1)).T

# We start with estimating the RTOP illustration.
fig = plt.figure(figsize=(10, 3))
ax = plt.subplot(1, 2, 1)
plt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,
             alpha=.5)

plot_mean_with_std(ax, taus * 1e3, rtops[0] ** (1 / 3.), 'r', '--',
                   label='RTOP:math:`^{1/3}` Test')
plot_mean_with_std(ax, taus * 1e3, rtops[1] ** (1 / 3.), 'g', '--',
                   label='RTOP:math:`^{1/3}` Retest')
ax.legend(fontsize=13)
ax.text(.0091 * 1e3, 162, 'D=3e-4', fontsize=12, rotation=-22)
ax.text(.0091 * 1e3, 140, 'D=4e-4', fontsize=12, rotation=-20)
ax.text(.0091 * 1e3, 113, 'D=6e-4', fontsize=12, rotation=-16)
ax.set_ylim(54, 170)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_title(r'Test-Retest RTOP(:math:`\tau`) Subject 1', fontsize=15)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_ylabel('RTOP:math:`^{1/3}` (1/mm)', fontsize=17)

ax = plt.subplot(1, 2, 2)
plt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,
             alpha=.5)
cb = plt.colorbar()
cb.set_label('Free Diffusivity (:math:`mm^2/s`)', fontsize=18)

plot_mean_with_std(ax, taus * 1e3, rtops[2] ** (1 / 3.), 'r', '--')
plot_mean_with_std(ax, taus * 1e3, rtops[3] ** (1 / 3.), 'g', '--')
ax.set_ylim(54, 170)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_title(r'Test-Retest RTOP(:math:`\tau`) Subject 2', fontsize=15)
plt.savefig('qt_indices_rtop.png')

# First, we estimate the RTAP trends.
fig = plt.figure(figsize=(10, 3))
ax = plt.subplot(1, 2, 1)
plt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,
             alpha=.5)

plot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[0]), 'r', '-',
                   label='RTAP:math:`^{1/2}` Test')
plot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[1]), 'g', '-',
                   label='RTAP:math:`^{1/2}` Retest')
ax.legend(fontsize=13)
ax.text(.0091 * 1e3, 162, 'D=3e-4', fontsize=12, rotation=-22)
ax.text(.0091 * 1e3, 140, 'D=4e-4', fontsize=12, rotation=-20)
ax.text(.0091 * 1e3, 113, 'D=6e-4', fontsize=12, rotation=-16)
ax.set_ylim(54, 170)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_title(r'Test-Retest RTAP(:math:`\tau`) Subject 1', fontsize=15)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_ylabel('RTAP:math:`^{1/2}` (1/mm)', fontsize=17)

ax = plt.subplot(1, 2, 2)
plt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,
             alpha=.5)
cb = plt.colorbar()
cb.set_label('Free Diffusivity (:math:`mm^2/s`)', fontsize=18)

plot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[2]), 'r', '-')
plot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[3]), 'g', '-')
ax.set_ylim(54, 170)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_title(r'Test-Retest RTAP(:math:`\tau`) Subject 2', fontsize=15)
plt.savefig('qt_indices_rtap.png')


# Finally the last one for RTPP.
fig = plt.figure(figsize=(10, 3))
ax = plt.subplot(1, 2, 1)
plt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,
             alpha=.5)

plot_mean_with_std(ax, taus * 1e3, rtpps[0], 'r', ':', label='RTPP Test')
plot_mean_with_std(ax, taus * 1e3, rtpps[1], 'g', ':', label='RTPP Retest')
ax.legend(fontsize=13)
ax.text(.0091 * 1e3, 113, 'D=6e-4', fontsize=12, rotation=-16)
ax.text(.0091 * 1e3, 91, 'D=9e-4', fontsize=12, rotation=-13)
ax.text(.0091 * 1e3, 69, 'D=15e-4', fontsize=12, rotation=-10)
ax.set_ylim(54, 170)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_title(r'Test-Retest RTPP(:math:`\tau`) Subject 1', fontsize=15)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_ylabel('RTPP (1/mm)', fontsize=17)

ax = plt.subplot(1, 2, 2)
plt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,
             alpha=.5)
cb = plt.colorbar()
cb.set_label('Free Diffusivity (:math:`mm^2/s`)', fontsize=18)

plot_mean_with_std(ax, taus * 1e3, rtpps[2], 'r', ':')
plot_mean_with_std(ax, taus * 1e3, rtpps[3], 'g', ':')
ax.set_ylim(54, 170)
ax.set_xlim(.009 * 1e3, 0.0185 * 1e3)
ax.set_xlabel('Diffusion Time (ms)', fontsize=17)
ax.set_title(r'Test-Retest RTPP(:math:`\tau`) Subject 2', fontsize=15)
plt.savefig('qt_indices_rtpp.png')