The GSoC Experience and Project Summary

Title

Model Fitting using Microstructure Imaging of Crossing (MIX): DIPY

Abstract

Diffusion MRI measures water diffusion in biological tissue, which can be used to probe its microstructure. The most common model for water diffusion in tissue is the diffusion tensor (DT), which assumes a Gaussian distribution. This assumption of Gaussian diffusion oversimplifies the diffusive behavior of water in complex media and is known experimentally to break down for relatively large b-values. DT derived indices, such as mean diffusivity or fractional anisotropy, can correlate with major tissue damage, but lack sensitivity and specificity to subtle pathological changes.

Microstructure Imaging of Crossing (MIX) is versatile and thus suitable to a broad range of generic multicompartment models, in particular for brain areas where axonal pathways cross.

Multicompartment models (assess the variability of diffusion in sub-voxel regions) enable the estimation of more specific indices, such as axon diameter, density, orientation, and permeability, and so potentially give much greater insight into tissue architecture and sensitivity to pathology.

The goal of Model Fitting: Identity which model compartments are essential to explain the data and parameters that are potentially estimable from a particular experiment.

As a part of GSoC, I worked a lot with Model Fitting using the Neurite Orientation Dispersion and Density Imaging (NODDI) model with its implementation using the MIX framework.

Achievements and Benchmarks

The MIX framework for Microstructure Imaging is a very novel and advanced technique, however, required a higher fitting time. As reported by the MIX paper the MATLAB implementation took 191.8 seconds to fit, whereas, the DIPY implementation that I worked on is significantly faster, taking only 10-14 seconds to fit.

This basically means that the current implementation is ~ 20x faster than the state-of-the-art.

Pull Requests with Detailed Descriptions:

https://github.com/nipy/dipy/pull/1600 -> Branch: https://github.com/ShreyasFadnavis/dipy/tree/noddix_speed

https://github.com/nipy/dipy/pull/1614 -> Branch: https://github.com/ShreyasFadnavis/dipy/tree/noddix_gsoc

The above mentioned PRs contain:
Main Code for the NODDIx model (contained in: dipy/dipy/reconst)
Simulation and Fitting Simulated Signal (contained in: dipy/dipy/sims)
Testing the model (contained in: dipy/dipy/reconst/tests)
Example for NODDIx using HCP example (contained in: dipy/docs/examples [only present in the noddix_speed branch])

Link to submission via the PSF GSoC blog : 
https://blogs.python-gsoc.org/shreyas-fadnavis/2018/08/10/the-gsoc-experience-and-project-summary/

Challenging aspects of the work

My project was particularly challenging as it was geared more towards research based on Mathematical Optimization and Model Fitting which I had not worked on before.

I feel that working with DIPY under PSF was one of the most amazing experiences as I really learnt how to overcome the following challenges and contribute quality code to DIPY:

  1. Optimization in Scientific Computing using Approximations and not losing out on model accuracy at the same time.
  2. Better coding practices and Software Design for new Models.
  3. Fitting the Model Parameters in Less time so that it can be used from an analytical standpoint.
  4. Understanding and implementing Cython modules to remove the Python Bottlenecks.
  5. Line-by-line profiling to understand and remove bottlenecks from the code.

As a part of this project, I have worked on [links to the Blogs written at each step]:

Implementing the NODDI model for fitting 2 Fiber Crossings using MIX

The NODDI model consists of normalized signals from intracellular and extracellular compartments.

The estimated dMRI signal Ŝ  comprises of the normalized signals from the following three compartments:

Ŝ  = (1−v_iso)(v_ic * S_ic (𝑂𝐷,𝜃,𝜑) + (1−v_ic)S_ec(d,𝜃,𝜑))+ v_iso * S_iso

where S_ic and S_ec are the normalized signals from intracellular and extracellular compartments respectively.

Parameters to be estimated: Six (v_icv_iso, 𝑂𝐷, d, 𝜃, 𝜑)

Noise: Rician noise needs is added to signal for each substrate for different noise realizations.

To visualize each of the compartments mentioned in the formulation above, please take a look at the following figure for the compartments:

Ref: https://users.fmrib.ox.ac.uk/~michielk/phdthesis.html

The above diagram is representative of a single fiber.

To further expand the above model to 2 fiber orientation crossings, the following formulation has been implemented:

Visual Representation of almost 0 ODs:

Parameters to be estimated: 13 but for implementation, we use only 11 since some of the parameters are reused.

The parameters of the model can be visualized as:

The signal was estimated with the model above and then fitted with SHORE (an analytical basis). We could then  visualize the fiber orientations using SHORE’s Orientation Distribution Functions (ODFs) as follows:

[Note: the steps to implement the above simulations and visualize the signals have been explained in detail below]

Implementation of NODDIx using MIX framework of Optimization

The Microstructure Imaging of Crossings is a novel and robust method using a 3 step optimization process. It enables to fit existing biophysical models with improved accuracy by utilizing the Variable Separation Method (VSM) to distinguish parameters that enter in both, linear and non-linear manner, in the model (Methods). The estimation of non-linear parameters is a non-convex problem and is handled first. This is done by using Differential Evolution since it is more effective in approximating exponential time series models.

The task to estimate linear parameters amounts to a convex problem and can be solved using standard least squares techniques. These parameter estimates provide a starting point for a Trust Region method in search for a refined solution.

4 Steps involved in Implementing MIX:

Step 1 – Variable Separation: The objective function has a separable structure which can be exploited to separate the variables by Variable Separation (VarPro) method. We can rewrite our objective function as a projection using the Moore-Penrose Inverse (Pseudoinverse) and get the variable projection functional.

This is a rather advanced and mathematically well-formed method which makes use of variable projections to transform the complicated computations between the variables of the NODDIx model into a space where they can be fit individually.

Taking advantage of this special structure of the model, the method of variable projections eliminates the linear variables obtaining a somewhat more complicated function that involves only the nonlinear parameters.

This procedure not only reduces the dimension of the parameter space but also results in a better-conditioned problem. The same optimization method applied to the original and reduced problems will always converge faster for the latter.

Further literature for this method can be found here:

Step 2 – Stochastic search for non-linear parameters ‘x’: The objective function is non-convex, particularly of non-linear least-square form. Any gradient based method employed to estimate the parameters will have critical dependence on a good starting point, which is unknown. Alternative approach can be regular grid search, which is time consuming and adds computational burden. This particular type of problem therefore points towards considering stochastic search methods like Differential Evolution (DE). In case of time series analysis, DE can be used efficiently for sum of exponentials functions. DE parameters can be varied for each selected biophysical model and time complexity may change with each choice.

I have written a different blog post for implementation of DE with a detailed explanation of its working and its nature for NODDIx.

Link: https://blogs.python-gsoc.org/shreyas-fadnavis/2018/07/23/differential-evolution-for-noddix

Step 3 – Constrained search for linear parameters ‘f ’: After estimating the parameters ‘x’, estimation of linear parameters ‘f ’ is a constrained linear least-squares estimation problem. I have made use of the cvxpy optimizer to perform the constrained search to find the f’s of the model compartments.

Step 4 – Non-Linear Least Squares Estimation (NLLS) using Trust Region Method: Step 2 and step 3 give a reliable initial guess of both ‘x’ and ‘f ’ as initial value for Trust Region method. This method has been implemented using the Levenberg Marquardt method to perform (NLLS) fitting from SciPy: Optimize module.

Testing and Running my code for NODDIx model in DIPY

Here is the link to the code for NODDIx implemented in DIPY: https://github.com/ShreyasFadnavis/dipy/tree/noddix_speed [1]

Steps to install DIPY:

  1. Install Anaconda: https://anaconda.org. Its free: you just need to create an account and you should be good to go!
  2. In the search of your Windows, search for: Anaconda Prompt.
  • This should open up a terminal where you will need to operate from. (try avoiding the normal windows command prompt as Anaconda uses a virtual environment to do so.)
  1. Install DIPY:
  •   conda install vtk
  •   You may need to install cvxpy (a convex optimization library) using:  
    • conda install -c cvxgrp cvxpy on Windows

Following are the steps you need to take to run the code:

I am not sure how well versed you are with Git, but all you need to do is:

  • Click on the above link in (1).
  • Clone the repository:
    • A green button will appear on the top right corner of the GitHub web-page, which will give you a link. Copy that link.
    • Open the command prompt and run the following 2 commands:
      • git clone <link copied from the above step>
    • This should create a folder in your local file-system which we need to compile using the following command:
      • cd path/to/dipy
      • Once inside the main dipy folder, run the following command to install the repository:
        • pip install –user -e .

Dabbling with the NODDIx code:

  • Anaconda will have already installed an IDE (editor) called Spyder . You can access it from the default Windows search.
  • Once inside spyder, click on open file and navigate to: ./dipy/dipy/reconst
  • There inside reconst is the file which contains the NODDIx code, namely:
    • NODDIx.py (has some code in cython: noddi_speed.pyx)

The sim_noddix.py file is what contains the simulation and fitting of the signal using the model written in NODDIx.py.

All you need to do is hit the RUN button in the navbar of the Spyder IDE and it will return:

  • Actual Parameters (11)
  • Fitted Parameters (11)
  • List of Errors in the Estimation (for each param) (11)
  • Sum of all errors

To run the simulations and visualize the fiber orientation crossings, please navigate to: ./dipy/dipy/sims

Within sims, you should find a file named sim_noddix.py which will simulate the signal and visualize with SHORE. The code contains the error functions and functions to visualize the data. They are as follows:

  • show_with_shore(gtab, reconst_signal) -> to visualize the simulated signal
  • sim_voxel_fit(reconst_signal) -> to fit the simulation of 1 voxel

Fitting the Real Data With NODDIx Model

Now it’s time to see how the model works on real data. To do so, I am working on an Example using a subject from the HCP (Human Connectome Project) data which is publicly available online.

The whole brain takes really long to fit, however, I have made use of a mask to fit only the region of the brain with the corpus callosum to see how the model works.

Here is the link you can use to test the model with your own data by just replacing the file paths: https://github.com/ShreyasFadnavis/dipy/blob/noddix_speed/doc/examples/example_noddi_hcp.py

Next Steps:

Approximately ~ 3000 lines of code have been written so far for the NODDIx model. Next steps include:

  • Creating an Example for the Gallery of DIPY to help new users use the MIX module.
  • Extend the 2 crossing model to 1 and 3 crossings.

Differential Evolution for NODDIx

In the previous post, we took a loot at how the NODDIx model could be made faster with Cython. In this post, we will take a look at how we need to set up the differential evolution with for optimizing and finding the minima for the parameters of the NODDIx model. The code to this function can be found at the branch here: https://github.com/ShreyasFadnavis/dipy/tree/noddix_speed .

In this post I will be explaining about how Differential Evolution works and in what context we have used it for optimization:

Setting the Parameters in SciPy for Differential Evolution Optimization:

On evaluating this, we see that the NODDIx model converges in 81 steps within 6 seconds approximately when f(x)  = 1.10912e-07

The accuracy we get is:

References: 

[1]. Differential Evolution: A Practical Approach to Global Optimization https://www.springer.com/us/book/9783540209508

[2]. Numerical optimization by differential evolution https://www.youtube.com/watch?v=UZGiapWcoA4&t=2151s

 

Cythonizing Bottlenecks in the NODDIx Model

+

Over the past few days, I have been trying to remove the bottlenecks in the Python code by vector implementation using Numpy extensions and a variety of other code optimization techniques.

However, I have come to realize that even Numpy at times does not solve problems like Global Interpreter Lock and Bounds Checking. Cython, on the other hand gives a very nice interface with the lower level C implementations by inserting an intermediate layer between Python and its interpreted code in C. The following diagram will help demonstrate what I mean:

Courtesy: https://www.software.ac.uk/blog/2017-07-25-speeding-python-cython-some-thoughts-python-programmers

Typed MemoryViews in Cython… Blissful interaction with NumPy Arrays!

So, it is often the case that we need to cythonize only a part of our code and not all of it. This requires smooth interaction of other functions with the cythonized functions. MemoryView provides a really good way of  doing this by allowing efficient access to memory buffers, such as those underlying NumPy arrays, without incurring any Python overhead. Memoryviews are similar to the current NumPy array buffer support (np.ndarray[np.float64_t, ndim=2]), but they have more features and cleaner syntax.

Memoryviews are more general than the old NumPy array buffer support, because they can handle a wider variety of sources of array data. For example, they can handle C arrays and the Cython array type (Cython arrays).

A memoryview can be used in any context (function parameters, module-level, cdef class attribute, etc) and can be obtained from nearly any object that exposes writable buffer through the `PEP 3118`_ buffer interface.

Following is an example of the Legendre Gauss Integral in Python and its Cythonized equivalent following it:

The Cythonized Version:

Need to Zoom in for this one..

The entire code for along with other cythonized functions can be found on this branch: noddix_speed !

References:

https://cython.readthedocs.io/en/latest/src/userguide/memoryviews.html

Visualizing the NODDI Signal with GQI and Shore 3D models

In the previous post, we looked at how we can improve the speedup of fitting the simulated signal with Cython and what tweaks did we perform with the Legendre-Gauss Integration.

In this post, we will look at how the simulated signal can be visualized using the SHORE-3D Model and the GQI Model from the DIPY Library:

SHORE-3D Model: Exploits the ability of Compressed Sensing (CS) to recover the whole 3D Diffusion MRI (dMRI) signal from a limited number of samples while efficiently recovering important diffusion features such as the Ensemble Average Propagator (EAP) and the Orientation Distribution Function (ODF). 

Following are the crossings generated from the NODDIx model with crossings. The model helps visualize in the form of lobes:

GQI (Generalized q-sampling imaging):

Based on the Fourier transform relation between diffusion magnetic resonance (MR) signals and the underlying diffusion displacement, a new relation is derived to estimate the spin distribution function (SDF) directly from diffusion MR signals. This relation leads to an imaging method called generalized -sampling imaging (GQI), which can obtain the SDF from the shell sampling scheme used in -ball imaging (QBI) or the grid sampling scheme used in diffusion spectrum imaging (DSI). The GQI method can be applied to grid or shell sampling schemes and can provide directional and quantitative information about the crossing fibers.

The code for these visualizations can be found here.

 

References:

  1. “Generalized q-Sampling Imaging”, IEEE Trans Med Imaging. 2010 Sep;29(9):1626-35. doi: 10.1109/TMI.2010.2045126. Epub 2010 Mar 18
  2. Continuous diffusion signal, EAP and ODF estimation via Compressive Sensing in diffusion MRI”, Merlet, Sylvain. L et al., Medical Image Analysis , Volume 17 , Issue 5 , 556 – 572

Speeding up the Legendre-Gauss Integration

As a part of the NODDI model, we were hitting a bottleneck while fitting due the Legendre integral required for computing the Watson Distribution. For those who are new to Numerical Methods, here are some of my notes that will help you understand how to Legendre-Gauss is computed using Legendre Polynomials:

We first start with evaluating the Polynomial between [-1 to 1] interval and the n extend it to any interval [a, b].

Trick to convert the above evaluation between any defined interval [a, b]...

Following is the link to the Cythonized code for Legendre Integrals: Branch Link

 

Simulating and Fitting the Signal using NODDIx

In the last post we took a loot at the NODDI model that provides neurite density and orientation dispersion estimates by disentangling two key contributing factors to FA. This enables the analysis of each factor individually.

The following image helps summarize the model components:

Source
Courtesy: Zhang, Hui et. al. Bingham-NODDI: Mapping anisotropic orientation dispersion of neurites using diffusion MRI. NeuroImage. 133. 10.1016/j.neuroimage.2016.01.046.

As mentioned above and in the previous post, NODDI model consists of 3 major sub-models for fitting the data in:

  • Intra-Cellular Region
  • Extra-Cellular Region
  • Cerebrospinal Fluid (CSF Area)

A major shortcoming of the NODDI model was that it could fit only one fiber at a time for neurite orientation patterns observed in brain tissue that include:

  • Highly coherently oriented white matter structures, such as the corpus
    callosum.
  • White matter structures composed of bending and fanning axons, such as the centrum semiovale.
  • The cerebral cortex and subcortical gray matter structures characterized by sprawling dendritic processes in all directions.

However, with the NODDIx model, we can fit 2 fibers that are in the form of crossings and visualize them using Microstructure Imaging. An example of the same would be as follows:

Courtesy: Hamza Farooq, et. al., Microstructure Imaging of Crossing (MIX) White Matter Fibers from diffusion MRI

Now that I have the 1st draft of the code written (which can bee found here and on my master), its time to simulate some data and try to fit it using the NODDIx model.

The code has been written in such a way that we can make use of the dame functions in the code to generate a signal per voxel and make sure that the model is fitting properly.

To do so, we fix the input parameters that we are giving to generate the signal and try to estimate them using the NODDIx model. Note, we will have 2 fiber crossings and we can explicitly test for different angles between the fibers with pre-defined volume fractions.

For simplicity and easy understanding of the estimates, we set the parameters as follows:

The volume fractions have been set to equal:

Intracellular Volume Fraction 1: 0.2
Intracellular Volume Fraction 2: 0.2
Extracellular Volume Fraction 1: 0.2
Extracellular Volume Fraction 2: 0.2
CSF Volume Fraction: 0.2

Lets now fix the Orientation Dispersions (OD) and the Thetas and Phis:

OD1: 0.2
Theta1: 0.72
Phi1: 1.57

OD2: 0.24
Theta2: 0.72
Phi2: 1.57

For both the fibers, the angle between the Theta and the Phi has been set to 41 degrees approximately.

These are the 11 parameters that the model is estimating… Lets see how the estimates look on the simulation data!

[Note: We expect the estimates to be almost the same as the input values…]

RESULTS:

The estimated volume fractions are:

Intracellular Volume Fraction 1: 0.19739277
Intracellular Volume Fraction 2: 0.1947075
Extracellular Volume Fraction 1: 0.20260892
Extracellular Volume Fraction 2: 0.2052905
CSF Volume Fraction: 0.20000128

The volume fractions seem to be estimated almost perfectly! Lets take a look at the ODs, Thetas and Phis;

OD1: 0.20045016
Theta1: 0.72000357
Phi1: 1.56999562

OD2: 0.24062461
Theta2: 0.71999553
Phi2: 1.5700052

Yep, they look almost the same!

So basically, the NODDIx Model is up and running!

[Note: Code for the above simulation and fit can be found here]

The speed of the fit needs to be taken care of though!

In the next post, I will try to speedup the code and have a rigorous code profiling!

End of 2 Weeks of GSoC!

GSoC with the Python Software foundation has been an incredible journey so far and has made me really productive a programmer. Here are the links to the work that I have done so far and is obviously work in progress!

I maintain a repo called dipy-tests where I perform all my experiments. Here is the link to it:
https://github.com/ShreyasFadnavis/dipy-tests
I have started working on implementing the NODDIx model for Microstructure Imaging which was mentioned in the previous post!
Here is the link to the code that I have written so far:
Following are the commits that I have made to the my master nipy/dipy repo over the past 2 weeks: https://github.com/ShreyasFadnavis/dipy

I have also made some commits and got some pull requests merged with the main nipy/dipy repo of my PSF sub-organization. Here is the link to it:

Watch out for my next post!
Will be working on simulating data using the NODDIx model and performing some model fitting on it (some rigorous testing expected!)
GitHub: https://github.com/ShreyasFadnavis

Neurite Orientation Dispersion and Density Imaging

We need to look at Neurite Orientation Dispersion and Density Imaging (NODDI) first!

NODDI is a method of quantifying the morphology of neurites (→ Projections of neurons and dendrites, collectively) using branching complexity of the dendritic trees in terms of dendritic density

Remember that we are looking at dMRI which measures the Displacement Pattern of water molecules undergoing diffusionNODDI is a Tissue Model which makes use of Orientation Dispersion Index as a Summary Statistic. This helps quantify Angular Variation of Neurite Orientation.

The NODDI tissue model:

So the above diagram has a tiny equation at the three types of microstructural environments. The water diffusion in each of the environments behaves differently and we will look at them in a little more depth below. (BTW: CSF == CerebroSpinal Fluid).

Intra-cellular Model

Intra-cellular compartment → space bounded by the membrane of Neurites: Model this as Set of Sticks. The normalized signal for this model  is given by:

Where,

q → Gradient Direction

b  → b-val 

f(n)dn → probability of finding sticks along orientation ‘n’

The exponential term   → Signal Attenuation due to the intrinsic diffusivity of the sticks

Now, 

The f(n): the Orientation Distribution function is modeled as a Watson Distribution. This distribution is the simplest distribution that can capture the dispersion in orientations.

where, 

M is a confluent hypergeometric function. 

𝛍 = mean dispersion 

қ = concentration parameter that measures the dispersion about 𝛍 
EXtra-cellular Model

This, as the name suggests takes care of the space around the neurites. The neurites hinder the diffusion, but does not restrict it in any way. Therefore we can model it using a ‘Gaussian Anisotropic Distribution‘.

The signal for this environment of the model looks something like:

 

 

Here, D(n) is a cylindrical symmetric tensor with the principal direction of diffusion ‘n’.  But, now we need to consider 2 directions of diffusion, perpendicular and parallel diffusivities: d⊥ and d∥

The parallel diffusivity is the same as the intrinsic
free diffusivity of the intra-cellular compartment; the perpendicular
diffusivity is set with a simple tortuosity model as
d⊥ = d∥ (1 − ν_ic), where ν_ic is the intra-cellular volume fraction.

where,

Csf Model

The CSF compartment models the space occupied by cerebrospinal fluid and is modeled as isotropic Gaussian diffusion with diffusivity d_iso. But in this model we make use of the Orientation Dispersion index instead of қ as follows:

OD = 2/pi * arctan(1/қ)

________________________________________________________

References:

[1] Zhang, H., Schneider, T., Wheeler-kingshott, C. A., & Alexander, D. C. (2012). NeuroImage NODDI : Practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage, 61(4), 1000–1016. https://doi.org/10.1016/j.neuroimage.2012.03.072

[2] Farooq, H., Xu, J., Nam, J. W., Keefe, D. F., Yacoub, E., Georgiou, T., & Lenglet, C. (2016). Microstructure Imaging of Crossing (MIX) White Matter Fibers from diffusion MRI. Scientific Reports, 6(September), 1–9. https://doi.org/10.1038/srep38927

[3] Ferizi, U., Schneider, T., Panagiotaki, E., Nedjati-Gilani, G., Zhang, H., Wheeler-Kingshott, C. A. M., & Alexander, D. C. (2014). A ranking of diffusion MRI compartment models with in vivo human brain data. Magnetic Resonance in Medicine, 72(6), 1785–1792. https://doi.org/10.1002/mrm.25080

Simulating dMRI Data using CAMINO

     +     

Data simulation forms a crucial component of any data-driven experiment which deals with model fitting as it is the ground-truth that we will be comparing against. In my project, I will be working with the following 2 tools for data simulation:

  • UCL Camino Diffusion MRI Toolkit
  • DIPY Simulations (… Obviously!)

This post will cover Camino 1st and I aim to get into DIPY in the next post!

The most confusing part about the Camino documentation  is understanding what the ‘scheme’ file is really made up of, because it needs to be passed as a parameter to the ‘datasynth’ command which we will look at for data simulation.

Scheme files accompany DWI data and describe imaging parameters that are used in image processing. For most users, we require the gradient directions and b-value of each measurement in the sequence.

Once you have this information, you can use the CAMINO commands described below to generate scheme files.

  • Comments are allowed, the line must start with ‘#’
  • The first non-comment line must be a header stating “VERSION: <version>”. In our case:
VERSION: BVECTOR
  • After removing comments and the header, measurements are described in order, one per line. The order must correspond to the order of the DWI data.
  • Entries on each line are separated by spaces or tabs.

The BVECTOR is the most common scheme format. Each line consists of four values: the (x, y, z) components of the gradient direction followed by the b-value. For example:

   # Standard 6 DTI gradient directions, [b] = s / mm^2
  VERSION: BVECTOR
   0.000000   0.000000   0.000000   0.0
   0.707107   0.000000   0.707107   1.000E03
  -0.707107   0.000000   0.707107   1.000E03
   0.000000   0.707107   0.707107   1.000E03
   0.000000   0.707107  -0.707107   1.000E03
   0.707107   0.707107   0.000000   1.000E03
  -0.707107   0.707107   0.000000   1.000E03

If the measurement is unweighted, its gradient direction should be zero. Otherwise, the gradient directions should be unit vectors, followed by a scalar b-value. The b-value can be in any units. Units are defined implicitly, in the above example we have used s / mm^2. The choice of units affects the scale of the output tensors, if we used this scheme file we would get tensors in units of mm^2 / s. We could change the units of b to s / m^2 by scaling the b-values by 1E6. Our reconstructed tensors would then be specified in units of m^2 / s.

Finding the information for the scheme file

The best way to find the information for your scheme file is to talk to the person who programmed your MRI sequence. There is software that can help you recover them from DICOM or other scanner-specific data formats. The dcm2nii program will attempt to recover b-values and vectors in FSL format.

Converting to Camino format

If you have a list of gradient directions, you can convert them to Camino format by hand or by using pointset2scheme. If you have FSL style bval and bvec files, you can use fsl2scheme. See the man pages for more information.

 

Simulating the Data

Finally! Now that we know what the scheme files are, lets look at how to simulate the voxels…

I will be making use of the 2 utilities which I feel are relevant to my project and will test the simulation functionalities using the 59.scheme file which is present on the Camino website tutorial.

1. Synthesis Using Analytic Models

 

This uses Camino to synthesize diffusion-weighted MRI data with the white matter analytic models.

The method is explained in detail in (Panagiotaki et al NeuroImage 2011, doi:10.1016/j.neuroimage.2011.09.081).

The following example synthesizes data using the three-compartment model “ZeppelinCylinderDot“, which has an intra-axonal compartment of single radius, a cylindrically symmetric tensor for the extra-axonal space and a stationary third compartment.

Example:

datasynth -synthmodel  compartment 3 CYLINDERGPD 0.6 1.7E-9 0.0  0.0  4E-6 zeppelin 0.1 1.7E-9 0.0 0.0 2E-10  Dot -schemefile 59.scheme -voxels 1 -outputfile ZCD.Bfloat
2. Crossing cylinders using Monte Carlo Diffusion Simulator

This simulator allows the simulation of diffusion from simple to extremely complex diffusion environments, called “substrates“. We will be looking at the Crossing fibres substrates as of now.

A substrate is envisaged to sit inside a single voxel, with spins diffusing across it. The boundaries of the voxel are usually periodic so that the substrate defines an environment made up of an infinite, 3D array of whatever you specify. The measurement model in the simulation does not capture the trade-off between voxel size and SNR and hence simulation "voxels" can be quite a bit smaller than those in actual scans. This simulation is, and has always been, intended as a tool to simulate signals due to sub-voxel structure, rather than large spatially-extended structures. [- UCL Camino Docs]
Crossing Cylinders

A situation that is often of interest in diffusion MR research is where we have more than one principle fibre direction. The simulation is able to model crossing fibres with a specified crossing angle. This substrate contains two populations of fibres in interleaved planes. One population is parallel to the z-axis and another is rotated about the y-axis by a given angle with respect to the first.

Cylinders on this substrate are arranged in parallel in the xz-plane in parallel layers one cylinder thick. i.e. a plane of cylinders parallel to the z-axis, with a rotated with respect the first, then another parallel z-axis and so on. Cylinders are all of a constant radius.

An example command to use here is:
datasynth -walkers 100000 -tmax 1000 -voxels 1 -p 0.0 -schemefile 59.scheme -initial uniform -substrate crossing -crossangle 0.7854 -cylinderrad 1E-6 -cylindersep 2.1E-6 > crossingcyls45.bfloat


Here we’ve specified a crossing substrate. The crossing angle is specified in radians (NOT degrees) using the -crossangle: 0.7854. This is approximately pi/4, or 45 degrees. The crossing angle can take on any value, just make sure you use radians!

 

REFERENCES:

[1] http://camino.cs.ucl.ac.uk/index.php

[2] Panagiotaki et al NeuroImage 2011, doi:10.1016/j.neuroimage.2011.09.081

Understanding Neuroimaging Data

        

I realized that before jumping into any kind of data analyses or modeling, I should devote a small post about what ‘Brain Data’ looks like.

In the future posts, we will be looking at a lot of Neuro data and I felt that this post would be a crucial component for the same. An NIPY-DIPY (and Nibabel) image is the association of three things:

  • The image data array: a 3D or 4D array of image data
  • An affine array that tells you the position of the image array data in a reference space.
  • Image metadata (data about the data) describing the image, usually in the form of an image header.
import nibabel as nib
import matplotlib.pyplot as plt

epi_img = nib.load('C:/Users/Shreyas/Desktop/DATA_understanding/someones_epi.nii')
epi_img_data = epi_img.get_data()
epi_img_data.shape


""" Function to display row of image slices """
def show_slices(slices):
 fig, axes = plt.subplots(1, len(slices))
 for i, slice in enumerate(slices):
 axes[i].imshow(slice.T, cmap="gray", origin="lower")
 
slice_0 = epi_img_data[26, :, :]
slice_1 = epi_img_data[:, 30, :]
slice_2 = epi_img_data[:, :, 16]
show_slices([slice_0, slice_1, slice_2])
plt.suptitle("Center slices for EPI image")

In the above grayscale image that looks like the cross-sections of the brain, each pixel is a ‘Voxel’, i.e. a pixel with Volume: Where black areas denote a minimum and white areas denote maximum values. A 3D array of such voxels is forms a voxel array. The following will help us understand how to get the voxel coordinates of the center slices and the value at that voxel!

n_i, n_j, n_k = epi_img_data.shape
center_i = (n_i - 1) // 2
center_j = (n_j - 1) // 2
center_k = (n_k - 1) // 2

centers = [center_i, center_j, center_k]
print("Co-ordinates in the voxel array: ", centers)

# O/P: Co-ordinates in the voxel array: [26, 30, 16]

center_vox_value = epi_img_data[center_i, center_j, center_k]
print(center_vox_value)
# O/P: 81.54928779602051

Note that the voxel coordiantes do not tell us anything about where the data came from, as in, the point of view of the scanner/ position of the subject/ left-right brain imaging,… etc. Therefore we need to perform some kind of affine transformations to these coordinates to get them into the subject-scanner reference space. 

(x, y, z) = f(i, j, k)

The above example, though of EPI (Echo Planar Imaging), gives a brief overview of what these data mean.

In Diffusion MRI (dMRI) usually we use three types of files, a Nifti file with the diffusion weighted data, and two text files: one with b-values and one with the b-vectors. In DIPY we provide tools to load and process these files and we also provide access to publicly available datasets for those who haven’t acquired yet their own datasets.

There are some brilliant tutorials which provide an abundance of high quality information:

http://nipy.org/dipy/examples_index.html#examples

http://nipy.org/nibabel/tutorials.html#tutorials