Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

MRtrix3Tissue

https://doi.org/10.5281/zenodo.18740151
Run this notebook

Authors: Thuy Dao, Monika Doerig

Date: 18 Nov 2024

License:

MIT License

Note: If this notebook uses neuroimaging tools from Neurocontainers, those tools retain their original licenses. Please see Neurodesk citation guidelines for details.

Citation and Resources:

Tools included in this workflow

MRtrix3Tissue:

  • MRtrix3Tissue, a fork of MRtrix3 (Tournier et al. 2019)

  • Tournier, J.-D., Smith, R., Raffelt, D., Tabbara, R., Dhollander, T., Pietsch, M., Christiaens, D., Jeurissen, B., Yeh, C.-H., & Connelly, A. (2019). MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 202, 116137. Tournier et al. (2019)

Dataset

Dataset from OpenNeuro: BTC_postop

Introduction

This notebook demonstrates advanced multi-tissue modelling using MRtrix3Tissue, an extension of the MRtrix3 framework designed for enhanced 3-tissue constrained spherical deconvolution (MSMT-CSD) and quantitative tissue analysis.

While the MRtrix series introduces preprocessing (Part 1), multi-tissue CSD and tissue boundary estimation (Part 2), and registration and streamline generation (Part 3), this notebook focuses on extended tissue modelling capabilities beyond tractography.

In particular, this tutorial demonstrates:

  • 3-tissue response function estimation using the dhollander algorithm

  • Multi-tissue CSD reconstruction

  • Intensity normalisation across tissue compartments (mtnormalise)

  • Tissue fraction–based masking

  • Derivation of quantitative tissue maps (e.g., ADC-based metrics)

Although the core modelling approach (3-tissue response estimation and MSMT-CSD) is conceptually similar to Part 2 of the MRtrix series, the goal here is not streamline generation but rather advanced tissue characterisation and quantitative analysis.

Load software tools and import python libraries

%%capture
!pip install numpy nibabel
import subprocess
import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
import tempfile
import os
import module
await module.load('mrtrix3tissue/5.2.8')
await module.list()
['mrtrix3tissue/5.2.8']

Data preparation

# download data
!datalad install https://github.com/OpenNeuroDatasets/ds002080.git
!cd ds002080 && datalad get sub-CON02
Cloning:   0%|                             | 0.00/2.00 [00:00<?, ? candidates/s]
Enumerating: 0.00 Objects [00:00, ? Objects/s]
                                              
Counting:   0%|                              | 0.00/39.9k [00:00<?, ? Objects/s]
                                                                                
Compressing:   0%|                           | 0.00/15.2k [00:00<?, ? Objects/s]
                                                                                
Receiving:   0%|                             | 0.00/41.9k [00:00<?, ? Objects/s]
Receiving:  19%|███▊                | 7.95k/41.9k [00:00<00:00, 78.0k Objects/s]
Receiving:  59%|████████████▍        | 24.7k/41.9k [00:00<00:00, 128k Objects/s]
                                                                                
Resolving:   0%|                              | 0.00/19.5k [00:00<?, ? Deltas/s]
Resolving:  40%|████████▍            | 7.80k/19.5k [00:00<00:00, 77.2k Deltas/s]
Resolving:  92%|███████████████████▎ | 17.9k/19.5k [00:00<00:00, 91.0k Deltas/s]
[INFO   ] Remote origin not usable by git-annex; setting annex-ignore 
[INFO   ] https://github.com/OpenNeuroDatasets/ds002080.git/config download failed: Not Found 
[INFO   ] access to 1 dataset sibling s3-BACKUP not auto-enabled, enable with:
| 		datalad siblings -d "/home/jovyan/workspace/books/examples/diffusion_imaging/ds002080" enable -s s3-BACKUP 
install(ok): /home/jovyan/workspace/books/examples/diffusion_imaging/ds002080 (dataset)
Total:   0%|                                   | 0.00/97.7M [00:00<?, ? Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   0%|            | 0.00/8.06M [00:00<?, ? Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   0%|    | 33.3k/8.06M [00:00<01:05, 122k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   1%|    | 68.2k/8.06M [00:00<00:58, 137k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   2%|     | 138k/8.06M [00:00<00:36, 217k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   2%|     | 190k/8.06M [00:00<00:34, 229k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   3%|▏    | 260k/8.06M [00:01<00:29, 264k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   4%|▏    | 312k/8.06M [00:01<00:24, 312k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   5%|▏    | 399k/8.06M [00:01<00:26, 293k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   6%|▎    | 469k/8.06M [00:01<00:24, 306k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   7%|▎    | 538k/8.06M [00:01<00:23, 314k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   8%|▍    | 625k/8.06M [00:02<00:21, 344k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:   9%|▍    | 712k/8.06M [00:02<00:20, 365k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  10%|▍    | 799k/8.06M [00:02<00:19, 380k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  11%|▌    | 886k/8.06M [00:02<00:18, 390k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  12%|▌    | 991k/8.06M [00:03<00:16, 422k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  14%|▌   | 1.10M/8.06M [00:03<00:15, 444k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  15%|▌   | 1.20M/8.06M [00:03<00:14, 460k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  16%|▋   | 1.30M/8.06M [00:03<00:14, 471k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  18%|▋   | 1.43M/8.06M [00:03<00:13, 503k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  19%|▊   | 1.55M/8.06M [00:04<00:12, 527k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  21%|▊   | 1.69M/8.06M [00:04<00:11, 566k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  23%|▉   | 1.83M/8.06M [00:04<00:10, 594k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  24%|▉   | 1.97M/8.06M [00:04<00:09, 615k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  26%|█   | 2.12M/8.06M [00:04<00:09, 653k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  28%|█▏  | 2.28M/8.06M [00:05<00:08, 681k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  30%|█▏  | 2.45M/8.06M [00:05<00:07, 724k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  33%|█▎  | 2.63M/8.06M [00:05<00:07, 755k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  35%|█▍  | 2.80M/8.06M [00:05<00:06, 780k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  37%|█▍  | 2.99M/8.06M [00:05<00:06, 819k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  40%|█▌  | 3.20M/8.06M [00:06<00:05, 868k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  42%|█▋  | 3.41M/8.06M [00:06<00:05, 911k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  45%|█▊  | 3.64M/8.06M [00:06<00:04, 958k Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  48%|█▍ | 3.88M/8.06M [00:06<00:04, 1.01M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  51%|█▌ | 4.12M/8.06M [00:07<00:03, 1.06M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  54%|█▋ | 4.39M/8.06M [00:07<00:03, 1.11M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  58%|█▋ | 4.66M/8.06M [00:07<00:02, 1.18M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  61%|█▊ | 4.94M/8.06M [00:07<00:02, 1.22M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  65%|█▉ | 5.24M/8.06M [00:07<00:02, 1.28M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  69%|██ | 5.57M/8.06M [00:08<00:01, 1.36M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  73%|██▏| 5.90M/8.06M [00:08<00:01, 1.43M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  78%|██▎| 6.25M/8.06M [00:08<00:01, 1.50M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  82%|██▍| 6.61M/8.06M [00:08<00:00, 1.57M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  87%|██▌| 7.00M/8.06M [00:08<00:00, 1.65M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  92%|██▊| 7.41M/8.06M [00:09<00:00, 1.75M Bytes/s]
Get sub-CON0 .. p_T1w.nii.gz:  97%|██▉| 7.85M/8.06M [00:09<00:00, 2.11M Bytes/s]
                                                                                
Get sub-CON0 .. p_T1w.nii.gz:   0%|            | 0.00/8.06M [00:00<?, ? Bytes/s]
Total:   8%|██▏                        | 8.06M/97.7M [00:11<02:06, 710k Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   0%|            | 0.00/48.5M [00:00<?, ? Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   1%|    | 483k/48.5M [00:00<00:16, 2.96M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   2%|   | 1.01M/48.5M [00:00<00:20, 2.34M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   3%|   | 1.55M/48.5M [00:00<00:19, 2.44M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   4%|▏  | 2.13M/48.5M [00:00<00:18, 2.56M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   6%|▏  | 2.72M/48.5M [00:01<00:17, 2.65M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   7%|▏  | 3.34M/48.5M [00:01<00:13, 3.26M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:   8%|▎  | 4.05M/48.5M [00:01<00:15, 2.81M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  10%|▎  | 4.75M/48.5M [00:01<00:14, 3.06M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  11%|▎  | 5.52M/48.5M [00:01<00:13, 3.16M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  13%|▍  | 6.32M/48.5M [00:02<00:12, 3.36M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  15%|▍  | 7.18M/48.5M [00:02<00:11, 3.56M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  17%|▍  | 8.06M/48.5M [00:02<00:10, 3.76M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  19%|▌  | 9.00M/48.5M [00:02<00:08, 4.39M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  20%|▌  | 9.58M/48.5M [00:02<00:08, 4.67M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  23%|▋  | 11.0M/48.5M [00:03<00:08, 4.20M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  25%|▊  | 12.1M/48.5M [00:03<00:08, 4.49M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  27%|▊  | 13.3M/48.5M [00:03<00:07, 4.76M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  30%|▉  | 14.5M/48.5M [00:03<00:06, 5.04M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  33%|▉  | 15.8M/48.5M [00:03<00:06, 5.33M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  35%|█  | 17.1M/48.5M [00:04<00:05, 5.92M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  37%|█  | 18.0M/48.5M [00:04<00:04, 6.46M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  41%|█▏ | 20.0M/48.5M [00:04<00:04, 6.03M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  44%|█▎ | 21.4M/48.5M [00:04<00:03, 7.29M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  47%|█▍ | 23.0M/48.5M [00:04<00:03, 7.45M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  51%|█▌ | 24.8M/48.5M [00:05<00:03, 7.55M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  55%|█▋ | 26.6M/48.5M [00:05<00:03, 6.98M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  58%|█▋ | 28.0M/48.5M [00:05<00:02, 8.00M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  63%|█▉ | 30.5M/48.5M [00:05<00:02, 7.89M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  67%|██ | 32.5M/48.5M [00:06<00:01, 8.43M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  70%|██ | 33.9M/48.5M [00:06<00:01, 9.32M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  74%|██▏| 36.1M/48.5M [00:06<00:01, 9.68M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  79%|██▎| 38.2M/48.5M [00:06<00:01, 9.85M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  83%|██▌| 40.4M/48.5M [00:06<00:00, 10.0M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  88%|██▋| 42.5M/48.5M [00:07<00:00, 10.0M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  92%|██▊| 44.7M/48.5M [00:07<00:00, 10.1M Bytes/s]
Get sub-CON0 .. P_dwi.nii.gz:  97%|██▉| 47.0M/48.5M [00:07<00:00, 10.3M Bytes/s]
                                                                                
Get sub-CON0 .. P_dwi.nii.gz:   0%|            | 0.00/48.5M [00:00<?, ? Bytes/s]
Total:  58%|███████████████           | 56.5M/97.7M [00:19<00:14, 2.90M Bytes/s]
Get sub-CON0 .. A_dwi.nii.gz:   0%|            | 0.00/1.18M [00:00<?, ? Bytes/s]
                                                                                
Get sub-CON0 .. A_dwi.nii.gz:   0%|            | 0.00/1.18M [00:00<?, ? Bytes/s]
                                                                                
Get sub-CON0 .. _bold.nii.gz:   0%|            | 0.00/40.0M [00:00<?, ? Bytes/s]
Get sub-CON0 .. _bold.nii.gz:   6%|▏  | 2.36M/40.0M [00:00<00:01, 23.6M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  12%|▎  | 4.77M/40.0M [00:00<00:03, 10.4M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  18%|▌  | 7.20M/40.0M [00:00<00:03, 10.9M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  23%|▋  | 9.31M/40.0M [00:00<00:02, 10.6M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  29%|▊  | 11.6M/40.0M [00:01<00:02, 10.7M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  34%|█  | 13.7M/40.0M [00:01<00:02, 12.4M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  38%|█▏ | 15.3M/40.0M [00:01<00:02, 10.7M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  44%|█▎ | 17.5M/40.0M [00:01<00:02, 10.8M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  49%|█▍ | 19.6M/40.0M [00:01<00:01, 10.5M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  56%|█▋ | 22.2M/40.0M [00:02<00:01, 9.50M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  61%|█▊ | 24.6M/40.0M [00:02<00:01, 9.95M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  66%|█▉ | 26.6M/40.0M [00:02<00:01, 9.85M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  72%|██▏| 28.8M/40.0M [00:02<00:01, 11.1M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  77%|██▎| 30.9M/40.0M [00:02<00:00, 9.74M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  83%|██▍| 33.0M/40.0M [00:03<00:00, 9.85M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  87%|██▌| 34.6M/40.0M [00:03<00:00, 10.8M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  91%|██▋| 36.6M/40.0M [00:03<00:00, 10.3M Bytes/s]
Get sub-CON0 .. _bold.nii.gz:  98%|██▉| 39.3M/40.0M [00:03<00:00, 11.1M Bytes/s]
                                                                                
Get sub-CON0 .. _bold.nii.gz:   0%|            | 0.00/40.0M [00:00<?, ? Bytes/s]
get(ok): sub-CON02/ses-postop/anat/sub-CON02_ses-postop_T1w.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-AP_dwi.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-PA_dwi.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-CON02/ses-postop/func/sub-CON02_ses-postop_task-rest_bold.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-CON02 (directory)
action summary:
  get (ok: 5)

Convert to .mif file

Once data is downloaded, it is recommended to import and store it as .mif file(s), the so-called “MRtrix Image Format”.

mrconvert is the tool to convert between all sorts of image formats (or extract parts of images or change properties) and enables you to read from e.g. DICOM folders or NIfTI (.nii) images and store an image as a .mif file.

If your diffusion MRI data comes as a NIfTI (.nii) image instead, the gradient orientations and b-values will typically be stored in 2 separate files (a “bvecs” and “bvals” file).

%%bash
mrconvert ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-AP_dwi.nii.gz AP.mif -fslgrad ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-AP_dwi.bvec ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-AP_dwi.bval
mrconvert ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-PA_dwi.nii.gz PA.mif
mrconvert PA.mif -fslgrad ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-PA_dwi.bvec ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-PA_dwi.bval - | mrmath - mean mean_b0_PA.mif -axis 3
mrconvert: [100%] uncompressing image "ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-AP_dwi.nii.gz"
mrconvert: [100%] copying from "ds002080/s...es-postop_acq-AP_dwi.nii.gz" to "AP.mif"
mrconvert: [100%] uncompressing image "ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-PA_dwi.nii.gz"
mrconvert: [100%] copying from "ds002080/s...es-postop_acq-PA_dwi.nii.gz" to "PA.mif"
mrconvert: [100%] copying from "PA.mif" to "/tmp/mrtrix-tmp-sNlnBg.mif"
mrmath: [100%] preloading data for "/tmp/mrtrix-tmp-sNlnBg.mif"
mrmath: [100%] computing mean along axis 3...

Preprocessing

The following are common preprocessing steps done with MRtrix.

Denoising and Gibbs ringing removal (“unringing”)

The first preprocessing step is denoise the data by using dwidenoise command. This requires an input and an output argument, and you also have the option to output the noise map with the -noise option.

The next step is to run mrdegibbs, which removes Gibbs’ ringing artifacts from the data.

%%bash
dwidenoise AP.mif dwi_denoised.mif -noise noise.mif
mrdegibbs dwi_denoised.mif dwi_denoised_unringed.mif
dwidenoise: [100%] preloading data for "AP.mif"
dwidenoise: [100%] running MP-PCA denoising
mrdegibbs: [100%] performing Gibbs ringing removal

Next, we extract the b-values from the primary phase-encoded image using dwiextract.

If you have also scanned a pair of reverse phase-encoded b=0 images, you can correct for susceptibility-induced EPI distortions by concatenating the reverse phase-encoded b=0 images using mrcat to create a new image b0_pair.mif

!dwiextract dwi_denoised.mif - -bzero | mrmath - mean mean_b0_AP.mif -axis 3
!mrcat mean_b0_AP.mif mean_b0_PA.mif -axis 3 b0_pair.mif
dwiextract: [100%] extracting volumes
mrmath: [100%] preloading data for "/tmp/mrtrix-tmp-RKriQO.mif"
mrmath: [100%] computing mean along axis 3...
mrcat: [100%] concatenating "mean_b0_AP.mif"
mrcat: [100%] concatenating "mean_b0_PA.mif"

Motion and distortion correction

Motion and distortion correction are performed by the dwifslpreproc command, most of the heavy lifting is done automatically by FSL’s topup and eddy tools.

-pe_dir option is used to tell dwifslpreproc the phase-encoding direction of the scan.

-se_epi is an additional image series consisting of spin-echo EPI images used by topup for estimating the inhomogeneity field

-eddy_options is an additional command-line options to the eddy

!dwifslpreproc dwi_denoised_unringed.mif dwi_denoised_unringed_preproc.mif -pe_dir AP -rpe_pair -se_epi b0_pair.mif  -eddy_options " --slm=linear --data_is_shelled"
dwifslpreproc: 
dwifslpreproc: Note that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.
dwifslpreproc: 
dwifslpreproc: Generated scratch directory: /home/jovyan/workspace/books/examples/diffusion_imaging/dwifslpreproc-tmp-U83OQG/
Command:  mrconvert /home/jovyan/workspace/books/examples/diffusion_imaging/dwi_denoised_unringed.mif /home/jovyan/workspace/books/examples/diffusion_imaging/dwifslpreproc-tmp-U83OQG/dwi.mif -json_export /home/jovyan/workspace/books/examples/diffusion_imaging/dwifslpreproc-tmp-U83OQG/dwi.json
Command:  mrconvert /home/jovyan/workspace/books/examples/diffusion_imaging/b0_pair.mif /home/jovyan/workspace/books/examples/diffusion_imaging/dwifslpreproc-tmp-U83OQG/se_epi.mif
dwifslpreproc: Changing to scratch directory (/home/jovyan/workspace/books/examples/diffusion_imaging/dwifslpreproc-tmp-U83OQG/)
dwifslpreproc: Total readout time not provided at command-line; assuming sane default of 0.1
Command:  mrinfo dwi.mif -export_grad_mrtrix grad.b
Command:  mrconvert se_epi.mif topup_in.nii -import_pe_table se_epi_manual_pe_scheme.txt -strides -1,+2,+3,+4 -export_pe_table topup_datain.txt
Command:  topup --imain=topup_in.nii --datain=topup_datain.txt --out=field --fout=field_map.nii.gz --config=/opt/fsl-6.0.3/etc/flirtsch/b02b0.cnf --verbose
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt - | mrinfo - -export_pe_eddy applytopup_config.txt applytopup_indices.txt
Command:  dwiextract dwi.mif -import_pe_table dwi_manual_pe_scheme.txt -pe 0.0,-1.0,0.0,0.1 - | mrconvert - dwi_pe_1.nii -json_export dwi_pe_1.json
Command:  applytopup --imain=dwi_pe_1.nii --datain=applytopup_config.txt --inindex=1 --topup=field --out=dwi_pe_1_applytopup.nii --method=jac
Command:  mrconvert dwi_pe_1_applytopup.nii.gz dwi_pe_1_applytopup.mif -json_import dwi_pe_1.json
Command:  dwi2mask dwi_pe_1_applytopup.mif - | maskfilter - dilate - | mrconvert - eddy_mask.nii -datatype float32 -strides -1,+2,+3
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt eddy_in.nii -strides -1,+2,+3,+4 -export_grad_fsl bvecs bvals -export_pe_eddy eddy_config.txt eddy_indices.txt
Command:  eddy_cuda9.1 --imain=eddy_in.nii --mask=eddy_mask.nii --acqp=eddy_config.txt --index=eddy_indices.txt --bvecs=bvecs --bvals=bvals --topup=field --slm=linear --data_is_shelled --out=dwi_post_eddy --verbose
dwifslpreproc: CUDA version of 'eddy' was not successful; attempting OpenMP version
Command:  eddy_openmp --imain=eddy_in.nii --mask=eddy_mask.nii --acqp=eddy_config.txt --index=eddy_indices.txt --bvecs=bvecs --bvals=bvals --topup=field --slm=linear --data_is_shelled --out=dwi_post_eddy --verbose
Command:  mrconvert dwi_post_eddy.nii.gz result.mif -strides -1,2,3,4 -fslgrad dwi_post_eddy.eddy_rotated_bvecs bvals
Command:  mrconvert result.mif /home/jovyan/workspace/books/examples/diffusion_imaging/dwi_denoised_unringed_preproc.mif
dwifslpreproc: Changing back to original directory (/home/jovyan/workspace/books/examples/diffusion_imaging)
dwifslpreproc: Deleting scratch directory (/home/jovyan/workspace/books/examples/diffusion_imaging/dwifslpreproc-tmp-U83OQG/)

Bias field correction

You can correct for bias fields at this point in the pipeline with the dwibiascorrect ants command

!dwibiascorrect ants dwi_denoised_unringed_preproc.mif dwi_denoised_unringed_preproc_unbiased.mif
dwibiascorrect: 
dwibiascorrect: Note that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.
dwibiascorrect: 
dwibiascorrect: Generated scratch directory: /home/jovyan/workspace/books/examples/diffusion_imaging/dwibiascorrect-tmp-AFN8NH/
Command:  mrconvert /home/jovyan/workspace/books/examples/diffusion_imaging/dwi_denoised_unringed_preproc.mif /home/jovyan/workspace/books/examples/diffusion_imaging/dwibiascorrect-tmp-AFN8NH/in.mif
dwibiascorrect: Changing to scratch directory (/home/jovyan/workspace/books/examples/diffusion_imaging/dwibiascorrect-tmp-AFN8NH/)
Command:  dwi2mask in.mif mask.mif
Command:  dwiextract in.mif - -bzero | mrmath - mean mean_bzero.mif -axis 3
Command:  mrconvert mean_bzero.mif mean_bzero.nii -strides +1,+2,+3
Command:  mrconvert mask.mif mask.nii -strides +1,+2,+3
Command:  N4BiasFieldCorrection -d 3 -i mean_bzero.nii -w mask.nii -o [corrected.nii,init_bias.nii] -s 4 -b [100,3] -c [1000,0.0]
Command:  mrcalc mean_bzero.mif mask.mif -mult - | mrmath - sum - -axis 0 | mrmath - sum - -axis 1 | mrmath - sum - -axis 2 | mrdump -
Command:  mrcalc corrected.nii mask.mif -mult - | mrmath - sum - -axis 0 | mrmath - sum - -axis 1 | mrmath - sum - -axis 2 | mrdump -
Command:  mrcalc init_bias.nii 0.569643924646 -mult bias.mif
Command:  mrcalc in.mif bias.mif -div result.mif
Command:  mrconvert result.mif /home/jovyan/workspace/books/examples/diffusion_imaging/dwi_denoised_unringed_preproc_unbiased.mif
dwibiascorrect: Changing back to original directory (/home/jovyan/workspace/books/examples/diffusion_imaging)
dwibiascorrect: Deleting scratch directory (/home/jovyan/workspace/books/examples/diffusion_imaging/dwibiascorrect-tmp-AFN8NH/)

Brain mask estimation

A brain mask can be estimated automatically from the diffusion MRI data using dwi2mask as follows:

!dwi2mask dwi_denoised_unringed_preproc_unbiased.mif dwi_mask.mif
dwi2mask: [100%] preloading data for "dwi_denoised_unringed_preproc_unbiased.mif"
dwi2mask: [100%] finding min/max of "mean b=0 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=700 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=1200 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=2800 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [done] computing dwi brain mask
dwi2mask: [done] applying mask cleaning filter

Decomposition

3-tissue response function estimation

An method to obtain 3-tissue response functions representing single-fibre white matter (WM), grey matter (GM) and CSF from the data itself, is available as the so-called dwi2response dhollander command from MRtrix3Tissue package.

!dwi2response dhollander dwi_denoised_unringed_preproc_unbiased.mif response_wm.txt response_gm.txt response_csf.txt
dwi2response: 
dwi2response: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
dwi2response: 
dwi2response: Generated scratch directory: /home/jovyan/workspace/books/examples/diffusion_imaging/dwi2response-tmp-DXIKY2/
dwi2response: Importing DWI data (/home/jovyan/workspace/books/examples/diffusion_imaging/dwi_denoised_unringed_preproc_unbiased.mif)...
dwi2response: Changing to scratch directory (/home/jovyan/workspace/books/examples/diffusion_imaging/dwi2response-tmp-DXIKY2/)
dwi2response: Computing brain mask (dwi2mask)...
dwi2response: -------
dwi2response: 4 unique b-value(s) detected: 0,700,1200,2800 with 6,16,30,50 volumes
dwi2response: -------
dwi2response: Preparation:
dwi2response: * Eroding brain mask by 3 pass(es)...
dwi2response:   [ mask: 96494 -> 69370 ]
dwi2response: * Computing signal decay metric (SDM):
dwi2response:  * b=0...
dwi2response:  * b=700...
dwi2response:  * b=1200...
dwi2response:  * b=2800...
dwi2response: * Removing erroneous voxels from mask and correcting SDM...
dwi2response:   [ mask: 69370 -> 69370 ]
dwi2response: -------
dwi2response: Crude segmentation:
dwi2response: * Crude WM versus GM-CSF separation (at FA=0.2)...
dwi2response:   [ 69370 -> 37155 (WM) & 32215 (GM-CSF) ]
dwi2response: * Crude GM versus CSF separation...
dwi2response:   [ 32215 -> 21449 (GM) & 10766 (CSF) ]
dwi2response: -------
dwi2response: Refined segmentation:
dwi2response: * Refining WM...
dwi2response:   [ WM: 37155 -> 33624 ]
dwi2response: * Refining GM...
dwi2response:   [ GM: 21449 -> 11941 ]
dwi2response: * Refining CSF...
dwi2response:   [ CSF: 10766 -> 5365 ]
dwi2response: -------
dwi2response: Final voxel selection and response function estimation:
dwi2response: * CSF:
dwi2response:  * Selecting final voxels (10.0% of refined CSF)...
dwi2response:    [ CSF: 5365 -> 537 ]
dwi2response:  * Estimating response function...
dwi2response: * GM:
dwi2response:  * Selecting final voxels (2.0% of refined GM)...
dwi2response:    [ GM: 11941 -> 239 ]
dwi2response:  * Estimating response function...
dwi2response: * single-fibre WM:
dwi2response:  * Selecting final voxels (0.5% of refined WM)...
dwi2response:    [ WM: 33624 -> 168 (single-fibre) ]
dwi2response:  * Estimating response function...
dwi2response: -------
dwi2response: Generating outputs...
dwi2response: -------
dwi2response: Changing back to original directory (/home/jovyan/workspace/books/examples/diffusion_imaging)
dwi2response: Deleting scratch directory (/home/jovyan/workspace/books/examples/diffusion_imaging/dwi2response-tmp-DXIKY2/)

Upsampling

Upsampling your preprocessed diffusion MRI data before the 3-tissue constrained spherical deconvolution (3-tissue CSD) step to get a higher quality result. Upsampling is useful if you’re after a higher quality visualisation, which reveals finer details using mrgrid.

!mrgrid dwi_denoised_unringed_preproc_unbiased.mif regrid dwi_denoised_unringed_preproc_unbiased_upsampled.mif -voxel 1.5
mrgrid: [100%] reslicing "dwi_denoised_unringed_preproc_unbiased.mif"

After upsampling the image, you need to regrid the brain mask to the exact same voxel grid as the image via the -template option to mrgrid. If the regridding is effectively upsampling, you can avoid jagged edges around masks by performing linear interpolation -interp linear. Then a median filter the result using the maskfilter command. To perform all of this at once (without storing an intermediate image), you an directly “pipe” the output of mrgrid into the input of maskfilter, by replacing the output in the one command and the input in the other command by a dash (“-”), and typing a “pipe” (“|”) between both commands.

!mrgrid dwi_mask.mif regrid - -template dwi_denoised_unringed_preproc_unbiased_upsampled.mif -interp linear -datatype bit | maskfilter - median dwi_mask_upsampled.mif
mrgrid: [100%] reslicing "dwi_mask.mif"
maskfilter: [100%] applying median filter to image -

3-tissue CSD modelling

Your data is now ready for 3-tissue CSD modelling with the previously obtained 3-tissue response functions, which will result in modelling the diffusion MRI data using WM-like (FOD), GM-like and CSF-like compartments. There are 2 different methods (or algorithms) available to perform 3-tissue CSD. The choice between both depends on what (part of your) data you intend to perform 3-tissue CSD modelling for.

If you want to perform 3-tissue CSD modelling for multi-shell data, this can be achieved using the multi-shell multi-tissue CSD (MSMT-CSD) method or algorithm, as follows:

!dwi2fod msmt_csd dwi_denoised_unringed_preproc_unbiased_upsampled.mif response_wm.txt wmfod.mif response_gm.txt gm.mif response_csf.txt csf.mif -mask dwi_mask_upsampled.mif
dwi2fod: [100%] preloading data for "dwi_denoised_unringed_preproc_unbiased_upsampled.mif"
dwi2fod: [100%] performing multi-shell, multi-tissue CSD

3-tissue bias field and intensity normalisation

Bias field correction and global intensity normalisation can be performed using the mtnormalise command, as follows:

!mtnormalise wmfod.mif wmfod_norm.mif gm.mif gm_norm.mif csf.mif csf_norm.mif -mask dwi_mask_upsampled.mif
mtnormalise: [done] loading input images
mtnormalise: [100%] performing log-domain intensity normalisation
mtnormalise: [done] writing output images
mtnormalise: [100%] performing log-domain intensity normalisation
mtnormalise: [done] loading input images

Temperature computation

Get CSF mask

You first need to get 2 CSF masks from 2 criteria:

  • csf_mask1.mif from CSF spherical harmonic coefficients applying the first criteria (aka the csf.mif file)

  • csf_mask2.mif from second criteria (note: 0.141047 might be half of max response function)

Then find the overlap of both masks resulting in the final mask csf_maskcombined.mif.

%%bash
mrconvert wmfod_norm.mif -coord 3 0 -axes 0,1,2 - | mrcalc csf_norm.mif gm_norm.mif - -add 5 -mult -gt csf_mask1.mif -datatype bit
mrthreshold csf_norm.mif csf_mask2.mif -abs 0.141047
mrcalc csf_mask1.mif csf_mask2.mif -mult csf_maskcombined.mif -datatype bit
mrconvert: [100%] copying from "wmfod_norm.mif" to "/tmp/mrtrix-tmp-xz7Ty7.mif"
mrcalc: [100%] computing: (csf_norm.mif > ((gm_norm.mif + /tmp/mrtrix-tmp-xz7Ty7.mif) * 5))
mrthreshold: [100%] Determining and applying per-volume thresholds
mrcalc: [100%] computing: (csf_mask1.mif * csf_mask2.mif)

Calculate diffusion coefficient in CSF

First, we filter the processed image to exclude b-values over 1000 into a new image reduced.mif.

Next, estimate the Apparent Diffusion Coefficient (ADC) using dwi2adc.

Then we extract the coefficient from the second array in dwi_adc.mif to get D values.

Finally, we filter the D values within CSF using generated mask.

# Get the index of b-values < 1000 from table and write to txt file for next step
bvals = np.loadtxt("ds002080/sub-CON02/ses-postop/dwi/sub-CON02_ses-postop_acq-AP_dwi.bval")
bvals_filtered = np.where(bvals < 1000)

with open("bval_less_1000.txt", "w") as txt_file:
    txt_file.write(",".join((str(x) for x in bvals_filtered[0])))
%%bash
mrconvert dwi_denoised_unringed_preproc_unbiased_upsampled.mif reduced.mif -coord 3 $(cat bval_less_1000.txt)
dwi2adc reduced.mif dwi_adc.mif
mrconvert dwi_adc.mif extracted_adc.mif -coord 3 1
mrcalc extracted_adc.mif csf_maskcombined.mif -mult masked_adc.mif
mrconvert: [100%] copying from "dwi_denois...proc_unbiased_upsampled.mif" to "reduced.mif"
dwi2adc: [100%] computing ADC values
mrconvert: [100%] copying from "dwi_adc.mif" to "extracted_adc.mif"
mrcalc: [100%] computing: (extracted_adc.mif * csf_maskcombined.mif)

Temperature estimation

We calculate the temperature of cerebrospinal fluid (CSF) using the following equation from this study:

T=2256.74ln⁡(4.39221D)−273T = \frac{2256.74}{\ln\left(\frac{4.39221}{D}\right)} - 273
(1)

where ( D ) is the diffusion constant.

!mrcalc 2256.74 4.39221 masked_adc.mif -divide -log -divide 273.15 -subtract masked_temperature.mif
mrcalc: [100%] computing: ((2256.74 / log ((4.39221 / masked_adc.mif))) - 273.15)

Visualization

# File paths for both images
mif_file_path1 = 'dwi_denoised_unringed_preproc_unbiased_upsampled.mif'  
mif_file_path2 = 'masked_temperature.mif'          

# Convert both .mif files to temporary .nii.gz files with -force flag
with tempfile.NamedTemporaryFile(suffix=".nii.gz") as temp_file1, \
     tempfile.NamedTemporaryFile(suffix=".nii.gz") as temp_file2:
         
    # Convert main image
    result1 = subprocess.run(["mrconvert", mif_file_path1, temp_file1.name, "-force"], capture_output=True, text=True)
    if result1.returncode != 0:
        print(f"Error in mrconvert for preprocessed image: {result1.stderr}")
        exit()
    
    # Convert overlay image
    result2 = subprocess.run(["mrconvert", mif_file_path2, temp_file2.name, "-force"], capture_output=True, text=True)
    if result2.returncode != 0:
        print(f"Error in mrconvert for original image: {result2.stderr}")
        exit()
         
    # Load the converted images
    nii_image1 = nib.load(temp_file1.name)
    nii_image2 = nib.load(temp_file2.name)

    data1 = nii_image1.get_fdata()
    data2 = nii_image2.get_fdata()

    # Select middle slice 
    slice_index = data1.shape[2] // 3

    fig, axes = plt.subplots(1, 1, figsize=(6, 6), constrained_layout=True)

    # Plot original diffusion data as overlay
    axes.imshow(np.rot90(data1[:, :, slice_index, 0]), cmap="Greys_r", vmin=0, vmax=1000) 
    overlay = axes.imshow(np.rot90(data2[:, :, slice_index, 0]), cmap="hot", alpha=0.5, vmin=-10, vmax=90) 
    axes.set_title("Preprocessed data with original diffusion data as overlay")
    axes.axis("off")
    fig.colorbar(overlay, ax=axes, location='right')
    plt.show()
<Figure size 600x600 with 2 Axes>

Dependencies in Jupyter/Python

  • Using the package watermark to document system environment and software versions used in this notebook, alongside the Neurodesktop version extracted from the JUPYTER_IMAGE or NEURODESKTOP_VERSION environment variables.

import os

%load_ext watermark

%watermark
%watermark --iversions

neurodesktop_version = (
    os.environ.get('JUPYTER_IMAGE', '').split(':')[-1] or
    os.environ.get('NEURODESKTOP_VERSION', 'unknown')
)

print(f"Neurodesktop version: {neurodesktop_version}")
Last updated: 2026-04-09T03:34:02.546828+00:00

Python implementation: CPython
Python version       : 3.13.9
IPython version      : 9.7.0

Compiler    : GCC 14.3.0
OS          : Linux
Release     : 5.15.0-171-generic
Machine     : x86_64
Processor   : x86_64
CPU cores   : 32
Architecture: 64bit

matplotlib: 3.10.8
nibabel   : 5.3.3
numpy     : 2.3.5

Neurodesktop version: 2025-12-20
References
  1. Tournier, J.-D., Smith, R., Raffelt, D., Tabbara, R., Dhollander, T., Pietsch, M., Christiaens, D., Jeurissen, B., Yeh, C.-H., & Connelly, A. (2019). MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 202, 116137. 10.1016/j.neuroimage.2019.116137
Neurodesk
MRtrix: Part 3
Neurodesk
QSIPrep