From fc2d6ffe354f054c544f957e0c59fd5f6ae04d7d Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Mon, 28 Oct 2019 15:26:57 +0100 Subject: [PATCH 1/2] add GTM parameters in set_tortuous_parameter This commit adds the `s0_responses` parameter to the `set_tortuous_parameter` function in the `DistributeModel` class. This allows a correct estimation of the volume fraction of the IC compartment which is not affected by the signal fraction bias. --- dmipy/distributions/distribute_models.py | 10 ++++++++-- dmipy/utils/utils.py | 18 ++++++++++++------ 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/dmipy/distributions/distribute_models.py b/dmipy/distributions/distribute_models.py index a1d162ef..1f769527 100644 --- a/dmipy/distributions/distribute_models.py +++ b/dmipy/distributions/distribute_models.py @@ -275,7 +275,8 @@ def set_fixed_parameter(self, parameter_name, value): def set_tortuous_parameter(self, lambda_perp, lambda_par, - volume_fraction_intra): + volume_fraction_intra, + s0_responses=None): """ Allows the user to set a tortuosity constraint on the perpendicular diffusivity of the extra-axonal compartment, which depends on the @@ -295,6 +296,8 @@ def set_tortuous_parameter(self, lambda_perp, volume_fraction_intra: string name of the intra-axonal volume fraction parameter, see self.parameter_names. + s0_responses: list + list of the s0 responses of the tissues associated to IC and EC. """ params = [lambda_perp, lambda_par, volume_fraction_intra] for param in params: @@ -308,7 +311,10 @@ def set_tortuous_parameter(self, lambda_perp, model, name = self._parameter_map[lambda_perp] self.parameter_links.append([model, name, T1_tortuosity, [ self._parameter_map[lambda_par], - self._parameter_map[volume_fraction_intra]] + self._parameter_map[volume_fraction_intra], + None, + s0_responses + ] ]) del self.parameter_ranges[lambda_perp] del self.parameter_scales[lambda_perp] diff --git a/dmipy/utils/utils.py b/dmipy/utils/utils.py index 41c6ef28..874cbdac 100644 --- a/dmipy/utils/utils.py +++ b/dmipy/utils/utils.py @@ -170,7 +170,7 @@ def rotation_matrix_100_to_theta_phi_psi(theta, phi, psi): return np.dot(R_100_to_theta_phi, R_around_100) -def T1_tortuosity(lambda_par, vf_intra, vf_extra=None): +def T1_tortuosity(lambda_par, vf_intra, vf_extra=None, s0_responses=None): """Tortuosity model for perpendicular extra-axonal diffusivity [1, 2, 3]. If vf_extra=None, then vf_intra must be a nested volume fraction, in the sense that E_bundle = vf_intra * E_intra + (1 - vf_intra) * E_extra, with @@ -186,6 +186,8 @@ def T1_tortuosity(lambda_par, vf_intra, vf_extra=None): intra-axonal volume fraction [0, 1]. vf_extra : float, (optional) extra-axonal volume fraction [0, 1]. + s0_responses : list, (optional) + s0 response of the tissues associated to IC and EC. Returns ------- @@ -204,11 +206,15 @@ def T1_tortuosity(lambda_par, vf_intra, vf_extra=None): .. [3] Szafer et al. "Theoretical model for water diffusion in tissues." Magnetic resonance in medicine 33.5 (1995): 697-712. """ - if vf_extra is None: - lambda_perp = (1 - vf_intra) * lambda_par - else: - fraction_intra = vf_intra / (vf_intra + vf_extra) - lambda_perp = (1 - fraction_intra) * lambda_par + if s0_responses is not None: + vf_intra /= s0_responses[0] + + fraction_intra = vf_intra + if vf_extra is not None: + if s0_responses is not None: + vf_extra /= s0_responses[1] + fraction_intra /= (vf_intra + vf_extra) + lambda_perp = (1. - fraction_intra) * lambda_par return lambda_perp From 0748eaa1d35007734411f7739fcbea721aa2554e Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Thu, 31 Oct 2019 16:41:41 +0100 Subject: [PATCH 2/2] change behaviour of S0 response assignment in GTM This commit changes the way in which the S0 of every tissue is passed to the model constructors. Now it is allowed to pass None, if the corresponding model has been defined using the `S0_tissue_responses` parameter. All of this happens in the `MultiCompartmentModel.fit` function. Other changes involve the proper definition of the tortuosity constraint that have been addressed in the parent commits. --- dmipy/core/modeling_framework.py | 26 +++++++++++++++++++++--- dmipy/distributions/distribute_models.py | 6 +++--- dmipy/utils/utils.py | 22 ++++++++++---------- 3 files changed, 37 insertions(+), 17 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index c9f114a8..c1c95e21 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -558,7 +558,8 @@ def _add_fixed_parameter_array(self, parameter_name, parameter_array): def set_tortuous_parameter(self, lambda_perp_parameter_name, lambda_par_parameter_name, volume_fraction_intra_parameter_name, - volume_fraction_extra_parameter_name): + volume_fraction_extra_parameter_name, + S0_responses=None): """ Allows the user to set a tortuosity constraint on the perpendicular diffusivity of the extra-axonal compartment, which depends on the @@ -581,6 +582,8 @@ def set_tortuous_parameter(self, lambda_perp_parameter_name, volume_fraction_extra_parameter_name: string name of the extra-axonal volume fraction parameter, see self.parameter_names. + S0_responses : list, (optional) + s0 response of the tissues associated to IC and EC. """ params = [lambda_perp_parameter_name, lambda_par_parameter_name, volume_fraction_intra_parameter_name, @@ -597,7 +600,8 @@ def set_tortuous_parameter(self, lambda_perp_parameter_name, self.parameter_links.append([model, name, T1_tortuosity, [ self._parameter_map[lambda_par_parameter_name], self._parameter_map[volume_fraction_intra_parameter_name], - self._parameter_map[volume_fraction_extra_parameter_name]] + self._parameter_map[volume_fraction_extra_parameter_name], + S0_responses] ]) del self.parameter_ranges[lambda_perp_parameter_name] del self.parameter_cardinality[lambda_perp_parameter_name] @@ -1212,8 +1216,24 @@ def fit(self, acquisition_scheme, data, start = time() mt_fractions = np.empty( np.r_[N_voxels, self.N_models], dtype=float) + + S0_responses = [] + for i, (m, s) in enumerate(zip(self.models, + self.S0_tissue_responses)): + if s is None: + # If the S0 of the tissue is None, the model is composite + # and the corresponding S0 is treated as the weighted + # average of the S0 tissues modelled by the "submodels". + s0 = 0.0 + for k, mm in enumerate(m.models): + s0 += (mm['partial_volume_{}'.format(k)] * + mm.S0_tissue_responses[k]) + S0_responses.append(s0) + else: + S0_responses.append(s) + fit_func = MultiTissueConvexOptimizer( - acquisition_scheme, self, self.S0_tissue_responses) + acquisition_scheme, self, S0_responses) for idx, pos in enumerate(zip(*mask_pos)): voxel_S = data_[pos] parameters = fitted_parameters_lin[idx] diff --git a/dmipy/distributions/distribute_models.py b/dmipy/distributions/distribute_models.py index 1f769527..51b0227f 100644 --- a/dmipy/distributions/distribute_models.py +++ b/dmipy/distributions/distribute_models.py @@ -276,7 +276,7 @@ def set_fixed_parameter(self, parameter_name, value): def set_tortuous_parameter(self, lambda_perp, lambda_par, volume_fraction_intra, - s0_responses=None): + S0_responses=None): """ Allows the user to set a tortuosity constraint on the perpendicular diffusivity of the extra-axonal compartment, which depends on the @@ -296,7 +296,7 @@ def set_tortuous_parameter(self, lambda_perp, volume_fraction_intra: string name of the intra-axonal volume fraction parameter, see self.parameter_names. - s0_responses: list + S0_responses: list list of the s0 responses of the tissues associated to IC and EC. """ params = [lambda_perp, lambda_par, volume_fraction_intra] @@ -313,7 +313,7 @@ def set_tortuous_parameter(self, lambda_perp, self._parameter_map[lambda_par], self._parameter_map[volume_fraction_intra], None, - s0_responses + S0_responses ] ]) del self.parameter_ranges[lambda_perp] diff --git a/dmipy/utils/utils.py b/dmipy/utils/utils.py index 874cbdac..c42b6e36 100644 --- a/dmipy/utils/utils.py +++ b/dmipy/utils/utils.py @@ -170,7 +170,7 @@ def rotation_matrix_100_to_theta_phi_psi(theta, phi, psi): return np.dot(R_100_to_theta_phi, R_around_100) -def T1_tortuosity(lambda_par, vf_intra, vf_extra=None, s0_responses=None): +def T1_tortuosity(lambda_par, vf_intra, vf_extra=None, S0_responses=None): """Tortuosity model for perpendicular extra-axonal diffusivity [1, 2, 3]. If vf_extra=None, then vf_intra must be a nested volume fraction, in the sense that E_bundle = vf_intra * E_intra + (1 - vf_intra) * E_extra, with @@ -186,7 +186,7 @@ def T1_tortuosity(lambda_par, vf_intra, vf_extra=None, s0_responses=None): intra-axonal volume fraction [0, 1]. vf_extra : float, (optional) extra-axonal volume fraction [0, 1]. - s0_responses : list, (optional) + S0_responses : list, (optional) s0 response of the tissues associated to IC and EC. Returns @@ -206,16 +206,16 @@ def T1_tortuosity(lambda_par, vf_intra, vf_extra=None, s0_responses=None): .. [3] Szafer et al. "Theoretical model for water diffusion in tissues." Magnetic resonance in medicine 33.5 (1995): 697-712. """ - if s0_responses is not None: - vf_intra /= s0_responses[0] - fraction_intra = vf_intra - if vf_extra is not None: - if s0_responses is not None: - vf_extra /= s0_responses[1] - fraction_intra /= (vf_intra + vf_extra) - lambda_perp = (1. - fraction_intra) * lambda_par - return lambda_perp + fraction_extra = 1.0 - fraction_intra if vf_extra is None else vf_extra + + if S0_responses is not None: + fraction_intra /= S0_responses[0] + fraction_extra /= S0_responses[1] + + fraction_intra /= fraction_extra + fraction_intra + + return (1. - fraction_intra) * lambda_par def parameter_equality(param):