examples/tutorial_setting_up_acquisition_scheme.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading acquisition parameters and setting up an acquisition scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "dmipy saves all acquisition parameter information inside an DmipyAcquisitionScheme instance. Acquisition schemes can be set up directly from b-values, q-values or gradient strength; can be converted from a dipy GradientTable; or be imported from a Camino-style schemefile."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up a Dmipy acquisition scheme by loading acquisition parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To set up an acquisition scheme directly from b-values one can load them directly from text files. The b-values and gradient directions, along with the used pulse duration time $\\delta$ and pulse separation time $\\Delta$ of the acquisition must be known. The dmipy toolbox uses SI units, so be careful, as bvalues are typically saved in $s/mm^2$, but as a rule we need them in $s/m^2$.\n",
    "\n",
    "As an example we load the acquisition parameters of the WU-MINN Human Connectome Project:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the necessary modules\n",
    "from dmipy.core import modeling_framework\n",
    "from dmipy.core.acquisition_scheme import acquisition_scheme_from_bvalues\n",
    "from os.path import join\n",
    "import numpy as np\n",
    "\n",
    "# the HCP acquisition parameters are saved in the following toolbox path:\n",
    "acquisition_path = modeling_framework.GRADIENT_TABLES_PATH\n",
    "\n",
    "# we can then load the parameters themselves and convert them to SI units:\n",
    "bvalues = np.loadtxt(join(acquisition_path, 'bvals_hcp_wu_minn.txt'))  # given in s/mm^2\n",
    "bvalues_SI = bvalues * 1e6  # now given in SI units as s/m^2\n",
    "gradient_directions = np.loadtxt(join(acquisition_path, 'bvecs_hcp_wu_minn.txt'))  # on the unit sphere\n",
    "\n",
    "# The delta and Delta times we know from the HCP documentation in seconds\n",
    "delta = 0.0106  \n",
    "Delta = 0.0431 \n",
    "\n",
    "# The acquisition scheme we use in the toolbox is then created as follows:\n",
    "acq_scheme = acquisition_scheme_from_bvalues(bvalues_SI, gradient_directions, delta, Delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the acquisition scheme, some basic information can be called using the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acquisition scheme summary\n",
      "\n",
      "total number of measurements: 288\n",
      "number of b0 measurements: 18\n",
      "number of DWI shells: 3\n",
      "\n",
      "shell_index |# of DWIs |bvalue [s/mm^2] |gradient strength [mT/m] |delta [ms] |Delta[ms] |TE[ms]\n",
      "0           |18        |0               |0                        |10.6       |43.1      |N/A  \n",
      "1           |90        |1000            |56                       |10.6       |43.1      |N/A  \n",
      "2           |90        |2000            |79                       |10.6       |43.1      |N/A  \n",
      "3           |90        |3000            |97                       |10.6       |43.1      |N/A  \n"
     ]
    }
   ],
   "source": [
    "acq_scheme.print_acquisition_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the summary shows, the function automatically separates the different acquisition shells and detects b0 measurements.  \n",
    "Other metrics such as qvalues, gradient strengths, diffusion times and the shell_indices can also be explicitly called as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "acq_scheme.bvalues;  # bvalues in s/m^2\n",
    "acq_scheme.gradient_directions;  # gradient directions on the unit sphere\n",
    "acq_scheme.gradient_strengths;  # the gradient strength in T/m\n",
    "acq_scheme.qvalues;  # describes the diffusion sensitization in 1/m\n",
    "acq_scheme.tau;  # diffusion time as Delta - delta / 3. in seconds\n",
    "acq_scheme.shell_indices;  # index assigned to each shell. 0 is assigned to b0 measurements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up a Dmipy acquisition scheme when delta and Delta are unknown\n",
    "It is possible we are using clinical data where the delta and Delta parameters are unknown or hard to find out for the acquisition parameters.\n",
    "\n",
    "Dmipy allows you to set up an acquisition scheme using only the bvalues and gradient directions (so just omitting delta and Delta), which will still allow you to use most of Dmipy's models. If a model is used that actually needs this information, this will be made clear with an appropriate error message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acquisition scheme summary\n",
      "\n",
      "total number of measurements: 288\n",
      "number of b0 measurements: 18\n",
      "number of DWI shells: 3\n",
      "\n",
      "shell_index |# of DWIs |bvalue [s/mm^2] |gradient strength [mT/m] |delta [ms] |Delta[ms] |TE[ms]\n",
      "0           |18        |0               |N/A                      |N/A        |N/A       |N/A  \n",
      "1           |90        |1000            |N/A                      |N/A        |N/A       |N/A  \n",
      "2           |90        |2000            |N/A                      |N/A        |N/A       |N/A  \n",
      "3           |90        |3000            |N/A                      |N/A        |N/A       |N/A  \n"
     ]
    }
   ],
   "source": [
    "acq_scheme = acquisition_scheme_from_bvalues(bvalues_SI, gradient_directions)\n",
    "acq_scheme.print_acquisition_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that now the gradient strength cannot be calculated from the bvalues and delta/Delta are unknown, but this has no effect on how shells are separated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up a Dmipy acquisition scheme by converting a Dipy GradientTable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dipy's GradientTable and Dmipy's AcquisitionScheme both store and verify a dataset's acquisition parameters to later be used in model fitting. For convenience, we provide the possibility to convert a dipy table to a dmipy scheme using a conversion function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acquisition scheme summary\n",
      "\n",
      "total number of measurements: 288\n",
      "number of b0 measurements: 18\n",
      "number of DWI shells: 3\n",
      "\n",
      "shell_index |# of DWIs |bvalue [s/mm^2] |gradient strength [mT/m] |delta [ms] |Delta[ms] |TE[ms]\n",
      "0           |18        |0               |0                        |10.6       |43.1      |N/A  \n",
      "1           |90        |1000            |56                       |10.6       |43.1      |N/A  \n",
      "2           |90        |2000            |79                       |10.6       |43.1      |N/A  \n",
      "3           |90        |3000            |97                       |10.6       |43.1      |N/A  \n"
     ]
    }
   ],
   "source": [
    "from dipy.core.gradients import gradient_table\n",
    "from dmipy.core.acquisition_scheme import gtab_dipy2dmipy\n",
    "gtab_dipy = gradient_table(bvalues, gradient_directions, big_delta=Delta, small_delta=delta)\n",
    "acq_scheme_mipy = gtab_dipy2dmipy(gtab_dipy)\n",
    "acq_scheme_mipy.print_acquisition_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the result is exactly the same. Note that the conversion requires that big_delta and small_delta are both set in the gradient_table!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up a Dmipy acquisition scheme from a Camino schemefile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to create a dmipy acquisition scheme directly from a Camino schemefile. Below we load a Camino schemefile containing the same wu-minn HCP acquisition parameters directly from the file itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acquisition scheme summary\n",
      "\n",
      "total number of measurements: 288\n",
      "number of b0 measurements: 18\n",
      "number of DWI shells: 3\n",
      "\n",
      "shell_index |# of DWIs |bvalue [s/mm^2] |gradient strength [mT/m] |delta [ms] |Delta[ms] |TE[ms]\n",
      "0           |18        |0               |0                        |10.6       |43.1      |65.3 \n",
      "1           |90        |1000            |56                       |10.6       |43.1      |65.3 \n",
      "2           |90        |2000            |79                       |10.6       |43.1      |65.3 \n",
      "3           |90        |3000            |97                       |10.6       |43.1      |65.3 \n"
     ]
    }
   ],
   "source": [
    "from dmipy.core.acquisition_scheme import acquisition_scheme_from_schemefile\n",
    "acq_scheme_mipy = acquisition_scheme_from_schemefile(\n",
    "    join(acquisition_path, \"schemefile_hcp_wu_minn.scheme1\")\n",
    ")\n",
    "acq_scheme_mipy.print_acquisition_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible so save an AcquisitionScheme to a schemefile file as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "acq_scheme_mipy.to_schemefile(join(acquisition_path, \"schemefile_hcp_wu_minn.scheme1\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the acquisition scheme has no TE indicated, then in the schemefile it will default to TE=Delta + 2delta + 0.001 seconds"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}