Note
Go to the end to download the full example code
Reconstruction with the Sparse Fascicle Model (SFM)#
In this example, we will use the Sparse Fascicle Model (SFM) [1], to reconstruct the fiber Orientation Distribution Function (fODF) in every voxel.
First, we import the modules we will use in this example:
from dipy.core.gradients import gradient_table
import dipy.data as dpd
import dipy.direction.peaks as dpp
from dipy.io.gradients import read_bvals_bvecs
from dipy.io.image import load_nifti
from dipy.reconst.csdeconv import auto_response_ssst
import dipy.reconst.sfm as sfm
from dipy.viz import actor, window
For the purpose of this example, we will use the Stanford HARDI dataset (150
directions, single b-value of 2000
hardi_fname, hardi_bval_fname, hardi_bvec_fname = dpd.get_fnames(name="stanford_hardi")
data, affine = load_nifti(hardi_fname)
bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
gtab = gradient_table(bvals, bvecs=bvecs)
# Enables/disables interactive visualization
interactive = False
Reconstruction of the fiber ODF in each voxel guides subsequent tracking steps. Here, the model is the Sparse Fascicle Model, described in [1]. This model reconstructs the diffusion signal as a combination of the signals from different fascicles. This model can be written as:
Where
Sparsity constraints on the fiber ODF (
Elastic Net optimizes the following cost function:
where
Just like Constrained Spherical Deconvolution (see Reconstruction with Constrained Spherical Deconvolution model (CSD)), the
SFM requires the definition of a response function. We’ll take advantage of
the automated algorithm in the csdeconv
module to find this response
function:
response, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)
The response
return value contains two entries. The first is an array
with the eigenvalues of the response function and the second is the average
S0 for this response.
It is a very good practice to always validate the result of
auto_response_ssst
. For, this purpose we can print it and have a look
at its values.
print(response)
(array([0.00139919, 0.0003007 , 0.0003007 ]), 416.7372408293461)
We initialize an SFM model object, using these values. We will use the default sphere (362 vertices, symmetrically distributed on the surface of the sphere), as a set of putative fascicle directions that are considered in the model
sphere = dpd.get_sphere()
sf_model = sfm.SparseFascicleModel(
gtab, sphere=sphere, l1_ratio=0.5, alpha=0.001, response=response[0]
)
For the purpose of the example, we will consider a small volume of data containing parts of the corpus callosum and of the centrum semiovale
data_small = data[20:50, 55:85, 38:39]
Fitting the model to this small volume of data, we calculate the ODF of this model on the sphere, and plot it.
sf_fit = sf_model.fit(data_small)
sf_odf = sf_fit.odf(sphere)
fodf_spheres = actor.odf_slicer(sf_odf, sphere=sphere, scale=0.8, colormap="plasma")
scene = window.Scene()
scene.add(fodf_spheres)
window.record(scene=scene, out_path="sf_odfs.png", size=(1000, 1000))
if interactive:
window.show(scene)

We can extract the peaks from the ODF, and plot these as well
sf_peaks = dpp.peaks_from_model(
sf_model,
data_small,
sphere,
relative_peak_threshold=0.5,
min_separation_angle=25,
return_sh=False,
)
scene.clear()
fodf_peaks = actor.peak_slicer(sf_peaks.peak_dirs, peaks_values=sf_peaks.peak_values)
scene.add(fodf_peaks)
window.record(scene=scene, out_path="sf_peaks.png", size=(1000, 1000))
if interactive:
window.show(scene)

Finally, we plot both the peaks and the ODFs, overlaid:
fodf_spheres.GetProperty().SetOpacity(0.4)
scene.add(fodf_spheres)
window.record(scene=scene, out_path="sf_both.png", size=(1000, 1000))
if interactive:
window.show(scene)

SFM Peaks and ODFs.
To see how to use this information in tracking, proceed to Tracking with the Sparse Fascicle Model.
References#
Total running time of the script: (5 minutes 46.220 seconds)