We show how to model the diffusion signal as a linear combination
of continuous functions from the SHORE basis [Merlet2013], [Özarslan2008],
[Özarslan2009]. We also compute the analytical Orientation Distribution
Function (ODF). First import the necessary modules: Download and read the data for this tutorial. The six parameters of these two functions define the ROI where to reconstruct
the data. They respectively correspond to Instantiate the SHORE Model. For details regarding these four parameters see [Cheng2011] and [Merlet2013]. Fit the SHORE model to the data Load an odf reconstruction sphere Compute the ODFs Display the ODFs Özarslan E. et al., “Simple harmonic oscillator based
estimation and reconstruction for one-dimensional q-space MR,” in Proc Intl
Soc Mag Reson Med, 16, Toronto, Canada, 2008 Özarslan E. et al., “Simple harmonic oscillator based
reconstruction and estimation for three-dimensional q-space MRI,” in Proc
Intl Soc Mag Reson Med, 17, Honolulu, HI, 2009 Merlet S. et al., “Continuous diffusion signal, EAP and ODF
estimation via Compressive Sensing in diffusion MRI”, Medical Image
Analysis, 2013. Cheng J. et al., “Theoretical Analysis and Pratical Insights on
EAP Estimation via Unified HARDI Framework”, MICCAI workshop on
Computational Diffusion MRI, 2011. Example source code You can download Continuous and analytical diffusion signal modelling with 3D-SHORE
from dipy.reconst.shore import ShoreModel
from dipy.viz import window, actor
from dipy.core.gradients import gradient_table
from dipy.data import get_fnames, get_sphere
from dipy.io.gradients import read_bvals_bvecs
from dipy.io.image import load_nifti
fetch_isbi2013_2shell()
provides data from the ISBI HARDI contest 2013 acquired for two shells at
b-values 1500 \(s/mm^2\) and 2500 \(s/mm^2\).(xmin,xmax,ymin,ymax,zmin,zmax)
with x, y, z and the three axis defining the spatial positions of the voxels.fraw, fbval, fbvec = get_fnames('isbi2013_2shell')
data, affine = load_nifti(fraw)
bvals, bvecs = read_bvals_bvecs(fbval, fbvec)
gtab = gradient_table(bvals, bvecs)
data_small = data[10:40, 22, 10:40]
print('data.shape (%d, %d, %d, %d)' % data.shape)
data
contains the voxel data and gtab
contains a GradientTable
object (gradient information e.g. b-values). For example, to show the b-values
it is possible to write:print(gtab.bvals)
radial_order
is the radial order of the SHORE basis.zeta
is the scale factor of the SHORE basis.lambdaN
and lambdaL
are the radial and angular regularization
constants, respectively.radial_order = 6
zeta = 700
lambdaN = 1e-8
lambdaL = 1e-8
asm = ShoreModel(gtab, radial_order=radial_order,
zeta=zeta, lambdaN=lambdaN, lambdaL=lambdaL)
asmfit = asm.fit(data_small)
sphere = get_sphere('repulsion724')
odf = asmfit.odf(sphere)
print('odf.shape (%d, %d, %d)' % odf.shape)
# Enables/disables interactive visualization
interactive = False
scene = window.Scene()
sfu = actor.odf_slicer(odf[:, None, :], sphere=sphere, colormap='plasma',
scale=0.5)
sfu.RotateX(-90)
sfu.display(y=0)
scene.add(sfu)
window.record(scene, out_path='odfs.png', size=(600, 600))
if interactive:
window.show(scene)
References
the full source code of this example
. This same script is also included in the dipy source distribution under the doc/examples/
directory.