BlochSim.jl

BlochSim is a Julia package for performing Bloch simulations for MRI.

GAMBAR

Gyromagnetic ratio for ¹H with units Hz/G.

GAMMA

Gyromagnetic ratio for ¹H with units rad/s/G.

bSSFPtuple1

Named tuple tissue parameter keys: (:Mz0, :T1_ms, :T2_ms, :Δf_Hz)

AbstractRF

Abstract type for representing radiofrequency (RF) pulses.

AbstractSpin

Abstract type for representing individual spins.

AbstractSpoiling

Abstract type for representing spoiling.

BalancedRF

Tuple type representing prephasing gradient, RF, and rephasing gradient as used in slice-selective bSSFP excitation.

Bloch3ExpmWorkspace{T<:Real}

Workspace for computing the matrix exponential of a 3×3 Bloch matrix.

BlochDynamicsMatrix(R1, R2, Δω)
BlochDynamicsMatrix{T}()
BlochDynamicsMatrix()

Create a mutable BlochDynamicsMatrix object.

Properties

• R1::Real: Spin-lattice relaxation rate
• R2::Real: Spin-spin relaxation rate
• Δω::Real: Off-resonance frequency

BlochMatrix(a11, a21, a31, a12, a22, a32, a13, a23, a33) # from scalars
BlochMatrix(A::AbstractMatrix) # construct from 3×3 matrix
BlochMatrix{T}() # zeros of type `T`
BlochMatrix() # zeros of type `Float64`

Create a mutable BlochMatrix object representing a fixed-size 3×3 matrix.

Properties

• aij::Real: Matrix entry (i,j) ∈ {1, 2, 3}²

BlochMcConnellDynamicsMatrix(A, E)
BlochMcConnellDynamicsMatrix{T}(N)
BlochMcConnellDynamicsMatrix(N)

Create a BlochMcConnellDynamicsMatrix object with N compartments.

Properties

• A::NTuple{N,BlochDynamicsMatrix{<:Real}}: List of BlochDynamicsMatrixes that make up the main block diagonal of the BlochMcConnellDynamicsMatrix
• E::NTuple{M,ExchangeDynamicsMatrix{<:Real}}: List of ExchangeDynamicsMatrixes that describe exchange between the different compartments; these matrices make up the remaining M = N^2 - N blocks of the BlochMcConnellDynamicsMatrix, sorted by column-major ordering

BlochMcConnellMatrix(A)
BlochMcConnellMatrix{T}(N)
BlochMcConnellMatrix(N)

Create a BlochMcConnellMatrix object with N compartments and representing a fixed-size 3N×3N matrix.

Properties

• A::NTuple{N,NTuple{N,BlochMatrix{<:Real}}}: List of 3×3 matrices that comprise the blocks of the BlochMcConnellMatrix; A[i][j] is the (i,j)th block

Examples

julia> B = BlochMcConnellMatrix(3)
BlochMcConnellMatrix{Float64,3}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> fill!(B, 0.5)

julia> B
BlochMcConnellMatrix{Float64,3}:
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5

ExchangeDynamicsMatrix(r)
ExchangeDynamicsMatrix{T}()
ExchangeDynamicsMatrix()

Create a mutable ExchangeDynamicsMatrix object.

Properties

• r::Real: Exchange rate from one compartment to another

ExcitationMatrix(A)
ExcitationMatrix{T}()
ExcitationMatrix()

Create an ExcitationMatrix object. Multiplying by a MagnetizationMC object has the effect of multiplying each component of the multi-compartment magnetization by the ExcitationMatrix.

Properties

• A::BlochMatrix{<:Real}: Matrix used to describe the excitation

Examples

julia> E = ExcitationMatrix(BlochMatrix(0, 1, 0, 1, 0, 0, 0, 0, 1))
ExcitationMatrix{Int64}:
 0  1  0
 1  0  0
 0  0  1

julia> ExcitationMatrix(reshape(1:9,3,3))
ExcitationMatrix{Int64}:
 1  4  7
 2  5  8
 3  6  9

julia> M = MagnetizationMC((1, 2, 3), (4, 5, 6))
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 1
  My = 2
  Mz = 3
 Compartment 2:
  Mx = 4
  My = 5
  Mz = 6

julia> E * M
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 2
  My = 1
  Mz = 3
 Compartment 2:
  Mx = 5
  My = 4
  Mz = 6

ExcitationWorkspace

Struct for in-place excitation operations.

FreePrecessionMatrix(E1, E2, θ)
FreePrecessionMatrix{T}()
FreePrecessionMatrix()

Create a mutable FreePrecessionMatrix object encoding the effects of relaxation and off-resonance precession.

Properties

• E1::Real: T1 relaxation
• E2::Real: T2 relaxation
• θ::Real: Angle of off-resonance precession (rad)

Examples

julia> T1 = 1000; T2 = 100; Δω = π/600; t = 100;

julia> F = FreePrecessionMatrix(exp(-t / T1), exp(-t / T2) * cos(Δω * t), exp(-t / T2) * sin(Δω * t))
FreePrecessionMatrix{Float64}:
 E1 = 0.9048374180359595
 E2 = 0.31859294158449203
 θ = 0.18393972058572114 rad

julia> Matrix(F)
3×3 Matrix{Float64}:
  0.313219  0.058272  0.0
 -0.058272  0.313219  0.0
  0.0       0.0       0.904837

Gradient(x, y, z)

Create a Gradient object representing x, y, and z B0 gradients. Units are G/cm.

Properties

• x::Real: x gradient
• y::Real: y gradient
• z::Real: z gradient

GradientSpoiling(grad, Tg) <: AbstractSpoiling
GradientSpoiling(gx, gy, gz, Tg)

Represents gradient spoiling, e.g., applying a gradient grad = Gradient(gx, gy, gz) for time Tg (ms). grad can be a Gradient (or gx, gy, and gz can be scalars), representing a constant gradient, or grad can be a collection of Gradients (or gx, gy, and gz can be collections of values), representing a gradient waveform with a constant time step.

IdealSpoiling() <: AbstractSpoiling

Represents ideal spoiling, i.e., setting the transverse (x and y) components of a spin's magnetization to 0.

idealspoiling = IdealSpoilingMatrix()

Matrix representing ideal spoiling. Multiplying by a Magnetization or MagnetizationMC has the effect of setting the x and y components to 0.

Examples

julia> idealspoiling * MagnetizationMC((1, 1, 1), (2, 2, 2))
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 0
  My = 0
  Mz = 1
 Compartment 2:
  Mx = 0
  My = 0
  Mz = 2

InstantaneousRF(α, θ = 0) <: AbstractRF

Represent an idealized instantaneous RF pulse with flip angle α and phase θ.

mese! = MESEBlochSim(TR, TE, nechoes, [rfex, rfref, rephaser, crusher, spoiling])
mese!(spin, [workspace])

Simulate a multi-echo spin echo (MESE) scan on spin, overwriting the spin's magnetization vector. Returns a Vector with the magnetization vectors at each echo.

Arguments

• TR::Real: Repetition time (ms)
• TE::Real: First echo time, and echo spacing (ms); the first echo time is measured from the middle of the excitation pulse
• nechoes::Integer: Number of echoes to readout
• rfex::AbstractRF = InstantaneousRF(π/2): Excitation RF pulse
• rfref::AbstractRF = InstantaneousRF(π, -π/2): Refocussing RF pulse
• rephaser::Union{<:GradientSpoiling,Nothing} = nothing: Slice-select excitation rephasing gradient
• crusher::Union{<:GradientSpoiling,Nothing} = nothing: Crusher gradient (placed on either side of each refocussing pulse)
• spoiling::Union{IdealSpoiling,<:GradientSpoiling,Nothing} = IdealSpoiling(): Type of spoiling to apply

workspace isa MESEBlochSimWorkspace.

Magnetization(x, y, z)
Magnetization{T}()
Magnetization()

Create a mutable Magnetization object representing a 3D magnetization vector.

Properties

• x::Real: x component of magnetization vector
• y::Real: y component of magnetization vector
• z::Real: z component of magnetization vector

Examples

julia> M = Magnetization()
Magnetization vector with eltype Float64:
 Mx = 0.0
 My = 0.0
 Mz = 0.0

julia> M.x = 1; M.y = 2; M.z = 3; M
Magnetization vector with eltype Float64:
 Mx = 1.0
 My = 2.0
 Mz = 3.0

julia> Magnetization(1:3)
Magnetization vector with eltype Int64:
 Mx = 1
 My = 2
 Mz = 3

MagnetizationMC((x1, y1, z1), ...)
MagnetizationMC{T}(N)
MagnetizationMC(N)

Create a MagnetizationMC object representing N 3D magnetization vectors in the same physical location.

One can access the ith component magnetization vector by indexing the MagnetizationMC object. Furthermore, iterating the MagnetizationMC object iterates through each of the component magnetization vectors.

Properties

• M::NTuple{N,Magnetization{<:Real}}: List of component magnetization vectors

Examples

julia> M = MagnetizationMC((1, 2, 3), (4, 5, 6))
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 1
  My = 2
  Mz = 3
 Compartment 2:
  Mx = 4
  My = 5
  Mz = 6

julia> M[2]
Magnetization vector with eltype Int64:
 Mx = 4
 My = 5
 Mz = 6

julia> foreach(m -> (m.x = 0; m.y = 1; m.z = 2), M)

julia> M
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 0
  My = 1
  Mz = 2
 Compartment 2:
  Mx = 0
  My = 1
  Mz = 2

Position(x, y, z)

Create a mutable Position object representing a 3D location. Units are cm.

Properties

• x::Real: x position
• y::Real: y position
• z::Real: z position

RF(waveform, Δt, [Δθ], [grad]) <: AbstractRF

Represent an RF pulse with the given (possibly complex-valued) waveform (G) and time step Δt (ms).

Options

• Δθ additional phase added to the waveform (defaults to 0 radians)
• grad B0 gradient that is turned on during the RF pulse (defaults to Gradient(0, 0, 0), i.e., turned off).

Properties

• α::Vector{<:Real}: Instantaneous flip angles (rad) at each time point; computed from the magnitude of waveform
• θ::Vector{<:Real}: Instantaneous phase (rad) at each time point; computed from the phase of waveform
• Δt::Real: Time step (ms)
• Δθ_initial::Real: Phase added to θ before any phase-cycling increment has been applied
• Δθ::Ref{<:Real}: Phase to be added to θ; can be updated to simulate phase-cycling/RF spoiling
• grad: Gradient applied during the RF pulse
  • ::Gradient: Constant gradient
  • ::Vector{<:Gradient}: Time-varying gradient

RFSpoiling(Δθ = 117°) <: AbstractSpoiling

Represents RF spoiling, i.e., quadratically incrementing the phase of the RF pulses from TR to TR.

RFandGradientSpoiling(grad_spoiling, rf_spoiling) <: AbstractSpoiling

Represents both RF and gradient spoiling.

RectRF(duration_ms, [α = π/2], [θ = 0], [grad = zero(Gradient)]) <: AbstractRF

Represent a rectangular RF pulse with

• duration_ms duration (ms)
• α flip angle (radians)
• θ phase (radians)
• grad constant B0 gradient 3-vector (G/cm) that is turned on during the RF pulse (defaults to Gradient(0, 0, 0), i.e., turned off).

spgr! = SPGRBlochSim(TR, TE, rf, [spoiling], [nTR], [save_transients])
spgr!(spin, [workspace])

Simulate a spoiled gradient-recalled echo (SPGR) scan on spin, overwriting the spin's magnetization vector. The resultant magnetization is stored in spin.M. If nTR > 0 and save_transients === true, then spgr!(...) returns a Vector with the magnetization vectors at the echo time for each of the nTR simulated TRs.

Arguments

• TR::Real: Repetition time (ms)
• TE::Real: Echo time (ms)
• rf:
  • ::Real: Excitation flip angle (rad) (assumes an instantaneous RF pulse)
  • ::AbstractRF: Excitation RF pulse
• spoiling::AbstractSpoiling = IdealSpoiling(): Type of spoiling to apply
• nTR::Val = Val(0): Number of TRs to simulate; Val(0) indicates to simulate a steady-state scan
• save_transients::Val = Val(false): Whether or not to return the magnetization vectors at the TE for each of the nTR simulated TRs; does nothing if nTR == 0

workspace isa SPGRBlochSimWorkspace.

Spin([M], M0, T1, T2, Δf, [pos]) <: AbstractSpin

Create an object that represents a single spin.

Properties

• M::Magnetization = Magnetization(0, 0, M0): Magnetization vector
• M0::Real: Equilibrium magnetization
• T1::Real: Spin-lattice recovery time constant (ms)
• T2::Real: Spin-spin recovery time constant (ms)
• Δf::Real: Off-resonance frequency (Hz)
• pos::Position = Position(0, 0, 0): Spatial location (cm)
• N::Int = 1: Number of compartments

Examples

julia> Spin(1, 1000, 100, 0, Position(1, 2, 3))
Spin{Float64}:
 M = Magnetization(0.0, 0.0, 1.0)
 M0 = 1.0
 T1 = 1000.0 ms
 T2 = 100.0 ms
 Δf = 0.0 Hz
 pos = Position(1.0, 2.0, 3.0) cm

julia> Spin(Magnetization(1, 2, 3), 1, 1000, 100, 0)
Spin{Float64}:
 M = Magnetization(1.0, 2.0, 3.0)
 M0 = 1.0
 T1 = 1000.0 ms
 T2 = 100.0 ms
 Δf = 0.0 Hz
 pos = Position(0.0, 0.0, 0.0) cm

SpinMC([M], M0, frac, T1, T2, Δf, τ, [pos]) <: AbstractSpin

Create an object that represents a single spin with multiple compartments.

Properties

• M::MagnetizationMC = Meq: Magnetization vector
• Meq::MagnetizationMC = MagnetizationMC((0, 0, frac[1] * M0), ...): Equilibrium magnetization vector
• M0::Real: Equilibrium magnetization
• frac::Tuple{<:Real}: Water fraction of each compartment
• T1::Tuple{<:Real}: Spin-lattice recovery time constants (ms)
• T2::Tuple{<:Real}: Spin-spin recovery time constants (ms)
• Δf::Tuple{<:Real}: Off-resonance frequencies (Hz)
• r::Tuple{Tuple{<:Real}}: Exchange rates (1/ms); r[i][j] is the exchange rate from compartment i to compartment j
• pos::Position = Position(0, 0, 0): Spatial location (cm)
• N::Int = length(frac): Number of compartments

Note

The SpinMC constructor takes τ (inverse exchange rate, or residence time) as input, not r. Furthermore, τ is given as a Tuple with N^2 - N elements, arranged like (τ12, τ13, ..., τ1N, τ21, τ23, ..., τ2N, ...).

Examples

julia> SpinMC(1, (0.2, 0.8), (400, 1000), (20, 100), (15, 0), (100, 25))
SpinMC{Float64,2}:
 M = MagnetizationMC((0.0, 0.0, 0.2), (0.0, 0.0, 0.8))
 M0 = 1.0
 frac = (0.2, 0.8)
 T1 = (400.0, 1000.0) ms
 T2 = (20.0, 100.0) ms
 Δf = (15.0, 0.0) Hz
 r = ((0.0, 0.01), (0.04, 0.0)) 1/ms
 pos = Position(0.0, 0.0, 0.0) cm

bSSFPmode

A type used to control how bSSFP signal is calculated.

• bSSFPbloch use BlochSim matrix computations (default)
• bSSFPbloch3 use analytical 3×3 matrix computations
• bSSFPellipse use ellipse model for 1-pool

PadeApproximantOfDegree!(expA, A, workspace, m)

Pade approximant to exponential.

Based on PADEAPPROXIMANTOFDEGREE F = PADEAPPROXIMANTOFDEGREE(M) is the degree M diagonal Pade approximant to EXP(A), where M = 3, 5, 7, 9 or 13. Series are evaluated in decreasing order of powers, which is in approx. increasing order of maximum norms of the terms.

rf = RF1(α_rad, tRF_ms, θ = 0, args...; nsamp=1, kwargs...)

RF "rectangular" pulse of duration tRF_ms for flip angle α_rad and phase θ (in radians), represented by a nsamp-sample constant "waveform" with phase θ.

_TE_ms(TE_ms, TR_ms, [rf])

Helper to convert TE symbols to a number

• Val{:midTR} for the usual TR/2 choice
• Val{:postRF} immediately after RF pulse ends

absolutesum(A)

Compute the sum of the absolute values of the elements of A.

add!(C, A, B)

Compute A + B, storing the result in C.

add!(A, B)

Compute A + B, storing the result in A.

applydynamics!(spin, BtoM, A, [B])

Apply dynamics to the given spin, overwriting the spin's magnetization vector. BtoM is used to store intermediate results (and is thus overwritten).

Examples

julia> s = Spin(1, 1000, 100, 3.75); s.M
Magnetization vector with eltype Float64:
 Mx = 0.0
 My = 0.0
 Mz = 1.0

julia> BtoM = Magnetization();

julia> (A,) = excite(s, InstantaneousRF(π/2)); applydynamics!(s, BtoM, A); s.M
Magnetization vector with eltype Float64:
 Mx = 1.0
 My = 0.0
 Mz = 6.123233995736766e-17

julia> (A, B) = freeprecess(s, 100); applydynamics!(s, BtoM, A, B); s.M
Magnetization vector with eltype Float64:
 Mx = -0.2601300475114444
 My = -0.2601300475114445
 Mz = 0.09516258196404054

b1 = b1_gauss(α_rad, tRF_ms)

Return finite-duration (rectangular) RF pulse amplitude

In

• α_rad tip angle (radians)
• tRF_ms pulse length (ms)

Notes:

• GAMMA has units rad/s/G
• Tip angle for constant pulse: α_rad = GAMMA * b1_gauss * tRF_s
• so b1_gauss = α_rad / GAMMA / tRF_s

 bssfp(bSSFPbloch3, tRF_ms, args...)

1-pool version for a constant RF of duration tRF_ms. (Account for gradient effects in Δf_Hz.)

bssfp(bSSFPellipse, Mz0, T1_ms, T2_ms, Δf_Hz,
   TR_ms, TE_ms, Δϕ_rad, α_rad, θ_rf_rad=0)

Elliptical signal model for bSSFP. This is the analytical solution to the 1-pool bSSFP signal with instantaneous RF.

See other bssfp methods for argument definitions.

Xiang et al. MRM 2014; https://doi.org/10.1002/mrm.25098

Keskin et al. IEEE T-MI 2022; https://doi.org/10.1109/TMI.2021.3102852

bssfp(spin, spin_no_rf_phase_fact, TR_ms, TE_ms, α_rad, θ_rf_rad=0)
bssfp(spin, spin_no_rf_phase_fact, TR_ms, TE_ms, rf)

Return steady-state magnetization signal value at the echo time for a phase-cycled bSSFP sequence using BlochSim.

Ref: Murthy, N., Nielsen, J. F., Whitaker, S. T., Haskell, M. W., Swanson, S. D., Seiberlich, N., & Fessler, J. A. (2022). Quantifying myelin water exchange using optimized bSSFP sequences. In Proc. Intl. Soc. Mag. Res. Med (p. 2068). [2]

In (tissue)

• 'spin' multi-compartment spin object with RF phase cycling factor
• 'spinnorfphasefact' multi-compartment spin object without RF phase cycling factor

In (scan)

• TR_ms repetition time (ms)
• TE_ms echo time (ms)
• α_rad flip angle of RF pulse (radians)
• θ_rf_rad phase angle of RF pulse (radians) [default 0]

Or instead provide

• rf::AbstractRF

Out

• signal steady-state magnetization (as a complex number)

bssfp(Mz0, T1_ms, T2_ms, Δf_Hz, TR_ms, TE_ms, Δϕ_rad, α_rad, θ_rf_rad=0)
bssfp(Mz0, T1_ms, T2_ms, Δf_Hz, TR_ms, TE_ms, Δϕ_rad, rf, [pos])
bssfp(spin, Δf_Hz, TR_ms, TE_ms | Val(:midTR) | Val(:postRF), Δϕ_rad, rf)

Return steady-state magnetization signal value at the given echo time for a single-pool tissue for a bSSFP sequence using BlochSim.

This generalizes Hargreaves et al., MRM 2001 by accounting for possibly finite RF duration.

In (tissue parameters):

• Mz0 initial condition for magnetization in the z-direction (constant)
• T1_ms MRI tissue parameter for T1 relaxation (ms)
• T2_ms MRI tissue parameter for T2 relaxation (ms)
• Δf_Hz off-resonance value (Hz) Use Val(:midTR) for TE = TR/2; Val(:postRF) for TE = tRF/2

In (scan parameters):

• TR_ms repetition time (ms)
• TE_ms echo time (ms), measured from middle of RF pulse
• Δϕ_rad RF phase cycling increment (radians)
• α_rad flip angle of RF pulse (radians)
• θ_rf_rad phase of RF pulse (radians)

Or, instead of α_rad and θ_rf_rad, provide:

• rf::AbstractRF, e.g., InstantaneousRF(α_rad, θ_rf_rad)

Option:

• pos::Position = Position(0.0, 0.0, 0.0) option if rf specified

Out

• signal steady-state transverse magnetization (as a complex number)

bssfp(...)

Compute bSSFP 2-pool signal value. Version with scalar arguments (convenient for autodiff).

In tissue:

• M0_phase # radians
• Mz0
• f_f
• T1_f_ms
• T1_s_ms
• T2_f_ms
• T2_s_ms
• τ_fs_ms
• Δff_Hz fast component frequency shift

In: system parameters (sometimes known):

• Δf0_Hz # B0 off resonance

In: scan parameters (always known):

• Δϕ_rad # RF phase cycling value (radians)
• TR_ms
• args.. remaining scan arguments, e.g., TE_ms...

bssfp(spin, TR_ms, TE_ms, rf::AbstractRF)

Classic version with no phase cycling increment, for InstantaneousRF only.

bssfp(spin, TR_ms, TE_ms, Δϕ_rad, rf::AbstractRF | rf::Tuple)

Signal accounting for phase cycling increment Δϕ_rad, allowing for finite duration rf pulse, or a 3-tuple combination of a

• prephasing gradient
• (typically slice selective) RF pulse
• rephasing gradient, provided as

Tuple{<:GradientSpoiling, <:AbstractRF, <:GradientSpoiling}. Typically the prephasing and rephasing are the same gradients.

Echo time is measured from the center of the RF pulse.

combine!(A, B, A1, B1, A2, B2)
combine!(A, A1, A2)
A, B = combine(A1, B1, A2, B2)

Combine the matrices and vectors that describe the dynamics of a spin into one matrix and one vector. The mutating versions combine! overwrite A and B (when relevant).

The dynamics described by A1 and B1 apply first, then those described by A2 and B2. In other words: A = A2 * A1 and B = A2 * B1 + B2.

Examples

julia> s = Spin(1, 1000, 100, 3.75);

julia> A = BlochMatrix(); B = Magnetization();

julia> (A1, B1) = excite(s, InstantaneousRF(π/2));

julia> (A2, B2) = freeprecess(s, 100);

julia> combine!(A, B, A1, B1, A2, B2); A * s.M + B
Magnetization vector with eltype Float64:
 Mx = -0.2601300475114444
 My = -0.2601300475114445
 Mz = 0.09516258196404054

A, B = combine(A1, B1, A2, B2)

See combine!.

crb(signal, x, σ)

Compute Cramer-Rao Bound (CRB) for estimating parameter vector x ∈ ℝᴺ from the noisy signal y = signal(x) + ϵ ∈ ℝᴹ for M ≥ N where ϵ denotes additive white Gaussian noise with standard deviation σ.

The output of signal(x) must be real! For complex signals, use the helper function Blochsim.real_imag.

cross!(v::Vector, a::AbstractVector, b::AbstractVector)

Compute cross-product v = a ⊗ b for 3-vectors; avoids allocations by mutating input v.

div!(A, a)

Compute A / a, storing the result in A.

duration(::GradientSpoiling) Duration of a GradientSpoiling "waveform" in ms.

duration(rf)

Return the duration (ms) of the RF pulse.

Caution. If the RF waveform has a long tail of zeros due to the presence of a built-in rephasing gradient, then this duration value includes that tail.

duration(rf)

Return the duration (ms) of the RF pulse.

(λ, V) = eigen_bloch3!(λ, V, row1, row2, row3, r1, r2, w, s, c)

Return eigendecomposition (λ, V) of 3×3 Bloch matrix A = [-r2 w s; -w -r2 c; -s -c -r1] Mutates λ, V and workspace row1 row2 row3, so allocation free.

eigvals_bloch3(r1, r2, w, s, c)

Return tuple of Complex-typed eigenvalues of 3×3 Bloch matrix A = [-r2 w s; -w -r2 c; -s -c -r1] using Cardano/trig formulas, withcbrt` called only on real arguments.

eigvec_3by3!(v, row1, row2, row3, [rtol = 1e-12])

Compute eigenvector v by crossing two rows; switch row pairs if needed. Mutates v and workspace row1 row2 row3

eigvec_bloch3!(v, row1, row2, row3, r1, r2, w, s, c, λ)

Compute eigenvector v for eigenvalue λ by crossing two rows; switch row pairs if needed. Mutates v and workspace row1 row2 row3

excite(spin, rf::InstantaneousRF, [nothing])
excite(spin, rf::RF, [workspace])

Simulate excitation for the given spin. Returns (A, B) such that A * M + B applies excitation to the magnetization M. If isnothing(B) (as is the case for InstantaneousRFs), then A * M applies excitation to M.

For RF objects, workspace isa ExcitationWorkspace. For SpinMC objects, use workspace = ExcitationWorkspace(spin, nothing) to use an approximate matrix exponential to solve the Bloch-McConnell equation.

For an in-place version, see excite!.

Examples

julia> s = Spin(1, 1000, 100, 3.75); s.M
Magnetization vector with eltype Float64:
 Mx = 0.0
 My = 0.0
 Mz = 1.0

julia> (A,) = excite(s, InstantaneousRF(π/2, π/4)); A * s.M
Magnetization vector with eltype Float64:
 Mx = 0.7071067811865476
 My = -0.7071067811865475
 Mz = 6.123233995736766e-17

excite!(spin, ...)

Apply excitation to the given spin, overwriting the spin's magnetization vector.

excite!(A, [nothing], spin, rf::InstantaneousRF, [nothing])
excite!(A, B, spin, rf::RF, [workspace])

Simulate excitation, overwriting A and B (in-place version of excite).

excite!(A, B, spin, rf, workspace)

Mutate A and B for an RF pulse described by a vector of samples with time step rf.Δt.

The method used here treats each sample using a "cascade" of three simpler steps:

• free precession for Δt/2
• instantaneous RF rotation
• free precession for Δt/2

See Gras et al. MRM Jul. 2018 https://doi.org/10.1002/mrm.27001. The accuracy of this "cascade" approximation depends on Δt, as explored in one of the repo demos.

(A1, b1) = excite(spin::Spin, rf::RectRF)

Version for a RectRF pulse of duration rf.duration.

excite_bloch3!(...)

Mutating version of excite_bloch3

(A1, b1) = excite_bloch3(spin::Spin, rf::InstantaneousRF)

Caution: used only for testing.

(A1, b1) = excite_bloch3(spin::Spin, rf::RF)

Version for a constant (rect) RF pulse of duration rf.Δt, currently represented by a single sample (often complex) RF waveform.

(A1, b1) = excite_bloch3(r1, r2, w, s, c, t)

Return solution to Bloch equation A1 = exp(A*t) and b1 = (A1 - I) * (A^{-1} * d) for the 3×3 Bloch matrix A = [-r2 w s; -w -r2 c; -s -c -r1] and d = [0, 0, r1] (i.e., for M₀ = 1).

expm!(expA, A, [workspace])

Compute the matrix exponential of A, storing it in expA.

workspace isa MatrixExponentialWorkspace.

expA = expm(A, [workspace])

Return the matrix exponential of BlochMcConnellDynamicsMatrix A, where workspace isa MatrixExponentialWorkspace.

expm_bloch3!(expAt, work, r1, r2, w, s, c, t)

Return exp(A*t) for the 3×3 Bloch matrix A = [-r2 w s; -w -r2 c; -s -c -r1]. Uses explicit eigendecomposition based on analytical roots of cubic characteristic polynomial. Mutates expAt and work::Bloch3ExpmWorkspace. Allocates very little memory; just some lu! overhead.

Caution. There is a π/2 difference between this code and excite! as seen in one of the demos.

expm_bloch3(r1, r2, w, s, c, t)

Return exp(A*t) for the 3×3 Bloch matrix A = [-r2 w s; -w -r2 c; -s -c -r1]. Uses explicit eigendecomposition based on analytical roots of cubic characteristic polynomial.

freeprecess(spin, t, [nothing])
freeprecess(spinmc, t, [workspace])
freeprecess(spin, t, grad, [nothing])
freeprecess(spinmc, t, grad, [workspace])

Simulate free-precession for the given spin for time t ms, optionally in the presence of a B0 gradient. Returns (A, B) such that A * M + B applies free-precession to the magnetization M.

For SpinMC objects, workspace isa BlochMcConnellWorkspace. Pass in nothing instead to use an approximate matrix exponential to solve the Bloch-McConnell equation.

For an in-place version, see freeprecess!.

Examples

julia> s = Spin(Magnetization(1, 0, 0), 1, 1000, 100, 3.75)
Spin{Float64}:
 M = Magnetization(1.0, 0.0, 0.0)
 M0 = 1.0
 T1 = 1000.0 ms
 T2 = 100.0 ms
 Δf = 3.75 Hz
 pos = Position(0.0, 0.0, 0.0) cm

julia> (A, B) = freeprecess(s, 100); A * s.M + B
Magnetization vector with eltype Float64:
 Mx = -0.2601300475114444
 My = -0.2601300475114445
 Mz = 0.09516258196404048

julia> s = Spin(Magnetization(1, 0, 0), 1, 1000, 100, 0, Position(0, 0, 3.75))
Spin{Float64}:
 M = Magnetization(1.0, 0.0, 0.0)
 M0 = 1.0
 T1 = 1000.0 ms
 T2 = 100.0 ms
 Δf = 0.0 Hz
 pos = Position(0.0, 0.0, 3.75) cm

julia> (A, B) = freeprecess(s, 100, Gradient(0, 0, 1/GAMBAR)); A * s.M + B
Magnetization vector with eltype Float64:
 Mx = -0.2601300475114444
 My = -0.2601300475114445
 Mz = 0.09516258196404048

freeprecess!(A, B, t, M0, T1, T2, Δf)
freeprecess!(A, B, spin, t, [nothing])
freeprecess!(A, B, spinmc, t, [workspace])
freeprecess!(A, B, spin, t, grad, [nothing])
freeprecess!(A, B, spinmc, t, grad, [workspace])

Simulate free-precession, overwriting A and B (in-place version of freeprecess).

freeprecess!(spin, ...)

Apply free-precession to the given spin, overwriting the spin's magnetization vector.

freeprecess(t, M0, T1, T2, Δf)

Simulate free-precession, i.e., relaxation and off-resonance precession. Returns (A, B) such that A * M + B applies free-precession to the magnetization M.

For an in-place version, see freeprecess!.

Arguments

• t::Real: Duration of free-precession (ms)
• M0::Real: Equilibrium magnetization
• T1::Real: Spin-lattice recovery time constant (ms)
• T2::Real: Spin-spin recovery time constant (ms)
• Δf::Real: Off-resonance frequency (Hz)

Examples

julia> (A, B) = freeprecess(100, 1, 1000, 100, 3.75); A * Magnetization(1, 0, 0) + B
Magnetization vector with eltype Float64:
 Mx = -0.2601300475114444
 My = -0.2601300475114445
 Mz = 0.09516258196404048

get_Δf_tuple(Δϕ_rad, Δf0_Hz, Δf_myelin_Hz, TR_ms)

Account for phase cycling increments as an effective frequency shift.

In:

• Δϕ_rad RF phase cycling value (radians)
• Δf0_Hz off-resonance value (Hz)
• Δf_myelin_Hz # additional off-resonance value experienced only by fast (myelin) water (Hz)
• TR_ms repetition time (ms)

Out:

• Δf_tuple_Hz tuple with off-resonance values for fast and slow compartments

Hinshaw, J. Appl. Phys. 1976.

get_τ_tuple(τ_fs_ms, f_f)

In:

• τ_fs_ms residence time for exchange from fast T2 pool (myelin water) to slow T2 pool (non-myelin water) (ms)
• f_f fast fraction (myelin fraction)

Out:

• τ_tuple_ms tuple with fast-to-slow and slow-to-fast residence times

gradient_frequency(grad, pos)

Compute the off-resonance frequency in Hz induced by the given B0 gradient grad at position pos.

gz = gz_sinc(tRF_ms, nlobe, slice_width)

Slice-selective gradient amplitude (in G/cm)

From Fourier analysis (small tip-angle approximation): 2π / (tRF_ms/2nlobe) = GAMMA * gz * slice_width

matrix_bloch3!(mat, r1, r2, w, s, c)

Mutating version of matrix_bloch3()

matrix_bloch3(r1, r2, w, s, c)

Return 3×3 Bloch matrix A = [-r2 w s; -w -r2 c; -s -c -r1]

muladd!(C, A, B)

Compute A * B + C, storing the result in C.

real_imag(y)

Return the stack [real(y); imag(y)] that is needed for CRB because ForwardDiff cannot handle complex signal values.

rf_gauss(rf::AbstractRF)

Return complex RF waveform in Gauss from rf object.

rf_sinc(α_rad, tRF_ms, Δt_ms, shape, nlobe)

Make truncated sinc RF waveform.

rf, rephasing = rf_slice(tRF_ms = 1 ; kwargs...)
rf = rf_slice( Val(:built_in_rephasing), tRF_ms = 1 ; kwargs...)

Make RF pulse for slice selection (with constant gradient along z direction), and rephasing gradient of length tRF_ms/2.

The version with Val(:built_in_rephasing) pads the RF waveform with zeros and includes the rephasing gradient; it is just for testing! It is not recommended because it is less efficient computationally and its duration(rf) = (3/2) * tRF_ms,

Caution: the bipolar gradient waveform used here ignores slew-rate constraints.

In

• Val(:built_in_rephasing) (see above)
• tRF_ms duration in ms of RF portion of the pulse; default 1 [ms] (excluding rephasing gradient)

Option

• α_rad flip angle; default π/2 radian
• Δt_ms dwell time (sampling period); default 1e-3 for 1μs
• shape pulse shape; default :sinc
• nlobe number of sinc lobes; default 5
• slice_width default 1 cm
• rephasing::Bool include rephasing gradient? default !isinf(slice_width)
• Δθ phase passed to RF; default 0 radian

Out

• rf a BlochSim.RF object

rfspoiling_increment(spoiling)

Return the quadratic phase increment used for RF spoiling.

rotatetheta(α::Real = π/2, θ::Real = 0)

Return a 3 × 3 BlochMatrix for (left-handed) flip angle α about an axis in the x-y plane that makes (left-handed) angle θ with the negative y-axis.

For an in-place version, see rotatetheta!.

Examples

julia> BlochSim.rotatetheta() * Magnetization(0, 0, 1) # Rotate towards positive x-axis
Magnetization vector with eltype Float64:
 Mx = 1.0
 My = 0.0
 Mz = 6.123233995736766e-17

julia> BlochSim.rotatetheta(π/2, π/2) * Magnetization(0, 0, 1) # Rotate towards negative y-axis
Magnetization vector with eltype Float64:
 Mx = 6.123233995736766e-17
 My = -1.0
 Mz = 6.123233995736766e-17

rotatetheta!(A, α, θ)

Simulate left-handed rotation by angle α about an axis in the x-y plane that makes left-handed angle θ with the negative y-axis, overwriting A.

This function is an in-place version of rotatetheta.

signal(spin)
signal(M)

Return the signal detected from the given spin or magnetization vector.

Examples

julia> signal(Spin(Magnetization(1, 2, 3), 1, 1000, 100, 0))
1.0 + 2.0im

julia> signal(MagnetizationMC((1, 2, 3), (1, 1, 1)))
2 + 3im

snr2sigma(snr_db, s::AbstractArray)

Convert SNR in dB to noise standard deviation σ for a complex-valued signal s to which complex-valued noise will be added.

spoil(spin, spoiling, [nothing])
spoil(spinmc, spoiling, [workspace])

Simulate gradient or ideal spoiling for the given spin. Returns (A, B), such that A * M + B applies spoiling to the magnetization M. If B is nothing (as is the case for IdealSpoiling), then A * M applies spoiling, and if both A and B are nothing (as is the case for RFSpoiling) then there is no spoiling.

For SpinMC objects and for GradientSpoiling and RFandGradientSpoiling, workspace isa BlochMcConnellWorkspace. Pass in nothing instead to use an approximate matrix exponential to solve the Bloch-McConnell equation.

Note

This function only simulates gradient or ideal spoiling, not RF spoiling. RF spoiling must be implemented by updating the phase of the RF pulse(s) in your sequence from TR to TR.

For an in-place version, see spoil!. ```

spoil!(A, B, spin, spoiling, [nothing])
spoil!(A, B, spinmc, spoiling, [workspace])

Simulate gradient or ideal spoiling, overwriting A and B (in-place version of spoil).

spoil!(spin, spoiling; A::FreePrecessionMatrix, B::Magnetization)

Mutates spin by applying effects of spoiling (which could be a refocusing gradient, for example). Provide optional (mutated) arguments A and B to avoid allocations.

spoil!(spin)

Apply ideal spoiling to the given spin.

spoiler_gradient(spoiling)

Get the Gradient object used for gradient spoiling.

spoiler_gradient_duration(spoiling)

Return the duration of the spoiler gradient (ms).

subtract!(A, B)

Compute A - B, storing the result in A.

subtract!(C, A, B)

Compute A - B, storing the result in C.

subtractmul!(C, X, A, B)

Compute (X - A) * B, storing the result in C.

subtractmuladd!(C, X, A, B)

Compute (X - A) * B + C, storing the result in C.

mul!(A, a)

Compute A * a, storing the result in A.

mul!(C, A, B)

Compute A * B, storing the result in C.