[Neuroimaging] [Dipy] discrepancies of SH coeffs with MRtrix

Tue Oct 25 11:42:33 EDT 2016

Regarding that last comment, it is unfortunately unrealistic to get
*exactly* the same result form two different software in general, but for
this specific case there is a difference in the implementation from mrtrix
and dipy (and everything else they use under the hood of course).

For the recursive calibration computation, the one in mrtrix is differently
implemented from Chantal Tax's original paper (they mention the differences
and why here [1] in the relevant section) and the one in dipy was made by
Chantal herself, so it should fit the theoretical description and the
implementation used in explore dti (which I assume is also made by Chantal).

[1]
http://mrtrix.readthedocs.io/en/latest/concepts/response_function_estimation.html

2016-10-25 11:28 GMT+02:00 Eric Moulton <eric.moulton.jr at gmail.com>:

> Hi everyone,
>
> Thanks to all those who answered me. It cleared a lot of stuff up for me.
> I spent some time playing with these functions to get the hang of it.
> Indeed, I performed the ConstrainedSphericalDeconvModel model and it does
> give me a spatially varying Y(0,0) volume, just like in MRtrix.
>
> That being said, I tried computing the SH_coeff volumes using similar
> algorithms in MRtrix and Dipy, and although I get close results, I do not
> get the same. From my understanding, Dipy seems to employ the Tax algorithm
> for calculating the response function, so that's what I tried to do in
> MRtrix. I started out with a 4D Dwi file (eddy corrected, etc.), a brain
> mask, and the bvecs and bvals.
>
> The parts of my code in bold are the parts where I tried to put in the
> exact same parameters for each program. The only extra options that could
> be responsible for the differences are in Dipy, where I have to give an
> initial FA and MD value.
>
> *************************
> *MRtrix Code*
> mrconvert dwi.nii.gz dwi.mif
> mrconvert dwi_mask.nii.gz dwi_mask.mif
> dwi2response tax *-peak_ratio 0.01 -max_iters 8 -convergence 0.001 -lmax
> 4 -mask dwi_mask.mif* -fslgrad subj.bvecs subj.bvals -force dwi.mif
> response_tax.txt
> dwi2fod -fslgrad subj.bvecs subj.bvals -lmax 4 -mask mask.mif -force csd
> dwi.mif response_tax.txt odf_tax.mif
>
> *Dipy Code*
> import numpy as np
> import nibabel as nib
> import os
> from dipy.data import get_sphere
> from dipy.reconst.shm import CsaOdfModel,
> from dipy.direction import peaks_from_model
> from dipy.tracking import utils
> from dipy.reconst.csdeconv import recursive_response
> from dipy.io import read_bvals_bvecs
> from dipy.core.gradients import gradient_table
> from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel
>
> base_rep='/path/to/files'
> fimg = os.path.join(base_rep,'dwi.nii.gz')
> fbvecs = os.path.join(base_rep,'subj.bvecs')
> fbvals = os.path.join(base_rep,'subj.bvals')
> img = nib.load(fimg)
> dwi = img.get_data()
> fdwi_mask = os.path.join(base_rep,'dwi_mask.nii.gz')
> dwi_mask = nib.load(fdwi_mask).get_data()
> dwi_mask = dwi_mask > 0
> affine = img.affine
> header = img.header
> from dipy.io import read_bvals_bvecs
> bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
> from dipy.core.gradients import gradient_table
> gtab = gradient_table(bvals, bvecs)
> sphere = get_sphere('symmetric724')
>
> response = recursive_response(gtab,dwi,*mask=dwi_mask*,sh_order=4,
>                               *peak_thr=0.01,iter=8,convergence = 0.001,*
>                               init_fa=0.05,init_trace=0.0021,
>                               parallel=True,sphere=sphere)
>
> csdmodel = ConstrainedSphericalDeconvModel(gtab,response,sh_order=4)
> SH_coeff = csdmodel.fit(dwi,mask=dwi_mask).shm_coeff
>
> #some trickery to put the coefficients in the same order as MRtrix
> SH_coeff_mrtrix = np.zeros(SH_coeff.shape)
> SH_coeff_mrtrix[...,0] = SH_coeff[...,0]
> SH_coeff_mrtrix[...,1:6] = SH_coeff[...,5:0:-1]
> SH_coeff_mrtrix[...,6:15] = SH_coeff[...,14:5:-1]
>
> nib.save(nib.Nifti1Image(SH_coeff_mrtrix,affine,header),os.
> path.join(base_rep,'odf_tax_dipy.nii.gz'))
> *************************
>
>
> When I open odf_tax.mif and odf_tax_dipy.nii.gz in mrview, I don't get the
> exact same thing, but both volumes are fairly similar upon visual
> inspection of each component and also the ODFs (although in areas of
> crossing fibers, they can sometimes be more or less "sharper", but
> sometimes MRtrix OR Dipy produces the sharper one depending on the
> location).
>
> Attached is a histogram of the differences of the images (MRtrix - Dipy),
> masked by the brain_mask so we don't get any 0 voxels outside the brain.
> Each component is the Y(l,m) volumes in MRtrix3 basis (i.e. 0th = Y(0,0),
> 1st = Y(2,-2), 2nd = Y(2,-1), etc.).
>
> I was wondering what you guys thought about my code to generate the
> spherical deconvolutions and if these differences are worrying. Also, is it
> unrealistic to hope to get the exact same result using two different
> softwares with lots of behind the scenes calculations?
>
> Best regards,
> Eric
>
>
>
> _______________________________________________
> Neuroimaging mailing list
> Neuroimaging at python.org
> https://mail.python.org/mailman/listinfo/neuroimaging
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/neuroimaging/attachments/20161025/1e1d6038/attachment.html>