Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

dti fit #103

Open
fayemckenna opened this issue Sep 24, 2020 · 8 comments
Open

dti fit #103

fayemckenna opened this issue Sep 24, 2020 · 8 comments

Comments

@fayemckenna
Copy link

the dipy dtifit is often called to fit the tensor and calculate the classical FA, AD and RD measures in various dmipy models. I would like to incorporate an anisotropic diffusion compartment into a multicompartmentmodel which I can then calculate FA, MD, RD, AD. Would you suggest calling dtifit within the multicomparmentmodel function? Thank you!

@rutgerfick
Copy link
Collaborator

rutgerfick commented Sep 25, 2020

HI Faye,

So if I understand correctly, you want to calculate the DTI metrics on part of a multi-compartment model?
Can you give me a little more information on exactly what model you are using and on which part you want these metrics?

@fayemckenna
Copy link
Author

fayemckenna commented Sep 25, 2020 via email

@rutgerfick
Copy link
Collaborator

So, the way to do this is quite simple but not directly "built into" the multi compartment model itself.

What you'll have to do is get the fitted parameters of your multi-compartment model and use those parameters to generate the signal of only the compartment you're interested in.

So as a simple example for a dummy ball and cylinder model

from dmipy.signal_models import cylinder_models, gaussian_models
from dmipy.core import modeling_framework
from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme
import numpy as np

scheme = wu_minn_hcp_acquisition_scheme()

# create some ball and cylinder model
ball = gaussian_models.G1Ball()
cyl = cylinder_models.C2CylinderStejskalTannerApproximation()
ballcyl = modeling_framework.MultiCompartmentModel([ball, cyl])

# some fake dummy parameters and data
params = {'G1Ball_1_lambda_iso': 2e-9,
         'C2CylinderStejskalTannerApproximation_1_lambda_par': 1.7e-9,
         'C2CylinderStejskalTannerApproximation_1_diameter': 6e-6,
         'C2CylinderStejskalTannerApproximation_1_mu':[0, 0],
         'partial_volume_0': .2,
         'partial_volume_1': .8}
dummy_data = ballcyl(scheme, **params)

# fit the ballcyl model to get parameters
fitted_ballcyl = ballcyl.fit(scheme, dummy_data)
fitted_params = fitted_ballcyl.fitted_parameters

# get the cylinder parameters and put them directly in the cylinder model to get
# the signal of that compartment only, the volume fractions don't matter in this case
# since DTI metrics only care about the "shape" of the data and not the amplitude
signal_cyl_only = cyl(acquisition_scheme=scheme,
                      lambda_par=fitted_params['C2CylinderStejskalTannerApproximation_1_lambda_par'],
                      mu=fitted_params['C2CylinderStejskalTannerApproximation_1_mu'][0],  # note the [0]
                      diameter=fitted_params['C2CylinderStejskalTannerApproximation_1_diameter'])

# put the cylinder signal into the dti tensor fit and get DTI metrics as usual
from dmipy.core.acquisition_scheme import gtab_dmipy2dipy
gtab_dipy = gtab_dmipy2dipy(scheme)

from dipy.reconst.dti import TensorModel
tenmod = TensorModel(gtab_dipy)
tenfit = tenmod.fit(cyl_only)
print(tenfit.fa, tenfit.md)
>> 0.9730780509802013 0.000600906351945028

Is this what you wanted to do?

@fayemckenna
Copy link
Author

fayemckenna commented Sep 27, 2020 via email

@fayemckenna
Copy link
Author

fayemckenna commented Sep 28, 2020 via email

@fayemckenna
Copy link
Author

fayemckenna commented Oct 4, 2020 via email

@rutgerfick
Copy link
Collaborator

Hi Faye,

To make it work the way (I think) you want, I'll slightly adjust the script.
If you make the multicompartment model with only the cylinder you can just pass all the multi-dimensional output parameters of the previous model as parameter in the cyl_model.simulate_signal function.

The reason why it was failing before is because the base cylinder model (outside of the multi-compartment model class) is not really built for multi-dimensional input. Put it in an MC-model and everything works fine :)

Let me know if this is what you actually wanted in the end!

from dmipy.signal_models import cylinder_models, gaussian_models
from dmipy.core import modeling_framework
from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme
import numpy as np

scheme = wu_minn_hcp_acquisition_scheme()

# create some ball and cylinder model
ball = gaussian_models.G1Ball()
cyl = cylinder_models.C2CylinderStejskalTannerApproximation()
ballcyl = modeling_framework.MultiCompartmentModel([ball, cyl])

# some fake dummy parameters and data
params = {'G1Ball_1_lambda_iso': 2e-9,
         'C2CylinderStejskalTannerApproximation_1_lambda_par': 1.7e-9,
         'C2CylinderStejskalTannerApproximation_1_diameter': 6e-6,
         'C2CylinderStejskalTannerApproximation_1_mu':[0, 0],
         'partial_volume_0': .2,
         'partial_volume_1': .8}
dummy_data = ballcyl(scheme, **params)

# make multi_dimensional data
dummy_data = np.tile(dummy_data, [10, 1])

# fit the ballcyl model to get parameters
fitted_ballcyl = ballcyl.fit(scheme, dummy_data)
fitted_params = fitted_ballcyl.fitted_parameters

# get the cylinder parameters and put them directly in the cylinder model to get
# the signal of that compartment only, the volume fractions don't matter in this case
# since DTI metrics only care about the "shape" of the data and not the amplitude
cyl_only_model = modeling_framework.MultiCompartmentModel([cyl])
parameters_for_cyl_only = {
 'C2CylinderStejskalTannerApproximation_1_mu': fitted_params['C2CylinderStejskalTannerApproximation_1_mu'],
 'C2CylinderStejskalTannerApproximation_1_diameter': fitted_params['C2CylinderStejskalTannerApproximation_1_diameter'],
 'C2CylinderStejskalTannerApproximation_1_lambda_par': fitted_params['C2CylinderStejskalTannerApproximation_1_lambda_par']
}
signal_cyl_only = cyl_only_model.simulate_signal(
    acquisition_scheme=scheme,
    parameters_array_or_dict=parameters_for_cyl_only)

# put the cylinder signal into the dti tensor fit and get DTI metrics as usual
from dmipy.core.acquisition_scheme import gtab_dmipy2dipy
gtab_dipy = gtab_dmipy2dipy(scheme)

from dipy.reconst.dti import TensorModel
tenmod = TensorModel(gtab_dipy)
tenfit = tenmod.fit(signal_cyl_only)
print(tenfit.fa, tenfit.md)
>> [0.97307805 0.97307805 0.97307805 0.97307805 0.97307805 0.97307805
 0.97307805 0.97307805 0.97307805 0.97307805] [0.00060091 0.00060091 0.00060091 0.00060091 0.00060091 0.00060091
 0.00060091 0.00060091 0.00060091 0.00060091]

@fayemckenna
Copy link
Author

fayemckenna commented Oct 26, 2020 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants