Tournier13 White Matter Response Function Estimation (#19)

* added tournier13 algorithm and peak finding fod
* added solver error catch to cvxpy fod solver
* optional inclusion of S0-response for MT-CSD fitting using 'fit_s0_reponse'
* update mtcsd example with and without S0-response signal fit illustration
* added and update tissue wm and 3 tissue response test

README.md

@@ -60,7 +60,7 @@ To get a feeling for how to use Dmipy, we provide a few tutorial notebooks:
 - [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)
-- Tissue Response Models (WM/GM/CSF, e.g. [Jeurissen et al. 2014])
+- 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)
 ### Global Multi-Compartment Optimizers

dmipy/core/fitted_modeling_framework.py

@@ -550,6 +550,12 @@ def __init__(self, model, S0, mask, fitted_parameters_vector):
         self.S0 = S0
         self.mask = mask
         self.fitted_parameters_vector = fitted_parameters_vector
+        self.S0_responses = model.S0_responses.copy()
+        self.max_S0_response = model.max_S0_response
+        self.fit_S0_response = model.fit_S0_response
+        self._fit_parameters = {
+            'fit_S0_response': self.fit_S0_response,
+            'S0_responses': self.S0_responses}
 
     @property
     def fitted_parameters(self):
@@ -618,6 +624,53 @@ def fod_sh(self):
         """
         return self.fitted_parameters['sh_coeff']
 
+    def peaks_directions(self, sphere, max_peaks=5,
+                         relative_peak_threshold=0.5,
+                         min_separation_angle=25):
+        """
+        Returns peak directions for estimated FODs. Uses dipy's peak_directions
+        function to get the local maximum on a sphere's tesselation.
+
+        Parameters
+        ----------
+        sphere : Sphere
+            The Sphere providing discrete directions for evaluation.
+        max_peaks : integer
+            The maximum number of peaks that is returned per fod.
+        relative_peak_threshold : float in [0., 1.]
+            Only peaks greater than min + relative_peak_threshold * scale are
+            kept, where min = max(0, odf.min()) and scale = odf.max() - min.
+        min_separation_angle : float in [0, 90]
+            The minimum distance between directions. If two peaks are too close
+            only the larger of the two is returned.
+
+        Returns
+        -------
+        directions : (Ndata, Npeaks, 3,) ndarray
+            N vertices for sphere, one for each peak
+        values : (Ndata, Npeaks,) ndarray
+            peak values
+        indices : (Ndata, Npeaks,) ndarray
+            peak indices of the directions on the sphere
+        """
+        peaks_shape = np.r_[self.mask.shape, max_peaks, 3]
+        peaks = np.zeros(peaks_shape)
+        values = np.zeros(peaks_shape[:-1])
+        indices = np.zeros(peaks_shape[:-1])
+
+        fods = self.fod(sphere.vertices)
+        mask_pos = np.where(self.mask)
+        for pos in zip(*mask_pos):
+            fod = fods[pos]
+            peaks_, values_, indices_ = peak_directions(
+                fod, sphere, relative_peak_threshold, min_separation_angle)
+            # if less peaks than max_peaks are found, only take those.
+            Npeaks = np.min([len(indices_), max_peaks])
+            peaks[pos, :Npeaks] = peaks_[:Npeaks]
+            values[pos, :Npeaks] = values_[:Npeaks]
+            indices[pos, :Npeaks] = indices_[:Npeaks]
+        return peaks, values, indices
+
     def anisotropy_index(self):
         """
         Estimates anisotropy index of spherical harmonics [1]_.
@@ -686,9 +739,15 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
             acquisition_scheme = self.model.scheme
         dataset_shape = self.fitted_parameters_vector.shape[:-1]
         if S0 is None:
-            S0 = self.S0
+            if self.fit_S0_response:
+                S0 = np.ones(dataset_shape) * self.max_S0_response
+            else:
+                S0 = self.S0
         elif isinstance(S0, float):
-            S0 = np.ones(dataset_shape) * S0
+            if self.fit_S0_response:
+                S0 = S0 * self.max_S0_response / self.S0
+            else:
+                S0 = np.ones(dataset_shape) * S0
         if mask is None:
             mask = self.mask
 
@@ -699,13 +758,18 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
         for pos in zip(*mask_pos):
             parameters = self.model.parameter_vector_to_parameters(
                 self.fitted_parameters_vector[pos])
-            predicted_signal[pos] = self.model(
-                acquisition_scheme, **parameters) * S0[pos]
+            predicted_signal[pos] = (
+                self.model(acquisition_scheme,
+                           **dict(parameters,
+                                  **self._fit_parameters)) * S0[pos])
         return predicted_signal
 
     def R2_coefficient_of_determination(self, data):
         "Calculates the R-squared of the model fit."
-        data_ = data / self.S0[..., None]
+        if self.model.scheme.TE is None:
+            data_ = data / self.S0[..., None]
+        else:
+            data_ = data / self.S0
 
         y_hat = self.predict(S0=1.)
         y_bar = np.mean(data_, axis=-1)

dmipy/core/modeling_framework.py

@@ -58,6 +58,9 @@
 
 class ModelProperties:
     "Contains various properties for CompartmentModels."
+
+    S0_response = 1.
+
     @property
     def parameter_ranges(self):
         """Returns the optimization ranges of the model parameters.
@@ -1642,7 +1645,8 @@ def _check_that_one_anisotropic_kernel_is_present(self):
 
     def fit(self, acquisition_scheme, data, mask=None, solver='csd',
             lambda_lb=1e-5, unity_constraint='kernel_dependent',
-            use_parallel_processing=have_pathos, number_of_processors=None):
+            fit_S0_response=False, use_parallel_processing=have_pathos,
+            number_of_processors=None):
         """ The main data fitting function of a
         MultiCompartmentSphericalHarmonicsModel.
 
@@ -1687,6 +1691,12 @@ def fit(self, acquisition_scheme, data, mask=None, solver='csd',
             enforce unity if the kernel is voxel-varying or when volume
             fractions are estimated. Otherwise unity_constraint is set to
             False.
+        fit_S0_response: bool,
+            whether or not to fit the raw signal or signal attenuation.
+            default: False, the signal is automatically divided by S0-value.
+            if True, the raw signal is fitted and the S0 intensities of the
+            biophysical models are used in the signal generation. This is
+            useful when using tissue_response_models for example.
         use_parallel_processing : bool,
             Whether or not to use parallel processing using pathos.
         number_of_processors : integer,
@@ -1730,8 +1740,14 @@ def fit(self, acquisition_scheme, data, mask=None, solver='csd',
         else:
             self.unity_constraint = unity_constraint
 
-        # S0_responses = np.array([model.S0_response for model in self.models])
-        # self.relative_responses = S0_responses / np.max(S0_responses)
+        self.fit_S0_response = fit_S0_response
+        if self.fit_S0_response:
+            S0_responses = np.r_[[model.S0_response for model in self.models]]
+            self.max_S0_response = S0_responses.max()
+            self.S0_responses = S0_responses / self.max_S0_response
+        else:
+            self.S0_responses = np.ones(len(self.models), dtype=float)
+            self.max_S0_response = 1.
 
         # estimate S0
         self.scheme = acquisition_scheme
@@ -1824,9 +1840,12 @@ def fit(self, acquisition_scheme, data, mask=None, solver='csd',
 
         start = time()
         for idx, pos in enumerate(zip(*mask_pos)):
-            voxel_E = data_[pos] / S0[pos]
+            if fit_S0_response:
+                data_to_fit = data_[pos] / self.max_S0_response
+            else:
+                data_to_fit = data_[pos] / S0[pos]
             voxel_x0_vector = x0_[pos]
-            fit_args = (voxel_E, voxel_x0_vector)
+            fit_args = (data_to_fit, voxel_x0_vector)
 
             if use_parallel_processing:
                 fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args)
@@ -1916,6 +1935,9 @@ def __call__(self, acquisition_scheme, **kwargs):
         kwargs = self.add_linked_parameters_to_parameters(
             kwargs
         )
+        self.S0_responses = kwargs.get('S0_responses', self.S0_responses)
+        self.fit_S0_response = kwargs.get(
+            'fit_S0_response', self.fit_S0_response)
 
         A = self._construct_convolution_kernel(
             self.parameters_to_parameter_vector(**kwargs))
@@ -1976,8 +1998,9 @@ def _construct_convolution_kernel(self, x0_vector):
             else:
                 partial_volumes = [1.]
             kernel = 0.
-            for model, partial_volume in zip(self.models,
-                                             partial_volumes):
+            for model, partial_volume, S0 in zip(self.models,
+                                                 partial_volumes,
+                                                 self.S0_responses):
                 parameters = {}
                 for parameter in model.parameter_ranges:
                     parameter_name = self._inverted_parameter_map[
@@ -1989,11 +2012,11 @@ def _construct_convolution_kernel(self, x0_vector):
                 model_rh = (
                     model.rotational_harmonics_representation(
                         self.scheme, **parameters))
-                kernel += partial_volume * construct_model_based_A_matrix(
+                kernel += S0 * partial_volume * construct_model_based_A_matrix(
                     self.scheme, model_rh, self.sh_order)
         else:
             kernel = []
-            for model in self.models:
+            for model, S0 in zip(self.models, self.S0_responses):
                 parameters = {}
                 for parameter in model.parameter_ranges:
                     parameter_name = self._inverted_parameter_map[
@@ -2006,10 +2029,10 @@ def _construct_convolution_kernel(self, x0_vector):
                     model.rotational_harmonics_representation(
                         self.scheme, **parameters))
                 if 'orientation' in model.parameter_types.values():
-                    kernel.append(construct_model_based_A_matrix(
+                    kernel.append(S0 * construct_model_based_A_matrix(
                         self.scheme, model_rh, self.sh_order))
                 else:
-                    kernel.append(construct_model_based_A_matrix(
+                    kernel.append(S0 * construct_model_based_A_matrix(
                         self.scheme, model_rh, 0))
 
             kernel = np.hstack(kernel)

dmipy/optimizers_fod/csd_cvxpy.py

@@ -163,7 +163,17 @@ def __call__(self, data, x0_vector):
             cost += (
                 self.lambda_lb * cvxpy.quad_form(sh_coef, self.R_smoothness))
         problem = cvxpy.Problem(cvxpy.Minimize(cost), constraints)
-        problem.solve()
+        try:
+            problem.solve()
+        except cvxpy.error.SolverError:
+            msg = 'cvxpy solver failed'
+            print(msg)
+            return np.zeros_like(x0_vector)
+
+        if problem.status in ["infeasible", "unbounded"]:
+            msg = 'cvxpy found {} problem'.format(problem.status)
+            print(msg)
+            return np.zeros_like(x0_vector)
 
         # return optimized fod sh coefficients
         fitted_params = self.model.parameter_vector_to_parameters(x0_vector)
@@ -172,6 +182,7 @@ def __call__(self, data, x0_vector):
         if not self.model.volume_fractions_fixed:  # if vf was estimated
             fractions_array = np.array(
                 sh_coef[self.vf_indices].value).squeeze() * 2 * np.sqrt(np.pi)
+            fractions_array /= np.sum(fractions_array)  # for small deviations
             for i, name in enumerate(self.model.partial_volume_names):
                 fitted_params[name] = fractions_array[i]
         fitted_parameter_vector = self.model.parameters_to_parameter_vector(

dmipy/tissue_response/tests/test_three_tissue_response.py

@@ -0,0 +1,21 @@
+import numpy as np
+from dmipy.tissue_response.three_tissue_response import (
+    optimal_threshold, signal_decay_metric)
+from dmipy.signal_models.gaussian_models import G1Ball
+from dmipy.data.saved_acquisition_schemes import (
+    wu_minn_hcp_acquisition_scheme)
+scheme = wu_minn_hcp_acquisition_scheme()
+
+
+def test_optimal_threshold():
+    data = np.linspace(0, 8, 101)
+    opt = optimal_threshold(data)
+    np.testing.assert_(np.round(opt) == 8 // 2)
+
+
+def test_signal_decay_metric():
+    data_iso1 = G1Ball(lambda_iso=2.5e-9)(scheme)
+    data_iso2 = G1Ball(lambda_iso=1.5e-9)(scheme)
+    data = np.array([data_iso1, data_iso2])
+    sdm = signal_decay_metric(scheme, data)
+    np.testing.assert_(sdm[0] > sdm[1])

dmipy/tissue_response/tests/test_tissue_response_csd.py

@@ -0,0 +1,37 @@
+from dmipy.core.modeling_framework import (
+    MultiCompartmentSphericalHarmonicsModel)
+from dmipy.signal_models.gaussian_models import G1Ball, G2Zeppelin
+from dmipy.signal_models.tissue_response_models import (
+    IsotropicTissueResponseModel,
+    AnisotropicTissueResponseModel)
+from dmipy.data.saved_acquisition_schemes import (
+    wu_minn_hcp_acquisition_scheme,)
+import numpy as np
+scheme = wu_minn_hcp_acquisition_scheme()
+
+
+def test_fit_S0_response(S0_iso=10., S0_aniso=1.):
+    ball = G1Ball(lambda_iso=3e-9)
+    data_iso = S0_iso * ball(scheme)
+    iso_model = IsotropicTissueResponseModel(scheme, np.atleast_2d(data_iso))
+
+    zeppelin = G2Zeppelin(
+        lambda_par=2.2e-9, lambda_perp=1e-9, mu=[np.pi / 2, np.pi / 2])
+    data_aniso = S0_aniso * zeppelin(scheme)
+    aniso_model = AnisotropicTissueResponseModel(
+        scheme, np.atleast_2d(data_aniso))
+
+    mtcsd = MultiCompartmentSphericalHarmonicsModel([iso_model, aniso_model])
+
+    data_to_fit = 0.3 * data_iso + 0.7 * data_aniso
+
+    csd_fit_no_S0 = mtcsd.fit(scheme, data_to_fit, fit_S0_response=False)
+    csd_fit_S0 = mtcsd.fit(scheme, data_to_fit, fit_S0_response=True)
+    np.testing.assert_almost_equal(
+        0.3, csd_fit_S0.fitted_parameters['partial_volume_0'], 3)
+    np.testing.assert_almost_equal(
+        0.7, csd_fit_S0.fitted_parameters['partial_volume_1'], 3)
+
+    # test iso volume fraction overestimated without S0 response
+    np.testing.assert_(
+        csd_fit_no_S0.fitted_parameters['partial_volume_0'] > 0.3)

dmipy/tissue_response/tests/test_tournier_wm_response.py

@@ -0,0 +1,64 @@
+from dmipy.signal_models.gaussian_models import G2Zeppelin
+from dmipy.distributions import distribute_models
+from dmipy.data.saved_acquisition_schemes import (
+    wu_minn_hcp_acquisition_scheme)
+from dmipy.signal_models.tissue_response_models import (
+    AnisotropicTissueResponseModel)
+from dmipy.tissue_response.white_matter_response import (
+    white_matter_response_tournier07,
+    white_matter_response_tournier13)
+import numpy as np
+from numpy.testing import assert_array_almost_equal, assert_raises
+scheme = wu_minn_hcp_acquisition_scheme()
+
+
+def _make_single_and_crossing():
+    zeppelin = G2Zeppelin(
+        lambda_par=1.7e-9, lambda_perp=1e-9, mu=[0., 0.])
+    data_aniso = zeppelin(scheme)
+    aniso_model = AnisotropicTissueResponseModel(
+        scheme, np.atleast_2d(data_aniso))
+
+    watson_mod = distribute_models.SD1WatsonDistributed(
+        [aniso_model])
+    watson_params_par = {
+        'SD1Watson_1_mu': np.array(
+            [0., 0.]),
+        'SD1Watson_1_odi': .2
+    }
+    watson_params_perp = {
+        'SD1Watson_1_mu': np.array(
+            [np.pi / 2., np.pi / 2.]),
+        'SD1Watson_1_odi': .2
+    }
+    data_watson_par = watson_mod(scheme, **watson_params_par)
+    data_watson_perp = watson_mod(scheme, **watson_params_perp)
+
+    data_cross = np.array(
+        [data_watson_par, data_watson_par + data_watson_perp])
+    return data_cross
+
+
+def test_tournier13_picks_single_peak():
+    data_cross = _make_single_and_crossing()
+
+    assert_raises(ValueError,
+                  white_matter_response_tournier13,
+                  scheme, data_cross, peak_ratio_setting='bla')
+
+    wm, _ = white_matter_response_tournier13(
+        scheme, data_cross, peak_ratio_setting='mrtrix',
+        N_candidate_voxels=1)
+    assert_array_almost_equal(
+        wm(scheme, mu=[0., 0.]),
+        data_cross[0], 4)
+
+
+def test_tournier07_picks_single_peak():
+    data_cross = _make_single_and_crossing()
+    wm, _ = white_matter_response_tournier07(
+        scheme, data_cross, peak_ratio_setting='mrtrix',
+        N_candidate_voxels=1)
+    assert_array_almost_equal(
+        wm(scheme, mu=[0., 0.]),
+        data_cross[0], 4)