Skip to content

Commit

Permalink
Merge pull request #120 from kaitj/tractography_exercises
Browse files Browse the repository at this point in the history 
ENH: Update deterministic tractography
  • Loading branch information
@kaitj
kaitj authored Apr 12, 2021
2 parents 6c3f5b9 + 18d790f commit 9179bcd
Show file tree
Hide file tree
Showing 7 changed files with 724 additions and 25 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/build-test.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,5 +61,5 @@ jobs:
pytest --nbval-lax -v code/introduction/solutions/introduction_solutions.ipynb
pytest --nbval-lax -v code/diffusion_tensor_imaging/solutions/diffusion_tensor_imaging_solutions.ipynb
pytest --nbval-lax -v code/constrained_spherical_deconvolution/constrained_spherical_deconvolution.ipynb
pytest --nbval-lax -v code/deterministic_tractography/deterministic_tractography.ipynb
pytest --nbval-lax -v code/deterministic_tractography/solutions/deterministic_tractography_solutions.ipynb
pytest --nbval-lax -v code/probabilistic_tractography/probabilistic_tractography.ipynb
120 changes: 116 additions & 4 deletions _episodes/deterministic_tractography.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,7 @@ keypoints:
## Deterministic tractography

Deterministic tractography algorithms perform tracking of streamlines by
following a predictable path, such as following the primary diffusion direction
($\lambda_1$).
following a predictable path, such as following the primary diffusion direction.

In order to demonstrate how to perform deterministic tracking on a diffusion MRI dataset, we will
build from the preprocessing presented in a previous episode and compute the diffusion tensor.
Expand Down Expand Up @@ -160,8 +159,8 @@ thresholding above our stopping criterion.
from dipy.tracking import utils
seed_mask = fa_img.copy()
seed_mask[seed_mask>=0.2] = 1
seed_mask[seed_mask<0.2] = 0
seed_mask[seed_mask >= 0.2] = 1
seed_mask[seed_mask < 0.2] = 0
seeds = utils.seeds_from_mask(seed_mask, affine=affine, density=1)
~~~
Expand Down Expand Up @@ -228,4 +227,117 @@ plt.show()

![EuDX Determinsitic Tractography]({{ relative_root_path }}/fig/deterministic_tractography/tractogram_deterministic_EuDX.png){:class="img-responsive"}

> ## Exercise 1
>
> In this episode, we applied a threshold stopping criteria
> to stop tracking when we reach a voxel where FA is below 0.2.
> There are also other stopping criteria available. We encourage
> you to read the `DIPY` documentation about the others. For this
> exercise, repeat the tractography, but apply a binary stopping
> criteria (`BinaryStoppingCriterion`) using the seed mask.
> Visualize the tractogram!
>
> > ## Solution
> >
> > ~~~
> > import os
> >
> > import nibabel as nib
> > import numpy as np
> >
> > from bids.layout import BIDSLayout
> >
> > from dipy.io.gradients import read_bvals_bvecs
> > from dipy.core.gradients import gradient_table
> > from dipy.data import get_sphere
> > from dipy.direction import peaks_from_model
> > import dipy.reconst.dti as dti
> > from dipy.segment.mask import median_otsu
> > from dipy.tracking import utils
> > from dipy.tracking.local_tracking import LocalTracking
> > from dipy.tracking.streamline import Streamlines
> >
> > from utils.visualization_utils import generate_anatomical_volume_figure
> > from fury import actor, colormap
> >
> > dwi_layout = BIDSLayout("../../data/ds000221/derivatives/uncorrected_topup_eddy", validate=False)
> > gradient_layout = BIDSLayout("../../data/ds000221/", > > validate=False)
> >
> > # Get subject data
> > subj = '010006'
> > dwi_fname = dwi_layout.get(subject=subj, suffix='dwi', extension='nii.gz', return_type='file')[0]
> > bvec_fname = dwi_layout.get(subject=subj, extension='eddy_rotated_bvecs', return_type='file')[0]
> > bval_fname = gradient_layout.get(subject=subj, suffix='dwi', extension='bval', return_type='file')[0]
> >
> > dwi_img = nib.load(dwi_fname)
> > affine = dwi_img.affine
> >
> > bvals, bvecs = read_bvals_bvecs(bval_fname, bvec_fname)
> > gtab = gradient_table(bvals, bvecs)
> >
> > dwi_data = dwi_img.get_fdata()
> > dwi_data, dwi_mask = median_otsu(dwi_data, vol_idx=[0], numpass=1) # Specify the volume index to the b0 volumes
> >
> > # Fit tensor and compute FA map
> > dti_model = dti.TensorModel(gtab)
> > dti_fit = dti_model.fit(dwi_data, mask=dwi_mask)
> > fa_img = dti_fit.fa
> > evecs_img = dti_fit.evecs
> >
> > sphere = get_sphere('symmetric362')
> > peak_indices = peaks_from_model(model=dti_model, data=dwi_data, sphere=sphere, relative_peak_threshold=.2, min_separation_angle=25, mask=dwi_mask, npeaks=2)
> >
> > # Create a binary seed mask
> > seed_mask = fa_img.copy()
> > seed_mask[seed_mask >= 0.2] = 1
> > seed_mask[seed_mask < 0.2] = 0
> >
> > seeds = utils.seeds_from_mask(seed_mask, affine=affine, density=1)
> >
> > # Set stopping criteria
> > stopping_criterion = BinaryStoppingCriterion(seed_mask==1)
> >
> > # Perform tracking
> > streamlines_generator = LocalTracking(peak_indices, stopping_criterion, seeds, affine=affine, step_size=.5)
> > streamlines = Streamlines(streamlines_generator)
> >
> > # Plot the tractogram
> > # Build the representation of the data
streamlines_actor = actor.line(streamlines, colormap.line_colors(streamlines))
> >
> > # Generate the figure
> > fig = generate_anatomical_volume_figure(streamlines_actor)
> > plt.show()
> > ~~~
> > {: .language-python}
> > ![Binary Stopping Criterion Tractography]({{ relative_root_path }}/fig/deterministic_tractography/tractogram_deterministic_ex1.png){:class="img-responsive"}
> {: .solution}
>
> ## Exercise 2
>
> As an additional challenge, set the color of the streamlines to display the values of the
> FA map and change the opacity to `0.05`. You may need to transform
> the streamlines from world coordinates to the subject's native space
> using `transform_streamlines` from `dipy.tracking.streamline`.
>
> > ## Solution
> >
> > ~~~
> > import numpy as np
> > from fury import actor, window
> >
> > from dipy.tracking.streamline import transform_streamlines
> >
> > streamlines_native = transform_streamlines(streamlines, np.linalg.inv(affine))
> > streamlines_actor = actor.line(streamlines_native, fa_img, opacity=0.05)
> >
> > fig = generate_anatomical_volume_figure(streamlines_actor)
> > plt.show()
> > ~~~
> > {: .language-python}
> >
> > ![FA Mapped Tractography]({{ relative_root_path }}/fig/deterministic_tractography/tractogram_deterministic_fa.png){:class="img-responsive"}
> {: .solution}
{: .challenge}
{% include links.md %}
89 changes: 69 additions & 20 deletions code/deterministic_tractography/deterministic_tractography.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
"source": [
"## Deterministic tractography\n",
"\n",
"Deterministic tractography algorithms perform tracking of streamlines by following a predictable path, such as following the primary diffusion direction ($\\lambda_1$).\n",
"Deterministic tractography algorithms perform tracking of streamlines by following a predictable path, such as following the primary diffusion direction.\n",
"\n",
"In order to demonstrate how to perform deterministic tracking on a diffusion MRI dataset, we will build from the preprocessing presented in a previous episode and compute the diffusion tensor."
]
Expand All @@ -28,14 +28,18 @@
"from dipy.core.gradients import gradient_table\n",
"\n",
"\n",
"dwi_layout = BIDSLayout(\"../../data/ds000221/derivatives/uncorrected_topup_eddy\", validate=False)\n",
"dwi_layout = BIDSLayout(\n",
" \"../../data/ds000221/derivatives/uncorrected_topup_eddy\", validate=False)\n",
"gradient_layout = BIDSLayout(\"../../data/ds000221/\", validate=False)\n",
"\n",
"subj = '010006'\n",
"\n",
"dwi_fname = dwi_layout.get(subject=subj, suffix='dwi', extension='nii.gz', return_type='file')[0]\n",
"bvec_fname = dwi_layout.get(subject=subj, extension='eddy_rotated_bvecs', return_type='file')[0]\n",
"bval_fname = gradient_layout.get(subject=subj, suffix='dwi', extension='bval', return_type='file')[0]\n",
"dwi_fname = dwi_layout.get(subject=subj, suffix='dwi',\n",
" extension='nii.gz', return_type='file')[0]\n",
"bvec_fname = dwi_layout.get(\n",
" subject=subj, extension='eddy_rotated_bvecs', return_type='file')[0]\n",
"bval_fname = gradient_layout.get(\n",
" subject=subj, suffix='dwi', extension='bval', return_type='file')[0]\n",
"\n",
"dwi_img = nib.load(dwi_fname)\n",
"affine = dwi_img.affine\n",
Expand All @@ -61,7 +65,8 @@
"from dipy.segment.mask import median_otsu\n",
"\n",
"dwi_data = dwi_img.get_fdata()\n",
"dwi_data, dwi_mask = median_otsu(dwi_data, vol_idx=[0], numpass=1) # Specify the volume index to the b0 volumes\n",
"# Specify the volume index to the b0 volumes\n",
"dwi_data, dwi_mask = median_otsu(dwi_data, vol_idx=[0], numpass=1)\n",
"\n",
"dti_model = dti.TensorModel(gtab)\n",
"dti_fit = dti_model.fit(dwi_data, mask=dwi_mask) # This step may take a while"
Expand All @@ -84,6 +89,8 @@
"source": [
"# Create the directory to save the results\n",
"\n",
"from scipy import ndimage # To rotate image for visualization purposes\n",
"import matplotlib.pyplot as plt\n",
"out_dir = \"../../data/ds000221/derivatives/dwi/tractography/sub-%s/ses-01/dwi/\" % subj\n",
"\n",
"if not os.path.exists(out_dir):\n",
Expand All @@ -99,15 +106,16 @@
"nib.save(fa_nii, os.path.join(out_dir, 'fa.nii.gz'))\n",
"\n",
"# Plot the FA\n",
"import matplotlib.pyplot as plt\n",
"from scipy import ndimage # To rotate image for visualization purposes\n",
"\n",
"%matplotlib inline\n",
"\n",
"fig, ax = plt.subplots(1, 3, figsize=(10, 10))\n",
"ax[0].imshow(ndimage.rotate(fa_img[:, fa_img.shape[1]//2, :], 90, reshape=False))\n",
"ax[1].imshow(ndimage.rotate(fa_img[fa_img.shape[0]//2, :, :], 90, reshape=False))\n",
"ax[2].imshow(ndimage.rotate(fa_img[:, :, fa_img.shape[-1]//2], 90, reshape=False))\n",
"ax[0].imshow(ndimage.rotate(\n",
" fa_img[:, fa_img.shape[1]//2, :], 90, reshape=False))\n",
"ax[1].imshow(ndimage.rotate(\n",
" fa_img[fa_img.shape[0]//2, :, :], 90, reshape=False))\n",
"ax[2].imshow(ndimage.rotate(\n",
" fa_img[:, :, fa_img.shape[-1]//2], 90, reshape=False))\n",
"fig.savefig(os.path.join(out_dir, \"fa.png\"), dpi=300, bbox_inches=\"tight\")\n",
"plt.show()"
]
Expand Down Expand Up @@ -147,7 +155,8 @@
"source": [
"from dipy.direction import peaks_from_model\n",
"\n",
"peak_indices = peaks_from_model(model=dti_model, data=dwi_data, sphere=sphere, relative_peak_threshold=.2, min_separation_angle=25, mask=dwi_mask, npeaks=2)"
"peak_indices = peaks_from_model(model=dti_model, data=dwi_data, sphere=sphere,\n",
" relative_peak_threshold=.2, min_separation_angle=25, mask=dwi_mask, npeaks=2)"
]
},
{
Expand Down Expand Up @@ -184,8 +193,8 @@
"from dipy.tracking import utils\n",
"\n",
"seed_mask = fa_img.copy()\n",
"seed_mask[seed_mask>=0.2] = 1\n",
"seed_mask[seed_mask<0.2] = 0\n",
"seed_mask[seed_mask >= 0.2] = 1\n",
"seed_mask[seed_mask < 0.2] = 0\n",
"\n",
"seeds = utils.seeds_from_mask(seed_mask, affine=affine, density=1)"
]
Expand All @@ -209,7 +218,8 @@
"from dipy.tracking.streamline import Streamlines\n",
"\n",
"# Initialize local tracking - computation happens in the next step.\n",
"streamlines_generator = LocalTracking(peak_indices, stopping_criterion, seeds, affine=affine, step_size=.5)\n",
"streamlines_generator = LocalTracking(\n",
" peak_indices, stopping_criterion, seeds, affine=affine, step_size=.5)\n",
"\n",
"# Generate streamlines object\n",
"streamlines = Streamlines(streamlines_generator)"
Expand All @@ -234,7 +244,8 @@
"sft = StatefulTractogram(streamlines, dwi_img, Space.RASMM)\n",
"\n",
"# Save the tractogram\n",
"save_tractogram(sft, os.path.join(out_dir, \"tractogram_deterministic_EuDX.trk\"))"
"save_tractogram(sft, os.path.join(\n",
" out_dir, \"tractogram_deterministic_EuDX.trk\"))"
]
},
{
Expand All @@ -251,10 +262,8 @@
"outputs": [],
"source": [
"# NBVAL_SKIP\n",
"\n",
"from fury import actor, colormap\n",
"\n",
"from utils.visualization_utils import generate_anatomical_volume_figure\n",
"from fury import actor, colormap\n",
"\n",
"# Plot the tractogram\n",
"\n",
Expand All @@ -268,6 +277,46 @@
" dpi=300, bbox_inches=\"tight\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 1\n",
"\n",
"In this episode, we applied a threshold stopping criteria\n",
"to stop tracking when we reach a voxel where FA is below 0.2.\n",
"There are also other stopping criteria available. We encourage\n",
"you to read the `DIPY` documentation about the others. For this\n",
"exercise, repeat the tractography, but apply a binary stopping \n",
"criteria (`BinaryStoppingCriterion`) using the seed mask. \n",
"Visualize the tractogram!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 2\n",
"\n",
"As an additional challenge, set the color of the streamlines to display the values of the FA map and change the opacity to `0.05`. You may need to transform the streamlines from world coordinates to the subject's native space using `transform_streamlines` from `dipy.tracking.streamline`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": []
}
],
"metadata": {
Expand All @@ -286,7 +335,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.10"
"version": "3.6.9"
}
},
"nbformat": 4,
Expand Down
Loading

0 comments on commit 9179bcd

Please sign in to comment.