🚀 Get started with with Tractography.jl

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

Basic use

In this example, we will sample Nmc streamlines from a Tractography Markov Chain (TMC).

using Tractography
const TG = Tractography

model = TMC(Δ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 TMC

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

model = TMC(Δt = 0.125f0,
            foddata = FODData((@__DIR__) * "/../../examples/fod-FC.nii.gz"),
            )
TMC with elype Float32
 ├─ Δt = 0.125
 ├─ minimal probability = 0.0
 ├─ cone                = Cone{Float32}(90.0f0)
 ├─ mollifier           = max_mollifier
 ├─ evaluation of SH    = PreComputeAllFOD()
 └─ data : (lmax = 8)

Step 2: Define the seeds

Nmc = 10 # Monte Carlo sample
seeds = rand(Float32, 6, Nmc)
6×10 Matrix{Float32}:
 0.189184  0.599411   0.369591  0.35136   …  0.436047   0.322606   0.323849
 0.873646  0.680399   0.537353  0.963496     0.103905   0.509582   0.274359
 0.794152  0.929741   0.834202  0.500795     0.0445515  0.190547   0.436875
 0.407167  0.0832273  0.20394   0.600732     0.706827   0.846928   0.0833818
 0.526066  0.150245   0.945063  0.628732     0.0758359  0.0747016  0.721972
 0.128673  0.220523   0.785185  0.4277    …  0.453732   0.152212   0.907935

Step 3: Chose a sample algorithm

alg = Probabilistic()
Probabilistic()

Step 4: Sample the streamlines

streamlines, tract_length = sample(model, alg, seeds);
(Float32[0.18918389 0.24677522 … 3.1390555 3.1390555; 0.87364626 0.95596635 … 1.006183 1.006183; 0.7941525 0.7197775 … -1.5477219 -1.5477219;;; 0.5994115 0.6584838 … -1.5639257 -1.5639257; 0.6803988 0.7252326 … -0.6427464 -0.6427464; 0.9297413 1.0303663 … 3.202866 3.202866;;; 0.36959058 0.40838584 … 6.53305 6.4994564; 0.5373532 0.645999 … 4.423181 4.5423937; 0.83420205 0.786077 … 1.5117011 1.4948261;;; 0.3513605 0.42365593 … -1.5495187 -1.5495187; 0.96349585 1.0505366 … 0.08961205 0.08961205; 0.5007953 0.4476703 … -0.14732972 -0.14732972;;; 0.04053074 0.11444362 … 2.4652789 2.4652789; 0.67750233 0.6572459 … -1.5026286 -1.5026286; 0.19839305 0.29714304 … 3.4571395 3.4571395;;; 0.2331658 0.27196106 … 9.787259 9.787259; 0.7703954 0.8790412 … 2.1000059 2.1000059; 0.61825603 0.57013106 … -1.58112 -1.58112;;; 0.6882717 0.8042444 … -1.5072757 -1.5072757; 0.21127158 0.25521687 … 3.0041265 3.0041265; 0.27018178 0.25455678 … -0.55731815 -0.55731815;;; 0.43604684 0.5129773 … 1.4604063 1.4604063; 0.103904605 0.12762308 … -1.5114424 -1.5114424; 0.04455149 0.14017649 … 0.8558017 0.8558017;;; 0.32260633 0.34604216 … -1.5287328 -1.5287328; 0.5095817 0.52448964 … -0.48310292 -0.48310292; 0.19054675 0.31242174 … -0.340703 -0.340703;;; 0.3238486 0.22100955 … 0.6504478 0.6504478; 0.2743587 0.21741238 … -1.5784667 -1.5784667; 0.43687546 0.47937545 … 1.4931253 1.4931253], UInt32[0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000])

Optimal use

It is often better to cache some data when computing batches of streamlines. This would be done as follows

model = TMC(Δ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)