from dipy.viz import window, actor
// Enables/disables interactive visualization
interactive = False
ren = window.Renderer()
evals = response[0]
evecs = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]).T
from dipy.data import get_sphere
sphere = get_sphere("symmetric724")
from dipy.sims.voxel import single_tensor_odf
response_odf = single_tensor_odf(sphere.vertices, evals, evecs)
// transform our data from 1D to 4D
response_odf = response_odf[None, None, None, :]
response_actor = actor.odf_slicer(response_odf, sphere=sphere, colormap="jet")
ren.add(response_actor)
print("Saving illustration as csd_response.png")
window.record(ren, out_path="csd_response.png", size=(200, 200))
if interactive:
window.show(ren)
.. figure:: csd_response.png
:align: center
Estimated response function.
ren.rm(response_actor)
Depending on the dataset, FA threshold may not be the best way to find the
best possible response function. For one, it depends on the diffusion tensor
(FA and first eigenvector), which has lower accuracy at high
b-values. Alternatively, the response function can be calibrated in a
data-driven manner [Tax2014]_.
First, the data is deconvolved with a "fat" response function. All voxels that
are considered to contain only one peak in this deconvolution (as determined by
the peak threshold which gives an upper limit of the ratio of the second peak
to the first peak) are maintained, and from these voxels a new response
function is determined. This process is repeated until convergence is
reached. Here we calibrate the response function on a small part of the data.
from dipy.reconst.csdeconv import recursive_response
A WM mask can shorten computation time for the whole dataset. Here it is created
based on the DTI fit.
import dipy.reconst.dti as dti
tenmodel = dti.TensorModel(gtab)
tenfit = tenmodel.fit(data, mask=data[..., 0] > 200)
from dipy.reconst.dti import fractional_anisotropy
FA = fractional_anisotropy(tenfit.evals)
MD = dti.mean_diffusivity(tenfit.evals)
wm_mask = (np.logical_or(FA >= 0.4, (np.logical_and(FA >= 0.15, MD >= 0.0011))))
response = recursive_response(gtab, data, mask=wm_mask, sh_order=8,
peak_thr=0.01, init_fa=0.08,
init_trace=0.0021, iter=8, convergence=0.001,
parallel=True)
We can check the shape of the signal of the response function, which should be
like a pancake:
response_signal = response.on_sphere(sphere)
response_actor = actor.odf_slicer(response_signal, sphere=sphere, colormap="jet")
ren = window.Renderer()
ren.add(response_actor)
print("Saving illustration as csd_recursive_response.png")
window.record(ren, out_path="csd_recursive_response.png", size=(200, 200))
if interactive:
window.show(ren)
.. figure:: csd_recursive_response.png
:align: center
Estimated response function using recursive calibration.
ren.rm(response_actor)
Now, that we have the response function, we are ready to start the deconvolution
process. Let"s import the CSD model and fit the datasets.
from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel
csd_model = ConstrainedSphericalDeconvModel(gtab, response)
For illustration purposes we will fit only a small portion of the data.
data_small = data[20:50, 55:85, 38:39]
csd_fit = csd_model.fit(data_small)
Show the CSD-based ODFs also known as FODFs (fiber ODFs).
csd_odf = csd_fit.odf(sphere)
Here we visualize only a 30x30 region.
fodf_spheres = actor.odf_slicer(csd_odf, sphere=sphere, scale=1.3, norm=False, colormap="jet")
ren.add(fodf_spheres)
print("Saving illustration as csd_odfs.png")
window.record(ren, out_path="csd_odfs.png", size=(600, 600))
if interactive:
window.show(ren)
.. figure:: csd_odfs.png
:align: center
CSD ODFs.
In Dipy we also provide tools for finding the peak directions (maxima) of the
ODFs. For this purpose we recommend using ``peaks_from_model``.
from dipy.direction import peaks_from_model
csd_peaks = peaks_from_model(model=csd_model,
data=data_small,
sphere=sphere,
relative_peak_threshold=.5,
min_separation_angle=25,
parallel=True)
window.clear(ren)
fodf_peaks = actor.peak_slicer(csd_peaks.peak_dirs, csd_peaks.peak_values, scale=1.3)
ren.add(fodf_peaks)
print("Saving illustration as csd_peaks.png")
window.record(ren, out_path="csd_peaks.png", size=(600, 600))
if interactive:
window.show(ren)
.. figure:: csd_peaks.png
:align: center
CSD Peaks.
We can finally visualize both the ODFs and peaks in the same space.
fodf_spheres.GetProperty().SetOpacity(0.4)
ren.add(fodf_spheres)
print("Saving illustration as csd_both.png")
window.record(ren, out_path="csd_both.png", size=(600, 600))
if interactive:
window.show(ren)
.. figure:: csd_both.png
:align: center
