{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Reconstruct with Q-space Trajectory Imaging (QTI)\n\nQ-space trajectory imaging (QTI) by Westin et al. [1]_ is a general framework\nfor analyzing diffusion-weighted MRI data acquired with tensor-valued diffusion\nencoding. This tutorial briefly summarizes the theory and provides an example\nof how to estimate the diffusion and covariance tensors using DIPY.\n\n## Theory\n\nIn QTI, the tissue microstructure is represented by a diffusion tensor\ndistribution (DTD). Here, DTD is denoted by $\\mathbf{D}$ and the voxel-level\ndiffusion tensor from DTI by $\\langle\\mathbf{D}\rangle$, where\n$\\langle \\ \rangle$ denotes averaging over the DTD. The covariance of\n$\\mathbf{D}$ is given by a fourth-order covariance tensor $\\mathbb{C}$ defined\nas\n\n\\begin{align}\\mathbb{C} = \\langle \\mathbf{D} \\otimes \\mathbf{D} \rangle - \\langle\n \\mathbf{D} \rangle \\otimes \\langle \\mathbf{D} \rangle ,\\end{align}\n\nwhere $\\otimes$ denotes a tensor outer product. $\\mathbb{C}$ has 21 unique\nelements and enables the calculation of several microstructural parameters.\n\nUsing the cumulant expansion, the diffusion-weighted signal can be approximated\nas\n\n\\begin{align}S \u0007pprox S_0 \\exp \\left(- \\mathbf{b} : \\langle \\mathbf{D} \rangle +\n \frac{1}{2}(\\mathbf{b} \\otimes \\mathbf{b}) : \\mathbb{C} \right) ,\\end{align}\n\nwhere $S_0$ is the signal without diffusion-weighting, $\\mathbf{b}$ is the\nb-tensor used in the acquisition, and $:$ denotes a tensor inner product.\n\nThe model parameters $S_0$, $\\langle\\mathbf{D}\rangle$, and $\\mathbb{C}$ can be\nestimated by solving the following equation:\n\n\\begin{align}S = \beta X,\\end{align}\n\nwhere\n\n\\begin{align}S = \begin{pmatrix} \\ln S_1 \\ \u000bdots \\ \\ln S_n \\end{pmatrix} ,\\end{align}\n\n\\begin{align}\beta = \begin{pmatrix} \\ln S_0 & \\langle \\mathbf{D} \rangle & \\mathbb{C}\n \\end{pmatrix}^ ext{T} ,\\end{align}\n\n\\begin{align}X =\n \begin{pmatrix}\n 1 & -\\mathbf{b}_1^ ext{T} & \frac{1}{2} (\\mathbf{b}_1 \\otimes \\mathbf{b}_1)\n ext{T} \\\n \u000bdots & \u000bdots & \u000bdots \\\n 1 & -\\mathbf{b}_n^ ext{T} & \frac{1}{2} (\\mathbf{b}_n \\otimes \\mathbf{b}_n)\n ^ ext{T}\n \\end{pmatrix} ,\\end{align}\n\nwhere $n$ is the number of acquisitions and $\\langle\\mathbf{D}\rangle$,\n$\\mathbb{C}$, $\\mathbf{b}_i$, and $(\\mathbf{b}_i \\otimes \\mathbf{b}_i)$ are\nrepresented by column vectors using Voigt notation. Estimation of the model\nparameters requires that $ext{rank}(\\mathbf{X}^ ext{T}\\mathbf{X}) = 28$.\nThis can be achieved by combining acquisitions with b-tensors with different\nshapes, sizes, and orientations.\n\nFor details, please see [1]_.\n\n## Usage example\n\nQTI can be fit to data using the module dipy.reconst.qti. Let's start by\nimporting the required modules and functions:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport dipy.reconst.qti as qti" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As QTI requires data with tensor-valued encoding, let's load an example dataset\nacquired with q-space trajectory encoding (QTE):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fdata, fbvals, fbvecs, fmask = get_fnames('qte_lte_pte')\ndata, affine = load_nifti(fdata)\nmask, _ = load_nifti(fmask)\nbvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)\nbtens = np.array(['LTE' for i in range(61)] + ['PTE' for i in range(61)])\ngtab = gradient_table(bvals, bvecs, btens=btens)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The dataset contains 122 volumes of which the first half were acquired with\nlinear tensor encoding (LTE) and the second half with planar tensor encoding\n(PTE). We can confirm this by calculating the ranks of the b-tensors in the\ngradient table.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ranks = np.array([np.linalg.matrix_rank(b) for b in gtab.btens])\nfor i, l in enumerate(['b = 0', 'LTE', 'PTE']):\n print('%s volumes with %s' % (np.sum(ranks == i), l))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have data acquired with tensor-valued diffusion encoding and the\ncorresponding gradient table containing the b-tensors, we can fit QTI to the\ndata as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "qtimodel = qti.QtiModel(gtab)\nqtifit = qtimodel.fit(data, mask)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "QTI parameter maps can accessed as the attributes of qtifit. For instance,\nfractional anisotropy (FA) and microscopic fractional anisotropy (\u03bcFA) maps can\nbe calculated as:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fa = qtifit.fa\nufa = qtifit.ufa" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, let's reproduce Figure 9 from [1]_ to visualize more QTI parameter\nmaps:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "z = 36\n\nfig, ax = plt.subplots(3, 4, figsize=(12, 9))\n\nbackground = np.zeros(data.shape[0:2]) # Black background for figures\nfor i in range(3):\n for j in range(4):\n ax[i, j].imshow(background, cmap='gray')\n ax[i, j].set_xticks([])\n ax[i, j].set_yticks([])\n\nax[0, 0].imshow(np.rot90(qtifit.md[:, :, z]), cmap='gray', vmin=0, vmax=3e-3)\nax[0, 0].set_title('MD')\nax[0, 1].imshow(np.rot90(qtifit.v_md[:, :, z]),\n cmap='gray', vmin=0, vmax=.5e-6)\nax[0, 1].set_title('V$_{MD}$')\nax[0, 2].imshow(np.rot90(qtifit.v_shear[:, :, z]), cmap='gray', vmin=0,\n vmax=.5e-6)\nax[0, 2].set_title('V$_{shear}$')\nax[0, 3].imshow(np.rot90(qtifit.v_iso[:, :, z]),\n cmap='gray', vmin=0, vmax=1e-6)\nax[0, 3].set_title('V$_{iso}$')\n\nax[1, 0].imshow(np.rot90(qtifit.c_md[:, :, z]), cmap='gray', vmin=0, vmax=.25)\nax[1, 0].set_title('C$_{MD}$')\nax[1, 1].imshow(np.rot90(qtifit.c_mu[:, :, z]), cmap='gray', vmin=0, vmax=1)\nax[1, 1].set_title('C$_{\u03bc}$ = \u03bcFA$^2$')\nax[1, 2].imshow(np.rot90(qtifit.c_m[:, :, z]), cmap='gray', vmin=0, vmax=1)\nax[1, 2].set_title('C$_{M}$ = FA$^2$')\nax[1, 3].imshow(np.rot90(qtifit.c_c[:, :, z]), cmap='gray', vmin=0, vmax=1)\nax[1, 3].set_title('C$_{c}$')\n\nax[2, 0].imshow(np.rot90(qtifit.mk[:, :, z]), cmap='gray', vmin=0, vmax=1.5)\nax[2, 0].set_title('MK')\nax[2, 1].imshow(np.rot90(qtifit.k_bulk[:, :, z]),\n cmap='gray', vmin=0, vmax=1.5)\nax[2, 1].set_title('K$_{bulk}$')\nax[2, 2].imshow(np.rot90(qtifit.k_shear[:, :, z]), cmap='gray', vmin=0,\n vmax=1.5)\nax[2, 2].set_title('K$_{shear}$')\nax[2, 3].imshow(np.rot90(qtifit.k_mu[:, :, z]), cmap='gray', vmin=0, vmax=1.5)\nax[2, 3].set_title('K$_{\u03bc}$')\n\nfig.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For more information about QTI, please read the article by Westin et al. [1]_.\n\n### References\n.. [1] Westin, Carl-Fredrik, et al. \"Q-space trajectory imaging for\n multidimensional diffusion MRI of the human brain.\" Neuroimage 135\n (2016): 345-362. https://doi.org/10.1016/j.neuroimage.2016.02.039.\n\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 0 }