🚀 Get started with with Tractography.jl

• 🚀 Get started with with Tractography.jl
• Basic use
• • Step 1: Define a Model
  • Step 2: Define the seeds
  • Step 3: Choose a sampling algorithm
  • Step 4: Sample the streamlines
• Optimal use

This tutorial will introduce you to the functionalities for computing streamlines.

Basic use

In this example, we sample Nmc streamlines from a Model model.

using Tractography
const TG = Tractography

model = Model(Δt = 0.125f0,
            foddata = FODData((@__DIR__) * "/../../examples/fod-FC.nii.gz"),
            )
Nmc = 10
seeds = rand(Float32, 6, Nmc)
alg = Probabilistic()
streamlines, tract_length = sample(model, alg, seeds);
size(streamlines)

(3, 1000, 10)

Step 1: Define a Model

We define a Tractography Markov Chain (Model) model as follows:

model = Model(Δt = 0.125f0,
            foddata = FODData((@__DIR__) * "/../../examples/fod-FC.nii.gz"),
            )

Model with elype Float32
 ├─ Δt = 0.125
 ├─ minimal probability     = 0.0
 ├─ cone                    = Cone{Float32}(90.0f0)
 ├─ mollifier               = max_mollifier
 ├─ evaluation of the basis = PreComputeAllFOD()
 └─ data : (lmax = 8)

Step 2: Define the seeds

seeds is a 6 x Nmc matrix. Each column stores:

• seeds[1:3, i]: initial position (x, y, z) in native space
• seeds[4:6, i]: initial direction (u1, u2, u3)

Nmc = 10 # Monte Carlo sample
seeds = rand(Float32, 6, Nmc)

6×10 Matrix{Float32}:
 0.958727   0.732877   0.972953   0.106595   …  0.647098  0.64341   0.486739
 0.284568   0.794337   0.0587938  0.834071      0.912675  0.149684  0.780145
 0.813771   0.54531    0.407041   0.779897      0.184815  0.890609  0.559893
 0.979626   0.14715    0.739527   0.889523      0.394945  0.421846  0.385374
 0.59004    0.0707152  0.567951   0.850859      0.876543  0.155006  0.318012
 0.0167293  0.88054    0.152792   0.0702036  …  0.268925  0.153168  0.402905

Step 3: Choose a sampling algorithm

alg = Probabilistic()

Probabilistic()

Step 4: Sample the streamlines

streamlines, tract_length = sample(model, alg, seeds);

(Float32[0.9587271 1.0367229 … 1.0760609 1.0760609; 0.28456777 0.35634938 … -1.5205965 -1.5205965; 0.8137715 0.74752146 … 1.6512716 1.6512716;;; 0.7328775 0.6464846 … -1.5124451 -1.5124451; 0.7943368 0.7637382 … -1.3866397 -1.3866397; 0.5453104 0.63031036 … 3.067185 3.067185;;; 0.9729532 0.9857268 … -0.6417298 -0.6417298; 0.058793843 0.11557364 … 2.9777167 2.9777167; 0.40704125 0.5176662 … -1.541709 -1.541709;;; 0.10659522 0.1498922 … -1.5378102 -1.5378102; 0.834071 0.8435897 … -1.2157539 -1.2157539; 0.7798966 0.8967716 … 3.6611466 3.6611466;;; 0.78961754 0.8703191 … 3.6800435 3.6800435; 0.07483721 -0.0086684 … -1.5166051 -1.5166051; 0.31199408 0.35824406 … 0.3351197 0.3351197;;; 0.497927 0.605739 … 0.13399333 0.13399333; 0.5730119 0.6244104 … -1.5001765 -1.5001765; 0.9473256 0.9104506 … 0.96732557 0.96732557;;; 0.091783345 0.19959533 … 0.52727395 0.52727395; 0.80801666 0.8594152 … 3.4576056 3.4576056; 0.588063 0.551188 … -1.5131872 -1.5131872;;; 0.6470983 0.5504283 … 2.175899 2.175899; 0.9126755 0.99060667 … 2.1313498 2.1313498; 0.18481529 0.19919029 … -1.5683095 -1.5683095;;; 0.64341044 0.7214062 … 1.5481051 1.5481051; 0.14968371 0.22146532 … -1.5661538 -1.5661538; 0.8906089 0.8243589 … 0.3487339 0.3487339;;; 0.48673928 0.5590347 … -1.6206806 -1.6206806; 0.7801454 0.86718625 … 3.4691792 3.4691792; 0.5598935 0.50676847 … 0.25176847 0.25176847], UInt32[0x00000020, 0x00000051, 0x0000008c, 0x000000f8, 0x0000011a, 0x0000004a, 0x000000a6, 0x0000007b, 0x00000034, 0x00000076])

Optimal use

When computing multiple batches for the same model, it is more efficient to precompute a cache once and reuse it.

model = Model(Δt = 0.125f0,
            foddata = FODData((@__DIR__) * "/../../examples/fod-FC.nii.gz"),
            )
Nmc = 10
seeds = rand(Float32, 6, Nmc)
streamlines = zeros(Float32, 6, 20, Nmc)
tract_length = zeros(UInt32, Nmc)
alg = Probabilistic()
cache = TG.init(model, alg)
# this can be called repeatedly after updating seeds for example
TG.sample!(streamlines, tract_length, model, cache, alg, seeds);
size(streamlines)

(6, 20, 10)