Library

from_fod(
    model::Model{𝒯},
    n_seeds::Int64;
    n_sphere,
    maxfod_start
) -> Any

Generate seeds from orientation distribution functions.

Keyword arguments

• maxfod_start = false The seeds are generated in voxels with non-zero average FOD with the orientations importance sampled.
• maxfod_start = true The seeds are generated in voxels with non-zero average FOD with the orientations corresponding to the maximum probability.

from_mask(
    model::Model{𝒯},
    mask::AbstractArray{Bool, 3},
    n_seeds::Int64
) -> Any

Generate seeds from a 3D mask. The seeds are generated randomly inside voxels with non-zero values in the mask. The orientations are random.

struct PlottingFOD <: Tractography.AbstractFODEvaluation

Set up for plotting ODF.

struct PreComputeAllFOD <: Tractography.AbstractFODEvaluation

Evaluation of the basis based on Fibonacci sampling. All ODF are pre-computed once and saved in a cache. Their positivity is enforced with a max(0,⋅) or a mollifier.

Details

If you have na angles for sampling the unit sphere and the data is of size (nx, ny, nz, nsph), it yields a matrix of dimensions (nx, ny, nz, na).

struct DirectFOD <: Tractography.AbstractFODEvaluation

The evaluation of the basis (for example spherical harmonics) is done on the fly. Requires little memory.

See also PreComputeAllFOD

struct FODData{𝒯, 𝒯d, 𝒯s, 𝒯t, 𝒯b, 𝒯i}

Structure to hold FOD data.

Internal fields

• filename::String: filename from which the (fod) data is read.

• data::Any: field which contains the FOD data.

• lmax::Int64: max l coordinate in of spherical harmonics.

• transform::Tractography.Transform: transform associated with data, see Transform.

• normalized::Bool: Are the data normalized? In this case fod[i,j,l,1] ∈ {0,1}.

• basis::Any: Basis for fODF. For example SphericalHarmonics()

• interpolation::Any: Interpolation of fODF

Methods

• get_lmax(::FODData) returns the max l coordinate in of spherical harmonics.
• size(::FODData) returns the size of the data.
• get_range(::FODData) returns the range of the data in the real world coordinates.
• is_normalized(::FODData) return whether the data is normalized.

Constructors

In the following, data has shape nx x ny x nz x nt.

• FODData(data::AbstractArray{𝒯, 4}, transform::Transform, normalize_it::Bool; file_name = "None").
• FODData(data::AbstractArray{𝒯, 4}, S, T, normalize_it::Bool; file_name = "None").
• FODData(data::AbstractArray{𝒯, 4}, normalize_it::Bool) uses a trivial transform.
• FODData(file_name::String; normalize_it::Bool = true, k...) the transform is extracted from the nii file.

struct Model{𝒯, 𝒯alg<:Tractography.AbstractFODEvaluation, 𝒯d, 𝒯C, 𝒯mol}

Model of streamlines.

Internal fields (with default values):

• Δt::Any: Step size of the Model. Default: 0.1

• evaluation_algo::Tractography.AbstractFODEvaluation: Spherical harmonics evaluation algorithm. Can be PreComputeAllFOD(), DirectFOD(). Default: PreComputeAllFOD()

• foddata::Any: ODF data, typically from nifti file but can be passed as an Array. Must be the list of ODF in the base of spherical harmonics. Hence, it should be an (abstract) 4d array. Default: nothing

• cone::Any: Cone function to restrict orientation sample. You can use a Cone or a custom function (d1, d2) -> return_a_boolean. Default: Cone(90.0f0)

• proba_min::Any: Probability below which we stop tracking. Default: 0.0

• mollifier::Any: Mollifier, used to make the fodf non negative. During odf evaluation, we effectively use mollifier(fodf[angle,i,j,k]). You can also use softplus using mollifier = Base.Fix2(softplus, 10) or remove the mollifier using mollifier = identity. Default: max_mollifier

Methods

• _apply_mask!(model, mask) apply a mask to the raw SH tensor. See its doc string.
• _getdata(model) return the fodf data associated with the model.
• size(model) return nx, ny, nz, nt.
• eltype(model) return the scalar type of the data (default Float64).
• get_lmax(model) return the max l coordinate in of spherical harmonics.

Constructors (pass the internal fields!)

• Model()
• Model(Δt = 0.1f0) for a Float32 Model
• Model(Δt = 0.1, proba_min = 0.) for a Float64 Model. You need to specify both fields Δt and proba_min
• Model(odfdata = rand(10,10,10,45)) for custom ODF

struct Deterministic <: Tractography.DeterministicSampler

Tractography sampling of the Tractography Markov Chain (TMC) performed with the argmax function. Can be used with the basis evalution strategy PreComputeAllFOD.

struct Probabilistic <: Tractography.AbstractNotPureRejectionSampler

Tractography sampling of the Tractography Markov Chain (TMC) performed with the cumulative sum distribution. Can be used with the basis evalution strategy PreComputeAllFOD.

Constructor

Probabilistic()

struct Diffusion{T, Tmol, Tdmol} <: Tractography.AbstractSDESampler{T}

Tractography sampling of the diffusive model performed with geometric Euler-Maruyama method [1]; its precision is weak order 1. The streamlines (Xₜ)ₜ are solution of the SDE

dXₜ = Uₜ⋅dt

dUₜ = γ⋅∇log f(Uₜ)⋅dt + √(2γ ⋅ γ_noise) ⋅ dnoiseₜ

Arguments (with default values):

• γ::Any: γ parameter of the diffusion process. It is related to the curvature of the streamline. Default: 1.0

• γ_noise::Any: parameter of the diffusion process to scale the variance. Default: 1.0

• mollifier::Any: mollifier. Default: Base.Fix2(softplus, 10)

• d_mollifier::Any: differential of mollifier. Default: Base.Fix2(∂softplus, 10)

• adaptive::Bool: Fixed time step? Default: false

Constructor

Example for Float32: Diffusion(γ = 1f0). If you want Float64, you have to pass the two scalars

```Diffusion(γ = 1.0, γ_noise = 1.0)```

Reference(s)

[1] Bharath, K., Lewis, A., Sharma, A., & Tretyakov, M. V. (n.d.). Sampling and Estimation on Manifolds using the Langevin Diffusion.

Define the transport algorithm. Options are the same as for Diffusion.

Arguments (with default values):

• γ::Any: γ parameter of the diffusion process. It is related to the curvature of the streamline. Default: 1.0

• mollifier::Any: mollifier. Default: Base.Fix2(softplus, 10)

• d_mollifier::Any: differential of mollifier. Default: Base.Fix2(∂softplus, 10)

• adaptive::Bool: Fixed time step? Default: false

struct Connectivity{Talg} <: Tractography.AbstractSampler

Tractography based sampling of structural connectivity. Do not compute the full streamline but only return the first/last points and the streamlines lengths. This allows to compute many more streamlines on GPU where memory is limited.

Constructor example

• Connectivity(Probabilistic())

init(
    model::Model{𝒯},
    alg;
    n_sphere,
    𝒯ₐ
) -> Tractography.ThreadedCache{_A, Nothing, _B, Nothing, Nothing} where {_A, _B}

Create a cache for computing streamlines in batches. This is useful for memory-limited environments (e.g. GPU).

Arguments

• alg sampling algorithm, Deterministic, Probabilistic, Diffusion, etc.
• n_sphere::Int = 400 number of points to discretize the sphere for spherical harmonics evaluation.
• 𝒯ₐ = Array{𝒯} array type for the cache. Pass a GPU array type like CuArray to run on GPU; leave as Array for CPU.

sample(
    model::Model{𝒯},
    alg::Tractography.AbstractSampler,
    seeds::AbstractArray{𝒯, 2};
    ...
) -> Tuple{Any, Any}
sample(
    model::Model{𝒯},
    alg::Tractography.AbstractSampler,
    seeds::AbstractArray{𝒯, 2},
    mask::Union{Nothing, AbstractArray{Bool}};
    nt,
    n_sphere,
    maxfod_start,
    reverse_direction,
    nthreads,
    gputhreads,
    saveat
) -> Tuple{Any, Any}

Sample the Model model.

Arguments

• model::Model
• alg sampling algorithm, Deterministic, Probabilistic, Diffusion, etc.
• seeds matrix of size 6 x Nmc where Nmc is the number of Monte-Carlo simulations to be performed.
• mask = nothing matrix of boolean where to stop computation. See also _apply_mask.

Optional arguments

• nt::Int maximal number of steps to compute each streamline.
• n_sphere::Int = 400 number of points to discretize the sphere on which we evaluate the spherical harmonics.
• maxfod_start::Bool for each locations, use direction provided by the argmax of the ODF.
• reverse_direction::Bool reverse initial direction.
• nthreads::Int = 8 number of threads on CPU.
• gputhreads::Int = 512 number of threads on GPU.
• saveat::Int record the streamline every saveat step. Only available for alg <: Diffusion.

Output

• streamlines with shape 3 x nt x Nmc

sample!(
    streamlines,
    streamlines_length,
    model::Model{𝒯},
    cache::Tractography.AbstractCache,
    alg,
    seeds;
    maxfod_start,
    reverse_direction,
    nthreads,
    gputhreads,
    nₜ,
    saveat,
    𝒯ₐ
) -> Any

Sample the model in place by overwriting streamlines. Uses minimal memory and can run indefinitely on GPU.

Arguments

• streamlines array with shape 3 x nt x Nmc. nt is the maximum length per streamline. Nmc is the number of Monte-Carlo simulations.
• streamlines_length lengths of the streamlines.
• model::Model model to sample from.
• alg sampling algorithm: Deterministic, Probabilistic, Diffusion, etc.
• seeds matrix of size 6 x Nmc with positions (x,y,z) and directions (u,v,w).

Optional arguments

• maxfod_start::Bool use the argmax direction of the ODF at each location.
• reverse_direction::Bool reverse the initial direction.
• nthreads::Int = 8 number of CPU threads.
• gputhreads::Int = 512 number of GPU threads.

plot_streamlines!(ax, streamlines; k...)

Plot the streamlines with lines colored by local orientation.

The optional parameters are passed to lines!.

plot_fod!(ax, model; n_sphere, radius, st, I, J, K, c0min)

Plot the the ODF with glyphs.

See also

• plot_fod(model; kwargs...) for an out-of-place method with same parameters.

Arguments

• ax GLMakie scene
• model::Model

Optional arguments

• I, J, K: range for displaying the fODFs. Some of them can be a single element, like K = 40:40.
• radius: radius of the glyph.
• st = 4: stride, only show one over st glyph in each direction.

plot_slice(model; k...)

Plot a slice of the data.

See also

• plot_slice!

Arguments

• model::Model

Optional arguments

• odf::Bool display the ODF with glyphs.
• slice::Bool display the brain slice from mean ODF value.
• interpolate::Bool interpolate the brain slice (reduces pixelization).
• alpha alpha value for brain slice.
• I,J,K range for displaying the ODF. One of them should be a single element, like K = 40:40.
• st::Int = 2
• c0min minimum mean ODF value for plotting the glyph.
• n_theta::Int number of points for showing the glyphs.
• radius = 0.1 radius of the glyphs.

struct Cone{𝒯<:Real}

Structure to encode a cone to limit sampling the direction. This ensures that the angle in degrees between to consecutive streamline directions is less than angle.

The implemented condition is for cn = Cone(angle).

(cn::Cone)(d1, d2) = dot(d1, d2) > cosd(cn.alpha)

Fields

• alpha::Real: half section angle in degrees

Constructor

Cone(angle)

_apply_mask!(model, mask)

Multiply the mask which is akin to a matrix of Bool with same size as the data stored in model. Basically, the mask mask[ix, iy, iz] ensures whether the voxel (ix, iy, iz) is discarded or not.

Arguments

• model::Model.
• mask can be a AbstractArray{3, Bool} or a NIVolume.

Exp𝕊²(p, X, t) -> Any

Exponential map on the sphere. Assumes t > 0.

See https://github.com/JuliaManifolds/ManifoldsBase.jl/blob/5c4a61ed3e5e44755a22f7872cb296a621905f87/test/ManifoldsBaseTestUtils.jl#L63

spherical_to_euclidean(θ, ϕ) -> Tuple{Any, Any, Any}

Transform spherical to cartesian coordinates, the chart on the sphere minus the north/south poles.

Its coordinates are (sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ)).

Recall that θ ∈ [0, π] and ϕ ∈ [0, 2π].

euclidean_to_spherical(x, y, z) -> Tuple{Any, Any}

Transform cartesian to spherical coordinates. Assume that the vector has norm one. Return (θ, ϕ) where θ ∈ [0, π] and ϕ ∈ [-π, π].

softplus(x) -> Any
softplus(x, k) -> Any

Smooth approximation of the ReLU function. Computed as log1p(exp(k*x)) / k. The parameter k controls the steepness: as k → ∞, softplus approaches ReLU. For k*x ≥ 30, returns x directly to avoid overflow.

orientation_probabilities(
    alg,
    fodf::AbstractArray{𝒯, 4},
    cone::AbstractArray{𝒯, 2},
    voxel_index₁,
    voxel_index₂,
    voxel_index₃,
    ind_u::UInt32,
    n_angles::UInt32
) -> Tuple{Any, Any, UInt32}

Compute orientation probabilities for a given voxel.

For each direction i, the conditioned probability is fodf[i] * cone[i, ind_u]. Returns a tuple (conditioned_proba, total_proba, ind_max) where:

• conditioned_proba: sum of fODF × cone over all directions
• total_proba: sum of raw fODF over all directions
• ind_max: index of maximum probability (only meaningful for DeterministicSampler)

next_orientation(
    _::Tractography.DeterministicSampler,
    ind_u::UInt32,
    ind_max::UInt32,
    conditioned_proba,
    fodf::AbstractArray{𝒯, 4},
    cone::AbstractArray{𝒯, 2},
    voxel_index₁,
    voxel_index₂,
    voxel_index₃,
    ind_u0::UInt32,
    n_angles::UInt32
) -> UInt32

Select the next orientation given the sampling algorithm.

• DeterministicSampler: returns the index of the direction with maximum probability (ind_max).
• Probabilistic: draws a random sample from the cumulative distribution of conditioned probabilities.