Multi-Shell Free Water Elimination Using MT-CSD

README.md

@@ -94,6 +94,10 @@ Constrained spherical deconvolution (CSD) models are primarily used for Fiber Or
 - [Multi-Shell Multi-Tissue CSD with Unsupervised 3-Tissue Response Function Estimation [Jeurissen 2014, Dhollander 2016a]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_multi_tissue_constrained_spherical_deconvolution.ipynb)
 - [Single-Shell (1 shell + b0) Multi-Tissue CSD [Dhollander et al. 2016b]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_single_shell_three_tissue_csd.ipynb)
 
+### Free Water Elimination
+Removing free water signal contributions can be useful for studying WM properties on the interface between WM and CSF.
+- Multi-Shell Free Water Elimination using MSMT-CSD (applied to MAP-MRI)
+
 ## How to contribute to Dmipy
 Dmipy's design is completely modular and can easily be extended with new models, distributions or optimizers. To contribute, view our [contribution guidelines](https://github.com/AthenaEPI/dmipy/blob/master/contribution_guidelines.rst).
 ## How to cite Dmipy

dmipy/core/fitted_modeling_framework.py

@@ -773,6 +773,43 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
                                   **self._fit_parameters)) * S0[pos])
         return predicted_signal
 
+    def return_filtered_signal(self, filtered_partial_fraction_names,
+                               acquisition_scheme=None, S0=None, mask=None):
+        """
+        Returns the fitted signal but omits the signal contributions of the
+        given partial volume names.
+
+        Parameters
+        ----------
+        filtered_partial_fraction_names: list of strings,
+            partial volume names to be omitted.
+        """
+        # replace with zero but save filtered volume fractions
+        fractions_tmp = {}
+        for fraction_name in filtered_partial_fraction_names:
+            if fraction_name not in self.model.partial_volume_names:
+                msg = "{} not a valid partial volume name".format(
+                    fraction_name)
+                raise ValueError(msg)
+            fractions_tmp[fraction_name] = self.fitted_parameters[
+                fraction_name].copy()
+            self.fitted_parameters[fraction_name] *= 0.
+        self.fitted_parameters_vector = (
+            self.model.parameters_to_parameter_vector(
+                **self.fitted_parameters))
+
+        # get filtered signal
+        filtered_signal = self.predict(acquisition_scheme, S0, mask)
+
+        # return original volume fractions
+        for fraction_name in filtered_partial_fraction_names:
+            self.fitted_parameters[fraction_name] += fractions_tmp[
+                fraction_name]
+        self.fitted_parameters_vector = (
+            self.model.parameters_to_parameter_vector(
+                **self.fitted_parameters))
+        return filtered_signal
+
     def R2_coefficient_of_determination(self, data):
         "Calculates the R-squared of the model fit."
         if self.model.scheme.TE is None:

dmipy/core/modeling_framework.py

@@ -2043,8 +2043,10 @@ def _construct_convolution_kernel(self, x0_vector):
                 model_rh = (
                     model.rotational_harmonics_representation(
                         self.scheme, **parameters))
-                kernel += S0 * partial_volume * construct_model_based_A_matrix(
+                A_matrix = construct_model_based_A_matrix(
                     self.scheme, model_rh, self.sh_order)
+                A_matrix = self._divide_by_positive_dirac_rh(A_matrix)
+                kernel += S0 * partial_volume * A_matrix
         else:
             kernel = []
             for model, S0 in zip(self.models, self.S0_responses):
@@ -2060,15 +2062,46 @@ def _construct_convolution_kernel(self, x0_vector):
                     model.rotational_harmonics_representation(
                         self.scheme, **parameters))
                 if 'orientation' in model.parameter_types.values():
-                    kernel.append(S0 * construct_model_based_A_matrix(
-                        self.scheme, model_rh, self.sh_order))
+                    A_matrix = construct_model_based_A_matrix(
+                        self.scheme, model_rh, self.sh_order)
+                    A_matrix = self._divide_by_positive_dirac_rh(A_matrix)
+                    kernel.append(S0 * A_matrix)
                 else:
                     kernel.append(S0 * construct_model_based_A_matrix(
                         self.scheme, model_rh, 0))
 
             kernel = np.hstack(kernel)
         return kernel
 
+    def _divide_by_positive_dirac_rh(self, A_matrix):
+        import cvxpy
+        from dipy.reconst.shm import gen_dirac, sph_harm_ind_list
+        from dipy.data import get_sphere, HemiSphere
+        from dipy.reconst.shm import real_sym_sh_mrtrix
+
+        sphere = get_sphere('symmetric724')
+        hemisphere = HemiSphere(phi=sphere.phi, theta=sphere.theta)
+        L_positivity = real_sym_sh_mrtrix(
+            self.sh_order, hemisphere.theta, hemisphere.phi)[0]
+
+        m, ll = sph_harm_ind_list(self.sh_order)
+        dirac = gen_dirac(m, ll, 0, 0)
+        dirac_data = np.dot(A_matrix, dirac)
+
+        sh = cvxpy.Variable(A_matrix.shape[1])
+        constraints = [L_positivity * sh >= 0.]
+        cost = cvxpy.sum_squares(A_matrix * sh - dirac_data)
+        problem = cvxpy.Problem(cvxpy.Minimize(cost), constraints)
+        problem.solve()
+
+        positive_dirac = sh.value[m == 0]
+        mult = positive_dirac / dirac[m == 0]
+
+        for i, ll_ in enumerate(ll):
+            A_matrix[:, i] /= mult[ll_ // 2]
+
+        return A_matrix
+
 
 def homogenize_x0_to_data(data, x0):
     """

dmipy/core/tests/test_return_filtered_signal.py

@@ -0,0 +1,32 @@
+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_free_water_elimination(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 = data_iso + data_aniso
+
+    csd_fit_S0 = mtcsd.fit(scheme, data_to_fit, fit_S0_response=True)
+    data_recovered_aniso = csd_fit_S0.return_filtered_signal(
+        ['partial_volume_0'])
+    np.testing.assert_array_almost_equal(
+        data_aniso, data_recovered_aniso[0], 1)