#!/usr/bin/env python
# coding: utf-8

# # Spherical Mean Representation of any Compartment Model

# The idea of Spherical Mean models is to fit the spherical mean of a model to the spherical mean of the data, rather than fitting the separate DWIs per shell. As such, spherical mean models have no concept of orientation, dispersion or crossings. *(Kaden et al. 2015)* proposed the spherical mean technique by providing the analytic expression of the Zeppelin model. Later, *(Kaden et al. 2016)* combined the spherical mean model of a Zeppelin and Stick model to create his multi-compartment spherical mean technique.
# 
# Taking a step forward, Dmipy allows to estimate the spherical mean representation of any Compartment Model, and combine them in a MultiCompartmentSphericalMeanModel in the same way as a standard MultiCompartmentModel. The spherical mean for a model is estimated as follows:
# - If available, the spherical mean for a model is estimated using an analytic expression (e.g. for C1Stick, G2Zeppelin and G3RestrictedZeppelin).
# - For round models such as G1Ball or restricted sphere models, the spherical mean is just a value for the given shell parameters.
# - For cylinder models without analytic expressions, we estimate the spherical mean using the its rotational_harmonics_representation. In short, the spherical mean of any spherical function (if sufficiently expanded) is equal to $c_{m=0}^{l=0} / (2 \sqrt{\pi})$, where $c_{m=0}^{l=0}$ is the zeroth order spherical harmonics coefficient.
# 
# A model's spherical mean for a given acquisition scheme can be called using model.spherical_mean(acquisition_scheme, **parameters). To illustrate, we show the composition of the well-known Wu-Minn HCP acquisition scheme

# In[1]:


# as an example we estimate the spherical mean of the Gaussian phase cylinder model.
from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme
scheme = wu_minn_hcp_acquisition_scheme()
scheme.print_acquisition_info


# It can be seen that this scheme has 4 measured shells (including the b0 measurements). This means that the spherical mean for any of dmipy's models will also only have 4 values. To illustrate, we estimate the spherical mean of the Gaussian phase cylinder approximation (which does not have an analytical representation):

# In[2]:


from dmipy.signal_models import cylinder_models
cylinder = cylinder_models.C4CylinderGaussianPhaseApproximation()

diameter = 2e-6  # diameter of 2 microns
lambda_par = 1.7e-9  # parallel intra-axonal diffusivity
cylinder.spherical_mean(scheme, diameter=diameter, lambda_par=lambda_par)


# Notice that for the spherical mean respresentation it is not necessary to give the model a value for its orientation 'mu'. This is because a spherical mean model has no concept of orientation (as it has been integrated out of the model representation). Currently, the dmipy models that have a spherical mean representation are the following:
# - cylinder_models.C1Stick
# - cylinder_models.C2CylinderStejskalTannerApproximation
# - cylinder_models.C3CylinderCallaghanApproximation
# - cylinder_models.C4CylinderGaussianPhaseApproximation
# - gaussian_models.G1Ball
# - gaussian_models.G2Zeppelin
# - gaussian_models.G3TemporalZeppelin
# - sphere_models.S1Dot
# - sphere_models.S2SphereStejskalTannerApproximation
# - sphere_models.S4SphereGaussianPhaseApproximation
# 
# In addition, the spherical mean is also available for any GammaDistributed model.

# # Constructing a MultiCompartmentSphericalMeanModel

# A MultiCompartmentModel for spherical mean representations of models is done by calling MultiCompartmentSphericalMeanModel with the same parameters as the MultiCompartmentModel. For example, making a spherical mean representation of the Multi Compartment Microscopic Diffusion Imaging (MC-MDI) *(Kaden et al. 2016)* model is done as follows:

# In[3]:


from dmipy.signal_models import cylinder_models, gaussian_models
from dmipy.core.modeling_framework import MultiCompartmentSphericalMeanModel
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()
mc_mdi_model = MultiCompartmentSphericalMeanModel(
    models=[stick, zeppelin])
mc_mdi_model.parameter_names


# Notice that the spherical mean representation has no orientation parameters 'mu'. Adding tortuosity and equality constraints on the model parameters is done the same as before:

# In[4]:


mc_mdi_model.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
                                    'C1Stick_1_lambda_par',
                                    'partial_volume_0',
                                    'partial_volume_1')
mc_mdi_model.set_equal_parameter('G2Zeppelin_1_lambda_par',
                                 'C1Stick_1_lambda_par')
mc_mdi_model.parameter_names


# With this the spherical mean model can be fitted to dMRI data as mc_mdi_model.fit(scheme, data). The application of this model on HCP data is given in the [MC-MDI model example](http://nbviewer.jupyter.org/github/AthenaEPI/mipy/blob/master/examples/example_multi_compartment_spherical_mean_technique.ipynb).

# ## References
# - Kaden, Enrico, Frithjof Kruggel, and Daniel C. Alexander. "Quantitative mapping of the per‐axon diffusion coefficients in brain white matter." Magnetic resonance in medicine 75.4 (2016): 1752-1763.
# - Kaden, Enrico, et al. "Multi-compartment microscopic diffusion imaging." NeuroImage 139 (2016): 346-359.