Skip to content

Commit

Permalink
Refactor model properties (#41)
Browse files Browse the repository at this point in the history 
* refactored spherical mean, rotational harmonics representation and convolution response matrix into signal properties
* working rotational harmonics, spherical mean and convolution kernel for distributed models
* spherical distributed models can now be used in MC-SM and MC-CSD models.
  • Loading branch information
@rutgerfick
rutgerfick committed Nov 18, 2018
1 parent e1ba0bd commit 0e49180
Show file tree
Hide file tree
Showing 16 changed files with 398 additions and 637 deletions.
2 changes: 1 addition & 1 deletion .travis.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
20 changes: 16 additions & 4 deletions dmipy/core/acquisition_scheme.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,9 +152,12 @@ def __init__(self, bvalues, gradient_directions, qvalues,

# calculates observation matrices to convert spherical harmonic
# coefficients to the positions on the sphere for every shell
self.unique_b0_indices = np.unique(self.shell_indices[self.b0_mask])
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), dtype=int)
self.shell_sh_orders = {}
for shell_index in self.unique_b0_indices:
self.shell_sh_orders[shell_index] = 0
for shell_index in self.unique_dwi_indices:
shell_mask = self.shell_indices == shell_index
bvecs_shell = self.gradient_directions[shell_mask]
Expand Down Expand Up @@ -409,6 +412,7 @@ def __init__(self, dmipy_acquisition_scheme, N_angular_samples=10):
Gdirs_all_shells = []
delta_all_shells = []
Delta_all_shells = []
shell_indices = []
for shell_index in scheme.unique_dwi_indices:
b = scheme.shell_bvalues[shell_index]
b_all_shells.append(np.tile(b, N_angular_samples))
Expand All @@ -419,7 +423,9 @@ def __init__(self, dmipy_acquisition_scheme, N_angular_samples=10):
Delta = scheme.shell_Delta[shell_index]
Delta_all_shells.append(np.tile(Delta, N_angular_samples))
Gdirs_all_shells.append(angles_cart)
shell_indices.append(np.tile(shell_index, N_angular_samples))

self.shell_indices = np.hstack(shell_indices)
self.bvalues = np.hstack(b_all_shells)
self.gradient_directions = np.vstack(Gdirs_all_shells)
self.delta = None
Expand All @@ -442,11 +448,17 @@ def __init__(self, dmipy_acquisition_scheme, N_angular_samples=10):
self.b0_mask = np.tile(False, len(self.bvalues))
self.shell_delta = scheme.shell_delta
self.shell_Delta = scheme.shell_Delta
self.shell_sh_orders = (
np.array(scheme.shell_sh_orders[scheme.unique_dwi_indices],
dtype=int))
self.unique_dwi_indices = scheme.unique_dwi_indices
self.number_of_measurements = len(self.bvalues)

self.shell_sh_matrices = {}
self.shell_sh_orders = {}
for shell_index in scheme.unique_dwi_indices:
self.shell_sh_orders[shell_index] = int(
scheme.shell_sh_orders[shell_index])
self.shell_sh_matrices[shell_index] = real_sym_sh_mrtrix(
self.shell_sh_orders[shell_index], thetas, phis)[0]

self.inverse_rh_matrix = {
rh_order: np.linalg.pinv(real_sym_rh_basis(
rh_order, thetas, phis
Expand Down
174 changes: 78 additions & 96 deletions dmipy/core/modeling_framework.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
FittedMultiCompartmentModel,
FittedMultiCompartmentSphericalMeanModel,
FittedMultiCompartmentSphericalHarmonicsModel)
from ..optimizers_fod.construct_observation_matrix import (
from ..utils.construct_observation_matrix import (
construct_model_based_A_matrix)
from ..optimizers.brute2fine import (
GlobalBruteOptimizer, Brute2FineOptimizer)
Expand Down Expand Up @@ -837,6 +837,80 @@ def _check_tissue_model_acquisition_scheme(self, acquisition_scheme):
msg += "model are not the same."
raise ValueError(msg)

def _construct_convolution_kernel(self, **kwargs):
"""
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.add_linked_parameters_to_parameters(kwargs)

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, S0 in zip(self.models,
partial_volumes,
self.S0_responses):
parameters = {}
for parameter in model.parameter_ranges:
parameter_name = self._inverted_parameter_map[
(model, parameter)
]
parameters[parameter] = parameters_dict.get(
parameter_name
)
kernel += S0 * partial_volume * (
model.convolution_kernel_matrix(
self.scheme, self.sh_order, **parameters))
else:
kernel = []
for model, S0 in zip(self.models, self.S0_responses):
parameters = {}
for parameter in model.parameter_ranges:
parameter_name = self._inverted_parameter_map[
(model, parameter)
]
parameters[parameter] = parameters_dict.get(
parameter_name
)
if 'orientation' in model.parameter_types.values():
kernel.append(S0 * model.convolution_kernel_matrix(
self.scheme, self.sh_order, **parameters))
else:
kernel.append(S0 * model.convolution_kernel_matrix(
self.scheme, 0, **parameters))

kernel = np.hstack(kernel)
return kernel


class MultiCompartmentModel(MultiCompartmentModelProperties):
r'''
Expand Down Expand Up @@ -1204,7 +1278,7 @@ def __init__(self, models, parameter_links=None):
if parameter_links is None:
self.parameter_links = []

self._check_for_dispersed_or_NMR_models()
self._check_for_NMR_models()
self._prepare_parameters()
self._delete_orientation_parameters()
self._prepare_partial_volumes()
Expand All @@ -1223,16 +1297,11 @@ def __init__(self, models, parameter_links=None):
msg += "multicore processing."
print(msg)

def _check_for_dispersed_or_NMR_models(self):
def _check_for_NMR_models(self):
for model in self.models:
if model._model_type is 'NMRModel':
msg = "Cannot estimate spherical mean of 1D-NMR models."
raise ValueError(msg)
if model._model_type is 'SphericalDistributedModel':
msg = "Cannot estimate spherical mean spherically distributed "
msg += "model. Please give the input models to the distributed"
msg += " model directly to MultiCompartmentSphericalMeanModel."
raise ValueError(msg)

def _delete_orientation_parameters(self):
"""
Expand Down Expand Up @@ -1963,16 +2032,11 @@ def __call__(self, acquisition_scheme, **kwargs):
kwargs: keyword arguments to the model parameter values,
Is internally given as **parameter_dictionary.
"""
kwargs = self.add_linked_parameters_to_parameters(
kwargs
)
A = self._construct_convolution_kernel(**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))

# if vf fixed then just multiply with sh_coeff
if self.volume_fractions_fixed:
E = np.dot(A, kwargs['sh_coeff'])
Expand All @@ -1987,88 +2051,6 @@ def __call__(self, acquisition_scheme, **kwargs):
E = np.dot(A, sh_coeff)
return E

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.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, S0 in zip(self.models,
partial_volumes,
self.S0_responses):
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 += S0 * partial_volume * construct_model_based_A_matrix(
self.scheme, model_rh, self.sh_order)
else:
kernel = []
for model, S0 in zip(self.models, self.S0_responses):
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():
kernel.append(S0 * construct_model_based_A_matrix(
self.scheme, model_rh, self.sh_order))
else:
kernel.append(S0 * 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

0 comments on commit 0e49180

Please sign in to comment.