Parallel Transport Tractography

Parallel Transport Tractography (PTT) [Aydogan2021]

from dipy.direction import peaks_from_model
from dipy.data import default_sphere
from dipy.io.streamline import save_trk
from dipy.io.stateful_tractogram import Space, StatefulTractogram
from dipy.data import get_sphere
from dipy.direction import PTTDirectionGetter
from dipy.reconst.shm import CsaOdfModel
from dipy.core.gradients import gradient_table
from dipy.data import get_fnames
from dipy.io.gradients import read_bvals_bvecs
from dipy.io.image import load_nifti, load_nifti_data
from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,
                                   auto_response_ssst)
from dipy.tracking import utils
from dipy.tracking.local_tracking import LocalTracking
from dipy.tracking.streamline import Streamlines
from dipy.tracking.stopping_criterion import ThresholdStoppingCriterion
from dipy.viz import window, actor, colormap, has_fury


# Enables/disables interactive visualization
interactive = False

hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')
label_fname = get_fnames('stanford_labels')

data, affine, hardi_img = load_nifti(hardi_fname, return_img=True)
labels = load_nifti_data(label_fname)
bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
gtab = gradient_table(bvals, bvecs)

seed_mask = (labels == 2)
white_matter = (labels == 1) | (labels == 2)
seeds = utils.seeds_from_mask(seed_mask, affine, density=1)

response, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)
csd_model = ConstrainedSphericalDeconvModel(gtab, response, sh_order=6)
csd_fit = csd_model.fit(data, mask=white_matter)
  0%|          | 0/58788 [00:00<?, ?it/s]
  1%|1         | 643/58788 [00:00<00:09, 6422.73it/s]
  3%|2         | 1530/58788 [00:00<00:07, 7853.65it/s]
  4%|4         | 2443/58788 [00:00<00:06, 8431.34it/s]
  6%|5         | 3348/58788 [00:00<00:06, 8671.79it/s]
  7%|7         | 4268/58788 [00:00<00:06, 8860.71it/s]
  9%|8         | 5193/58788 [00:00<00:05, 8992.74it/s]
 10%|#         | 6124/58788 [00:00<00:05, 9096.39it/s]
 12%|#2        | 7064/58788 [00:00<00:05, 9191.11it/s]
 14%|#3        | 8012/58788 [00:00<00:05, 9279.94it/s]
 15%|#5        | 8940/58788 [00:01<00:05, 9257.52it/s]
 17%|#6        | 9878/58788 [00:01<00:05, 9294.45it/s]
 18%|#8        | 10823/58788 [00:01<00:05, 9340.29it/s]
 20%|##        | 11758/58788 [00:01<00:05, 9335.06it/s]
 22%|##1       | 12702/58788 [00:01<00:04, 9364.68it/s]
 23%|##3       | 13666/58788 [00:01<00:04, 9445.47it/s]
 25%|##4       | 14627/58788 [00:01<00:04, 9491.27it/s]
 26%|##6       | 15577/58788 [00:01<00:04, 9467.71it/s]
 28%|##8       | 16543/58788 [00:01<00:04, 9522.29it/s]
 30%|##9       | 17502/58788 [00:01<00:04, 9542.36it/s]
 31%|###1      | 18457/58788 [00:02<00:04, 9489.56it/s]
 33%|###3      | 19407/58788 [00:02<00:04, 9452.42it/s]
 35%|###4      | 20353/58788 [00:02<00:04, 9440.66it/s]
 36%|###6      | 21298/58788 [00:02<00:03, 9436.38it/s]
 38%|###7      | 22266/58788 [00:02<00:03, 9507.31it/s]
 39%|###9      | 23217/58788 [00:02<00:03, 9437.38it/s]
 41%|####1     | 24161/58788 [00:02<00:03, 9411.71it/s]
 43%|####2     | 25114/58788 [00:02<00:03, 9446.75it/s]
 44%|####4     | 26059/58788 [00:02<00:03, 9440.33it/s]
 46%|####5     | 27004/58788 [00:02<00:03, 9382.84it/s]
 48%|####7     | 27943/58788 [00:03<00:03, 9349.28it/s]
 49%|####9     | 28879/58788 [00:03<00:03, 9234.36it/s]
 51%|#####     | 29803/58788 [00:03<00:03, 9104.07it/s]
 52%|#####2    | 30714/58788 [00:03<00:03, 9071.47it/s]
 54%|#####3    | 31626/58788 [00:03<00:02, 9085.19it/s]
 55%|#####5    | 32551/58788 [00:03<00:02, 9133.51it/s]
 57%|#####6    | 33504/58788 [00:03<00:02, 9248.82it/s]
 59%|#####8    | 34449/58788 [00:03<00:02, 9308.42it/s]
 60%|######    | 35394/58788 [00:03<00:02, 9350.24it/s]
 62%|######1   | 36350/58788 [00:03<00:02, 9412.81it/s]
 63%|######3   | 37304/58788 [00:04<00:02, 9449.86it/s]
 65%|######5   | 38257/58788 [00:04<00:02, 9472.59it/s]
 67%|######6   | 39225/58788 [00:04<00:02, 9533.16it/s]
 68%|######8   | 40182/58788 [00:04<00:01, 9541.69it/s]
 70%|######9   | 41137/58788 [00:04<00:01, 9532.98it/s]
 72%|#######1  | 42119/58788 [00:04<00:01, 9616.45it/s]
 73%|#######3  | 43081/58788 [00:04<00:01, 9565.64it/s]
 75%|#######4  | 44039/58788 [00:04<00:01, 9568.41it/s]
 77%|#######6  | 44996/58788 [00:04<00:01, 9545.14it/s]
 78%|#######8  | 45951/58788 [00:04<00:01, 9491.55it/s]
 80%|#######9  | 46901/58788 [00:05<00:01, 9449.36it/s]
 81%|########1 | 47847/58788 [00:05<00:01, 9446.39it/s]
 83%|########2 | 48792/58788 [00:05<00:01, 9432.49it/s]
 85%|########4 | 49739/58788 [00:05<00:00, 9441.22it/s]
 86%|########6 | 50708/58788 [00:05<00:00, 9513.17it/s]
 88%|########7 | 51660/58788 [00:05<00:00, 9483.05it/s]
 89%|########9 | 52609/58788 [00:05<00:00, 9338.64it/s]
 91%|#########1| 53544/58788 [00:05<00:00, 9258.99it/s]
 93%|#########2| 54471/58788 [00:05<00:00, 9164.29it/s]
 94%|#########4| 55388/58788 [00:05<00:00, 9145.18it/s]
 96%|#########5| 56303/58788 [00:06<00:00, 9055.36it/s]
 97%|#########7| 57209/58788 [00:06<00:00, 8953.44it/s]
 99%|#########8| 58105/58788 [00:06<00:00, 8835.14it/s]
100%|##########| 58788/58788 [00:06<00:00, 9243.01it/s]

We use the GFA of the CSA model to build a stopping criterion.

csa_model = CsaOdfModel(gtab, sh_order=6)
gfa = csa_model.fit(data, mask=white_matter).gfa
stopping_criterion = ThresholdStoppingCriterion(gfa, .25)

Prepare the PTT direction getter using the fiber ODF (FOD) obtain with CSD. Start the local tractography using PTT direction getter.

sphere = get_sphere(name='repulsion724')
fod = csd_fit.odf(sphere)
pmf = fod.clip(min=0)
ptt_dg = PTTDirectionGetter.from_pmf(pmf, max_angle=15, probe_length=0.5,
                                     sphere=sphere)

# Parallel Transport Tractography
streamline_generator = LocalTracking(direction_getter=ptt_dg,
                                     stopping_criterion=stopping_criterion,
                                     seeds=seeds,
                                     affine=affine,
                                     step_size=0.2)
streamlines = Streamlines(streamline_generator)
sft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)
save_trk(sft, "tractogram_ptt_dg_pmf.trk")

if has_fury:
    scene = window.Scene()
    scene.add(actor.line(streamlines, colormap.line_colors(streamlines)))
    window.record(scene, out_path='tractogram_ptt_dg_pmf.png',
                  size=(800, 800))
    if interactive:
        window.show(scene)
tracking ptt
examples_built/13_fiber_tracking/tractogram_ptt_dg_pmf.png

Corpus Callosum using ptt direction getter from PMF

References

[Aydogan2021]

Aydogan DB, Shi Y. Parallel Transport Tractography. IEEE Trans Med Imaging. 2021 Feb;40(2):635-647. doi: 10.1109/TMI.2020.3034038. Epub 2021 Feb 2. PMID: 33104507; PMCID: PMC7931442.

Total running time of the script: ( 1 minutes 30.201 seconds)

Gallery generated by Sphinx-Gallery