API Reference

GeneralGradient(g⃗, TE)

General gradient sequence $\vec{g}(t) \in \mathbb{R}^3$. The direction and amplitude may vary in time.

ScalarGradient(dir, profile, amplitude)

Gradient sequence with amplitude amplitude, normal direction dir and scalar time profile profile.

The direction is constant, while the amplitude is controlled by the time profile.

The direction is normalized upon construction.

f = PGSE(δ, Δ)

Pulsed Gradient Spin Echo sequence with pulse duration δ and time between pulses Δ-δ.

f = DoublePGSE(δ, Δ)

Double Pulsed Gradient Spin Echo sequence with four pulses of duration δ, separated by pauses of duration Δ-δ, 0, and Δ-δ repsectively.

f = CosOGSE(δ, Δ, nperiod)

Oscillating Gradient Spin Echo sequence with two cos-pulses of duration δ separated by a pause of duration Δ-δ for nperiod periods per pulse.

f = SinOGSE(δ, Δ, nperiod)

Oscillating Gradient Spin Echo sequence with two sin-pulses of duration δ separated by a pause of duration Δ-δ for nperiod periods per pulse.

integral(f, t = echotime(f))

Integral of time profile f between 0 and t. Unless specified, the echotime is used as the upper integral limit.

For the PGSE, SinOGSE, CosOGSE and DoublePGSE sequences, analytical expressions are available. Otherwise a numerical integral is computed.

int_F²(f)

Compute the time profile contribution to the $b$-value given the time profile f. The $b$-value is given by

b = γ^2 * g^2 * int_F²(f)

where γ is the gyromagnetic ratio of the water proton and g is the gradient amplitude.

intervals(f)

Get characteristic intervals of the time profile f.

isconstant(profile)

Return true if the time profile profile is intervalwise constant.

echotime(f::TimeProfile)

Get echo time $TE$ of the time profile f, which is the end of the last characteristic interval.

echotime(grad::AbstractGradient)

Get echo time $TE$ of gradient.

assemble_matrices(model)

Assemble finite element matrices.

split_mesh(mesh_all)

Split global mesh into compartments.

split_field(mesh, ξ)

Split the field ξ into a vector containing a view of each compartment.

get_cmpt_volumes(mesh)

Get compartment volumes.

Model(; mesh, ρ, D, T₂, κ, γ)

Finite element discretized biological model with initial spin densities ρ, diffusion tensors D, T₂-relaxation times T₂, wall permeabilities κ, and gyromacnetic ratio γ. The vectors ρ, D, T₂, are of length ncompartment, while κ is of length nboundary.

initial_conditions(model)

Get initial conditions at all degrees of freedom.

savefield(mesh, ξ, filename::String; fieldname = "Magnetization")

Save field ξ to a VTK file. It may then be visualized in Paraview.

compute_adc_sta(model)

Compute the ADC for each compartment in the short diffusion time regime.

fit_adc(bvalues, signals)

Fit the apparent diffusion coefficient (ADC) using a polynomial logfit of the normalized signals signals against the b-values bvalues.

fit_tensors(directions, adcs)

Fit effective diffusion tensors to directionalized ADCs. The six components of the symmetric diffusion tensors are fitted by least squares to the gradient directions and resulting ADCs.

plot_mesh(femesh, compartments = 1:ncompartment)

Plot finite element mesh, with a subset of the compartments.

plot_field(femesh, ξ)

Plot field ξ on the finite element mesh.

GeneralBTPDE(; model, matrices)

Bloch-Torrey PDE problem.

HADC(; model, matrices)

HADC problem.

Karger(; model, difftensors)

Karger problem.

Laplace(; model, matrices, neig_max)

Laplace eigenvalue problem.

MatrixFormalism(; model, matrices, lap_eig)

Matrix formalism problem. Given a Laplace eigendecomposition lap_eig, this problem consists of computing the MF magnetization.

AnalyticalLaplace(; ρ, r, D, W, T₂, γ, dim, eiglim, eigstep)

Analytical radial Laplace eigenvalue problem.

AnalyticalMatrixFormalism(; analytical_laplace, lap_mat, volumes)

Analytical Matrix formalism problem. Given a radial Laplace eigendecomposition analytical_laplace, this problem consists of computing the MF compartment signals.

IntervalConstantSolver(; θ = 0.5, timestep)

BTPDE solver specialized on intervalwise constant ScalarGradients, e.g PGSE, DoublePGSE. Raises an error for other gradients.

solve(problem::Laplace)

Compute the Laplace eigenvalues, eigenfunctions and first order moments of products of pairs of eigenfunctions.

solve(problem::MatrixFormalism, gradient; ninterval = 500)

Solve for magnetization using Matrix Formalism.

solve(problem::AnalyticalLaplace)

Solve analytical Laplace eigenvalue problem for symmetrical geometries.

solve(problem::AnalyticalMatrixFormalism, gradient)

Compute the signal in a multilayered cylinder or sphere using an analytical matrix formalism solution.

This function is based on the following articles and corresponding code: [1] D. S. Grebenkov, NMR Survey of Reflected Brownian Motion, Rev. Mod.Phys. 79, 1077 (2007) [2] D. S. Grebenkov, Pulsed-gradient spin-echo monitoring of restricted diffusion in multilayered structures, J. Magn. Reson. 205, 181-195 (2010).

solve(
    problem::GeneralBTPDE,
    gradient,
    odesolver = QNDF();
    abstol = 1e-6,
    reltol = 1e-4,
    callbacks = [].
)

Solve the Bloch-Torrey partial differential equation using P1 finite elements in space and odesolver in time.

function solve(
    problem::BTPDE,
    gradient::ScalarGradient,
    odesolver::IntervalConstantSolver;
    callbacks = [],
)

Solve the Bloch-Torrey partial differential equation using P1 finite elements in space and a theta-rule in time. This time stepping scheme requires a degree of implicitness θ and a time step Δt:

• θ = 0.5: Crank-Nicolson (second order)
• θ = 1.0: Implicit Euler (first order)

The function only works for interval-wise constant ScalarGradients, and errors otherwise.

solve(problem::Karger, gradient, odesolver = MagnusGL6(); timestep)

Solve the finite pulse Karger model (FPK) using precomputed effective diffusion tensors difftensors.

solve(problem::HADC, gradient::ScalarGradient, odesolver = QNDF();
    abstol = 1e-6,
    reltol = 1e-4,
)

Compute the ADC using a homogenized ADC model (HADC). This is currently only implemented for scalar gradients.

solve_multigrad(problem, gradients, args...; kwargs...)

Solve problem for multiple magnetic field gradients.

length2eig(length, D)

Convert length scale to Laplace eigenvalue.

eig2length(λ, D)

Convert Laplace eigenvalue to length scale.

limit_lengthscale(lap_eig, λ_max)

Only keep modes with length scales larger than minimum

compute_mf_diffusion_tensor(mesh, M, lap_eig, gradient::ScalarGradient)

Compute effective diffusion tensors using the matrix formalism.

unitcircle(ndirection; half = false, normal = [0, 0, 1])

Create ndirection directions unformly distributed on the unit circle defined by normal.

If half is true, the directions lie on a half-circle instead of the whole circle.

unitsphere(ndirection; half = false, normal = [0, 0, 1])

Create ndirection directions unformly distributed on the unit sphere.

If half is true, the points will be discributed on the hemisphere defined by normal.

compute_signal(M, ξ)

Compute signal from magnetization ξ, using the mass matrix M for integration.

Given a mesh mesh and a vector of compartment mass matrices M_cmpts, a vector of compartment signals may be obtained by compute_signal.(M_cmpts, split_field(mesh, ξ)).

Printer(; nupdate = 1)

Print time stepping information to the console.

VTKWriter(; nupdate = 1, dir = "output", filename = "solution")

Write magnetization field to a VTK file after every nupdate timestep. The files are stored in a ParaView data collection file (PVD). The magnetization time series may be visualized in ParaView by opening the file "$dir/$filename.pvd". The compartments are labeled.

Plotter(; nupdate = 1)

Plot the evolution of the BTPDE during time stepping. This requires loading the GLMakie plotting backend (]add GLMakie; using GLMakie). The plot is updated every nupdate time step. The resulting figure contains a plot of the time profile, total signal attenuation, and magnetization field (complex magnitude and phase shift).

Setup recipe for a set of stacked plates.

Setup recipe for a set of cylinders immersed in a ECS.

Setup recipe for a set of spheres immersed in a ECS.

Setup recipe for a neuron, possibly wrapped in a ECS.

coefficients(setup; D, T₂, ρ, κ, γ)

Order coefficients compartment arrays.

analytical_coefficients(setup, coeffs)

Get coefficients for the analytical module.

create_geometry(setup; recreate = true)

Create cells, surfaces and finite element mesh. If recreate = false, previous geometry will be reused.

This function does the following:

• Check geometry setup consistency
• Create or load cell configuration
• Create or load surface triangulation
• Call TetGen
• Deform domain
• Split mesh into compartments

For custom geometries with more than one compartment, call split_mesh directly instead. This requires facet and element labels.