General Comments:

• nM, number of spins, as magnetic spin vectors are often denoted as 𝑀.
• nT, number of steps/time-points.

Unitful.jl related:

• 𝐁 = 𝐌*𝐈^-1*𝐓^-2, dimension of magnetic field strength.
• 𝐅 = 𝐓^-1, dimension of temporal frequency.
• 𝐊 = 𝐋^-1, dimension of spatial frequency.
• 𝚪 = 𝐅/𝐁, dimension of gyro ratio.

Some utilities functions routinely used in MR simulations.

Some steady state properties of common sequences.

Analytical expressions of common steady states sequences signals.

Tools for simulating RF spoiling in gradient echo sequences.

const γ¹H = 4257.6u"Hz/Gauss"

Gyromagnetic ratio of water proton.

An abstract type for pulses.

See also: Pulse.

This type keeps the essentials of magnetic spins. Its instance struct must contain all fields listed listed in the exemplary struct SpinArray.

Misc

Might make AbstractSpinArray <: AbstractArray in a future version

See also: SpinArray, AbstractSpinCube.

UNDER CONSTRUCTION

AbstractSpinBolus <: AbstractSpinArray. This type inherits AbstractSpinArray as a field. Its instance struct must contain all fields listed in the exemplary struct SpinBolus.

See also: AbstractSpinArray, SpinBolus.

AbstractSpinCube <: AbstractSpinArray. This type inherits AbstractSpinArray as a field. Its instance struct must contain all fields listed in the exemplary struct SpinCube.

See also: AbstractSpinArray, SpinCube.

B0D{T<:Real} = Quantity{T, 𝐁}

Type of magetic field strength. Based on Unitful.Quantity.

Examples:

julia> (1u"Gauss")::B0D
1 Gauss

F0D{T<:Real} = Quantity{T, 𝐅}

Type of temporal frequency. Based on Unitful.Quantity.

Examples:

julia> (1u"s^-1")::F0D
1 s^-1

GR0D{T<:Real} = Quantity{T, 𝐁/𝐋}

Type of magnetic gradient. Based on Unitful.Quantity.

Examples:

julia> (1u"Gauss/cm")::GR0D
1 Gauss cm^-1

K0D{T<:Real} = Quantity{T, 𝐊}

Type of spatial frequency. Based on Unitful.Quantity.

Examples:

julia> (1u"cm^-1")::K0D
1 cm^-1

L0D{T<:Real} = Quantity{T, 𝐋}

Type of length. Based on Unitful.Quantity.

Examples:

julia> (1u"cm")::L0D
1 cm

A struct for typical MR pulses: Pulse <: AbstractPulse.

Fields:

Mutable:

• rf::TypeND(RF0D, [1,2]) (nT,) or (nT, nCoils).
• gr::TypeND(GR0D, [2]) (nT, 3), where 3 accounts for x-y-z channels.
• dt::T0D (1,), simulation temporal step size, i.e., dwell time.
• des::String, an description of the pulse to be constructed.

See also: AbstractPulse.

Pulse(rf, gr; dt=(4e-6)u"s", des="generic pulse")

Create Pulse object with prescribed parameters.

RF0D{T<:Complex} = Quantity{T, 𝐁}

Type of magnetic RF. Based on Unitful.Quantity.

Examples:

julia> ((1+1im)u"Gauss")::RF0D
1 + 1im Gauss

SL0D{T<:Real} = Quantity{T, 𝐁/𝐋/𝐓}

Type of magnetic gradient. Based on Unitful.Quantity.

Examples:

julia> (1u"Gauss/cm/s")::SL0D
1 Gauss cm^-1 s^-1

An exemplary struct instantiating AbstractSpinArray.

Fields:

Immutable:

• dim::Dims (nd,): nM ← prod(dim), dimension of the object.
• mask::BitArray (nx,(ny,(nz))): Mask for M, dim == (nx,(ny,(nz)))

Mutable:

• T1::TypeND(T0D, [0,1]) (1,) or (nM,): Longitudinal relaxation coeff.
• T2::TypeND(T0D, [0,1]) (1,) or (nM,): Transversal relaxation coeff.
• γ::TypeND(Γ0D, [0,1]) (1,) or (nM,): Gyromagnetic ratio.
• M::TypeND(Real, [2]) (count(mask), 3): Magnetic spins, (𝑀x,𝑀y,𝑀z).

Notes:

off-resonance, Δf, and locations, loc, are intentionally unincluded, as they are not intrinsic to spins, and can change over time. Unincluding them allows extensional subtypes specialized for, e.g., arterial spin labelling.

See also: AbstractSpinArray.

SpinArray(dim::Dims; T1=1.47u"s", T2=0.07u"s", γ=γ¹H, M=[0. 0. 1.])

Create SpinArray object with prescribed parameters, with mask = trues(dim).

SpinArray(mask::BitArray; T1=1.47u"s", T2=0.07u"s", γ=γ¹H, M=[0. 0. 1.])

Create SpinArray object with prescribed parameters, with dim = size(mask).

UNDER CONSTRUCTION

An exemplary struct instantiating AbstractSpinBolus, designed to model a set of moving spins, e.g., a blood bolus in ASL context.

See also: AbstractSpinBolus.

An exemplary struct instantiating AbstractSpinCube, designed to model a set of regularly spaced spins, e.g., a volume.

Fields:

Immutable:

• spinarray::AbstractSpinArray (1,): inherited AbstractSpinArray struct
• fov ::TypeND(L0D, [2]) (1, 3): field of view.
• ofst::TypeND(L0D, [2]) (1, 3): fov offset from magnetic field iso-center.
• loc ::TypeND(L0D, [2]) (nM, 3): location of spins.

Mutable:

• Δf::TypeND(F0D, [0,1]) (1,) or (nM,): off-resonance map.

See also: AbstractSpinCube.

spincube = SpinCube(mask::BitArray{3}, fov; ofst, Δf, T1, T2, γ)

dim, mask, T1, T2, and γ are passed to SpinArray constructors.

Create SpinCube object with prescribed parameters, with dim = size(mask).

spincube = SpinCube(dim::Dims{3}, fov; ofst, Δf, T1, T2, γ)

Create SpinCube object with prescribed parameters, with mask = trues(dim).

T0D{T<:Real} = Quantity{T, 𝐓}

Type of time. Based on Unitful.Quantity.

Examples:

julia> (1u"s")::T0D
1 s

Γ0D{T<:Real} = Quantity{T, 𝚪}

Type of gyro magnetic ratio. Based on Unitful.Quantity.

Examples:

julia> (1u"Hz/Gauss")::Γ0D
1 Hz Gauss^-1

TypeND(T,Ns) = Union{AbstractArray{<:T,Ns[1]}, AbstractArray{<:T,Ns[2]},...}

Sugar for creating Union{<:T typed array of different dimensions}.

Usage

INPUTS:

• T::Type (1,), the underlying type of the union.
• Ns::Array{Int64,1} (# diff dims,), an array of wanted dimensions.

TypeND(T, ::Colon) = AbstractArray{<:T}

Sugar for creating <:T typed array of arbitrary dimensions.

B = rfgr2B(rf, gr, loc=[0 0 0]u"cm"; Δf=0u"Hz", b1Map=1, γ=γ¹H)

Turn rf, rf, and gradient, gr, into 𝐵-effective magnetic field.

INPUTS:

• rf::TypeND(RF0D, [1,2]) (nT, (nCoil))
• gr::TypeND(GR0D, [2]) (nT, 3)
• loc::TypeND(L0D, [2]) (1,3) or (nM, 3), locations.

KEYWORDS:

• Δf::TypeND(F0D, [0,1,2]) (1,) or (nM,), off-resonance.
• b1Map::TypeND(Union{Real,Complex},[0,1,2]) (1,) or (nM,(nCoils)), transmit sensitivity.
• γ::TypeND(Γ0D, [0,1]) (1,) or (nM,), gyro-ratio

OUTPUS:

• B, generator of TypeND(B0D, [2]) (1,1,nT), 𝐵 field.

See also: Pulse2B, blochsim.

TODO:

Support loc, Δf, and b1Map being Base.Generators.

B = Pulse2B(pulse::Pulse, loc; Δf, b1Map, γ)

Create effective magnetic field, 𝐵, from input pulse.

See also: rfgr2B, B2UΦ, blochsim.

B = Pulse2B(pulse::Pulse, spa::AbstractSpinArray, loc; Δf, b1Map)

...with γ=spa.γ.

B = Pulse2B(pulse::Pulse, cb::AbstractSpinCube; b1Map)

...with loc, Δf, γ = cb.loc, cb.Δf, cb.γ.

blochsim(M, B; T1, T2, γ, dt, doHist)

Same as blochsim!(M, B; T1,T2,γ,dt,doHist), M will not be updated.

See also: blochsim!

blochsim(M, A, B)

Same as blochsim(M, A, B), M will not be updated.

blochsim!(M, B; T1=(Inf)u"s",T2=(Inf)u"s",γ=γ¹H,dt=(4e-6)u"s",doHist=false)

Old school 𝐵-effective magnetic field, B, based bloch simulation. Globally or spin-wisely apply B over spins, M. M will be updated by the results.

INPUTS:

• M::TypeND(Real, [2]) (nM, xyz): input spins' magnetizations.
• B::Union{TypeND(B0D, [2,3]), Base.Generator}: Global, (nT,xyz); Spin-wise, (nM,xyz,nT).

KEYWORDS:

• T1 & T2 ::TypeND(T0D, [0,1]): Global, (1,); Spin-wise, (nM,1).
• γ::TypeND(Γ0D, [0,1]): Global, (1,); Spin-wise, (nM, 1). gyro ratio
• dt::T0D (1,), simulation temporal step size, i.e., dwell time.
• doHist::Bool, whether to output spin history through out B.

OUTPUTS:

• M::TypeND(Real, [2]) (nM, xyz): spins after applying B.
• Mhst::TypeND(Real, [3]) (nM, xyz, nT): spins history during B.

See also: applyPulse, blochsim.

Notes:

1. Not much sanity check inside this function, user is responsible for matching up the dimensions.
2. Put relax at the end of each time step may still be inaccurate, since physically spins relax continuously, this noise/nuance may worth study for applications like fingerprinting simulations, etc.

blochsim!(M, A, B)

Hargreave's 𝐴/𝐵, mat/vec, based bloch simulation. Globally or spin-wisely apply matrix A and vector B over spins, M, described in doi:10.1002/mrm.1170

INPUTS:

• M::TypeND(Real, [2]) (nM, xyz): input spins' magnetizations.
• A::TypeND(AbstractFloat,[3]) (nM, 3,3), A[iM,:,:] is the iM-th 𝐴.
• B::TypeND(AbstractFloat,[2]) (nM, 3), B[iM,:] is the iM-th 𝐵.

OUTPUTS:

• M::TypeND(Real, [2]) (nM, xyz): spins after applying B.

B2AB(B; T1=(Inf)u"s", T2=(Inf)u"s", γ=γ¹H, dt=(4e-6)u"s")

Turn B-effective into Hargreave's 𝐴/𝐵, mat/vec, see: doi:10.1002/mrm.1170.

INPUTS:

• B::Union{TypeND(B0D, [2,3]), Base.Generator}: Global, (nT,xyz); Spin-wise, (nM,xyz,nT).

KEYWORDS:

• T1 & T2 ::TypeND(T0D, [0,1]): Global, (1,); Spin-wise, (nM,1).
• γ::TypeND(Γ0D, [0,1]): Global, (1,); Spin-wise, (nM, 1). gyro ratio
• dt::T0D (1,), simulation temporal step size, i.e., dwell time.

OUTPUTS:

• A::TypeND(AbstractFloat,[3]) (nM, 3,3), A[iM,:,:] is the iM-th 𝐴.
• B::TypeND(AbstractFloat,[2]) (nM, 3), B[iM,:] is the iM-th 𝐵.

See also: rfgr2B, Pulse2B.

B2UΦ(B::TypeND(B0D,[2,3]); γ::TypeND(Γ0D,[0,1]), dt::T0D=4e-6u"s")

Given 𝐵-effective, B, compute rotation axis/angle, U/Φ.

INPUTS:

• B::TypeND(B0D, [2,3]) (1,3,nT) or (nM, 3, nT), 𝐵 field.

KEYWORDS:

• γ::TypeND(Γ0D, [0,1]): Global, (1,); Spin-wise, (nM, 1). gyro ratio
• dt::T0D (1,), simulation temporal step size, i.e., dwell time.

OUTPUTS:

• U::TypeND(Real, [2,3]) (1,3,(nT)) or (nM,3,(nT)), axis.
• Φ::TypeND(Real, [2,3]) (1,1,(nT)) or (nM,1,(nT)), angle.

See also: B2UΦ!, UΦRot!.

Notes:

Somehow, in-place version, B2UΦ!(B,U,Φ; γ,dt), provokes more allocs in julia.

B2UΦ!(B, U; Φ, γ, dt=(4e-6)u"s")

In-place version of B2UΦ. Somehow, B2UΦ!, provokes more allocs in julia.

See also: B2UΦ, blochsim.

UΦRot(U, Φ, V)

Same as UΦRot!(U, Φ, V, R), except not in-place.

See also: UΦRot!.

UΦRot!(U, Φ, V, R)

Apply axis-angle, U-Φ based rotation on V. Rotation is broadcasted on V along its 3rd dimension. Results will overwrite into R.

INPUTS:

• U::TypeND(AbstractFloat,[2]) (nM, 3), rotation axes in 3D, assumed unitary;
• Φ::TypeND(AbstractFloat,[1]) (nM,), rotation angles;
• V::TypeND(AbstractFloat,[2,3]) (nM, 3, (3)), vectors to be rotated;
• R::TypeND(AbstractFloat,[2,3]) (nM, 3, (3)), vectors rotated, i.e., results;

OUTPUTS:

• R the input container R is also returned for convenience.

See also: UΦRot, B2UΦ!.

applyPulse(spa::AbstractSpinArray, p::Pulse, loc; Δf, b1Map, doHist)

Turn p into 𝐵-effective and apply it on spa.M, using its own M, T1, T2, γ.

See also: blochsim, freePrec.

applyPulse(cb::AbstractSpinCube, p::Pulse; b1Map, doHist)

Turn p into 𝐵-effective and apply it on cb.M, using its own M, T1, T2, γ.

applyPulse!(spa::AbstractSpinArray, p::Pulse, loc; Δf, b1Map, doHist)

Update spa.M before return.

applyPulse!(cb::AbstractSpinCube, p::Pulse; b1Map, doHist)

Update cb.M before return.

freePrec(M, t; Δf, T1, T2)

Same as freePrec!(M, t; Δf, T1, T2), M will not be updated.

See also: freePrec!.

freePrec(spa::AbstractSpinArray, t; Δf)

spa::AbstractSpinArray free precess by t. spa.M will not be updated.

See also: applyPulse, freePrec!.

freePrec(cb::AbstractSpinCube, t)

cb::AbstractSpinCube free precess by t. cb.M will not be updated.

See also: applyPulse, freePrec.

freePrec!(M, t; Δf=0u"Hz", T1=(Inf)u"s", T2=(Inf)u"s")

Spins, M, free precess by time t. M will be updated by the results.

INPUTS:

• M::TypeND(Real, [2]) (nM, xyz): input spins' magnetizations.
• t::T0D (1,): duration of free precession.

KEYWORDS:

• T1 & T2 ::TypeND(T0D, [0,1]): Global, (1,); Spin-wise, (nM,1).

OUTPUTS:

• M::TypeND(Real, [2]) (nM, xyz): output spins' magnetizations.

See also: freePrec.

freePrec!(spa::AbstractSpinArray, t; Δf)

...spa.M will updated by the results.

freePrec!(cb::AbstractSpinCube, t)

...cb.M will be updated by the results.

Throw ArgumentError when $x is an immutable field.

CartesianLocations(dim::Dims, doShift::Bool=true)

Retuns a (prod(dim),length(dim)) sized array of grid locations. doShift shifts the locations to be consistent with fftshift.

Examples

julia> loc = [CartesianLocations((2,2)), CartesianLocations((2,2),false)]
2-element Array{Array{Int64,2},1}:
 [-1 -1; 0 -1; -1 0; 0 0]
 [1 1; 2 1; 1 2; 2 2]

ctrInd(dim::Dims) = sum((dim.÷2) .* [1; cumprod([dim[1:end-1]...])])+1

As a separate fn, ensure consistent behariour of getting the linear index to the center of a Nd-array of size dim. This center should match fftshift's center.

See also: ctrSub.

ctrSub(dim::Dims) = CartesianIndex(dim .÷ 2 .+ 1)

As a separate function, ensure consistent behaviour of getting ::CartesianIndex to the center of a Nd-Array of size dim. This center should match fftshift's center.

See also: ctrInd.

Notes:

The function may be removed once julia FFT packages provides this functionality.

g2k(g::TypeND(GR0D,:); isTx::Bool=false, dt::T0D=4e-6u"s", γ::Γ0D=γ¹H)

Compute k-space from gradient.

Usage

INPUTS:

• g::TypeND(GR0D, :) (nSteps, Nd...), gradient
• isTx::Bool, if true, compute transmit k-space, k, ends at the origin.

KEYWORDS:

• dt::T0D (1,), gradient temporal step size, i.e., dwell time.
• γ::Γ0D (1,), gyro-ratio.

OUTPUTS:

• k::TypeND(K0D, :) (nSteps, Nd...), k-space, w/ unit u"cm^-1".

See also: k2g, g2s.

g2s(g::TypeND(GR0D,:); dt::T0D=4e-6u"s")

Slew rate sl, of the gradient, g.

Usage

INPUTS:

• g::TypeND(GR0D, :) (nSteps, Nd...)

KEYWORDS:

• dt::T0D (1,), gradient temporal step size, i.e., dwell time.

OUTPUTS:

• sl::TypeND(Quantity{<:Real, 𝐁/𝐋/𝐓, :) (nSteps, Nd...), slew rate

See also: g2k, k2g.

Note

No s2g is provided for the moment.

k2g(k::TypeND(K0D,:), isTx::Bool=false; dt::T0D=4e-6u"s", γ::Γ0D=γ¹H)

Gradient, g, of the TxRx k-space, (trasmit/receive, excitation/imaging).

Usage

INPUTS:

• k::TypeND(K0D, :) (nSteps, Nd...), Tx or Rx k-space, w/ unit u"cm^-1".
• isTx::Bool, if true, compute transmit k-space, k, ends at the origin.

KEYWORDS:

• dt::T0D (1,), gradient temporal step size, i.e., dwell time.
• γ::Γ0D (1,), gyro-ratio.

OUTPUTS:

• g::TypeND(GR0D, :) (nSteps, Nd...), gradient

Note

The function asserts if k ends at the origin for isTx==true.

See also: g2k, g2s.

SPGR(α; TR, T1)

TE=0 Steady state SPGR signal. 10.1002/mrm.1910130109, eq.(1), ideal spoiling.

INPUTS:

• α::Real (1,), tip angle in degree;

KEYWORDS:

• TR::T0D (1,), repetition time;
• T1::T0D (1,), longitudinal relaxation coefficient;

OUTPUTS:

• sig::Real (1,), steady-state signal.

See also: bSSFP, STFR.

STFR(α, β; ϕ, Δf, T1, T2, Tg, Tf)

TE=0 Steady state STFR signal. 10.1002/mrm.25146, eq.(2), ideal spoiling.

INPUTS:

• α::Real (1,), tip-down angle in degree;
• β::Real (1,), tip-up angle in degree;

KEYWORDS:

• ϕ::Real (1,), phase of the tip-up pulse in radians;
• Δf::F0D (1,), off-resonance in Hz;
• T1::T0D (1,), longitudinal relaxation coefficient;
• T2::T0D (1,), transverse relaxation coefficient;
• Tg::T0D (1,), duration of gradient crusher;
• Tf::T0D (1,), duration of free precession in each TR;

OUTPUTS:

• sig::Real (1,), steady-state signal.

See also: bSSFP, SPGR.

bSSFP(α; TR, Δf, T1, T2)

TE=0 Steady state bSSFP signal. 10.1002/jmri.24163, eq.(4), with ϕ=2π*Δf*TR.

INPUTS:

• α::Real (1,), tip angle in degree;

KEYWORDS:

• TR::T0D (1,), repetition time;
• Δf::F0D (1,), off-resonance in Hz;
• T1::T0D (1,), longitudinal relaxation coefficient;
• T2::T0D (1,), transverse relaxation coefficient;

OUTPUTS:

• sig::Complex (1,), steady-state signal.

See also: SPGR, STFR.

FZstates(Φ, α; TR, T1, T2, FZ)

𝐹, 𝑍 from: 10.1002/(SICI)1099-0534(1999)11:5<291::AID-CMR2>3.0.CO;2-J, eq.(7,8). eq.(7) refined to ÷√(2), instead of ÷2, as in 10.1002/mrm20736: eq.(2).

Assume constant gradient spoiling of m⋅2π dephasing in each TR, m∈ℤ. In practice, if dephase is a constant but not exactly m⋅2π, the resulting states can be computed by convolving a sinc with the m⋅2π dephased results.

INPUTS:

• Φ::TypeND(Real,[2]) (nC,nTR), nC: #C as C in QuadPhase. Typically, one simulates a range of Cs, picking a C yielding a signal intensity equals to that of ideal spgr spoiling.
• α::Real (1,), flip-angle.

KEYWORDS:

• TR::T0D (1,), repetition time;
• T1::T0D (1,), longitudinal relaxation coefficient;
• T2::T0D (1,), transverse relaxation coefficient;
• FZ::NamedTuple, (Fs,Fcs,Zs), simulate from prescribed states if given:
  Fs ::TypeND(Complex,[2]), transversal dephasing states, 𝐹ₙ;
  Fcs::TypeND(Complex,[2]), conjugate transversal dephasing states, 𝐹₋ₙ*;
  Zs ::TypeND(Complex,[2]), longitudinal states, 𝑍ₙ;

OUTPUTS:

• FZ::NamedTuple, (Fs,Fcs,Zs), simulation results.

QuadPhase(nTR::Integer, C::Real, B::Real, A::Real)

Quadratically cycling phases in (°): Φ(n) = mod.(C⋅n² + B⋅n + A, 360), n=0:nTR-1