AthenaEPI / dmipy

The open source toolbox for reproducible diffusion MRI-based microstructure estimation
MIT License
96 stars 30 forks source link

dti fit #103

Open fayemckenna opened 3 years ago

fayemckenna commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

Hi,

Yes that is correct! There is a chance I am just overlooking some simple way to do this ( perhaps include a cylinder as a tensor into the multicomparmental model). But since I don’t see an example of it, only calling upon dipy’s tensor fit in other models I want to double check.

Thank you for your help!

Faye

On Fri, Sep 25, 2020 at 10:19 AM Rutger Fick notifications@github.com wrote:

HI Faye,

So if I understand correctly, you want to calculate the DTI metrics on part of a multi-compartment model?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/AthenaEPI/dmipy/issues/103#issuecomment-698956780, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHTA6N42BUTPSNQQNJBKE6LSHSRHJANCNFSM4RX2EVNA .

rutgerfick commented 3 years ago

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 commented 3 years ago

Hi,

Thank you for this- it is helpful! I can now see how it is possible to take the 'cyl_only' signal from the multicompartmental and calculate FA, MD etc metrics. Specifically, I am looking to build a multicompartmental model that I can fit the ivim signal via the dstar fixed method but instead the diffusion of the brain tissue is modeled as a tensor rather than scalar. I believe in this case the cylinder compartment should be entered into the intra_voxel_incoherent_motion procedure/script and then the resultant cylinder parameters can be used to calculate DTI metrics. I hope that makes sense!

Best,

Faye

On Fri, Sep 25, 2020 at 10:51 AM Rutger Fick notifications@github.com wrote:

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?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/AthenaEPI/dmipy/issues/103#issuecomment-698974634, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHTA6N7JKGD5OXQN2FNTIHLSHSU6JANCNFSM4RX2EVNA .

fayemckenna commented 3 years ago

Hi,

Sorry one more comment. Why in this example does the tenfit.fa and tenfit.md only output one value rather when cyl_only is entered rather than the whole volume when for example original data is entered?

Thank you!

On Sun, Sep 27, 2020 at 5:49 PM Faye McKenna ffmckenna@gmail.com wrote:

Hi,

Thank you for this- it is helpful! I can now see how it is possible to take the 'cyl_only' signal from the multicompartmental and calculate FA, MD etc metrics. Specifically, I am looking to build a multicompartmental model that I can fit the ivim signal via the dstar fixed method but instead the diffusion of the brain tissue is modeled as a tensor rather than scalar. I believe in this case the cylinder compartment should be entered into the intra_voxel_incoherent_motion procedure/script and then the resultant cylinder parameters can be used to calculate DTI metrics. I hope that makes sense!

Best,

Faye

On Fri, Sep 25, 2020 at 10:51 AM Rutger Fick notifications@github.com wrote:

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?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/AthenaEPI/dmipy/issues/103#issuecomment-698974634, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHTA6N7JKGD5OXQN2FNTIHLSHSU6JANCNFSM4RX2EVNA .

fayemckenna commented 3 years ago

Just following up that there seems to be a bug where entering the mu variable back into the cyl_only signal will only allow 2 inputs (theta and phi) which works for your example code but not with a dataset with many voxels. It looks like if the code called upon unitsphere2cart_Nd(mu) instead of unitsphere2cart_1d(mu) to recreate the signal potentially it wouldn't give the error? Do you have any suggestions for this? Thanks again!

On Mon, Sep 28, 2020 at 10:56 AM Faye McKenna ffmckenna@gmail.com wrote:

Hi,

Sorry one more comment. Why in this example does the tenfit.fa and tenfit.md only output one value rather when cyl_only is entered rather than the whole volume when for example original data is entered?

Thank you!

On Sun, Sep 27, 2020 at 5:49 PM Faye McKenna ffmckenna@gmail.com wrote:

Hi,

Thank you for this- it is helpful! I can now see how it is possible to take the 'cyl_only' signal from the multicompartmental and calculate FA, MD etc metrics. Specifically, I am looking to build a multicompartmental model that I can fit the ivim signal via the dstar fixed method but instead the diffusion of the brain tissue is modeled as a tensor rather than scalar. I believe in this case the cylinder compartment should be entered into the intra_voxel_incoherent_motion procedure/script and then the resultant cylinder parameters can be used to calculate DTI metrics. I hope that makes sense!

Best,

Faye

On Fri, Sep 25, 2020 at 10:51 AM Rutger Fick notifications@github.com wrote:

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?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/AthenaEPI/dmipy/issues/103#issuecomment-698974634, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHTA6N7JKGD5OXQN2FNTIHLSHSU6JANCNFSM4RX2EVNA .

rutgerfick commented 3 years ago

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 commented 3 years ago

Hi,

Thank you for this! It worked well, I just found I needed to remove NaNs from signal_cyl_only before entering it into the DIPY tensor model in case anyone comes across that same issue.

Best,

Faye

On Sun, Oct 4, 2020 at 4:41 PM Rutger Fick notifications@github.com wrote:

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]

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/AthenaEPI/dmipy/issues/103#issuecomment-703313111, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHTA6NZKUFXHHJZ2GOJH4B3SJDMYTANCNFSM4RX2EVNA .