Note
Click here to download the full example code
Fiber to bundle coherence measures
This demo presents the fiber to bundle coherence (FBC) quantitative measure of the alignment of each fiber with the surrounding fiber bundles [Meesters2016]. These measures are useful in ‘cleaning’ the results of tractography algorithms, since low FBCs indicate which fibers are isolated and poorly aligned with their neighbors, as shown in the figure below.
On the left this figure illustrates (in 2D) the contribution of two fiber points to the kernel density estimator. The kernel density estimator is the sum over all such locally aligned kernels. The local fiber to bundle coherence, shown on the right, color-coded for each fiber, is obtained by evaluating the kernel density estimator along the fibers. One spurious fiber is present which is isolated and badly aligned with the other fibers, and can be identified by a low LFBC value in the region where it deviates from the bundle. Figure adapted from [Portegies2015].
Here we implement FBC measures based on kernel density estimation in the non-flat 5D position-orientation domain. First we compute the kernel density estimator induced by the full lifted output (defined in the space of positions and orientations) of the tractography. Then, the Local FBC (LFBC) is the result of evaluating the estimator along each element of the lifted fiber. A whole fiber measure, the relative FBC (RFBC), is calculated by the minimum of the moving average LFBC along the fiber. Details of the computation of FBC can be found in [Portegies2015].
The FBC measures are evaluated on the Stanford HARDI dataset (150 orientations, b=2000 \(s/mm^2\)) which is one of the standard example datasets in DIPY_.
# Enables/disables interactive visualization
interactive = False
import numpy as np
from dipy.core.gradients import gradient_table
from dipy.data import get_fnames
from dipy.io.image import load_nifti_data, load_nifti
from dipy.io.gradients import read_bvals_bvecs
# Fix seed
np.random.seed(1)
# Read data
hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')
label_fname = get_fnames('stanford_labels')
t1_fname = get_fnames('stanford_t1')
data, affine = load_nifti(hardi_fname)
labels = load_nifti_data(label_fname)
t1_data = load_nifti_data(t1_fname)
bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
gtab = gradient_table(bvals, bvecs)
# Select a relevant part of the data (left hemisphere)
# Coordinates given in x bounds, y bounds, z bounds
dshape = data.shape[:-1]
xa, xb, ya, yb, za, zb = [15, 42, 10, 65, 18, 65]
data_small = data[xa:xb, ya:yb, za:zb]
selectionmask = np.zeros(dshape, 'bool')
selectionmask[xa:xb, ya:yb, za:zb] = True
The data is first fitted to the Constant Solid Angle (CDA) ODF Model. CSA is a good choice to estimate general fractional anisotropy (GFA), which the stopping criterion can use to restrict fiber tracking to those areas where the ODF shows significant restricted diffusion, thus creating a region-of-interest in which the computations are done.
# Perform CSA
from dipy.reconst.shm import CsaOdfModel
from dipy.data import default_sphere
from dipy.direction import peaks_from_model
csa_model = CsaOdfModel(gtab, sh_order=6)
csa_peaks = peaks_from_model(csa_model, data, default_sphere,
relative_peak_threshold=.6,
min_separation_angle=45,
mask=selectionmask)
# Stopping Criterion
from dipy.tracking.stopping_criterion import ThresholdStoppingCriterion
stopping_criterion = ThresholdStoppingCriterion(csa_peaks.gfa, 0.25)
In order to perform probabilistic fiber tracking we first fit the data to the Constrained Spherical Deconvolution (CSD) model in DIPY. This model represents each voxel in the data set as a collection of small white matter fibers with different orientations. The density of fibers along each orientation is known as the Fiber Orientation Distribution (FOD), used in the fiber tracking.
# Perform CSD on the original data
from dipy.reconst.csdeconv import auto_response_ssst
from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel
response, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)
csd_model = ConstrainedSphericalDeconvModel(gtab, response)
csd_fit = csd_model.fit(data_small)
csd_fit_shm = np.lib.pad(csd_fit.shm_coeff, ((xa, dshape[0]-xb),
(ya, dshape[1]-yb),
(za, dshape[2]-zb),
(0, 0)), 'constant')
# Probabilistic direction getting for fiber tracking
from dipy.direction import ProbabilisticDirectionGetter
prob_dg = ProbabilisticDirectionGetter.from_shcoeff(csd_fit_shm,
max_angle=30.,
sphere=default_sphere)
0%| | 0/69795 [00:00<?, ?it/s]
1%| | 479/69795 [00:00<00:14, 4779.54it/s]
1%|1 | 957/69795 [00:00<00:15, 4544.46it/s]
2%|2 | 1413/69795 [00:00<00:15, 4438.25it/s]
3%|2 | 1858/69795 [00:00<00:15, 4355.78it/s]
3%|3 | 2301/69795 [00:00<00:15, 4379.30it/s]
4%|3 | 2762/69795 [00:00<00:15, 4453.91it/s]
5%|4 | 3208/69795 [00:00<00:14, 4451.16it/s]
5%|5 | 3654/69795 [00:00<00:15, 4383.20it/s]
6%|5 | 4100/69795 [00:00<00:14, 4405.99it/s]
7%|6 | 4541/69795 [00:01<00:14, 4381.12it/s]
7%|7 | 4985/69795 [00:01<00:14, 4397.20it/s]
8%|7 | 5459/69795 [00:01<00:14, 4499.11it/s]
8%|8 | 5910/69795 [00:01<00:14, 4456.73it/s]
9%|9 | 6356/69795 [00:01<00:14, 4427.16it/s]
10%|9 | 6799/69795 [00:01<00:14, 4425.70it/s]
10%|# | 7242/69795 [00:01<00:14, 4420.22it/s]
11%|#1 | 7686/69795 [00:01<00:14, 4421.94it/s]
12%|#1 | 8160/69795 [00:01<00:13, 4516.49it/s]
12%|#2 | 8612/69795 [00:01<00:13, 4436.24it/s]
13%|#2 | 9056/69795 [00:02<00:13, 4432.06it/s]
14%|#3 | 9503/69795 [00:02<00:13, 4441.03it/s]
14%|#4 | 9949/69795 [00:02<00:13, 4446.03it/s]
15%|#4 | 10400/69795 [00:02<00:13, 4461.89it/s]
16%|#5 | 10859/69795 [00:02<00:13, 4499.72it/s]
16%|#6 | 11310/69795 [00:02<00:13, 4429.07it/s]
17%|#6 | 11754/69795 [00:02<00:13, 4431.28it/s]
17%|#7 | 12202/69795 [00:02<00:12, 4443.18it/s]
18%|#8 | 12650/69795 [00:02<00:12, 4453.75it/s]
19%|#8 | 13114/69795 [00:02<00:12, 4507.13it/s]
19%|#9 | 13565/69795 [00:03<00:12, 4488.57it/s]
20%|## | 14014/69795 [00:03<00:12, 4455.55it/s]
21%|## | 14472/69795 [00:03<00:12, 4489.15it/s]
21%|##1 | 14922/69795 [00:03<00:12, 4463.63it/s]
22%|##2 | 15369/69795 [00:03<00:12, 4460.34it/s]
23%|##2 | 15838/69795 [00:03<00:11, 4527.96it/s]
23%|##3 | 16291/69795 [00:03<00:11, 4489.66it/s]
24%|##3 | 16745/69795 [00:03<00:11, 4501.63it/s]
25%|##4 | 17204/69795 [00:03<00:11, 4526.26it/s]
25%|##5 | 17657/69795 [00:03<00:11, 4490.05it/s]
26%|##5 | 18108/69795 [00:04<00:11, 4493.27it/s]
27%|##6 | 18564/69795 [00:04<00:11, 4512.97it/s]
27%|##7 | 19017/69795 [00:04<00:11, 4515.91it/s]
28%|##7 | 19469/69795 [00:04<00:11, 4514.50it/s]
29%|##8 | 19932/69795 [00:04<00:10, 4547.83it/s]
29%|##9 | 20387/69795 [00:04<00:10, 4534.53it/s]
30%|##9 | 20850/69795 [00:04<00:10, 4561.21it/s]
31%|### | 21307/69795 [00:04<00:10, 4545.02it/s]
31%|###1 | 21762/69795 [00:04<00:10, 4514.27it/s]
32%|###1 | 22217/69795 [00:04<00:10, 4523.82it/s]
32%|###2 | 22670/69795 [00:05<00:10, 4516.97it/s]
33%|###3 | 23122/69795 [00:05<00:10, 4486.38it/s]
34%|###3 | 23579/69795 [00:05<00:10, 4509.02it/s]
34%|###4 | 24030/69795 [00:05<00:10, 4457.67it/s]
35%|###5 | 24480/69795 [00:05<00:10, 4468.94it/s]
36%|###5 | 24943/69795 [00:05<00:09, 4515.10it/s]
36%|###6 | 25395/69795 [00:05<00:09, 4502.81it/s]
37%|###7 | 25862/69795 [00:05<00:09, 4551.95it/s]
38%|###7 | 26331/69795 [00:05<00:09, 4592.70it/s]
38%|###8 | 26791/69795 [00:05<00:09, 4553.02it/s]
39%|###9 | 27249/69795 [00:06<00:09, 4560.46it/s]
40%|###9 | 27714/69795 [00:06<00:09, 4584.53it/s]
40%|#### | 28186/69795 [00:06<00:09, 4623.17it/s]
41%|####1 | 28661/69795 [00:06<00:08, 4657.98it/s]
42%|####1 | 29127/69795 [00:06<00:08, 4603.62it/s]
42%|####2 | 29588/69795 [00:06<00:08, 4523.38it/s]
43%|####3 | 30057/69795 [00:06<00:08, 4570.41it/s]
44%|####3 | 30526/69795 [00:06<00:08, 4605.08it/s]
44%|####4 | 30992/69795 [00:06<00:08, 4620.47it/s]
45%|####5 | 31455/69795 [00:06<00:08, 4619.91it/s]
46%|####5 | 31918/69795 [00:07<00:08, 4602.56it/s]
46%|####6 | 32379/69795 [00:07<00:08, 4580.94it/s]
47%|####7 | 32861/69795 [00:07<00:07, 4651.66it/s]
48%|####7 | 33342/69795 [00:07<00:07, 4697.58it/s]
48%|####8 | 33812/69795 [00:07<00:07, 4696.62it/s]
49%|####9 | 34282/69795 [00:07<00:07, 4681.61it/s]
50%|####9 | 34751/69795 [00:07<00:07, 4681.01it/s]
50%|##### | 35239/69795 [00:07<00:07, 4739.21it/s]
51%|#####1 | 35734/69795 [00:07<00:07, 4800.90it/s]
52%|#####1 | 36215/69795 [00:08<00:07, 4721.10it/s]
53%|#####2 | 36698/69795 [00:08<00:06, 4752.59it/s]
53%|#####3 | 37174/69795 [00:08<00:06, 4682.81it/s]
54%|#####3 | 37658/69795 [00:08<00:06, 4727.30it/s]
55%|#####4 | 38142/69795 [00:08<00:06, 4760.36it/s]
55%|#####5 | 38619/69795 [00:08<00:06, 4703.45it/s]
56%|#####6 | 39090/69795 [00:08<00:06, 4701.09it/s]
57%|#####6 | 39561/69795 [00:08<00:06, 4656.27it/s]
57%|#####7 | 40027/69795 [00:08<00:06, 4600.45it/s]
58%|#####8 | 40507/69795 [00:08<00:06, 4656.71it/s]
59%|#####8 | 40973/69795 [00:09<00:06, 4632.74it/s]
59%|#####9 | 41437/69795 [00:09<00:06, 4575.76it/s]
60%|###### | 41907/69795 [00:09<00:06, 4611.43it/s]
61%|###### | 42369/69795 [00:09<00:06, 4530.03it/s]
61%|######1 | 42839/69795 [00:09<00:05, 4578.86it/s]
62%|######2 | 43315/69795 [00:09<00:05, 4631.59it/s]
63%|######2 | 43799/69795 [00:09<00:05, 4691.79it/s]
63%|######3 | 44276/69795 [00:09<00:05, 4714.57it/s]
64%|######4 | 44748/69795 [00:09<00:05, 4626.76it/s]
65%|######4 | 45212/69795 [00:09<00:05, 4596.32it/s]
65%|######5 | 45701/69795 [00:10<00:05, 4681.58it/s]
66%|######6 | 46183/69795 [00:10<00:05, 4722.05it/s]
67%|######6 | 46666/69795 [00:10<00:04, 4753.73it/s]
68%|######7 | 47142/69795 [00:10<00:04, 4700.70it/s]
68%|######8 | 47613/69795 [00:10<00:04, 4593.53it/s]
69%|######8 | 48091/69795 [00:10<00:04, 4646.30it/s]
70%|######9 | 48559/69795 [00:10<00:04, 4655.09it/s]
70%|####### | 49046/69795 [00:10<00:04, 4717.78it/s]
71%|####### | 49521/69795 [00:10<00:04, 4724.61it/s]
72%|#######1 | 49994/69795 [00:10<00:04, 4605.09it/s]
72%|#######2 | 50456/69795 [00:11<00:04, 4541.21it/s]
73%|#######2 | 50923/69795 [00:11<00:04, 4578.47it/s]
74%|#######3 | 51384/69795 [00:11<00:04, 4586.32it/s]
74%|#######4 | 51863/69795 [00:11<00:03, 4645.48it/s]
75%|#######4 | 52328/69795 [00:11<00:03, 4603.94it/s]
76%|#######5 | 52789/69795 [00:11<00:03, 4559.13it/s]
76%|#######6 | 53251/69795 [00:11<00:03, 4574.03it/s]
77%|#######6 | 53713/69795 [00:11<00:03, 4586.02it/s]
78%|#######7 | 54179/69795 [00:11<00:03, 4607.94it/s]
78%|#######8 | 54641/69795 [00:11<00:03, 4608.31it/s]
79%|#######8 | 55102/69795 [00:12<00:03, 4488.42it/s]
80%|#######9 | 55552/69795 [00:12<00:03, 4380.01it/s]
80%|######## | 55991/69795 [00:12<00:03, 4371.12it/s]
81%|######## | 56429/69795 [00:12<00:03, 4361.47it/s]
81%|########1 | 56871/69795 [00:12<00:02, 4377.77it/s]
82%|########2 | 57314/69795 [00:12<00:02, 4392.81it/s]
83%|########2 | 57754/69795 [00:12<00:02, 4281.68it/s]
83%|########3 | 58183/69795 [00:12<00:02, 4252.39it/s]
84%|########4 | 58643/69795 [00:12<00:02, 4352.58it/s]
85%|########4 | 59095/69795 [00:13<00:02, 4401.27it/s]
85%|########5 | 59562/69795 [00:13<00:02, 4480.86it/s]
86%|########6 | 60029/69795 [00:13<00:02, 4536.92it/s]
87%|########6 | 60484/69795 [00:13<00:02, 4436.46it/s]
87%|########7 | 60929/69795 [00:13<00:02, 4404.33it/s]
88%|########7 | 61381/69795 [00:13<00:01, 4437.43it/s]
89%|########8 | 61837/69795 [00:13<00:01, 4470.49it/s]
89%|########9 | 62309/69795 [00:13<00:01, 4543.48it/s]
90%|########9 | 62769/69795 [00:13<00:01, 4558.72it/s]
91%|######### | 63226/69795 [00:13<00:01, 4458.54it/s]
91%|#########1| 63673/69795 [00:14<00:01, 4418.26it/s]
92%|#########1| 64116/69795 [00:14<00:01, 4411.47it/s]
92%|#########2| 64558/69795 [00:14<00:01, 4376.26it/s]
93%|#########3| 65026/69795 [00:14<00:01, 4465.38it/s]
94%|#########3| 65473/69795 [00:14<00:00, 4363.32it/s]
94%|#########4| 65914/69795 [00:14<00:00, 4374.26it/s]
95%|#########5| 66368/69795 [00:14<00:00, 4421.08it/s]
96%|#########5| 66830/69795 [00:14<00:00, 4479.13it/s]
96%|#########6| 67295/69795 [00:14<00:00, 4529.63it/s]
97%|#########7| 67749/69795 [00:14<00:00, 4496.13it/s]
98%|#########7| 68199/69795 [00:15<00:00, 4376.79it/s]
98%|#########8| 68638/69795 [00:15<00:00, 4353.54it/s]
99%|#########8| 69079/69795 [00:15<00:00, 4368.07it/s]
100%|#########9| 69529/69795 [00:15<00:00, 4405.29it/s]
100%|##########| 69795/69795 [00:15<00:00, 4528.27it/s]
The optic radiation is reconstructed by tracking fibers from the calcarine sulcus (visual cortex V1) to the lateral geniculate nucleus (LGN). We seed from the calcarine sulcus by selecting a region-of-interest (ROI) cube of dimensions 3x3x3 voxels.
# Set a seed region region for tractography.
from dipy.tracking import utils
mask = np.zeros(data.shape[:-1], 'bool')
rad = 3
mask[26-rad:26+rad, 29-rad:29+rad, 31-rad:31+rad] = True
seeds = utils.seeds_from_mask(mask, affine, density=[4, 4, 4])
Local Tracking is used for probabilistic tractography which takes the direction getter along with the stopping criterion and seeds as input.
# Perform tracking using Local Tracking
from dipy.tracking.local_tracking import LocalTracking
streamlines_generator = LocalTracking(prob_dg, stopping_criterion, seeds,
affine, step_size=.5)
# Compute streamlines.
from dipy.tracking.streamline import Streamlines
streamlines = Streamlines(streamlines_generator)
In order to select only the fibers that enter into the LGN, another ROI is created from a cube of size 5x5x5 voxels. The near_roi command is used to find the fibers that traverse through this ROI.
# Set a mask for the lateral geniculate nucleus (LGN)
mask_lgn = np.zeros(data.shape[:-1], 'bool')
rad = 5
mask_lgn[35-rad:35+rad, 42-rad:42+rad, 28-rad:28+rad] = True
# Select all the fibers that enter the LGN and discard all others
filtered_fibers2 = utils.near_roi(streamlines, affine, mask_lgn, tol=1.8)
sfil = []
for i in range(len(streamlines)):
if filtered_fibers2[i]:
sfil.append(streamlines[i])
streamlines = Streamlines(sfil)
Inspired by [Rodrigues2010], a lookup-table is created, containing rotated versions of the fiber propagation kernel \(P_t\) [DuitsAndFranken2011] rotated over a discrete set of orientations. See the Contextual enhancement example for more details regarding the kernel. In order to ensure rotationally invariant processing, the discrete orientations are required to be equally distributed over a sphere. By default, a sphere with 100 directions is obtained from electrostatic repulsion in DIPY.
# Compute lookup table
from dipy.denoise.enhancement_kernel import EnhancementKernel
D33 = 1.0
D44 = 0.02
t = 1
k = EnhancementKernel(D33, D44, t)
The FBC measures are now computed, taking the tractography results and the lookup tables as input.
# Apply FBC measures
from dipy.tracking.fbcmeasures import FBCMeasures
fbc = FBCMeasures(streamlines, k)
After calculating the FBC measures, a threshold can be chosen on the relative FBC (RFBC) in order to remove spurious fibers. Recall that the relative FBC (RFBC) is calculated by the minimum of the moving average LFBC along the fiber. In this example we show the results for threshold 0 (i.e. all fibers are included) and 0.2 (removing the 20 percent most spurious fibers).
# Calculate LFBC for original fibers
fbc_sl_orig, clrs_orig, rfbc_orig = \
fbc.get_points_rfbc_thresholded(0, emphasis=0.01)
# Apply a threshold on the RFBC to remove spurious fibers
fbc_sl_thres, clrs_thres, rfbc_thres = \
fbc.get_points_rfbc_thresholded(0.125, emphasis=0.01)
The results of FBC measures are visualized, showing the original fibers colored by LFBC (see The optic radiation obtained through probabilistic tractography colored by local fiber to bundle coherence.), and the fibers after the cleaning procedure via RFBC thresholding (see The tractography result is cleaned (shown in bottom) by removing fibers with a relative FBC (RFBC) lower than the threshold \tau = 0.2.).
# Visualize the results
from dipy.viz import window, actor
# Create scene
scene = window.Scene()
# Original lines colored by LFBC
lineactor = actor.line(fbc_sl_orig, np.vstack(clrs_orig), linewidth=0.2)
scene.add(lineactor)
# Horizontal (axial) slice of T1 data
vol_actor1 = actor.slicer(t1_data, affine=affine)
vol_actor1.display(z=20)
scene.add(vol_actor1)
# Vertical (sagittal) slice of T1 data
vol_actor2 = actor.slicer(t1_data, affine=affine)
vol_actor2.display(x=35)
scene.add(vol_actor2)
# Show original fibers
scene.set_camera(position=(-264, 285, 155),
focal_point=(0, -14, 9),
view_up=(0, 0, 1))
window.record(scene, n_frames=1, out_path='OR_before.png', size=(900, 900))
if interactive:
window.show(scene)
# Show thresholded fibers
scene.rm(lineactor)
scene.add(actor.line(fbc_sl_thres, np.vstack(clrs_thres), linewidth=0.2))
window.record(scene, n_frames=1, out_path='OR_after.png', size=(900, 900))
if interactive:
window.show(scene)
The optic radiation obtained through probabilistic tractography colored by local fiber to bundle coherence.
The tractography result is cleaned (shown in bottom) by removing fibers with a relative FBC (RFBC) lower than the threshold \(\tau = 0.2\).
Acknowledgments
The techniques are developed in close collaboration with Pauly Ossenblok of the Academic Center of Epileptology Kempenhaeghe & Maastricht UMC+.
References
S. Meesters, G. Sanguinetti, E. Garyfallidis, J. Portegies, P. Ossenblok, R. Duits. (2016) Cleaning output of tractography via fiber to bundle coherence, a new open source implementation. Human Brain Mapping Conference 2016.
J. Portegies, R. Fick, G. Sanguinetti, S. Meesters, G.Girard, and R. Duits. (2015) Improving Fiber Alignment in HARDI by Combining Contextual PDE flow with Constrained Spherical Deconvolution. PLoS One.
R. Duits and E. Franken (2011) Left-invariant diffusions on the space of positions and orientations and their application to crossing-preserving smoothing of HARDI images. International Journal of Computer Vision, 92:231-264.
P. Rodrigues, R. Duits, B. Romeny, A. Vilanova (2010). Accelerated Diffusion Operators for Enhancing DW-MRI. Eurographics Workshop on Visual Computing for Biology and Medicine. The Eurographics Association.
Total running time of the script: ( 3 minutes 10.237 seconds)