IntraVoxel Incoherent Motion (IVIM) and S0 optimization

.travis.yml

@@ -55,7 +55,7 @@ install:
 - python setup.py install --record installed_files.txt
 
 script:
-- flake8 --ignore N802,N806,E731,F401,W504 `find . -name \*.py | grep -v setup.py | grep -v /doc/`
+- flake8 --ignore N802,N806,E731,F401,W504,W605 `find . -name \*.py | grep -v setup.py | grep -v /doc/`
 - nosetests --with-coverage -v
   #- travis-sphinx build -s doc
   #- for a in examples/*ipynb; do 

README.md

@@ -34,7 +34,7 @@ Dmipy allows the user to do Microstructure Imaging research at the highest level
 - clone repository
 - python setup.py install
 
-See solutions to [common issues](https://github.com/AthenaEPI/mipy/blob/master/common_issues.md)
+See solutions to [common issues](https://github.com/AthenaEPI/dmipy/blob/master/common_issues.md)
 ## Dependencies
 Recommended to use Anaconda Python distribution.
 - numpy >= 1.13
@@ -47,47 +47,48 @@ Recommended to use Anaconda Python distribution.
 
 ## Getting Started
 To get a feeling for how to use Dmipy, we provide a few tutorial notebooks:
-- [Setting up an acquisition scheme](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/tutorial_setting_up_acquisition_scheme.ipynb)
-- [Simulating and fitting data using a simple Stick model](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/tutorial_simulating_and_fitting_using_a_simple_model.ipynb)
-- [Combining biophysical models into a Microstructure model](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/tutorial_combining_biophysical_models_into_microstructure_model.ipynb)
+- [Setting up an acquisition scheme](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/tutorial_setting_up_acquisition_scheme.ipynb)
+- [Simulating and fitting data using a simple Stick model](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/tutorial_simulating_and_fitting_using_a_simple_model.ipynb)
+- [Combining biophysical models into a Microstructure model](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/tutorial_combining_biophysical_models_into_microstructure_model.ipynb)
 - [Creating a dispersed and/or distributed axon bundle representation](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/tutorial_distributed_model_representations.ipynb)
-- [Imposing parameter links and constraints](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/tutorial_imposing_parameter_links.ipynb)
-- [Parameter Cascading: Using a simple model to initialize a complex one](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/tutorial_parameter_cascading_and_simulating_nd_datasets.ipynb)
+- [Imposing parameter links and constraints](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/tutorial_imposing_parameter_links.ipynb)
+- [Parameter Cascading: Using a simple model to initialize a complex one](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/tutorial_parameter_cascading_and_simulating_nd_datasets.ipynb)
 
 ## Explanations and Illustrations of Dmipy Contents
 ### Biophysical Models and Distributions
-- [Cylinder Models (Axons, e.g. [Assaf et al. 2004])](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_cylinder_models.ipynb)
-- [Sphere Models (Tumor cells, e.g. [Panagiotaki et al. 2014])](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_sphere_models.ipynb)
-- [Parameter Distribution Models (Axon Diameter Distribution, e.g. [Assaf et al. 2008])](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_diameter_distributions.ipynb)
-- [Gaussian Models (Extra-axonal, e.g. [Behrens et al. 2003])](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_gaussian_models.ipynb)
+- [Cylinder Models (Axons, e.g. [Assaf et al. 2004])](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_cylinder_models.ipynb)
+- [Sphere Models (Tumor cells, e.g. [Panagiotaki et al. 2014])](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_sphere_models.ipynb)
+- [Parameter Distribution Models (Axon Diameter Distribution, e.g. [Assaf et al. 2008])](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_diameter_distributions.ipynb)
+- [Gaussian Models (Extra-axonal, e.g. [Behrens et al. 2003])](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_gaussian_models.ipynb)
 - Tissue Response Function Models and Estimation (WM/GM/CSF, e.g. [Jeurissen et al. 2014])
-- [Spherical Distribution Models (Axon Dispersion, e.g. [Kaden et al. 2007])](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_watson_bingham.ipynb)
-- [Spherical Mean of any Compartment Model](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_spherical_mean_models.ipynb)
+- [Spherical Distribution Models (Axon Dispersion, e.g. [Kaden et al. 2007])](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_watson_bingham.ipynb)
+- [Spherical Mean of any Compartment Model](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_spherical_mean_models.ipynb)
 ### Global Multi-Compartment Optimizers
 - [Brute Force to Gradient Descent (Brute2Fine)](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_brute_force_optimization.ipynb)
 - [Stochastic (MIX) [Farooq et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_stochastic_mix_optimization.ipynb)
 ### Constrained Spherical Deconvolution Optimizers
 - [Generalized Multi-Shell Multi-Compartment CSD [Tournier et al. 2007, Jeurissen et al. 2014]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_generalized_csd_optimizer.ipynb)
 ## Dmipy implementations of Microstructure Models in Literature
-Dmipy uses HCP data to illustrate microstructure model examples. To reproduce these examples, dmipy provides a direct way to download HCP data (using your own AWS credentials) in the [HCP tutorial](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/tutorial_human_connectome_project_aws.ipynb).
+Dmipy uses HCP data to illustrate microstructure model examples. To reproduce these examples, dmipy provides a direct way to download HCP data (using your own AWS credentials) in the [HCP tutorial](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/tutorial_human_connectome_project_aws.ipynb).
 ### Single Bundle Models
-- [Ball and Stick [Behrens et al. 2003]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_ball_and_stick.ipynb)
-- [Ball and Racket [Sotiropoulos et al. 2012]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_ball_and_racket.ipynb)
-- [NODDI-Watson [Zhang et al. 2012]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_noddi_watson.ipynb)
-- [NODDI-Bingham [Tariq et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_noddi_bingham.ipynb)
-- [AxCaliber [Assaf et al. 2008]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_axcaliber.ipynb)
+- [Ball and Stick [Behrens et al. 2003]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_ball_and_stick.ipynb)
+- [Ball and Racket [Sotiropoulos et al. 2012]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_ball_and_racket.ipynb)
+- [NODDI-Watson [Zhang et al. 2012]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_noddi_watson.ipynb)
+- [NODDI-Bingham [Tariq et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_noddi_bingham.ipynb)
+- [AxCaliber [Assaf et al. 2008]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_axcaliber.ipynb)
 - [AxCaliber with Extra-Axonal Diffusion Time-Dependence [De Santis et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_axcaliber_temporal_zeppeline.ipynb)
 
 ### Crossing Bundle Models
-- [Ball and Sticks [Behrens et al. 2003]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_ball_and_sticks.ipynb)
-- [NODDI in Crossings (NODDIx) [Farooq et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_mix_microstructure_imaging_in_crossings.ipynb)
+- [Ball and Sticks [Behrens et al. 2003]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_ball_and_sticks.ipynb)
+- [NODDI in Crossings (NODDIx) [Farooq et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_mix_microstructure_imaging_in_crossings.ipynb)
 
 ### Tumor Models
-- [VERDICT [Panagiotaki et al. 2014]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_verdict.ipynb)
+- [IVIM [Le Bihan et al. 1988, Gurney-Champion et al. 2018]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_ivim.ipynb)
+- [VERDICT [Panagiotaki et al. 2014]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_verdict.ipynb)
 ### Spherical Mean-Based Models
 Any Spherical Mean model can also estimate parametric FODs.
 - [Spherical Mean Technique [Kaden et al. 2015]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_spherical_mean_technique.ipynb)
-- [Multi-Compartment Microscopic Diffusion Imaging [Kaden et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_multi_compartment_spherical_mean_technique.ipynb)
+- [Multi-Compartment Microscopic Diffusion Imaging [Kaden et al. 2016]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_multi_compartment_spherical_mean_technique.ipynb)
 ### Spherical Deconvolution-Based Models
 Constrained spherical deconvolution (CSD) models are primarily used for Fiber Orientation Distribution (FOD) estimation. Multi-Tissue CSD models improve FOD estimation by separating WM/GM/CSF signal contributions using multiple tissue response functions.
 - [Multi-Shell Multi-Compartment CSD [model-based version of Jeurissen et al. 2014]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_multi_compartment_constrained_spherical_deconvolution.ipynb)

dmipy/core/acquisition_scheme.py

@@ -165,9 +165,11 @@ def __init__(self, bvalues, gradient_directions, qvalues,
                 self.shell_sh_orders[shell_index], theta_, phi_)[0]
         # warning in case there are no b0 measurements
         if sum(self.b0_mask) == 0:
-            msg = "No b0 measurements were detected. Check if the b0_threshold"
-            msg += " option is high enough, or if there is a mistake in the "
-            msg += "acquisition design."
+            msg = "No b0 measurements were detected. If this is the case, you "
+            msg += "need to estimate it from the data by setting "
+            msg += "optimize_S0=True when fitting. If not, check if the "
+            msg += "b0_threshold option is high enough, or if there is a "
+            msg += "mistake in the acquisition design."
             warn(msg)
 
         self.spherical_mean_scheme = SphericalMeanAcquisitionScheme(
@@ -362,8 +364,12 @@ def return_pruned_acquisition_scheme(self, shell_indices, data=None):
 
         bvals = self.bvalues[mask]
         gradient_directions = self.gradient_directions[mask]
-        delta = self.delta[mask]
-        Delta = self.Delta[mask]
+        delta = None
+        Delta = None
+        if self.delta is not None:
+            delta = self.delta[mask]
+        if self.Delta is not None:
+            Delta = self.Delta[mask]
         if self.TE is not None:
             TE = self.TE[mask]
         else:

dmipy/core/fitted_modeling_framework.py

@@ -35,11 +35,11 @@ class FittedMultiCompartmentModel:
         fitted model parameters array.
     """
 
-    def __init__(self, model, S0, mask, fitted_parameters_vector):
+    def __init__(self, model, mask, fitted_parameters_vector):
         self.model = model
-        self.S0 = S0
         self.mask = mask
         self.fitted_parameters_vector = fitted_parameters_vector
+        self.S0 = self.fitted_and_linked_parameters['S0']
 
     @property
     def fitted_parameters(self):
@@ -179,9 +179,9 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
 
         dataset_shape = self.fitted_parameters_vector.shape[:-1]
         if S0 is None:
-            S0 = self.S0
+            S0_mult = np.ones_like(self.S0)
         elif isinstance(S0, float):
-            S0 = np.ones(dataset_shape) * S0
+            S0_mult = S0 / self.S0
         if mask is None:
             mask = self.mask
 
@@ -193,7 +193,8 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
             parameters = self.model.parameter_vector_to_parameters(
                 self.fitted_parameters_vector[pos])
             predicted_signal[pos] = self.model(
-                acquisition_scheme, **parameters) * S0[pos]
+                acquisition_scheme, **parameters) * np.dot(
+                    self.model.S0_mapping, S0_mult[pos])
         return predicted_signal
 
     def R2_coefficient_of_determination(self, data):
@@ -241,11 +242,11 @@ class FittedMultiCompartmentSphericalMeanModel:
         fitted model parameters array.
     """
 
-    def __init__(self, model, S0, mask, fitted_parameters_vector):
+    def __init__(self, model, mask, fitted_parameters_vector):
         self.model = model
-        self.S0 = S0
         self.mask = mask
         self.fitted_parameters_vector = fitted_parameters_vector
+        self.S0 = self.fitted_and_linked_parameters['S0']
 
     @property
     def fitted_parameters(self):
@@ -292,9 +293,9 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
 
         dataset_shape = self.fitted_parameters_vector.shape[:-1]
         if S0 is None:
-            S0 = self.S0
+            S0_mult = np.ones_like(self.S0)
         elif isinstance(S0, float):
-            S0 = np.ones(dataset_shape) * S0
+            S0_mult = S0 / self.S0
         if mask is None:
             mask = self.mask
 
@@ -306,7 +307,8 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
             parameters = self.model.parameter_vector_to_parameters(
                 self.fitted_parameters_vector[pos])
             predicted_signal[pos] = self.model(
-                acquisition_scheme, **parameters) * S0[pos]
+                acquisition_scheme, **parameters) * np.dot(
+                    self.model.S0_mapping_sm, S0_mult[pos])
         return predicted_signal
 
     def R2_coefficient_of_determination(self, data):
@@ -447,6 +449,8 @@ def return_parametric_fod_model(
                                                      parameters[i])
 
         for parameters in parameter_pairs:
+            if smt_par_name == 'S0':  # temporary solution
+                continue
             smt_parameter_name = parameters[0]
             mc_parameter_name = parameters[1]
             mc_bundles_model.set_fixed_parameter(
@@ -526,6 +530,8 @@ def return_spherical_harmonics_fod_model(self, sh_order=8):
             del sh_model.parameter_optimization_flags[param_to_delete]
 
         for smt_par_name in self.model.parameter_names:
+            if smt_par_name == 'S0':  # temporary solution
+                continue
             sh_model.set_fixed_parameter(
                 smt_par_name, self.fitted_parameters[smt_par_name])
         return sh_model