Wrapping up : Google Summer of Code 2016

Software development is challenging and developing a robust package for use in the real world is a completely different experience than writing code for yourself. Google Summer of Code 2016 was a valueable learning experience which introduced me to Free and Open Source Software. I learnt the importance of writing robust tested code which is easy to read and understand for other developers. GSoC also gave me an opportunity to interact with experienced mentors who were very patient in answering my doubts and suggesting good coding practices. With thorough code reviews and regular hangout sessions, my mentors went out of their way to help me with my project. The entire process of writing a proposal, getting involved with a new community, developing a package which will be used by a large number of users, testing the code and writing examples for it has given me the confidence to take up any project and work independently. I am sure this is the first of many more projects and contributions that I will get involved with in the FOSS world.

One of the most interesting part of GSoC for me was to get to know and work with people from different parts of the world. This introduced me to a completely new way of working in a team remotely and I am now looking to get involved in more such projects.

As someone who had very little previous experience in software developement one of the most difficut task was getting selected for the project but I am very grateful that I got a chance. I already had some experience with Mathematical Modelling and as a pre-final year student in Physics this project has helped me explore the topic thorougly. I am sure this will be a great help in deciding my final year thesis and whenever I develop code from now on, I will have the mindset of an open source developer and try to write code such that it can be used freely and developed further.

Search mailing list archive

Search this site

Intravoxel incoherent motion¶

The intravoxel incoherent motion (IVIM) model describes diffusion and perfusion in the signal acquired with a diffusion MRI sequence that contains multiple low b-values. The IVIM model can be understood as an adaptation of the work of Stejskal and Tanner [Stejskal65] in biological tissue, and was proposed by Le Bihan [LeBihan84]. The model assumes two compartments: a slow moving compartment, where particles diffuse in a Brownian fashion as a consequence of thermal energy, and a fast moving compartment (the vascular compartment), where blood moves as a consequence of a pressure gradient. In the first compartment, the diffusion coefficient is \(\mathbf{D}\) while in the second compartment, a pseudo diffusion term \(\mathbf{D^*}\) is introduced that describes the displacement of the blood elements in an assumed randomly laid out vascular network, at the macroscopic level. According to [LeBihan84], \(\mathbf{D^*}\) is greater than \(\mathbf{D}\).
The IVIM model expresses the MRI signal as follows:

\[S(b)=S_0(fe^{-bD^*}+(1-f)e^{-bD})\]

where \(\mathbf{b}\) is the diffusion gradient weighing value (which is dependent on the measurement parameters), \(\mathbf{S_{0}}\) is the signal in the absence of diffusion gradient sensitization, \(\mathbf{f}\) is the perfusion fraction, \(\mathbf{D}\) is the diffusion coefficient and \(\mathbf{D^*}\) is the pseudo-diffusion constant, due to vascular contributions.
In the following example we show how to fit the IVIM model on a diffusion-weighteddataset and visualize the diffusion and pseudo diffusion coefficients. First, we import all relevant modules:

import matplotlib.pyplot as plt

from dipy.reconst.ivim import IvimModel
from dipy.data.fetcher import read_ivim

We get an IVIM dataset using Dipy’s data fetcher read_ivim. This dataset was acquired with 21 b-values in 3 different directions. Volumes corresponding to different directions were registered to each other, and averaged across directions. Thus, this dataset has 4 dimensions, with the length of the last dimension corresponding to the number of b-values. In order to use this model the data should contain signals measured at 0 bvalue.

img, gtab = read_ivim()

The variable img contains a nibabel NIfTI image object (with the data) and gtab contains a GradientTable object (information about the gradients e.g. b-values and b-vectors). We get the data from img using read_data.

data = img.get_data()
print('data.shape (%d, %d, %d, %d)' % data.shape)

The data has 54 slices, with 256-by-256 voxels in each slice. The fourth dimension corresponds to the b-values in the gtab. Let us visualize the data by taking a slice midway(z=27) at \(\mathbf{b} = 0\).

z = 27
b = 20

plt.imshow(data[:, :, z, b].T, origin='lower', cmap='gray',
           interpolation='nearest')
plt.axhline(y=100)
plt.axvline(x=170)
plt.savefig("ivim_data_slice.png")
plt.close()

Heat map of a slice of data

The region around the intersection of the cross-hairs in the figure contains cerebral spinal fluid (CSF), so it so it should have a very high \(\mathbf{f}\) and \(\mathbf{D^*}\), the area between the right and left is white matter so that should be lower, and the region on the right is gray matter and CSF. That should give us some contrast to see the values varying across the regions.

x1, x2 = 160, 180
y1, y2 = 90, 110
data_slice = data[x1:x2, y1:y2, z, :]

plt.imshow(data[x1:x2, y1:y2, z, b].T, origin='lower',
           cmap="gray", interpolation='nearest')
plt.savefig("CSF_slice.png")
plt.close()

Heat map of the CSF slice selected.

Now that we have prepared the datasets we can go forward with the ivim fit. Instead of fitting the entire volume, we focus on a small section of the slice as selected aboove, to fit the IVIM model. First, we instantiate the Ivim model. Using a two-stage approach: first, a tensor is fit to the data, and then initial guesses for the parameters \(\mathbf{S_{0}}\) and \(\mathbf{D}\) obtained from this this tensor by _estimate_S0_D is used as the starting point for the non-linear fit of IVIM parameters using Scipy’s leastsq or least_square function depending on which Scipy version you are using. All initializations for the model such as split_b are passed while creating the IvimModel. If you are using Scipy 0.17, you can also set bounds by setting bounds=([0., 0., 0., 0.], [np.inf, 1., 1., 1.])) while initializing the IvimModel. It is recommeded that you upgrade to Scipy 0.17 since the fitting results might at times return values which do not make sense physically. (For example a negative \(\mathbf{f}\))

ivimmodel = IvimModel(gtab)

To fit the model, call the fit method and pass the data for fitting.

ivimfit = ivimmodel.fit(data_slice)

The fit method creates a IvimFit object which contains the parameters of the model obtained after fitting. These are accessible through the model_params attribute of the IvimFit object. The parameters are arranged as a 4D array, corresponding to the spatial dimensions of the data, and the last dimension (of length 4) corresponding to the model parameters according to the following order : \(\mathbf{S_{0}, f, D^*, D}\).

ivimparams = ivimfit.model_params
print("ivimparams.shape : {}".format(ivimparams.shape))

As we see, we have a 20x20 slice at the height z = 27. Thus we have 400 voxels. We will now plot the parameters obtained from the fit for a voxel and also various maps for the entire slice. This will give us an idea about the diffusion and perfusion in that section. Let(i, j) denote the coordinate of the voxel. We have already fixed the z component as 27 and hence we will get a slice which is 27 units above.

i, j = 10, 10
estimated_params = ivimfit.model_params[i, j, :]
print(estimated_params)

Next, we plot the results relative to the model fit. For this we will use the predict method of the IvimFit object to get the estimated signal.

estimated_signal = ivimfit.predict(gtab)[i, j, :]

plt.scatter(gtab.bvals, data_slice[i, j, :],
            color="green", label="Actual signal")
plt.plot(gtab.bvals, estimated_signal, color="red", label="Estimated Signal")
plt.xlabel("bvalues")
plt.ylabel("Signals")

S0_est, f_est, D_star_est, D_est = estimated_params
text_fit = """Estimated \n S0={:06.3f} f={:06.4f}\n
            D*={:06.5f} D={:06.5f}""".format(S0_est, f_est, D_star_est, D_est)

plt.text(0.65, 0.50, text_fit, horizontalalignment='center',
         verticalalignment='center', transform=plt.gca().transAxes)
plt.legend(loc='upper right')
plt.savefig("ivim_voxel_plot.png")
plt.close()

Plot of the signal from one voxel.

Now we can map the perfusion and diffusion maps for the slice. We will plot a heatmap showing the values using a colormap. It will be useful to define a plotting function for the heatmap here since we will use it to plot for all the IVIM parameters. We will need to specify the lower and upper limits for our data. For example, the perfusion fractions should be in the range (0,1). Similarly, the diffusion and pseudo-diffusion constants are much smaller than 1. We pass an argument called variable to out plotting function which gives the label for the plot.

def plot_map(raw_data, variable, limits, filename):
    lower, upper = limits
    plt.title('Map for {}'.format(variable))
    plt.imshow(raw_data.T, origin='lower', clim=(lower, upper),
               cmap="gray", interpolation='nearest')
    plt.colorbar()
    plt.savefig(filename)
    plt.close()

Let us get the various plots so that we can visualize them in one page

plot_map(ivimfit.S0_predicted, "Predicted S0", (0, 10000), "predicted_S0.png")
plot_map(data_slice[:, :, 0], "Measured S0", (0, 10000), "measured_S0.png")
plot_map(ivimfit.perfusion_fraction, "f", (0, 1), "perfusion_fraction.png")
plot_map(ivimfit.D_star, "D*", (0, 0.01), "perfusion_coeff.png")
plot_map(ivimfit.D, "D", (0, 0.001), "diffusion_coeff.png")

Heatmap of S0 predicted from the fit

Heatmap of measured signal at bvalue = 0.

Heatmap of perfusion fraction values predicted from the fit

Heatmap of perfusion coefficients predicted from the fit.

Heatmap of diffusion coefficients predicted from the fit

References:

[Stejskal65]

Stejskal, E. O.; Tanner, J. E. (1 January 1965). “Spin Diffusion Measurements: Spin Echoes in the Presence of a Time-Dependent Field Gradient”. The Journal of Chemical Physics 42 (1): 288. Bibcode: 1965JChPh..42..288S. doi:10.1063/1.1695690.

[LeBihan84]

(1, 2) Le Bihan, Denis, et al. “Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging.” Radiology 168.2 (1988): 497-505.

Example source code

You can download the full source code of this example. This same script is also included in the dipy source distribution under the doc/examples/ directory.

© Copyright 2008-2016, dipy developers . Created using Sphinx 1.3.6.

Model fitting in Python - A basic tutorial on least square fitting, unit testing and plotting Magnetic Resonance images

Fitting a model using scipy.optimize.leastsq¶

I am halfway through my Google Summer of Code project with Dipy under the Python Software Foundation and I have published a few short posts about the project before but in this post I am going to walk through the entire project from the start. This can also be a good tutorial if you want to learn about curve fitting in python, unit testing and getting started with Dipy, a diffusion imaging library in Python.
The first part of this post will cover general curve fitting with leastsq, plotting and writing unit tests. The second part will focus on Dipy and the specific model I am working on - The Intravoxel Incoherent motion model.

Development branch for the IVIM model in Github : (https://github.com/nipy/dipy/pull/1058)¶

Curve fitting with leastsq¶

Mathematical models proposed to fit observed data can be characterized by several parameters. The simplest model is a straight line characterized by the slope (m) and intercept (C).

y = mx + C¶

Given a set of data points, we can determine the parameters m and C by using least squares regression which minimizes the error between the actual data points and those predicted by the model function. To demonstrate how to use leastsq for fitting let us generate some data points by importing the library numpy, defining a function and making a scatter plot.

In [140]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
# This is a magic function that renders the plots in the Ipython notebook itself
% matplotlib inline 

Let us define a model function which is an exponential decay curve and depends on two parameters.

y = A exp(-x * B)¶

Documentation is important and we should write what each function does and specify the input and output parameters for each function.

In [141]:

# Define a linear model function
def func(params, x):
    """
    Model function. The function is defined by the following equation.
                    y = A exp(-x * B)
    Parameters
    ----------
    params : array (2,)
        An array with two elements, A and B
        
    x      : array (N,)
        The independent variable
    
    Returns
    -------
    y      : array (N,)
        The values of the dependent variable
    """
    A, B = params[0], params[1]
    return A*np.exp(-x*B)

Use numpy's arange to generate an array of "x" and apply the function to get a set of data points. We can add some noise by the use of the random number generator in numpy.

In [143]:

x = np.arange(0,400,10)
N = x.shape[0]
params = 1., 0.01
Y = func(params, x)
noise = np.random.rand(N)

y = Y + noise/5.

Let us make a scatter plot of our data

In [144]:

plt.scatter(x, y, color='blue')
plt.title("Plot")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

Now that we have a set of data points, we can go ahead with the fitting using leastsq. Leastsq is a wrapped around the FORTRAN package MINPACK’s lmdif and lmder algorithms. It uses the Levenberg Marquardt algorithm for finding the best fit parameters.
Leastsq requires an error function which gives the difference between the observed target data (y) and a (non-linear) function of the parameters f(x, params). You can also specify initial guesses for the parameters using the variable x0. For more reference check out the documentation at http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

In [91]:

def error_fun(params, x, y):
    """
    Error function for model function.
    
    Parameters
    ----------
    params : array (2,)
        Parameters for the function
        
    x      : array (N,)
        The independent variable
    
    y      : array (N,)
        The observed variable
        
    Returns
    -------
    residual : array (N,)
        The difference between observed data and predicted data from the model function

    """
    residual = func(params, x) - y
    return residual

Now, we are ready to use leastsq for getting the fit. The call to leastsq returns an array whose first element is the set of parameters obtained and the other elements give more information about the fitting.

In [92]:

fit = leastsq(error_fun, x0=[0.5, 0.05], args=(x,y))
fitted_params = fit[0]
print ("Parameters after leastsq fit are :", fitted_params)

Parameters after leastsq fit are : [ 1.03770447  0.00783634]

Let us compare our fitting and plot the results.

In [93]:

y_predicted = func(fitted_params, x)
plt.scatter(y_predicted, x, color="red", label="Predicted signal")
plt.scatter(y, x, color="blue", label="Actual data")
plt.title("Plot results")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.plot()

Out[93]:

[]

Congratualations ! You have done your first model fitting with Python. You can now look at other fitting routines such as scipy.optimize.minimize, use different types of model functions and try to write code for more complicated functions and data. But at the heart of it, its as simple as this basic example.

Writing unit tests¶

Unit testing is a software testing method by which individual units of source code, sets of one or more modules together are tested to determine whether they are fit for use. Unit tests are usually very simple to write and involve very elementary tests on the expected output of your code. In the current context, it can be as simple as checking if your fit gives back the same parameters you used to generate the data. We will be using numpy's unit testing suite for our test.

In [145]:

from numpy.testing import (assert_array_equal, assert_array_almost_equal)

def test_fit():
    """
    Test the implementation of leastsq fitting
    """
    x = np.arange(0,400,10)
    N = x.shape[0]
    params = 1., 0.01
    y = func(params, x)
    
    fit = leastsq(error_fun, x0=[0.5, 0.05], args=(x,y))
    fitted_params = fit[0]
    predicted_signal = func(fitted_params, x)
    
    # Check if fit returns parameters for noiseless simulated data
    assert_array_equal(params, fitted_params)
    assert_array_almost_equal(y, predicted_signal)

test_fit()

# Usually, the package nosetests is used for unit testing. You can run the command nosetest test_file.py to run
# all the tests for a module.

So your tests pass for a noiseless fit and you get back your original parameters. But what if there is noise in the data ? Let us write a more extreme version of the test and try to see if it fails.

In [146]:

def test_fit_with_noise():
    """
    Test the fitting in presence of noise
    """
    x = np.arange(0,400,10)
    N = x.shape[0]
    params = 1., 0.01
    Y = func(params, x)
    
    noise = np.random.rand(N)

    y = Y + noise/5.
    
    fit = leastsq(error_fun, x0=[0.5, 0.05], args=(x,y))
    fitted_params = fit[0]
    predicted_signal = func(fitted_params, x)
    
    # Check if fit returns parameters for noiseless simulated data
    assert_array_equal(params, fitted_params)
    assert_array_almost_equal(y, predicted_signal)

test_fit_with_noise()

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-146-dd7240a6a422> in <module>()
     20     assert_array_almost_equal(y, predicted_signal)
     21 
---> 22 test_fit_with_noise()

<ipython-input-146-dd7240a6a422> in test_fit_with_noise()
     17 
     18     # Check if fit returns parameters for noiseless simulated data
---> 19     assert_array_equal(params, fitted_params)
     20     assert_array_almost_equal(y, predicted_signal)
     21 

/usr/local/lib/python3.4/dist-packages/numpy/testing/utils.py in assert_array_equal(x, y, err_msg, verbose)
    737     """
    738     assert_array_compare(operator.__eq__, x, y, err_msg=err_msg,
--> 739                          verbose=verbose, header='Arrays are not equal')
    740 
    741 def assert_array_almost_equal(x, y, decimal=6, err_msg='', verbose=True):

/usr/local/lib/python3.4/dist-packages/numpy/testing/utils.py in assert_array_compare(comparison, x, y, err_msg, verbose, header, precision)
    663                                 names=('x', 'y'), precision=precision)
    664             if not cond :
--> 665                 raise AssertionError(msg)
    666     except ValueError as e:
    667         import traceback

AssertionError: 
Arrays are not equal

(mismatch 100.0%)
 x: array([ 1.  ,  0.01])
 y: array([ 1.079017,  0.007812])

Thus, unit tests can make sure that your code is working properly and it behaves as it is expected to behave.¶

===================================================================================================================

The Intravoxel incoherent motion model¶

The Intravoxel incoherent motion model assumes that biological tissue includes a volume fraction of water flowing in perfused capillaries, with a perfusion coefficient D* and a fraction (1-f) of static (diffusion only), intra and extracellular water, with a diffusion coefficient D. Magnetic resonance imaging can be used to reconstruct a model of the brain (or any other tissue) using these parameters. Dipy is a python library for the analysis of diffusion MRI data. Let us get started by importing dipy, loading an image dataset and visualizing a slice of the brain. Dipy has "data fetchers" which can be used to download different datasets and read an image and load other attributes. We will use the Stanford HARDI dataset in this example.
If you want an IVIM dataset you have to pull the ivim_dev branch from dipy to get the fetcher read_ivim and replace read_stanford_hardi with read_ivim. The dataset is provided by Eric Peterson. For more details refer to https://figshare.com/articles/IVIM_dataset/3395704

In [147]:

from dipy.data import read_ivim
img, gtab = read_ivim()

Dataset is already in place. If you want to fetch it again please first remove the folder /home/shahnawaz/.dipy/ivim 

The function read_stanford_hardi downloads the dataset, reads it and returns two objects: an image and a gradient table. The gradient table has b-values and b-vectors which can be accessed by calling gtab.bvals and gtab.bvecs. The b-value is a factor that reflects the strength and timing of the gradients used to generate diffusion-weighted images. The higher the b-value, the stronger the diffusion effects.

In [148]:

print("b-values are :\n", gtab.bvals)

b-values are :
 [    0.    10.    20.    30.    40.    60.    80.   100.   120.   140.
   160.   180.   200.   300.   400.   500.   600.   700.   800.   900.
  1000.]

"img" is a NIFTI image which can be explored using the nibabel library.

In [149]:

print(img)

<class 'nibabel.nifti1.Nifti1Image'>
data shape (256, 256, 54, 21)
affine: 
[[  -0.9375       -0.            0.          120.03099823]
 [  -0.            0.9375       -0.          -91.56149292]
 [   0.            0.            2.5         -75.15000153]
 [   0.            0.            0.            1.        ]]
metadata:
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr      : 348
data_type       : b''
db_name         : b''
extents         : 0
session_error   : 0
regular         : b'r'
dim_info        : 0
dim             : [  4 256 256  54  21   1   1   1]
intent_p1       : 0.0
intent_p2       : 0.0
intent_p3       : 0.0
intent_code     : none
datatype        : float32
bitpix          : 32
slice_start     : 0
pixdim          : [-1.      0.9375  0.9375  2.5     1.      0.      0.      0.    ]
vox_offset      : 0.0
scl_slope       : nan
scl_inter       : nan
slice_end       : 0
slice_code      : unknown
xyzt_units      : 2
cal_max         : 0.0
cal_min         : 0.0
slice_duration  : 0.0
toffset         : 0.0
glmax           : 0
glmin           : 0
descrip         : b''
aux_file        : b''
qform_code      : aligned
sform_code      : scanner
quatern_b       : -0.0
quatern_c       : 1.0
quatern_d       : 0.0
qoffset_x       : 120.03099822998047
qoffset_y       : -91.56149291992188
qoffset_z       : -75.1500015258789
srow_x          : [  -0.9375       -0.            0.          120.03099823]
srow_y          : [ -0.           0.9375      -0.         -91.56149292]
srow_z          : [  0.           0.           2.5        -75.15000153]
intent_name     : b''
magic           : b'n+1'

For our use, we just need the data which can be read by calling the get_data() method of the img object.

In [150]:

data = img.get_data()

The data has 54 slices, with 256-by-256 voxels in each slice. The fourth dimension corresponds to the b-values in the gtab.

In [151]:

print (data.shape)

(256, 256, 54, 21)

Let us take a slice of the data at a particular b value and height (z) and visualize it

In [152]:

plt.imshow(data[:,:,31, 0], cmap='gray', origin ="lower")
plt.show()

Thus, for each b value we have a 3d matrix which gives us the signal value in a 256x256x54 box. To see how the signal in one of the points in this volume (one voxel) varies with bvalue, let us consider a random point from the data slice. Say we want the data in a 20x20 square at the height z=31 for all the bvalues, we specify the range of indices from our array using ":" as :

In [153]:

x1, x2 = 100, 120
y1, y2 = 150, 170
z = 31

data_slice = data[x1:x2, y1:y2, z, :]
print(data_slice.shape)

(20, 20, 21)

The data slice is now a 20x20 portion with the third dimension specifying the bvalues. We can visulaize this slice at a particular bvalue and extract the x and y data for the voxel located at (0,0) and for all bvals (:).

In [154]:

plt.imshow(data_slice[0], origin="lower", cmap="gray")
plt.show()

In [155]:

xdata = gtab.bvals
ydata = data_slice[0,0,:]

plt.scatter(xdata, ydata)
plt.xlabel("b-values")
plt.ylabel("Signal")
plt.plot()

Out[155]:

[]

The fitting method remains the same now. Let us try to fit the IVIM model function on this data and get IVIM parameters for this voxel. The IVIM model function is a bi exponential curve and has 4 parameters S0, f, D_star and D and is defined as follows :

S(b) = S0(f e^(- b D_star) + (1 - f) e^(-b D))¶

In [156]:

def ivim_function(params, bvals):
    """The Intravoxel incoherent motion (IVIM) model function.

        S(b) = S_0[f*e^{(-b*D\*)} + (1-f)e^{(-b*D)}]

        S_0, f, D\* and D are the IVIM parameters.

    Parameters
    ----------
        params : array
                parameters S0, f, D_star and D of the model

        bvals : array
                bvalues

    References
    ----------
    .. [1] Le Bihan, Denis, et al. "Separation of diffusion
               and perfusion in intravoxel incoherent motion MR
               imaging." Radiology 168.2 (1988): 497-505.
    .. [2] Federau, Christian, et al. "Quantitative measurement
               of brain perfusion with intravoxel incoherent motion
               MR imaging." Radiology 265.3 (2012): 874-881.
    """
    S0, f, D_star, D = params
    S = S0 * (f * np.exp(-bvals * D_star) + (1 - f) * np.exp(-bvals * D))
    return S


def _ivim_error(params, bvals, signal):
    """Error function to be used in fitting the IVIM model
    """
    return (signal - ivim_function(params, bvals))

Once we have our error function, we can use the same procedure to fit an IVIM model by using the leastsq function from scipy.optimize

In [157]:

x0 = [1000., 0.1, 0.01, 0.001]
fit = leastsq(_ivim_error,
              x0,
              args=(xdata, ydata),)

estimated_parameters = fit[0]
print('Parameters estimated using leastsq :', estimated_parameters)

Parameters estimated using leastsq : [  4.44684213e+03   1.41584992e-02  -1.88101978e-03   1.10032569e-03]

In [158]:

predicted_signal = ivim_function(estimated_parameters, xdata)

plt.plot(xdata, predicted_signal, label="Estimated signal")
plt.scatter(xdata, ydata, color="red", label="Observed data")
plt.xlabel("b-values")
plt.ylabel("Signal")
plt.legend()
plt.show()

This forms the core of the project and once you get the hang of this, you can easily implement any new model for fitting data.¶