FCUP

• FCUP
• Define the TMC
• Define the seeds
• Compute 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,
            odfdata = TG.ODFData((@__DIR__) * "/../../examples/fod-FC.nii.gz"),
            cone = TG.Cone(15),
            proba_min = 0.005f0,
            )

TMC with elype Float32
 ├─ Δt = 0.125
 ├─ minimal probability = 0.005
 ├─ cone                = Cone{Int64}(15)
 ├─ mollifier           = default_mollifier
 ├─ evaluation of SH    = Tractography.PreComputeAllODF()
 └─ data : ⋯

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

using CairoMakie

f, sc = TG.plot_odf(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

using LinearAlgebra

Nmc = 100_000
_ind1 = findall(mask .== 1)
seed_id = rand(1:length(_ind1), Nmc)
seeds = zeros(Float32, 6, Nmc)
for i=1:Nmc
    seeds[:,i] .= vcat(TG.transform(model.odfdata, _ind1[seed_id[i]])[1:3]..., normalize(randn(3)))
end

Compute the streamlines

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

kernel : 2.107880 seconds (952.58 k allocations: 48.703 MiB, 22.35% compilation time)
Dimension of computed streamlines = (3, 1000, 100000)

f, sc = @time TG.plot_odf(model; n_sphere = 500, radius = 0.3, st = 2);
TG.plot_streamlines!(sc, streamlines[1:3, 1:1:end, 1:10:Nmc])
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 : 1.930161 seconds (402.18 k allocations: 19.969 MiB, 19.83% compilation time)
Dimension of computed streamlines = (3, 2, 100000)