Skip to content

Commit

Permalink
Dhollander RF estimation (#17)
Browse files Browse the repository at this point in the history 
* dhollander16 produces 3 tissue_response_models representing WM/GM/CSF, useable with multi-shell data.
* dhollander16 uses tournier07 white matter response estimation.
* tissue response models can be used by any type of MC-model (MC, MC-SM, MC-CSD).
* added tests for tissue_response_models and their use in MC-models.
* added mt-csd example
* Does not use non-unity S0 intensities yet, only works on signal attenuation shape.
  • Loading branch information
@rutgerfick
rutgerfick committed Aug 29, 2018
1 parent f635d70 commit 4814485
Show file tree
Hide file tree
Showing 12 changed files with 1,230 additions and 182 deletions.
2 changes: 2 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +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])
- [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
Expand Down Expand Up @@ -90,6 +91,7 @@ Any Spherical Mean model can also estimate parametric FODs.
### Spherical Deconvolution-Based Models
Constrained spherical deconvolution (CSD) models are primarily used for Fiber Orientation Distribution (FOD) estimation. Multi-compartment CSD models improve FOD estimation by separating the signal contributions of white matter from CSF and grey matter.
- [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)
- [Multi-Tissue CSD [Jeurissen et al. 2014]](http://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/dhollander16/examples/example_multi_tissue_constrained_spherical_deconvolution.ipynb)

## 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).
Expand Down
2 changes: 1 addition & 1 deletion dmipy/core/acquisition_scheme.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ def __init__(self, bvalues, gradient_directions, qvalues,
# coefficients to the positions on the sphere for every shell
self.unique_dwi_indices = np.unique(self.shell_indices[~self.b0_mask])
self.shell_sh_matrices = {}
self.shell_sh_orders = np.zeros(len(self.shell_bvalues))
self.shell_sh_orders = np.zeros(len(self.shell_bvalues), dtype=int)
for shell_index in self.unique_dwi_indices:
shell_mask = self.shell_indices == shell_index
bvecs_shell = self.gradient_directions[shell_mask]
Expand Down
1 change: 1 addition & 0 deletions dmipy/core/fitted_modeling_framework.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from ..utils.spherical_mean import (
estimate_spherical_mean_multi_shell)
from dipy.reconst.shm import sph_harm_ind_list
from dipy.direction import peak_directions


__all__ = [
Expand Down
172 changes: 131 additions & 41 deletions dmipy/core/modeling_framework.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@
FittedMultiCompartmentModel,
FittedMultiCompartmentSphericalMeanModel,
FittedMultiCompartmentSphericalHarmonicsModel)
from ..optimizers_fod.construct_observation_matrix import (
construct_model_based_A_matrix)
from ..optimizers.brute2fine import (
GlobalBruteOptimizer, Brute2FineOptimizer)
from ..optimizers_fod.csd_tournier import CsdTournierOptimizer
Expand Down Expand Up @@ -792,6 +794,35 @@ def _add_parameter_nodes(self, graph_model, entry_uuid, entry_model):
graph_model.node(parameter_uuid, parameter_name)
graph_model.edge(parameter_uuid, entry_uuid)

def _check_tissue_model_acquisition_scheme(self, acquisition_scheme):
"""Tests if acquisition scheme between MC-model and tissue response
model are the same.
Parameters
----------
acquisition_scheme : DmipyAcquisitionScheme instance,
An acquisition scheme that has been instantiated using dMipy.
"""
for model in self.models:
if model._model_type == 'TissueResponseModel':
mc_scheme_params = [
acquisition_scheme.shell_bvalues,
acquisition_scheme.shell_delta,
acquisition_scheme.shell_Delta,
acquisition_scheme.shell_gradient_strengths]
tr_scheme_params = [
model.acquisition_scheme.shell_bvalues,
model.acquisition_scheme.shell_delta,
model.acquisition_scheme.shell_Delta,
model.acquisition_scheme.shell_gradient_strengths]
try:
np.testing.assert_array_almost_equal(
mc_scheme_params, tr_scheme_params)
except AssertionError:
msg = "Acquisition scheme of MC-model and tissue response "
msg += "model are not the same."
raise ValueError(msg)


class MultiCompartmentModel(MultiCompartmentModelProperties):
r'''
Expand Down Expand Up @@ -920,6 +951,7 @@ def fit(self, acquisition_scheme, data,
Can be used to recover parameters themselves or other useful
functions.
"""
self._check_tissue_model_acquisition_scheme(acquisition_scheme)

# estimate S0
self.scheme = acquisition_scheme
Expand Down Expand Up @@ -1281,6 +1313,7 @@ def fit(self, acquisition_scheme, data,
Can be used to recover parameters themselves or other useful
functions.
"""
self._check_tissue_model_acquisition_scheme(acquisition_scheme)

# estimate S0
self.scheme = acquisition_scheme
Expand Down Expand Up @@ -1526,6 +1559,7 @@ def __init__(self, models, sh_order=8):
self._prepare_parameters_to_optimize()
self._add_spherical_harmonics_parameters(sh_order)
self._check_that_one_anisotropic_kernel_is_present()
# self._check_for_tissue_response_models()

self.x0_parameters = {}
self.sh_order = sh_order
Expand Down Expand Up @@ -1682,8 +1716,8 @@ def fit(self, acquisition_scheme, data, mask=None, solver='csd',
Official Journal of the International Society for Magnetic
Resonance in Medicine 58.3 (2007): 497-510.
"""

self._check_if_kernel_parameters_are_fixed()
self._check_tissue_model_acquisition_scheme(acquisition_scheme)

self.voxel_varying_kernel = False
if bool(self.x0_parameters): # if the dictionary is not empty
Expand All @@ -1696,9 +1730,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)

# estimate S0
self.scheme = acquisition_scheme
data_ = np.atleast_2d(data)
# if self.tissue_response_kernels_present:
# S0 = np.ones(data_.shape[:-1], dtype=float) * S0_responses.max()
if self.scheme.TE is None or len(np.unique(self.scheme.TE)) == 1:
S0 = np.mean(data_[..., self.scheme.b0_mask], axis=-1)
else: # if multiple TE are in the data
Expand Down Expand Up @@ -1878,52 +1917,103 @@ def __call__(self, acquisition_scheme, **kwargs):
kwargs
)

if len(self.models) > 1:
partial_volumes = [
kwargs[p] for p in self.partial_volume_names
]
A = self._construct_convolution_kernel(
self.parameters_to_parameter_vector(**kwargs))

# if vf fixed then just multiply with sh_coeff
if self.volume_fractions_fixed:
E = np.dot(A, kwargs['sh_coeff'])
else:
partial_volumes = [1.]
sh_coeff = np.zeros(self.optimizer.Ncoef_total)
for i, name in enumerate(self.partial_volume_names):
sh_coeff[self.optimizer.vf_indices[i]] = (
kwargs[name] / (2 * np.sqrt(np.pi)))
sh_coeff[self.optimizer.sh_start:
self.optimizer.Ncoef + self.optimizer.sh_start] = kwargs[
'sh_coeff']
E = np.dot(A, sh_coeff)
return E

rh_models = []
sh_coeffs = []
for model, partial_volume in zip(self.models, partial_volumes):
parameters = {}
for parameter in model.parameter_ranges:
parameter_name = self._inverted_parameter_map[
(model, parameter)
]
parameters[parameter] = kwargs.get(
parameter_name
)
rh_model = model.rotational_harmonics_representation(
acquisition_scheme, **parameters)
rh_models.append(rh_model)
if 'orientation' in model.parameter_types.values():
sh_coeffs.append(kwargs['sh_coeff'])
else:
sh_coeffs.append(
np.atleast_1d(partial_volume / (2 * np.sqrt(np.pi))))
def _construct_convolution_kernel(self, x0_vector):
"""
Helper function that constructs the convolution kernel for the given
multi-compartment model and the initial condition x0_vector.
First the parameter vector is converted to a dictionary with the
corresponding parameter names. Then, the linked parameters are added to
the given ones. Finally, the rotational harmonics of the model is
passed to the construct_model_based_A_matrix, which constructs the
kernel for an arbitrary PGSE-acquisition scheme.
For multiple models with fixed volume fractions, the A-matrices
are combined to have a combined convolution kernel.
E = np.ones(acquisition_scheme.number_of_measurements)
for i, shell_index in enumerate(acquisition_scheme.unique_dwi_indices):
shell_mask = acquisition_scheme.shell_indices == shell_index
sh_mat = acquisition_scheme.shell_sh_matrices[shell_index]
sh_order = int(acquisition_scheme.shell_sh_orders[shell_index])
For multiple models without fixed volume fractions, the convolution
kernels for anisotropic and isotropic models are concatenated, with
the isotropic kernels always having a spherical harmonics order of 0.
Parameters
----------
x0_vector: array of size (N_parameters),
Contains the fixed parameters of the convolution kernel.
shell_E = 0.
for j, model in enumerate(self.models):
Returns
-------
kernel: array of size (N_coef, N_data),
Observation matrix that maps the FOD spherical harmonics
coefficients to the DWI signal values.
"""
parameters_dict = self.parameter_vector_to_parameters(
x0_vector)
parameters_dict = self.add_linked_parameters_to_parameters(
parameters_dict)

if self.volume_fractions_fixed:
if len(self.models) > 1:
partial_volumes = [
parameters_dict[p] for p in self.partial_volume_names
]
else:
partial_volumes = [1.]
kernel = 0.
for model, partial_volume in zip(self.models,
partial_volumes):
parameters = {}
for parameter in model.parameter_ranges:
parameter_name = self._inverted_parameter_map[
(model, parameter)
]
parameters[parameter] = parameters_dict.get(
parameter_name
)
model_rh = (
model.rotational_harmonics_representation(
self.scheme, **parameters))
kernel += partial_volume * construct_model_based_A_matrix(
self.scheme, model_rh, self.sh_order)
else:
kernel = []
for model in self.models:
parameters = {}
for parameter in model.parameter_ranges:
parameter_name = self._inverted_parameter_map[
(model, parameter)
]
parameters[parameter] = parameters_dict.get(
parameter_name
)
model_rh = (
model.rotational_harmonics_representation(
self.scheme, **parameters))
if 'orientation' in model.parameter_types.values():
E_dispersed_sh = sh_convolution(
sh_coeffs[j], rh_models[j][i, :sh_order // 2 + 1])
shell_E += np.dot(sh_mat[:, :len(E_dispersed_sh)],
E_dispersed_sh)
kernel.append(construct_model_based_A_matrix(
self.scheme, model_rh, self.sh_order))
else:
E_dispersed_sh = sh_convolution(
sh_coeffs[j], rh_models[j][i])
shell_E += np.dot(sh_mat[0, 0], E_dispersed_sh)
E[shell_mask] = shell_E
return E
kernel.append(construct_model_based_A_matrix(
self.scheme, model_rh, 0))

kernel = np.hstack(kernel)
return kernel


def homogenize_x0_to_data(data, x0):
Expand Down
2 changes: 1 addition & 1 deletion dmipy/optimizers_fod/construct_observation_matrix.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ def construct_model_based_A_matrix(acquisition_scheme, model_rh, lmax):
# prepare the rotational harmonics of the kernel
counter = 0
for n_ in range(0, lmax + 1, 2):
if n_ // 2 > model_rh.shape[1]:
if n_ // 2 >= model_rh.shape[1]:
# in case an isotropic kernel is given
break
coef_in_order = 2 * n_ + 1
Expand Down
86 changes: 3 additions & 83 deletions dmipy/optimizers_fod/csd_cvxpy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ def __init__(self, acquisition_scheme, model, x0_vector=None, sh_order=8,
x0_vector, (-1, x0_vector.shape[-1]))[0]
if np.all(np.isnan(x0_single_voxel)):
self.single_convolution_kernel = True
self.A = self._construct_convolution_kernel(
self.A = self.model._construct_convolution_kernel(
x0_single_voxel)
else:
self.single_convolution_kernel = False
Expand All @@ -109,7 +109,7 @@ def __init__(self, acquisition_scheme, model, x0_vector=None, sh_order=8,
self.vf_indices = np.where(np.hstack(vf_array))[0]

sh_l = sph_harm_ind_list(sh_order)[1]
lb_weights = sh_l ** 2 * (sh_l + 1) ** 2 # laplace-beltrami
lb_weights = sh_l ** 2 * (sh_l + 1) ** 2 # laplace-beltrami [3]
if self.model.volume_fractions_fixed:
self.R_smoothness = np.diag(lb_weights)
else:
Expand Down Expand Up @@ -145,7 +145,7 @@ def __call__(self, data, x0_vector):
if self.single_convolution_kernel:
A = self.A
else:
A = self._construct_convolution_kernel(x0_vector)
A = self.model._construct_convolution_kernel(x0_vector)

sh_coef = cvxpy.Variable(self.Ncoef_total)
sh_fod = sh_coef[self.sh_start: self.Ncoef + self.sh_start]
Expand Down Expand Up @@ -177,83 +177,3 @@ def __call__(self, data, x0_vector):
fitted_parameter_vector = self.model.parameters_to_parameter_vector(
**fitted_params)
return fitted_parameter_vector

def _construct_convolution_kernel(self, x0_vector):
"""
Helper function that constructs the convolution kernel for the given
multi-compartment model and the initial condition x0_vector.
First the parameter vector is converted to a dictionary with the
corresponding parameter names. Then, the linked parameters are added to
the given ones. Finally, the rotational harmonics of the model is
passed to the construct_model_based_A_matrix, which constructs the
kernel for an arbitrary PGSE-acquisition scheme.
For multiple models with fixed volume fractions, the A-matrices
are combined to have a combined convolution kernel.
For multiple models without fixed volume fractions, the convolution
kernels for anisotropic and isotropic models are concatenated, with
the isotropic kernels always having a spherical harmonics order of 0.
Parameters
----------
x0_vector: array of size (N_parameters),
Contains the fixed parameters of the convolution kernel.
Returns
-------
kernel: array of size (N_coef, N_data),
Observation matrix that maps the FOD spherical harmonics
coefficients to the DWI signal values.
"""
parameters_dict = self.model.parameter_vector_to_parameters(
x0_vector)
parameters_dict = self.model.add_linked_parameters_to_parameters(
parameters_dict)

if self.model.volume_fractions_fixed:
if len(self.model.models) > 1:
partial_volumes = [
parameters_dict[p] for p in self.model.partial_volume_names
]
else:
partial_volumes = [1.]
kernel = 0.
for model, partial_volume in zip(self.model.models,
partial_volumes):
parameters = {}
for parameter in model.parameter_ranges:
parameter_name = self.model._inverted_parameter_map[
(model, parameter)
]
parameters[parameter] = parameters_dict.get(
parameter_name
)
model_rh = (
model.rotational_harmonics_representation(
self.acquisition_scheme, **parameters))
kernel += partial_volume * construct_model_based_A_matrix(
self.acquisition_scheme, model_rh, self.sh_order)
else:
kernel = []
for model in self.model.models:
parameters = {}
for parameter in model.parameter_ranges:
parameter_name = self.model._inverted_parameter_map[
(model, parameter)
]
parameters[parameter] = parameters_dict.get(
parameter_name
)
model_rh = (
model.rotational_harmonics_representation(
self.acquisition_scheme, **parameters))
if 'orientation' in model.parameter_types.values():
kernel.append(construct_model_based_A_matrix(
self.acquisition_scheme, model_rh, self.sh_order))
else:
kernel.append(construct_model_based_A_matrix(
self.acquisition_scheme, model_rh, 0))
kernel = np.hstack(kernel)
return kernel

0 comments on commit 4814485

Please sign in to comment.