API

struct Isochromat{T<:Real} <: FieldVector{3,T}
    x::T
    y::T
    z::T
end

Holds the x,y,z components of a spin isochromat in a FieldVector, which is a StaticVector (from the package StaticArrays) with custom fieldnames.

decay(m::Isochromat{T}, E₁, E₂) where T

Apply T₂ decay to transverse component and T₁ decay to longitudinal component of Isochromat.

initialize_states(::AbstractResource, ::IsochromatSimulator{T}) where T

Initialize a spin isochromat to be used throughout a simulation of the sequence.

This may seem redundant but to is necessary to share the same programming interface with EPGSimulators.

invert(m::Isochromat{T}, p::AbstractTissueProperties) where T

Invert Isochromat with B₁ insenstive (i.e. adiabatic) inversion pulse

invert(m::Isochromat{T}, p::AbstractTissueProperties) where T

Invert z-component of Isochromat (assuming spoiled transverse magnetization so xy-component zero).

regrowth(m::Isochromat{T}, E₁) where T

Apply T₁ regrowth to longitudinal component of Isochromat.

rotate(m::Isochromat, γΔtGRz, z, Δt, p::AbstractTissueProperties)

Rotation of Isochromat without RF (so around z-axis only) due to gradients and B0 (i.e. refocussing slice select gradient).

rotate(m::Isochromat{T}, γΔtRF, γΔtGR, r, Δt, p::AbstractTissueProperties, ΔtΔω=zero(T)) where T

RF, gradient and/or ΔB₀ induced rotation of Isochromat computed using Rodrigues' rotation formula.

Arguments

• m: Input isochromat state.
• γΔtRF: Complex value representing RF pulse effect (radians). B₁ scaling from p is applied internally if hasB₁(p).
• γΔtGR: Tuple/Vector representing gradient effect γ * G * Δt (radians/centimeter).
• r: Position vector (x,y,z) [cm].
• Δt: Time step duration (seconds).
• p: Tissue properties (AbstractTissueProperties). B₀ effects from p are applied internally if hasB₀(p).
• ΔtΔω: Additional off-resonance phase accumulated during Δt (radians). Often Δt * Δω, where Δω is in rad/s.

Returns

• Rotated isochromat state.

sample!(output, index::Union{Integer,CartesianIndex}, m::Isochromat)

Sample transverse magnetization from Isochromat. The "+=" is needed for 2D sequences where slice profile is taken into account.

sample_xyz!(output, index::Union{Integer,CartesianIndex}, m::Isochromat)

Sample m.x, m.y and m.z components from Isochromat. The "+=" is needed for 2D sequences where slice profile is taken into account.

F̄₋(Ω)

View into the second row of the configuration state matrix Ω, corresponding to the F̄₋ states.

F₊(Ω)

View into the first row of the configuration state matrix Ω, corresponding to the F₊ states.

RF_rotation_matrix(RF, α)

The 3x3 rotation matrix applied to each configuration state by an RF pulse with flip angle α (radians) and phase angle(RF). Split out of excite! so that forward-sensitivity propagation can differentiate it w.r.t. the flip angle (see ∂RF_rotation_matrix∂B₁ in src/derivatives/forward_sensitivity.jl).

RF_rotation_matrix(RF::Real, α)

Version for a real-valued RF pulse (zero phase). See the complex version above.

Z(Ω)

View into the third row of the configuration state matrix Ω, corresponding to the Z states.

decay!(Ω::AbstractConfigurationStates, E₁, E₂)

Apply T₁ and T₂ relaxation effects to the EPG states Ω.

Arguments

• Ω: The configuration state matrix.
• E₁: T₁ relaxation factor, exp(-Δt/T₁), where Δt is the time duration (seconds) and T₁ is from the tissue properties (seconds).
• E₂: T₂ relaxation factor, exp(-Δt/T₂), where Δt is the time duration (seconds) and T₂ is from the tissue properties (seconds).

dephasing!(Ω::AbstractConfigurationStates)

Shift states around due to dephasing gradient: The F₊ go up one, the F̄₋ go down one and Z do not change

diffuse!(Ω::AbstractConfigurationStatestates, diffusion_decay)

Apply diffusion decay according to state number by element-wise multiplication with the pre-computed diffusion decay matrix.

• Ω: The configuration state matrix
• diffusion_decay: The pre-calculated diffusion decay terms

diffusion_decay_matrix(Ω::EPGSAbstractConfigurationStatestates, D)

Pre-calculate diffusion decay according to state number. Store in a matrix with the same size/type as Ω s.t. later on the decay can be applied by element-wise multiplication.

Arguments

• Ω: The configuration state matrix
• D: The diffusion coefficient

On GPU, when each thread only holds parts of the configuration state matrix, each thread only computes the corresponding columns of the diffusion decay matrix.

excite!(Ω::AbstractConfigurationStates, RF, p::AbstractTissueProperties)

Apply RF pulse rotation to the EPG states Ω.

Arguments

• Ω: The configuration state matrix.
• RF: Complex RF pulse value. abs(RF) is the flip angle (degrees), angle(RF) is the pulse phase (radians). B₁ scaling from p is applied internally if hasB₁(p).
• p: Tissue properties (AbstractTissueProperties).

excite!(Ω::AbstractConfigurationStates, RF, p::AbstractTissueProperties) where {T<:Union{Real, Quantity{<:Real}}}

Apply RF pulse rotation to the EPG states Ω (version for real-valued RF pulse, assuming zero phase).

Arguments

• Ω: The configuration state matrix.
• RF: Real RF pulse value representing the flip angle (degrees). B₁ scaling from p is applied internally if hasB₁(p).
• p: Tissue properties (AbstractTissueProperties).

flip_angle(RF, p::AbstractTissueProperties)

Flip angle in radians of an RF pulse, including B₁ scaling if hasB₁(p). For a complex RF the nominal flip angle in degrees is abs(RF) (its phase is handled separately); a real RF is the nominal flip angle in degrees, signed.

initial_conditions!(Ω::AbstractConfigurationStates)

Set all components of all states to 0, except the Z-component of the 0th state which is set to 1.

initialize_states(::AbstractResource, sequence::EPGSimulator{T,Ns}) where {T,Ns}

Initialize an MMatrix of EPG states on CPU to be used throughout the simulation.

initialize_states(::CUDALibs, sequence::EPGSimulator{T,Ns}) where {T,Ns}

Initialize an array of EPG states on a CUDA GPU to be used throughout the simulation.

invert!(Ω::AbstractConfigurationStates, p::AbstractTissueProperties)

Invert Z-component of states of all orders. Assumes fully spoiled transverse magnetization.

invert!(Ω::AbstractConfigurationStates)

Invert with B₁ insenstive (i.e. adiabatic) inversion pulse

regrowth!(Ω::AbstractConfigurationStates, E₁)

Apply T₁ regrowth to the longitudinal magnetization equilibrium state (Z₀).

Arguments

• Ω: The configuration state matrix.
• E₁: T₁ relaxation factor, exp(-Δt/T₁), where Δt is the time duration (seconds) and T₁ is from the tissue properties (seconds). The regrowth amount is (1 - E₁).

rotate!(Ω::AbstractConfigurationStates, eⁱᶿ)

Apply phase accrual due to off-resonance to the transverse EPG states (F₊, F̄₋).

Arguments

• Ω: The configuration state matrix.
• eⁱᶿ: Complex rotation factor, typically exp(-im * Δω * Δt), where Δω is the off-resonance frequency (rad/s, potentially derived from p.B₀) and Δt is the time duration (seconds).

rotate_decay!(Ω::AbstractConfigurationStates, E₁, E₂, eⁱᶿ)

Apply combined off-resonance rotation and T₁/T₂ relaxation to the EPG states Ω.

Arguments

• Ω: The configuration state matrix.
• E₁: T₁ relaxation factor (exp(-Δt/T₁)).
• E₂: T₂ relaxation factor (exp(-Δt/T₂)).
• eⁱᶿ: Complex off-resonance rotation factor (exp(-im * Δω * Δt)). (See rotate! and

decay! for details on arguments).

sample_transverse!(output, index::Union{Integer,CartesianIndex}, Ω::AbstractConfigurationStates)

Sample the measurable transverse magnetization, that is, the F₊ component of the 0th state. The += is needed for 2D sequences where slice profile is taken into account.

sample_Ω!(output, index::Union{Integer,CartesianIndex}, Ω::AbstractConfigurationStates)

Sample the entire configuration state matrix Ω. The += is needed for 2D sequences where slice profile is taken into account.

Ω_eltype(sequence::EPGSimulator{T,Ns}) where {T,Ns} = Complex{T}

By default, configuration states are complex. For some sequences, they will only ever be real (no RF phase, no complex slice profile correction) and for these sequences a method needs to be added to this function.

AbstractTissueProperties{N,T} <: FieldVector{N,T}

Abstract type for custom structs that hold tissue properties used for a simulation within one voxel. For simulations, SimulationParameterss are used that can be assembled with the @parameters macro.

Possible fields:

• T₁::T: Longitudinal relaxation time constant in seconds.
• T₂::T: Transverse relaxation time constant in seconds.
• B₁::T: Relative transmit B₁ field scaling factor (dimensionless, typically around 1.0).
• B₀::T: Off-resonance frequency in Hz (Hertz).
• D::T: Diffusion coefficient in m²/s. The decay time constant T_D is calculated from D, TR, and the spoiling gradient area Δk in the sequence implementation.
• ρˣ::T: Real part of proton density or equilibrium magnetization M₀ (arbitrary units, dimensionless scaling factor).
• ρʸ::T: Imaginary part of proton density or equilibrium magnetization M₀ (arbitrary units, dimensionless scaling factor).

Implementation details:

The structs are subtypes of FieldVector, which is a StaticVector with named fields (see the documentation of StaticArrays.jl). There are three reasons for letting the structs be subtypes of FieldVector:

1. FieldVectors/StaticVectors have sizes that are known at compile time. This is beneficial for performance reasons
2. The named fields improve readability of the code (e.g. p.B₁ vs p[3])
3. Linear algebra operations can be performed on instances of the structs. This allows, for example, subtraction (without having to manually define methods) and that is useful for comparing parameter maps.

T₁T₂{T} <: AbstractTissueProperties{2,T}

Tissue properties struct containing T₁ and T₂. Units are defined in AbstractTissueProperties.

T₁T₂B₀{T} <: AbstractTissueProperties{3,T}

Tissue properties struct containing T₁, T₂, and B₀. Units are defined in AbstractTissueProperties.

T₁T₂B₀D{T} <: AbstractTissueProperties{5,T}

Tissue properties struct containing T₁, T₂, B₀, and D. Units are defined in AbstractTissueProperties.

T₁T₂B₀Dρˣρʸ{T} <: AbstractTissueProperties{6,T}

Tissue properties struct containing T₁, T₂, B₀, D, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

T₁T₂B₀ρˣρʸ{T} <: AbstractTissueProperties{5,T}

Tissue properties struct containing T₁, T₂, B₀, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

T₁T₂B₁{T} <: AbstractTissueProperties{3,T}

Tissue properties struct containing T₁, T₂, and B₁. Units are defined in AbstractTissueProperties.

T₁T₂B₁B₀{T} <: AbstractTissueProperties{4,T}

Tissue properties struct containing T₁, T₂, B₁, and B₀. Units are defined in AbstractTissueProperties.

T₁T₂B₁B₀D{T} <: AbstractTissueProperties{5,T}

Tissue properties struct containing T₁, T₂, B₁, B₀, and D. Units are defined in AbstractTissueProperties.

T₁T₂B₁B₀Dρˣρʸ{T} <: AbstractTissueProperties{7,T}

Tissue properties struct containing T₁, T₂, B₁, B₀, D, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

T₁T₂B₁B₀ρˣρʸ{T} <: AbstractTissueProperties{6,T}

Tissue properties struct containing T₁, T₂, B₁, B₀, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

T₁T₂B₁D{T} <: AbstractTissueProperties{4,T}

Tissue properties struct containing T₁, T₂, B₁, and D. Units are defined in AbstractTissueProperties.

T₁T₂B₁Dρˣρʸ{T} <: AbstractTissueProperties{6,T}

Tissue properties struct containing T₁, T₂, B₁, D, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

T₁T₂B₁ρˣρʸ{T} <: AbstractTissueProperties{5,T}

Tissue properties struct containing T₁, T₂, B₁, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

T₁T₂D{T} <: AbstractTissueProperties{3,T}

Tissue properties struct containing T₁, T₂, and D. Units are defined in AbstractTissueProperties.

T₁T₂Dρˣρʸ{T} <: AbstractTissueProperties{5,T}

Tissue properties struct containing T₁, T₂, D, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

T₁T₂ρˣρʸ{T} <: AbstractTissueProperties{4,T}

Tissue properties struct containing T₁, T₂, ρˣ, ρʸ. Units are defined in AbstractTissueProperties.

macro parameters(args...)

Create a SimulationParameters with the actual struct type being determined by the arguments passed to the macro.

Examples

# Create a StructArray{T₁T₂} with T₁ and T₂ values
T₁, T₂ = rand(100), 0.1*rand(100)
parameters = @parameters T₁ T₂

# Create a StructArray{T₁T₂B₁} with T₁, T₂ and B₁ values
T₁, T₂, B₁ = rand(100), 0.1*rand(100), ones(100)
parameters = @parameters T₁ T₂ B₁

# Create a StructArray{T₁T₂B₀} with T₁, T₂ and B₀ values
# This time use the aliases that don't use unicode characters
T1, T2, B0 = rand(100), 0.1*rand(100), ones(100)
parameters = @parameters T1 T2 B0

BlochSimulator{T}

The abstract type of which all sequence simulators will be a subtype. The parameter T should be a number type (e.g. Float64, Float32) and the tissueparameters that are used as input to the simulator should have the same number type. By convention, a BlochSimulator will be used to simulate magnetization at echo times only without taking into account spatial encoding gradients (i.e. readout or phase encoding gradients). To simulate the magnetization at other readout times, including phase from spatial encoding gradients, an AbstractTrajectory will be needed as well.

To make a simulator for a particular pulse sequence:

1. Make a struct that's a subtype of either IsochromatSimulator or EPGSimulator. The struct will hold parameters that are necessary for performing the simulations.

2. Add a method to simulate_magnetization! that implements the pulse sequence. For both performance and GPU compatibility, make sure that simulate_magnetization! does not do any heap allocations. Examples for pSSFP and FISP sequences are found in src/sequences.

3. Add methods to output_eltype and output_size that are used to allocate an output array within the simulate function.

4. [Optional] Add a method to show for nicer printing of the sequence in the REPL

5. [Optional] Add a method to getindex to easily reduce the length of the sequence

6. [Optional] Add a constructor for the struct that takes in data from Matlab or something else and assembles the struct

IMPORTANT

The simulate_magnetization! functions (which dispatch on the provided sequence) are assumed to be type-stable and non-allocating Should be possible to achieve when using functions from operators/epg.jlandoperators/isochromat.jl` and a properly parametrized sequence struct.

IsochromatSimulator{T} <: BlochSimulator{T}

Abstract type of which all sequence simulators that are based on the isochromat model will be a subtype. The parameter T should be a number type (e.g. Float64, Float32) and the tissueparameters that are used as input to the simulator should have the same number type.

EPGSimulator{T,Ns} <: BlochSimulator{T}

Abstract type of which all sequence simulators that are based on the EPG model will be a subtype. The parameter T should be a number type (e.g. Float64, Float32) and the tissueparameters that are used as input to the simulator should have the same number type. The parameter Ns corresponds to the maximum order of configuration states that are tracked in the simulations.

output_eltype(::BlochSimulator)

For each <:BlochSimulator, a method should be added to this function that specifies the output type of the simulation. For MR signal simulation, this is typically a complex number representing the transverse magnetization. For other types of simulations, one may want to retrieve the x,y and z components of an isochromat as output (implemented as a FieldVector perhaps) or state configuration matrices Ω.

output_size(::BlochSimulator)

For each <:BlochSimulator, a method should be added to this function that specifies the output size of the simulation for a single ::AbstractTissueProperties.

simulate_magnetization!(magnetization, sequence::BlochSimulator, state, p::AbstractTissueProperties) end

For each <:BlochSimulator, a method should be added to this function that implements the actual pulse sequence using information contained in the sequence struct together with the operators from src/operators/{isochromat,epg}.jl. For performance reasons as well as GPU compatibility it is important that the implementation is type-stable and non-allocating.

Arguments

• magnetization: Pre-allocated array with size(magnetization) = output_size(sequence) and eltype(magnetization) = output_eltype(sequence) to store the output of the simulation.
• sequence: Sequence struct containing fields that are used to implement the actual pulse sequence.
• state: Either an Isochromat or EPGStates, depending on which model is used.
• p: Custom struct (<:AbstractTissueProperties) containing input parameters to the simulation (e.g. T₁T₂)

FISP2D{T, Ns, U<:AbstractVector, V<:AbstractMatrix} <: EPGSimulator{T,Ns}

This struct is used to simulate gradient-spoiled sequence with varying flip angle scheme and adiabatic inversion prepulse using the EPG model. The TR and TE are fixed throughout the sequence. Instantenous RF excitations are assumed. Slice profile correction is done using the partitioned EPG model, where for each flip angle a vector of RF scaling factors are determined prior to the simulation (using, for example, the small tip-angle approximation or Shinnar LeRoux forward model).

Within each TR, a single time steps is used to simulate the RF excitation. Then, in one time step we go from the end of the RF excitation to the echo time (applying T₁ and T₂ decay, T₁ regrowth and B₀ rotation), and again in one time step from the echo time to the start of the next RF excitation.

Fields

• RF_train::U: Vector with flip angle for each TR. abs.(RF_train) are RF flip angles in degrees and angle.(RF_train) are RF phases in radians.
• sliceprofiles::V: Matrix with RF scaling factors (dimensionless) to simulate slice profile effects. Each column represents the (flip angle dependent) scaling factors for one position along the slice direction.
• TR::T: Repetition time in seconds, assumed constant during the sequence
• TE::T: Echo time in seconds, assumed constant during the sequence
• max_state::Val{Ns}: Maximum number of states to keep track of in EPG simulation (dimensionless)
• TI::T: Inversion delay after the inversion prepulse in seconds
• Δk_spoil::T: Spoiling gradient area (in rad/m). Used to calculate diffusion decay time constant T_D = 1/(D * Δk²)

Generic2D{T,V,M,S} where {T<:AbstractFloat, V<:AbstractVector, M<:AbstractMatrix, S} <: IsochromatSimulator{T}

Simulate a generic 2D sequence defined by arrays containing RF and gradient waveforms. Contains a loop over z locations to take into account slice profile effects. The Δt vector stores the time intervals for the waveforms.

Fields

• RF::V{Complex{T}}: Vector with (complex) RF values during each time interval. Units: radians (γ·B₁·Δt)
• GR::M{T}: Matrix with GRx, GRy and GRz values during each time interval. Units: rad·s⁻¹·cm⁻¹ (γ·G, gyromagnetic ratio times gradient strength)
• sample::S: Vector with Bool's to indicate the sample points (dimensionless)
• Δt::V{T}: Vector with time intervals in seconds
• z::V{T}: Vector with different positions along the slice direction in cm (centimeters)

Generic3D{T,V<:AbstractVector{Complex{T}},W<:AbstractVector{T},M<:AbstractMatrix{T},S} <: IsochromatSimulator{T}

Simulate a generic sequence defined by arrays containing RF and gradient waveforms. Unlike the Generic2D sequence, it is assumed that the excitation is homogenous over the voxel and therefore no summation over a slice direction is applied. The Δt vector stores the time intervals for the waveforms.

Fields

• RF::V{Complex{T}}: Vector with (complex) RF values during each time interval. Units: radians (γ·B₁·Δt)
• GR::M{T}: Matrix with GRx, GRy and GRz values during each time interval. Units: rad·s⁻¹·cm⁻¹ (γ·G, gyromagnetic ratio times gradient strength)
• sample::S: Vector with Bool's to indicate the sample points (dimensionless)
• Δt::V{T}: Vector with time intervals in seconds

pSSFP2D{T<:AbstractFloat,N,M,U<:AbstractVector{Complex{T}},V<:Number} <: IsochromatSimulator{T}

This struct is used to simulate a inversion-recovery, gradient-balanced transient-state sequence with varying flip angle scheme based on the isochromat model. The TR and TE are fixed throughout the sequence. The TR and TE are fixed throughout the sequence. Slice profile correction is done by discretizing the RF excitation waveform in time and using multiple Isochromats with different positions along the slice direction (z) per voxel. The sequence also uses an 'α/2' prepulse after the inversion.

Within each TR, multiple time steps are used to simulate the RF excitation. Then, in one time step we go from the end of the RF excitation to the echo time (applying slice refocussing gradient, T₂ decay and B₀ rotation), and again in one time step from the echo time to the start of the next RF excitation.

Fields

• RF_train::U: Vector with flip angle for each TR. abs.(RF_train) are RF flip angles in degrees and angle.(RF_train) are RF phases in radians.
• TR::T: Repetition time in seconds, assumed constant during the sequence
• γΔtRF::SVector{N}{V}: Time-discretized RF waveform, normalized to flip angle of 1 degree. Units: radians
• Δt::NamedTuple{(:ex, :inv, :pr),NTuple{3,T}}: Time intervals in seconds for each sample of excitation pulse (ex), inversion delay (inv) and time between RF and TE (pr)
• γΔtGRz::NamedTuple{(:ex, :inv, :pr),NTuple{3,T}}: Slice select gradients for ex, inv and pr. Units: rad/cm (γ·G·Δt; multiplied by z position in cm gives phase in radians)
• z::SVector{M}{T}: Vector with different positions along the slice direction in cm (centimeters)

pSSFP3D{T<:AbstractFloat,N,U<:AbstractVector{Complex{T}},V<:Number} <: IsochromatSimulator{T}

This struct is used to simulate an inversion-recovery, gradient-balanced, transient-state sequence with varying flip angle scheme based on the isochromat model. The TR and TE are fixed throughout the sequence. The RF excitation waveform can be discretized in time but no slice profile mechanism is provided. The sequence also uses an 'α/2' prepulse after the inversion.

Within each TR, multiple time steps are used to simulate the RF excitation. Then, in one time step we go from the end of the RF excitation to the echo time (applying slice refocussing gradient, T₂ decay and B₀ rotation), and again in one time step from the echo time to the start of the next RF excitation.

Fields

• RF_train::U: Vector with flip angle for each TR. abs.(RF_train) are RF flip angles in degrees and angle.(RF_train) are RF phases in radians.
• TR::T: Repetition time in seconds, assumed constant during the sequence
• γΔtRF::SVector{N}{V}: Time-discretized RF waveform, normalized to flip angle of 1 degree. Units: radians
• Δt::NamedTuple{(:ex, :inv, :pr),NTuple{3,T}}: Time intervals in seconds for each sample of excitation pulse (ex), inversion delay (inv) and time between RF and TE (pr)

AdiabaticInversion{T<:Real, V<:AbstractVector} <: IsochromatSimulator{T}

This struct is used to simulate an adiabatic inversion pulse. This struct itself could be used as field in other sequence structs.

Fields

• γΔtA::V: Time-dependent amplitude modulation. Units: radians (γ·B₁·Δt)
• Δω::V: Time-dependent frequency modulation. Units: rad/s (angular frequency offset)
• Δt::T: Time discretization step in seconds, assumed constant

AbstractTrajectory{T}

The abstract type of which all gradient trajectories will be a subtype. The subtypes should contain fields that can describe the full trajectory during a sequence. The type T refers to the precision of the floating point values within the trajectory struct.

SpokesTrajectory{T} <: AbstractTrajectory{T}

Typical Cartesian and radial trajectories have a lot in common: a readout can be described by a starting point in k-space and a Δk per sample point. To avoid code repetition, both type of trajectories are made a subtype of SpokesTrajectory such that some methods that would be the same for both trajectories otherwise are written for SpokesTrajectory instead.

_get_readout_and_sample_idx(trajectory, t)

Given time index t, compute the associated readout and sample indices (r,s). For trajectories where each readout has the same length, this is equivalent to r,s = fld1(t,ns), mod1(t,ns) with ns being the (constant) number of samples per readout.

nreadouts(::AbstractTrajectory)

For each ::AbstractTrajectory, a method should be added to this function that specifies how many readouts the trajectory consists of.

nsamples(trajectory::AbstractTrajectory)

Determines the total number of samples acquired with the trajectory. Requires nreadouts and nsamplesperreadout to be implemented.

nsamplesperreadout(::AbstractTrajectory, readout_idx)

For each ::AbstractTrajectory, a method should be added to this function that specifies how many samples in total are acquired during the trajectory.

phase_encoding!(magnetization, trajectory::AbstractTrajectory, coordinates::StructArray{<:Coordinates})

For each ::AbstractTrajectory, a method should be added to this function if it does any kind of phase encoding (so far Cartesian only).

to_sample_point(m, trajectory, readout_idx, sample_idx, parameters)

For each ::AbstractTrajectory, a method should be added to this function that, given the magnetization m at the readout with index readout_idx, it computes the magnetization at the readout point with index sample_idx (by applying spatial encoding gradients, T₂ decay, B₀ rotation, etc...) based on the trajectory and parameters.

Arguments

• m: Magnetization at the echo with index readout_idx.
• trajectory: Trajectory struct containing fields that are used to compute the magnetization at other readout times, including the effects of spatial encoding gradients.
• readout_idx: Index that corresponds to the current readout.
• sample_idx: Index for the desired sample during this readout.
• parameters: Tissue parameters of current voxel, including spatial coordinates.

Output:

• mₛ: Magnetization at sample with index sample_idx

CartesianTrajectory2D{T,I,U,V} <: CartesianTrajectory{T}

Struct that is used to implement a typical Cartesian 2D gradient trajectory. The trajectory is described in a compact fashion by only storing the starting position in k-space (k_start_readout) for each readout as well as the step in k-space per readout point Δk_adc.

Note that CartesianTrajectory2D and RadialTrajectory2D are essentially the same when using this compact description. A SpokesTrajectory struct is therefore defined as a supertype of both and methods are defined for SpokesTrajectory instead to avoid code repetition.

The type parameters are intentionally left vague. The J, for example, may be an integer for sequences where each readout has the same number of samples, but for sequences with different numbers of samples per readout it may be a vector of integers.

Fields

• nreadouts::I: The total number of readouts for this trajectory (dimensionless)
• nsamplesperreadout::I: The total number of samples per readout (dimensionless)
• Δt::T: Time between sample points in seconds
• k_start_readout::U: Starting position in k-space for each readout. Units: rad/cm (k-space, where k = γ·G·t)
• Δk_adc::T: k-space step Δkₓ per sample point (same for all readouts). Units: rad/cm
• py::V: Phase encoding index for each readout (dimensionless)
• readout_oversampling::I: Readout oversampling factor (dimensionless)

kspace_coordinates(tr::CartesianTrajectory2D)

Return matrix (nrsamplesperreadout, nrreadouts) with kspace coordinates for the trajectory. Needed for nuFFT reconstructions.

magnetization_to_signal(::Union{CPU1,CPUThreads,CUDALibs}, magnetization, parameters, trajectory::CartesianTrajectory2D, coordinates, coil_sensitivities)

Arguments

• magnetization: Matrix{Complex} of size (# readouts, # voxels) with phase-encoded magnetization at echo times.
• parameters: Tissue parameters of all voxels, including spatial coordinates.
• trajectory: Cartesian trajectory struct.
• coordinates: Vector{Coordinates} with spatial coordinates for each voxel.
• coil_sensitivities: Matrix{Complex} of size (# voxels, # coils) with coil sensitivities.

Returns

• signal: Vector of length (# coils) with each element a Matrix{Complex} of size (# readouts, # samples per readout)

Extended help

As noted in the description of the simulate_signal function (see src/simulate/signal.jl), we simulate the MR signal at timepoint t from coil i as: signalᵢ[t] = sum(m[t,v] * cᵢ[v] * ρ[v] for v in 1:(# voxels)), where cᵢis the coil sensitivity profile of coil i, ρ is the proton density map and m the matrix with the magnetization at all timepoints for each voxel obtained through Bloch simulations.

The output (signalᵢ) for each coil is in principle a Vector{Complex}of length (# samples per readout) * (# readouts). If we reshape the output into aMatrix{Complex}of size (# samples per readout, # readouts) instead, and do something similar form`, then the signal value associated with the s-th sample point of the r-th readout can be expressed as signalᵢ[r,s] = sum( m[r,s,v]] * cᵢ[v] * ρ[v] for v in 1:(# voxels)).

The problem here is that we typically cannot store the full m. Instead, we compute the magnetization at echo times only. The reason is that, if mᵣ is the magnetization at the r-th echo time in some voxel, and E = exp(-ΔtR₂[v]) * exp(im(Δkₓ*x[v])) is the change per sample point (WHICH FOR CARTESIAN SEQUENCES IS THE SAME FOR ALL READOUTS AND SAMPLES), then the magnetization at the s-th sample relative the the echo time can can be computed as mₛ = mᵣ * E[v]^s

Therefore we can write

signalⱼ[r,s] = sum( magnetization[r,v] * E[v]^s * ρ[v] * cⱼ[v] for v in 1:(# voxels)) signalⱼ[r,s] = magnetization[r,:] * (E.^s .* ρ .* cⱼ)

Because the (E.^s .* ρ .* cⱼ)-part is the same for all readouts, we can simply perform this computation for all readouts simultaneously as signalⱼ[:,s] = magnetization * (E.^s .* ρ .* cⱼ)

If we define the matrix Eˢ as E .^ (-(ns÷2):(ns÷2)-1), then we can do the computation for all different sample points at the same time as well using a single matrix-matrix multiplication: signalⱼ = magnetization * (Eˢ .* (ρ .* cⱼ))

The signalⱼ array is of size (# readouts, # samples per readout). We prefer to have it transposed, therefore we compute signalⱼ = transpose(Eˢ .* (ρ .* cⱼ)) * transpose(magnetization) instead.

For the final output, we do this calculation for each coil j and get a vector of signal matrices (one matrix for each coil) as a result.

Note that this implementation relies entirely on vectorized code and works on both CPU and GPU. The matrix-matrix multiplications are - I think - already multi-threaded so a separate multi-threaded implementation is not needed.

sampling_mask(tr::CartesianTrajectory2D)

For undersampled Cartesian trajectories, the gradient trajectory can also be described by a sampling mask.

RadialTrajectory2D{T,I,U,V} <: SpokesTrajectory{T}

Struct that is used to implement a typical radial gradient trajectory. The trajectory can is described in a compact fashion by only storing the starting position in k-space (k_start_readout) for each readout as well as the step in k-space per readout point Δk_adc.

Note that CartesianTrajectory2D and RadialTrajectory2D are essentially the same when using this compact description. A SpokesTrajectory struct is therefore defined as a supertype of both and methods are defined for SpokesTrajectory instead to avoid code repetition.

The type parameters are intentionally left vague. The J, for example, may be an integer for sequences where each readout has the same number of samples, but for sequences with different numbers of samples per readout it may be a vector of integers.

Fields

• nreadouts::I: The total number of readouts for this trajectory (dimensionless)
• nsamplesperreadout::I: The total number of samples per readout (dimensionless)
• Δt::T: Time between sample points in seconds
• k_start_readout::U: Starting position in k-space for each readout. Units: rad/cm (k-space, where k = γ·G·t)
• Δk_adc::U: k-space step Δk between each readout. Units: rad/cm (complex valued: Δkₓ + i·Δkᵧ)
• φ::V: Radial angle for each readout in degrees
• readout_oversampling::I: Readout oversampling factor (dimensionless)

add_gradient_delay!(tr::RadialTrajectory2D, S)

Apply gradient delay to radial trajectory in in-place fashion. The delay is described by the 2x2 matrix S and is assumed to influence the start of the readout only, not the readout direction.

kspace_coordinates(tr::RadialTrajectory2D)

Return matrix (nrsamplesperreadout, nrreadouts) with kspace coordinates for the trajectory. Needed for nuFFT reconstructions.

_allocate_array_on_resource(resource, _eltype, _size)

Allocate an array on the specified resource with the given element type _eltype and size _size. If resource is CPU1() or CPUThreads(), the array is allocated on the CPU. If resource is CUDALibs(), the array is allocated on the GPU. For CPUProcesses(), the array is distributed in the "voxel"-dimension over multiple CPU workers.

This function is called by _allocate_magnetization_array and is not considered part of the public API.

_allocate_magnetization_array(resource, sequence, parameters)

Allocate an array to store the output of the Bloch simulations (per voxel, echo times only) to be performed with the sequence. For each BlochSimulator, methods should have been added to output_eltype and output_size for this function to work properly.

This function is called by simulate_magnetization and is not considered part of the public API.

Returns

• magnetization_array: An array allocated on the specified resource, formatted to store the simulation results for each voxel across the specified echo times.

simulate_magnetization!(magnetization, resource, sequence, parameters)

Simulate the magnetization response for all combinations of tissue properties contained in parameters and stores the results in the pre-allocated magnetization array. The actual implementation depends on the computational resource specified in resource.

This function is called by simulate_magnetization and is not considered part of the public API.

If the tissue properties for a single voxel are provided only, the simulation is performed on the CPU in a single-threaded fashion.

If no resource is provided, the simulation is performed on the CPU in a multi-threaded fashion by default.

However, if the parameters are a CuArray, the simulation is performed on the GPU.

If the parameters are a DArray, the simulation is performed on the multiple workers.

simulate_magnetization(sequence, parameters)

Convenience function to simulate magnetization without specifying the computational resource. The function automatically selects the appropriate resource based on the type of the sequence and parameters. The fallback case is to use multi-threaded CPU computations.

Automatic Resource Selection Rules

• Single AbstractTissueProperties → CPU1() (single-threaded)
• StructArray on CPU → CPUThreads() (multi-threaded, default)
• CuArray → CUDALibs() (GPU)
• DArray → CPUProcesses() (distributed across workers)

simulate_magnetization(resource, sequence, parameters)

Simulate the magnetization response (typically the transverse magnetization at echo times without any spatial encoding gradients applied) for all combinations of tissue properties contained in parameters.

This function can also be used to generate dictionaries for MR Fingerprinting purposes.

Arguments

• resource::AbstractResource: Either CPU1(), CPUThreads(), CPUProcesses() or CUDALibs()
• sequence::BlochSimulator: Custom sequence struct
• parameters::SimulationParameters: Array (typically a StructArray) containing AbstractTissueProperties for each voxel. Ensure tissue properties (like T₁, T₂, B₀) use the units specified in the AbstractTissueProperties docstring (e.g., seconds for relaxation times, Hz for off-resonance).

Note

• If resource == CUDALibs(), the sequence and parameters must have been moved to the GPU using gpu(sequence) and gpu(parameters) prior to calling this function.
• If resource == CPUProcesses(), the parameters must be a DArray with the first dimension corresponding to the number of workers. The function will distribute the simulation across the workers in the first dimension of the DArray.

Returns

• magnetization::AbstractArray: Array of size (output_size(sequence), length(parameters)) containing the magnetization response of the sequence for all combinations of input tissue properties.

_allocate_signal_array(resource, trajectory::AbstractTrajectory, coil_sensitivities)

Allocate an array to store the output of the signal simulation (all readout points, integrated over all voxels).

_signal_per_coil!(signal, resource, magnetization, parameters, trajectory, coordinates, coil_sensitivities)

Compute the signal for a given coil by calculating a volume integral of the transverse magnetization in each voxel for each time point separately (using the signal_at_time_point! function). Each time point is computed in parallel for multi-threaded CPU computation and on the GPU for CUDA computation.

magnetization_to_signal(resource, magnetization, parameters, trajectory, coordinates, coil_sensitivities)

Allocates memory for the signal and computes the signal for each coil separately using the _signal_per_coil! function.

Implementation details

The _signal_per_coil! function has different implementations depending on the computational resources (i.e. the type of resource). The default implementations loop over all time points and compute the volume integral of the transverse magnetization in each voxel for each time point separately. This loop order is not necessarily optimal (and performance may be) across all trajectories and computational resources. If a better implementation is available, add new methods to this function for those specific combinations of resources and trajectories.

The "voxels" are assumed to be distributed over the workers. Each worker computes performs a volume integral over the voxels that it owns only (for all time points) using the CPU1() code. The results are then summed up across all workers.

Note

When using multiple CPU's, the "voxels" are distributed over the workers. Each worker computes the signal for its own voxels in parallel and the results are summed up across all workers.

signal_at_time_point!(signal, time_point, magnetization, parameters, trajectory, coordinates, coil_sensitivities)

If to_sample_point has been defined for the provided trajectory, this (generic but not optimized) function computes the signal at timepoint time_point for all receive coils. It does so by computing the readout- and sample indices for the given time_point, reading in the magnetization at echo time of the r-th readout, using to_sample_point to compute the magnetization at the s-th sample index, and then it integrates over all voxels (while scaling the magnetization with the proper coil sensitivity and proton density).

Better performance can likely be achieved by incorporating more trajectory-specific information together with different loop orders.

simulate_signal(sequence, partitioned_parameters::AbstractVector{<:SimulationParameters})

In situations where the number of voxels is too large to store the intermediate magnetization array, the signal can be calculated in batches: the voxels are divided (by the user) into partitions and the signal is calculated for each partition separately. The final signal is the sum of the signals from all partitions.

When coil sensitivities are not provided, use a single coil with sensitivity = 1 everywhere

simulate_signal(resource, sequence, parameters, trajectory, coordinates, coil_sensitivities)

Simulate the MR signal at timepoint t from coil i as: sᵢ(t) = ∑ⱼ cᵢⱼρⱼmⱼ(t), where cᵢⱼ is the coil sensitivity of coil i at the position of voxel j, ρⱼ is the proton density of voxel j, and mⱼ(t) is the (normalized) transverse magnetization in voxel j obtained through Bloch simulations.

Note that it calls simulate_magnetization to compute the magnetization response at echo times in all voxels, phase_encoding! to apply phase encoding (typically only relevant for Cartesian trajectories), and magnetization_to_signal to compute the signal from the (phase-encoded) magnetization at echo times.

Arguments

• resource::AbstractResource: Either CPU1(), CPUThreads(), CPUProcesses() or CUDALibs()
• sequence::BlochSimulator: Custom sequence struct
• parameters::SimulationParameters: Array (typically a StructArray) containing AbstractTissueProperties for each voxel. Must include proton density fields (ρˣ, ρʸ). Ensure tissue properties (like T₁, T₂, B₀) use the units specified in the AbstractTissueProperties docstring (e.g., seconds for relaxation times, Hz for off-resonance).
• trajectory::AbstractTrajectory: Custom trajectory struct
• coordinates::StructArray{<:Coordinates}: Array with spatial coordinates for each voxel (in cm)
• coil_sensitivities::AbstractMatrix: Sensitivity of coil j in voxel v is given by coil_sensitivities[v,j]

Returns

• signal::AbstractArray{<:Complex}: Simulated MR signal for the sequence and trajectory. The array is of size (# samples per readout, # readouts, # coils).

f32(x)

Change precision of x to Float32. It uses Functors.fmap to recursively traverse the fields of the struct x. For custom structs (e.g. <:BlochSimulator or <:AbstractTrajectory), it is required that typeof(x) be made a Functors.@functor (e.g. @functor FISP).

It may be necessary to add new adapt rules (by adding new methods to adapt_storage) if new structs with complicated nested fields are introduced.

f64(x)

Change precision of x to Float64. It uses Functors.fmap to recursively traverse the fields of the struct x. For custom structs (e.g. <:BlochSimulator or <:AbstractTrajectory), it is required that typeof(x) be made a Functors.@functor (e.g. @functor FISP).

It may be necessary to add new adapt rules (by adding new methods to adapt_storage) if new structs with complicated nested fields are introduced.

_all_arrays_are_cuarrays(x)

Returns true if all AbstractArray fields in x are CuArrays and false otherwise. Will also return if x does not have any AbstractArray fields.

gpu(x)

Move x to CUDA device. It uses Functors.fmap to recursively traverse the fields of the struct x, converting <:AbstractArrays to CuArrays, and ignoring isbitsarrays. For custom structs (e.g. <:BlochSimulator or <:AbstractTrajectory), it is required that typeof(x) be made a Functors.@functor (e.g. @functor FISP).