Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: Update deterministic tractography #120

Merged
merged 18 commits into from
Apr 12, 2021
Merged

ENH: Update deterministic tractography #120

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
e2de540
add exercise to det. tractography
kaitj Apr 3, 2021
149c213
add figure to solutions
kaitj Apr 3, 2021
5ae7242
fix solution formatting
kaitj Apr 3, 2021
9b2b6bb
rm commented line in solutions
kaitj Apr 3, 2021
89a076f
lint deterministic notebook
kaitj Apr 5, 2021
d16e8eb
lint markdown code & fix sentence structure
kaitj Apr 6, 2021
eb10991
add solutions notebook for exercise, update workflow
kaitj Apr 6, 2021
58b3203
Update exercise, add solutions notebook and include in workflow
kaitj Apr 6, 2021
e1f4bd6
Update non-solutions notebook
kaitj Apr 6, 2021
e8e0d24
Update figures
kaitj Apr 6, 2021
c3cecf2
Update solutions figure for det tract
kaitj Apr 6, 2021
160a09f
rm extra fig
kaitj Apr 6, 2021
594f1b4
rm whitespace
kaitj Apr 6, 2021
34c821a
rm savefig from solutions cell
kaitj Apr 6, 2021
eca626d
fix figure white space final time
kaitj Apr 6, 2021
75e998c
rm latex call
kaitj Apr 6, 2021
b886a5a
fix line breaks
kaitj Apr 8, 2021
18d790f
update with fig for ex1, fix spelling
kaitj Apr 9, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/build-test.yaml
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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