FCUP

• FCUP
• Define the TMC
• Define the seeds
• Compute the streamlines
• plot the streamlines
• Compute the connections

This is a more advanced tutorial because we want to show how to apply a mask.

Define the TMC

import Tractography as TG

model = TG.TMC(Δt = 0.125f0,
            foddata = TG.FODData((@__DIR__) * "/../../examples/fod-FC.nii.gz"),
            cone = TG.Cone(45),
            proba_min = 0.015f0,
            )

TMC with elype Float32
 ├─ Δt = 0.125
 ├─ minimal probability = 0.015
 ├─ cone                = Cone{Int64}(45)
 ├─ mollifier           = max_mollifier
 ├─ evaluation of SH    = PreComputeAllFOD()
 └─ data : (lmax = 8)

Just for fun, we plot the FODs of the model.

using CairoMakie

f, sc = TG.plot_fod(model; n_sphere = 1500, radius = 0.3, st = 2);
cam3d = Makie.cameracontrols(sc)
cam3d.eyeposition[] = Vec3f(85, 95, -28)
cam3d.lookat[] = Vec3f(84, 95, 59)
rotate_cam!(sc.scene, 0, 0, -pi/2)
f

Define the seeds

We next apply a mask on the boundary of which the streamlines stop.

using NIfTI
mask = NIfTI.niread((@__DIR__) * "/../../examples/wm-FC.nii.gz");
TG._apply_mask!(model, mask);

┌ Warning: #= line 0 =#:
│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.
│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.
└ @ LoopVectorization ~/.julia/packages/LoopVectorization/GKxH5/src/condense_loopset.jl:1166

We compute Nmc streamlines, hence we need Nmc seeds

Nmc = 100_000
seeds = TG.from_fod(model, Nmc; maxfod_start = true)

6×100000 Matrix{Float32}:
 111.235     93.3188     74.7658     …  129.314      111.013     83.5252
 131.827     60.3483    119.558         111.176       74.931     33.3063
   6.02401    2.62542     9.02902         2.84324      9.25144    9.36183
   0.871696   0.647271    0.0316319      -0.0785332    0.350142  -0.618368
   0.407629  -0.753382    0.992977       -0.91915      0.923275  -0.778433
  -0.272     -0.116      -0.114      …    0.386        0.158      0.108

Compute the streamlines

streamlines, tract_length = TG.sample(model, TG.Deterministic(), seeds; nt = 1000);
println("Dimension of computed streamlines = ", size(streamlines))

kernel : 6.877422 seconds (985.81 k allocations: 48.376 MiB, 7.09% compilation time)
Dimension of computed streamlines = (5, 1000, 100000)

plot the streamlines

f, scene = @time TG.plot_fod(model; n_sphere = 500, radius = 0.3, st = 1);
ind_st = findall(tract_length .> 60)
TG.plot_streamlines!(scene, streamlines[:, :, ind_st[1:10:end]])
f

We can also add the seeds

scatter!(scene, seeds[1:3, ind_st[1:10:end]], color = :white)
f

Compute the connections

When computing structural connectivity, we don't need to record the entire streamline but only its extremities.

streamlines, tract_length = TG.sample(model, TG.Connectivity(TG.Deterministic()), seeds; nt = 1000);
println("Dimension of computed streamlines = ", size(streamlines))

kernel : 6.032587 seconds (405.27 k allocations: 19.388 MiB, 6.48% compilation time)
Dimension of computed streamlines = (5, 2, 100000)