API

In the following, you find the documentation of all exported functions of the MRIgeneralizedBloch.jl package:

MRIgeneralizedBloch.CRB_gradient_OCT — Method
CRB, grad_ω1, grad_TRF = CRB_gradient_OCT(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, weights; isInversionPulse = [true; falses(length(ω1)-1)])

Calculate the Cramer-Rao bound of a pulse sequence along with the derivatives wrt. ω1 and TRF.

Arguments

  • ω1::Vector{Real}: Control vector of length = Npulses
  • TRF::Vector{Real}: Control vector of length = Npulses
  • TR::Real: Repetition time in seconds
  • ω0::Real: Off-resonance frequency in rad/s
  • B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
  • m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
  • R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
  • R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
  • Rx::Real: Exchange rate between the two spin pools in 1/seconds
  • R1f::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
  • T2s::Real: Transversal relaxation time of the semi-solid pool in seconds
  • R2slT::NTuple{3, Function}: Tuple of three functions: R2sl(TRF, ω1, B1, T2s), dR2sldB1(TRF, ω1, B1, T2s), and R2sldT2s(TRF, ω1, B1, T2s). Can be generated with precompute_R2sl
  • grad_list::Tuple{<:grad_param}: Tuple that specifies the gradients that are calculated; the vector elements can either be any subset/order of grad_list=(grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
  • weights::transpose(Vector{Real}): Row vector of weights applied to the Cramer-Rao bounds (CRB) of the individual parameters. The first entry always refers to the CRB of M0, followed by the values defined in grad_list in the order defined therein. Hence, the vector weights has to have one more entry than grad_list

Optional Keyword Arguments:

  • isInversionPulse::Vector{Bool}: Indicates all inversion pulses; must have the same length as α; the default = [true; falses(length(ω1)-1)] indicates that the first pulse is an inversion pulse and all others are not

Examples

julia> CRB, grad_ω1, grad_TRF = MRIgeneralizedBloch.CRB_gradient_OCT(rand(100) .* π, rand(100) .* 400e-6 .+ 100e-6, 3.5e-3, 0, 1, 0.15, 0.5, 15, 30, 2, 10e-6, precompute_R2sl(), [grad_m0s(), grad_R1f()], transpose([0, 1, 1]); isInversionPulse = [true, falses(99)...])
(2.6266536440386683e20, [0.0, -8.357210433553662e19, 1.8062863407658156e20, -9.181952733568582e19, 2.0889419004304123e20, -1.0127412004909923e20, 1.1472963520187394e20, -6.048455202064828e19, 1.6635577264610125e20, -1.2997982001201938e20  …  -4.0462197701237735e19, 4.4051154836362985e19, -5.703747921741744e19, 1.1580676614266505e20, -1.2930234020298534e20, 1.4073548384507303e20, -9.192708958806614e19, 1.3584033382847213e20, -3.697066939905562e19, 6.313101282386484e19], [0.0, -7.51230331957692e23, 9.811932053428692e23, -1.0734285487552513e24, 6.675582483464475e23, -3.1051435697300785e23, 2.8969707405246626e23, -1.1612336440328984e24, 6.698477560905162e23, -1.8718360662340176e22  …  -2.7429211167215447e23, 2.5368127989367466e23, -6.640000159002342e23, 1.7977260470624765e23, -3.6616011555760077e23, 4.9307219096771845e23, -7.650701790011881e23, 4.5704084508410106e23, -1.0229952455676927e24, 9.526419421729279e23])

c.f. Optimal Control

source
MRIgeneralizedBloch.RF_power! — Method
F = RF_power!(grad_ω1, grad_TRF, ω1, TRF; idx=(ω1 .* TRF) .≉ π, λ=1, Pmax=3e6, TR=3.5e-3)

Calculate RF power penalty and over-write the gradients in place.

Arguments

  • grad_ω1::Vector{Real}: Gradient of control, which will be over-written in place
  • grad_TRF::Vector{Real}: Gradient of control, which will be over-written in place
  • ω1::Vector{Real}: Control vector
  • TRF::Vector{Real}: Control vector

Optional Keyword Arguments:

  • idx::Vector{Bool}: index of flip angles that are considered. Set individual individual pulses to false to exclude, e.g., inversion pulses
  • λ::Real: regularization parameter
  • Pmax::Real: Maximum average power deposition in (rad/s)²; everything above this value will be penalized and with an appropriate λ, the resulting power will be equal to or less than this value.
  • TR::Real: Repetition time of the pulse sequence

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

julia> grad_ω1 = similar(ω1);

julia> grad_TRF = similar(ω1);

julia> F = MRIgeneralizedBloch.RF_power!(grad_ω1, grad_TRF, ω1, TRF; λ=1e3, Pmax=2.85e5)
1.2099652735600044e16
source
MRIgeneralizedBloch.TRF_TV! — Method
F = TRF_TV!(grad_TRF, ω1, TRF; idx=(ω1 .* TRF) .≉ π, λ = 1)

Calculate the total variation penalty of TRF and over-write grad_TRF in place.

Arguments

  • grad_TRF::Vector{Real}: Gradient of control, which will be over-written in place
  • ω1::Vector{Real}: Control vector
  • TRF::Vector{Real}: Control vector

Optional Keyword Arguments:

  • idx::Vector{Bool}: index of flip angles that are considered. Set individual individual pulses to false to exclude, e.g., inversion pulses
  • λ::Real: regularization parameter

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

julia> grad_TRF = similar(ω1);

julia> F = MRIgeneralizedBloch.TRF_TV!(grad_TRF, ω1, TRF; λ = 1e-3)
1.5456176321183175e-5
source
MRIgeneralizedBloch.apply_hamiltonian_gbloch! — Method
apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p, t)

Apply the generalized Bloch Hamiltonian to m and write the resulting derivative wrt. time into ∂m∂t.

Arguments

  • ∂m∂t::Vector{Real}: Vector describing to derivative of m wrt. time; this vector has to be of the same size as m, but can contain any value, which is replaced by H * m
  • m::Vector{Real}: Vector the spin ensemble state of the form [xf, yf, zf, zs, 1] if now gradient is calculated or of the form [xf, yf, zf, zs, 1, ∂xf/∂θ1, ∂yf/∂θ1, ∂zf/∂θ1, ∂zs/∂θ1, 0, ..., ∂xf/∂θn, ∂yf/∂θn, ∂zf/∂θn, ∂zs/∂θn, 0] if n derivatives wrt. θn are calculated
  • mfun: History function; can be initialized with mfun(p, t; idxs=nothing) = typeof(idxs) <: Real ? 0.0 : zeros(5n + 5) for n gradients, and is then updated by the delay differential equation solvers
  • p::NTuple{6,Any}: (ω1, B1, ω0, R1s, T2s, g) or
  • p::NTuple{6,Any}: (ω1, B1, φ, R1s, T2s, g) or
  • p::NTuple{10,Any}: (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g) or
  • p::NTuple{10,Any}: (ω1, B1, φ, m0s, R1f, R2f, Rx, R1s, T2s, g) or
  • p::NTuple{12,Any}: (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s, grad_list) or
  • p::NTuple{12,Any}: (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s, grad_list) with the following entries
    • ω1::Real: Rabi frequency in rad/s (rotation about the y-axis) or
    • ω1(t)::Function: Rabi frequency in rad/s as a function of time for shaped RF-pulses
    • B1::Real: B1 scaling normalized so that B1=1 corresponds to a perfectly calibrated RF field
    • ω0::Real: Larmor or off-resonance frequency in rad/s or
    • φ::Function: RF-phase in rad as a function of time for frequency/phase-sweep pulses (works only in combination with ω1(t)::Function)
    • m0s::Real: Fractional semi-solid spin pool size in the range of 0 to 1
    • R1f::Real: Longitudinal spin relaxation rate of the free pool in 1/seconds
    • R2f::Real: Transversal spin relaxation rate of the free pool in 1/seconds
    • Rx::Real: Exchange rate between the two pools in 1/seconds
    • R1s::Real: Longitudinal spin relaxation rate of the semi-solid pool in 1/seconds
    • T2s::Real: Transversal spin relaxation time of the semi-solid pool in seconds
    • g::Function: Green's function of the form G(κ) = G((t-τ)/T2s)
    • dG_o_dT2s_x_T2s::Function: Derivative of the Green's function wrt. T2s, multiplied by T2s; of the form dG_o_dT2s_x_T2s(κ) = dG_o_dT2s_x_T2s((t-τ)/T2s)
    • grad_list::Vector{grad_param}: List of gradients to be calculated, i.e., any subset of [grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]; length of the vector must be n (cf. arguments m and ∂m∂t); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
  • t::Real: Time in seconds

Optional:

  • pulsetype=:normal: Use default for a regular RF-pulse; the option pulsetype=:inversion should be handled with care as it is only intended to calculate the saturation of the semi-solid pool and its derivative.

Examples

julia> using DifferentialEquations


julia> α = π/2;

julia> TRF = 100e-6;

julia> ω1 = α/TRF;

julia> B1 = 1;

julia> ω0 = 0;

julia> m0s = 0.2;

julia> R1f = 1/3;

julia> R2f = 15;

julia> R1s = 2;

julia> T2s = 10e-6;

julia> Rx = 30;

julia> G = interpolate_greens_function(greens_superlorentzian, 0, TRF / T2s);


julia> m0 = [0; 0; 1-m0s; m0s; 1];

julia> mfun(p, t; idxs=nothing) = typeof(idxs) <: Real ? 0.0 : zeros(5);

julia> sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G)))
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 9-element Vector{Float64}:
 0.0
 1.375006182301112e-7
 1.512506800531223e-6
 8.042561696923577e-6
 2.107848894861101e-5
 3.9114182159866e-5
 6.26879358261189e-5
 9.147711414688425e-5
 0.0001
u: 9-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.8, 0.2, 1.0]
 [0.0017278806030763402, 0.0, 0.7999981340131751, 0.19999953350448, 1.0]
 [0.019004717382235078, 0.0, 0.7997742277135814, 0.19994357804868362, 1.0]
 [0.10079111348917136, 0.0, 0.7936248122939504, 0.19842287240368398, 1.0]
 [0.2600257867257624, 0.0, 0.7565529666157949, 0.1898191304278861, 1.0]
 [0.46104237829774064, 0.0, 0.6537239462232086, 0.16937683398576228, 1.0]
 [0.6661740376622253, 0.0, 0.44261209248221817, 0.13589311206074786, 1.0]
 [0.7923117772809817, 0.0, 0.10713073823030607, 0.09390260581965477, 1.0]
 [0.7994211188442756, 0.0, 0.0004403374305009168, 0.08214809659226184, 1.0]

julia> using Plots

julia> plot(sol, labels=["xf" "yf" "zf" "zs" "1"], xlabel="t (s)", ylabel="m(t)");




julia> dG_o_dT2s_x_T2s = interpolate_greens_function(dG_o_dT2s_x_T2s_superlorentzian, 0, TRF / T2s);


julia> grad_list = (grad_R2f(), grad_m0s());


julia> m0 = [0; 0; 1-m0s; m0s; 1; zeros(5*length(grad_list))];


julia> mfun(p, t; idxs=nothing) = typeof(idxs) <: Real ? 0.0 : zeros(5 + 5*length(grad_list));

julia> sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G, dG_o_dT2s_x_T2s, grad_list)));




julia> plot(sol);
source
MRIgeneralizedBloch.apply_hamiltonian_sled! — Method
apply_hamiltonian_sled!(∂m∂t, m, p, t)

Apply Sled's Hamiltonian to m and write the resulting derivative wrt. time into ∂m∂t.

Arguments

  • ∂m∂t::Vector{<:Real}: Vector of length 1 describing to derivative of m wrt. time; this vector can contain any value, which is replaced by H * m
  • m::Vector{<:Real}: Vector of length 1 describing the zs magnetization
  • p::NTuple{6 or 10, Any}: (ω1, B1, ω0, R1s, T2s, g) for a simulating an isolated semi-solid pool or (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g) for simulating a coupled spin system; with
    • ω1::Real: Rabi frequency in rad/s (rotation about the y-axis) or
    • ω1(t)::Function: Rabi frequency in rad/s as a function of time for shaped RF-pulses
    • B1::Real: B1 scaling normalized so that B1=1 corresponds to a perfectly calibrated RF field
    • ω0::Real: Larmor or off-resonance frequency in rad/s (is only used for the free spin pool)
    • R1f::Real: Longitudinal spin relaxation rate of the free pool in 1/seconds
    • R2f::Real: Transversal spin relaxation rate of the free pool in 1/seconds
    • R1s::Real: Longitudinal spin relaxation rate of the semi-solid in 1/seconds
    • Rx::Real: Exchange rate between the two pools in 1/seconds
    • T2s::Real: Transversal spin relaxation time in seconds
    • g::Function: Green's function of the form G(κ) = G((t-τ)/T2s)
  • t::Real: Time in seconds

Examples

julia> using DifferentialEquations


julia> α = π/2;

julia> TRF = 100e-6;

julia> ω1 = α/TRF;

julia> B1 = 1;

julia> ω0 = 0;

julia> R1s = 2;

julia> T2s = 10e-6;

julia> G = interpolate_greens_function(greens_superlorentzian, 0, TRF / T2s);

julia> m0 = [1];

julia> sol = solve(ODEProblem(apply_hamiltonian_sled!, m0, (0, TRF), (ω1, 1, ω0, R1s, T2s, G)), Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 3-element Vector{Float64}:
 0.0
 7.475414666720001e-5
 0.0001
u: 3-element Vector{Vector{Float64}}:
 [1.0]
 [0.6313928231811967]
 [0.4895365449661915]

julia> using Plots

julia> plot(sol, labels=["zs"], xlabel="t (s)", ylabel="m(t)");
source
MRIgeneralizedBloch.bound_ω1_TRF! — Method
x = bound_ω1_TRF!(ω1, TRF; ω1_min = 0, ω1_max = 2e3π, TRF_min = 100e-6, TRF_max = 500e-6)

Bound the controls ω1 and TRF (over-written in place) and return a vector of length 2Npulses with values in the range [-Inf, Inf] that relate to the bounded ω1 and TRF via the tanh function.

Arguments

  • ω1::Vector{Real}: Control vector of length Npulses
  • TRF::Vector{Real}: Control vector of length Npulses

Optional Keyword Arguments (see above for defaults):

  • ω1_min::Real: lower bound for ω1 in rad/s
  • ω1_max::Real: upper bound for ω1 in rad/s
  • TRF_min::Real: lower bound for TRF in s
  • TRF_max::Real: upper bound for TRF in s

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

julia> x = MRIgeneralizedBloch.bound_ω1_TRF!(ω1, TRF)
200-element Vector{Float64}:
  Inf
   0.3084393995336273
  -0.222010414402936
   0.1981882670715111
  Inf
  -0.3590504865635206
  Inf
  Inf
  Inf
  Inf
   ⋮
  -1.3069326972020265
 -Inf
   4.449542946745821
   0.06378289522574039
   0.2699178281115771
   1.271412541641303
  -0.46991502137668045
  -1.4377324062335655
  -0.9746661706753423
source
MRIgeneralizedBloch.calculatesignal_gbloch_ide — Method
calculatesignal_gbloch_ide(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s[; grad_list, Ncyc=2, output=:complexsignal])

Calculate the signal or magnetization evolution with the full generalized Bloch model assuming a super-Lorentzian lineshape (slow).

The simulation assumes a sequence of rectangular RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

Arguments

  • α::Vector{Real}: Array of flip angles in radians
  • TRF::Vector{Real}: Array of the RF-pulse durations in seconds
  • TR::Real: Repetition time in seconds
  • ω0::Real: Off-resonance frequency in rad/s
  • B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
  • m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
  • R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
  • R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
  • Rx::Real: Exchange rate between the two spin pools in 1/seconds
  • R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
  • T2s::Real: Transversal relaxation time of the semi-solid pool in seconds

Optional:

  • grad_list=(): Tuple that specifies the gradients that are calculated; the Tuple can either be empty () for no gradient, or contain any subset/order of grad_list=(grad_m0s(), grad_R1s(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
  • Ncyc=2: The magnetization is initialized with thermal equilibrium and then performed Ncyc times and only the last cycle is stored. The default value is usually a good approximation for antiperiodic boundary conditions. Increase the number for higher precision at the cost of computation time.
  • output=:complexsignal: The default keywords triggers the function to output a complex-valued signal (xf + 1im yf); the keyword output=:realmagnetization triggers an output of the entire (real valued) vector [xf, yf, zf, xs, zs]
  • greens=(greens_superlorentzian, dG_o_dT2s_x_T2s_superlorentzian): Tuple of a Greens function G(κ) = G((t-τ)/T2s) and its partial derivative wrt. T2s, multiplied by T2s ∂G((t-τ)/T2s)/∂T2s * T2s. This package supplies the three Greens functions greens=(greens_superlorentzian, dG_o_dT2s_x_T2s_superlorentzian) (default), greens=(greens_lorentzian, dG_o_dT2s_x_T2s_lorentzian), and greens=(greens_gaussian, dG_o_dT2s_x_T2s_gaussian)

Examples

julia> calculatesignal_gbloch_ide(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6)
100×1 Matrix{ComplexF64}:
  -0.007966316445310092 + 0.0im
  0.0012590590420428192 - 0.0im
  -0.006088855588249714 + 0.0im
   0.002418738940913121 - 0.0im
  -0.004361339394954344 + 0.0im
  0.0034891358222515403 - 0.0im
 -0.0027633710605856443 + 0.0im
   0.004483217942995644 - 0.0im
 -0.0012812573504847747 + 0.0im
    0.00540885403652597 - 0.0im
                        ⋮
   0.017760808273166742 - 0.0im
    0.01757611897551337 + 0.0im
    0.01781395094611674 - 0.0im
    0.01764385633643124 + 0.0im
   0.017863575856212052 - 0.0im
   0.017706926034820797 + 0.0im
   0.017909914934264796 - 0.0im
    0.01776565037417232 + 0.0im
   0.017953184893722073 - 0.0im

julia> calculatesignal_gbloch_ide(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6; grad_list=(grad_R1f(), grad_T2s()), output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.00796627   0.0   0.000637773  …   0.0   -10.8757  -335.26   0.0
  0.00125903  -0.0  -0.00700671      -0.0   125.882   -326.977  0.0
 -0.00608882   0.0   0.00185086       0.0   -30.4187  -317.56   0.0
  0.00241873  -0.0  -0.00520622      -0.0    96.1776  -309.906  0.0
 -0.00436133   0.0   0.00296471       0.0   -47.5803  -302.948  0.0
  0.003489    -0.0  -0.00354518   …  -0.0    69.5148  -298.697  0.0
 -0.00276366   0.0   0.00399588       0.0   -62.8453  -294.886  0.0
  0.00448273  -0.0  -0.00200673      -0.0    45.3179  -292.783  0.0
 -0.00128187   0.0   0.00495478       0.0   -76.6573  -290.321  0.0
  0.00540814  -0.0  -0.000578836     -0.0    23.1756  -289.245  0.0
  ⋮                               ⋱
  0.0177563   -0.0   0.0175372       -0.0  -290.779   -349.855  0.0
  0.0175716    0.0   0.0177845        0.0  -295.347   -350.002  0.0
  0.0178094   -0.0   0.0176073       -0.0  -292.44    -350.163  0.0
  0.0176393    0.0   0.0178359        0.0  -296.668   -350.3    0.0
  0.017859    -0.0   0.0176727    …  -0.0  -294.001   -350.451  0.0
  0.0177024    0.0   0.0178838        0.0  -297.914   -350.579  0.0
  0.0179053   -0.0   0.0177335       -0.0  -295.467   -350.72   0.0
  0.0177611    0.0   0.0179286        0.0  -299.09    -350.84   0.0
  0.0179486   -0.0   0.0177902       -0.0  -296.845   -350.972  0.0
source
MRIgeneralizedBloch.calculatesignal_graham_ode — Method
calculatesignal_graham_ode(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s[; grad_list, Ncyc=2, output=:complexsignal])

Calculate the signal or magnetization evolution with Graham's spectral model assuming a super-Lorentzian lineshape.

The simulation assumes a sequence of rectangular RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

Arguments

  • α::Vector{Real}: Array of flip angles in radians
  • TRF::Vector{Real}: Array of the RF-pulse durations in seconds
  • TR::Real: Repetition time in seconds
  • ω0::Real: Off-resonance frequency in rad/s
  • B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
  • m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
  • R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
  • R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
  • Rx::Real: Exchange rate between the two spin pools in 1/seconds
  • R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
  • T2s::Real: Transversal relaxation time of the semi-solid pool in seconds

Optional:

  • grad_list=(): Tuple that specifies the gradients that are calculated; the Tuple can either be empty () for no gradient, or contain any subset/order of grad_list=(grad_m0s(), grad_R1s(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
  • Ncyc=2: The magnetization is initialized with thermal equilibrium and then performed Ncyc times and only the last cycle is stored. The default value is usually a good approximation for antiperiodic boundary conditions. Increase the number for higher precision at the cost of computation time.
  • output=:complexsignal: The default keywords triggers the function to output a complex-valued signal (xf + 1im yf); the keyword output=:realmagnetization triggers an output of the entire (real valued) vector [xf, yf, zf, xs, zs]

Examples

julia> calculatesignal_graham_ode(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6)
100×1 Matrix{ComplexF64}:
   -0.00807345119181352 + 0.0im
   0.001268643290482942 - 0.0im
 -0.0061786944372399285 + 0.0im
  0.0024358658178670984 - 0.0im
  -0.004437476277201395 + 0.0im
   0.003516464649840319 - 0.0im
 -0.0028315542520438736 + 0.0im
   0.004523902250379144 - 0.0im
 -0.0013422299483389865 + 0.0im
   0.005454562034842185 - 0.0im
                        ⋮
   0.018148222040953822 - 0.0im
    0.01795769661414946 + 0.0im
   0.018204860876661903 - 0.0im
    0.01802936310261858 + 0.0im
    0.01825782014256292 - 0.0im
   0.018096168912467753 + 0.0im
   0.018307337729977783 - 0.0im
    0.01815844450117701 + 0.0im
   0.018353636229654767 - 0.0im

julia> calculatesignal_graham_ode(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6; grad_list=(grad_R1f(), grad_T2s()), output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.0080756    0.0   0.000643162  …   0.0   -10.4986  -323.634  0.0
  0.00126867  -0.0  -0.00710692      -0.0   123.078   -316.358  0.0
 -0.00618067   0.0   0.00186482       0.0   -29.4458  -307.862  0.0
  0.00243634  -0.0  -0.0052899       -0.0    94.1692  -300.821  0.0
 -0.00443746   0.0   0.00298646       0.0   -46.1422  -294.116  0.0
  0.00351386  -0.0  -0.00361718   …  -0.0    68.262   -289.421  0.0
 -0.00282882   0.0   0.00402793       0.0   -61.0236  -285.242  0.0
  0.00451808  -0.0  -0.00206741      -0.0    44.7888  -282.867  0.0
 -0.00133608   0.0   0.00499316       0.0   -74.4075  -280.79   0.0
  0.00544896  -0.0  -0.00062663      -0.0    23.2536  -280.148  0.0
  ⋮                               ⋱
  0.0181755   -0.0   0.0179431       -0.0  -285.071   -337.48   0.0
  0.0179911    0.0   0.0181963        0.0  -289.591   -337.634  0.0
  0.0182343   -0.0   0.018019        -0.0  -286.777   -337.804  0.0
  0.0180645    0.0   0.0182531        0.0  -290.964   -337.948  0.0
  0.0182893   -0.0   0.0180898    …  -0.0  -288.383   -338.106  0.0
  0.0181329    0.0   0.0183062        0.0  -292.261   -338.241  0.0
  0.0183407   -0.0   0.0181559       -0.0  -289.895   -338.39   0.0
  0.0181968    0.0   0.0183559        0.0  -293.488   -338.516  0.0
  0.0183888   -0.0   0.0182175       -0.0  -291.318   -338.656  0.0
source
MRIgeneralizedBloch.calculatesignal_linearapprox — Method
calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT[; grad_list=(undef,), rfphase_increment=[π], m0=:periodic, output=:complexsignal, isInversionPulse = [true; falses(length(α)-1)]])

Calculate the signal or magnetization evolution with the linear approximation of the generalized Bloch model assuming a super-Lorentzian lineshape.

The simulation assumes a sequence of rectangular RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

Arguments

  • α::Vector{Real}: Array of flip angles in radians
  • TRF::Vector{Real}: Array of the RF-pulse durations in seconds
  • TR::Real: Repetition time in seconds
  • ω0::Real: Off-resonance frequency in rad/s
  • B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
  • m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
  • R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
  • R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
  • Rx::Real: Exchange rate between the two spin pools in 1/seconds
  • R1f::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
  • T2s::Real: Transversal relaxation time of the semi-solid pool in seconds
  • R2slT::NTuple{3, Function}: Tuple of three functions: R2sl(TRF, ω1, B1, T2s), dR2sldB1(TRF, ω1, B1, T2s), and R2sldT2s(TRF, ω1, B1, T2s). Can be generated with precompute_R2sl

Optional:

  • grad_list=(undef,): Tuple that specifies the gradients that are calculated; the vector elements can either be undef for no gradient, or any subset/order of grad_list=(grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
  • rfphase_increment=[π]::Vector{Real}: Increment of the RF phase between consecutive pulses. The default value π, together with $ω0=0$ corresponds to the on-resonance condition. When more than one value is supplied, their resulting signal is stored along the second dimension of the output array.
  • m0=:periodic: With the default keyword :periodic, the signal and their derivatives are calculated assuming $m(0) = -m(T)$, where T is the duration of the RF-train. With the keyword :thermal, the magnetization $m(0)$ is initialized with thermal equilibrium [xf, yf, zf, xs, zs] = [0, 0, 1-m0s, 0, m0s], followed by a α[1]/2 - TR/2 prep pulse; and with the keyword :IR, this initialization is followed an inversion pulse of duration TRF[1], (set α[1]=π) and a α[2]/2 - TR/2 prep pulse.
  • preppulse=false: if true, a α/2 - TR/2 preparation is applied. In the case of m0=:IR, it is applied after the inversion pulse based on α[2], otherwise it is based on α[1]
  • output=:complexsignal: The default keywords triggers the function to output a complex-valued signal (xf + 1im yf); the keyword output=:realmagnetization triggers an output of the entire (real valued) vector [xf, yf, zf, xs, zs, 1]
  • isInversionPulse::Vector{Bool}: Indicates all inversion pulses; must have the same length as α; the default = [true; falses(length(α)-1)] indicates that the first pulse is an inversion pulse and all others are not

Examples

julia> R2slT = precompute_R2sl();


julia> calculatesignal_linearapprox(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.1, 1, 15, 30, 6.5, 10e-6, R2slT)
100×1×1 Array{ComplexF64, 3}:
[:, :, 1] =
  -0.029305987774458163 - 0.0im
   0.004329424678273618 + 0.0im
  -0.022761409373843695 + 0.0im
   0.008280330224850314 + 0.0im
  -0.016751305727868086 + 0.0im
   0.011921649143708494 + 0.0im
   -0.01119057989041819 + 0.0im
   0.015305608440952356 + 0.0im
 -0.0060267918180664974 + 0.0im
   0.018463027499697613 + 0.0im
                        ⋮
   0.061531473509660414 + 0.0im
    0.06081400768357244 + 0.0im
     0.0617257515931871 + 0.0im
   0.061064216222629225 + 0.0im
     0.0619074570531536 + 0.0im
    0.06129750535994419 + 0.0im
   0.062077399758002534 + 0.0im
   0.061515021919442275 + 0.0im
    0.06223633770306246 + 0.0im
source
MRIgeneralizedBloch.dG_o_dT2s_x_T2s_gaussian — Method
dG_o_dT2s_x_T2s_gaussian(κ)

Evaluate the derivative of Green's function, corresponding to a Gaussian lineshape, wrt. T2s at κ = (t-τ)/T2s and multiply it by T2s.

The multiplication is added so that the function merely depends on κ = (t-τ)/T2s. The actual derivative is given by dG_o_dT2s_x_T2s_gaussian((t-τ)/T2s)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> dGdT2s = dG_o_dT2s_x_T2s_gaussian((t-τ)/T2s)/T2s
1.9287498479639177e-15
source
MRIgeneralizedBloch.dG_o_dT2s_x_T2s_lorentzian — Method
dG_o_dT2s_x_T2s_lorentzian(κ)

Evaluate the derivative of Green's function, corresponding to a Lorentzian lineshape, wrt. T2s at κ = (t-τ)/T2s and multiply it by T2s.

The multiplication is added so that the function merely depends on κ = (t-τ)/T2s. The actual derivative is given by dG_o_dT2s_x_T2s_lorentzian((t-τ)/T2s)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> dGdT2s = dG_o_dT2s_x_T2s_lorentzian((t-τ)/T2s)/T2s
45.39992976248485
source
MRIgeneralizedBloch.dG_o_dT2s_x_T2s_superlorentzian — Method
dG_o_dT2s_x_T2s_superlorentzian(κ)

Evaluate the derivative of Green's function, corresponding to a super-Lorentzian lineshape, wrt. T2s at κ = (t-τ)/T2s and multiply it by T2s.

The multiplication is added so that the function merely depends on κ = (t-τ)/T2s. The actual derivative is given by dG_o_dT2s_x_T2s_superlorentzian((t-τ)/T2s)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> dGdT2s = dG_o_dT2s_x_T2s_superlorentzian((t-τ)/T2s)/T2s
15253.095033670965
source
MRIgeneralizedBloch.fit_gBloch — Method
qM = fit_gBloch(data, α, TRF, TR;
    reM0 = (-Inf,   1,  Inf),
    imM0 = (-Inf,   0,  Inf),
    m0s  = (   0, 0.2,    1),
    R1f  = (   0, 0.3,  Inf),
    R2f  = (   0,  15,  Inf),
    Rx   = (   0,  20,  Inf),
    R1s  = (   0,   3,  Inf),
    T2s  = (8e-6,1e-5,12e-6),
    ω0   = (-Inf,   0,  Inf),
    B1   = (   0,   1,  1.5),
    R1a  = (   0, 0.7,  Inf),
    u=1,
    fit_apparentR1=false,
    show_trace=false,
    maxIter=100,
    R2slT = precompute_R2sl(TRF_min=minimum(TRF), TRF_max=maximum(TRF), T2s_min=minimum(T2s), T2s_max=maximum(T2s), ω1_max=maximum(α ./ TRF), B1_max=maximum(B1)),
    )

Fit the generalized Bloch model for a train of RF pulses and balanced gradient moments to data.

Arguments

  • data::Vector{Number}: Array of measured data points, either in the time or a compressed domain (cf. u)
  • α::Vector{Real}: Array of flip angles in radians; can also be a Vector{Vector{Real}} which simulates each RF pattern and concatenates the signals of each simulation
  • TRF::Vector{Real}: Array of the RF-pulse durations in seconds (or Vector{Vector{Real}} if α::Vector{Vector{Real}}`)
  • TR::Real: Repetition time in seconds
  • ω0::Real: Off-resonance frequency in rad/s

Optional Keyword Arguments:

  • reM0::Union{Real, Tuple{Real, Real, Real}}: Real part of M0; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • imM0::Union{Real, Tuple{Real, Real, Real}}: Imaginary part of M0; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • m0s::Union{Real, Tuple{Real, Real, Real}}: Fractional size of the semi-solid pool (should be in range of 0 to 1); either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • R1f::Union{Real, Tuple{Real, Real, Real}}: Longitudinal relaxation rate of the free pool in 1/s; only used in combination with fit_apparentR1=false; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • R2f::Union{Real, Tuple{Real, Real, Real}}: Transversal relaxation rate of the free pool in 1/s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • Rx::Union{Real, Tuple{Real, Real, Real}}: Exchange rate between the two spin pools in 1/s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • R1s::Union{Real, Tuple{Real, Real, Real}}: Longitudinal relaxation rate of the semi-solid pool in 1/s; only used in combination with fit_apparentR1=false; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • T2s::Union{Real, Tuple{Real, Real, Real}}: Transversal relaxation time of the semi-solid pool in s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • ω0::Union{Real, Tuple{Real, Real, Real}}: Off-resonance frequency in rad/s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • B1::Union{Real, Tuple{Real, Real, Real}}: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • R1a::Union{Real, Tuple{Real, Real, Real}}: Apparent longitudinal relaxation rate in 1/s; only used in combination with fit_apparentR1=true; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
  • u::Union{Number, Matrix}: Compression matrix that transform the simulated time series to a series of coefficients. Set to 1 by default to enable the fitting in the time domain
  • fit_apparentR1::Bool: Switch between fitting R1f and R1s separately (false; default) and an apparent R1a = R1f = R1s (true)
  • show_trace::Bool: print output during the optimization; default=false
  • maxIter::Int: Maximum number of iteration; default=100
  • R2slT::NTuple{3, Function}: Tuple of three functions: R2sl(TRF, ω1, B1, T2s), dR2sldB1(TRF, ω1, B1, T2s), and R2sldT2s(TRF, ω1, B1, T2s). By default generated with precompute_R2sl

Examples

c.f. Non-Linear Least Square Fitting

source
MRIgeneralizedBloch.get_bounded_ω1_TRF — Method
ω1, TRF = get_bounded_ω1_TRF(x; ω1_min = 0, ω1_max = 2e3π, TRF_min = 100e-6, TRF_max = 500e-6)

Transform a vector of length 2Npulses with values in the range [-Inf, Inf] into two vectors of length Npulses, which describe the bounded controls ω1 and TRF.

Arguments

  • x::Vector{Real}: Control vector of length = 2Npulses with values in the range [-Inf, Inf]

Optional Keyword Arguments:

  • ω1_min::Real: lower bound for ω1 in rad/s
  • ω1_max::Real: upper bound for ω1 in rad/s
  • TRF_min::Real: lower bound for TRF in s
  • TRF_max::Real: upper bound for TRF in s

Examples

julia> x = 1000 * randn(2 * 100);

julia> ω1, TRF = MRIgeneralizedBloch.get_bounded_ω1_TRF(x)
([0.0, 6283.185307179586, 0.0, 0.0, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 0.0, 0.0  …  0.0, 0.0, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 0.0, 0.0, 1.4115403811440933e-9, 0.0], [0.0005, 0.0001, 0.0001, 0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.0001, 0.0005  …  0.0005, 0.0005, 0.0001, 0.0001, 0.0001, 0.0004999999277453363, 0.0005, 0.0001, 0.0001, 0.0001])
source
MRIgeneralizedBloch.graham_saturation_rate_single_frequency — Method
graham_saturation_rate_single_frequency(lineshape, ω1, TRF, Δω)

Calculate saturation rate (in units of 1/s) according to Graham's single frequency approximation.

Arguments

  • lineshape::Function: as a function of ω₀ (in rad/s). Supply, e.g., the anonymous function ω₀ -> lineshape_superlorentzian(ω₀, T2s). Note that the integral over the lineshape has to be 1.
  • ω1::Function: ω1 in rad/s as a function of time (in units of s) where the puls shape is defined for t ∈ [0,TRF]
  • TRF::Real: duration of the RF pulse in s
  • Δω::Real: offset frequency in rad/s

Examples

julia> using SpecialFunctions

julia> T2s = 10e-6;

julia> α = π;

julia> TRF = 100e-6;

julia> NSideLobes = 1;

julia> ω1(t) = sinc(2(NSideLobes+1) * t/TRF - (NSideLobes+1)) * α / (sinint((NSideLobes+1)π) * TRF/π / (NSideLobes+1));

julia> Δω = 200;

julia> graham_saturation_rate_single_frequency(ω₀ -> lineshape_superlorentzian(ω₀, T2s), ω1, TRF, Δω)
419969.3376658947
source
MRIgeneralizedBloch.graham_saturation_rate_spectral — Method
graham_saturation_rate_spectral(lineshape, ω1, TRF, Δω)

Calculate saturation rate (in units of 1/s) according to Graham's spectral model.

Arguments

  • lineshape::Function: as a function of ω₀ (in rad/s). Supply, e.g., the anonymous function ω₀ -> lineshape_superlorentzian(ω₀, T2s). Note that the integral over the lineshape has to be 1.
  • ω1::Function: ω1 in rad/s as a function of time (in units of s) where the puls shape is defined for t ∈ [0,TRF]
  • TRF::Real: duration of the RF pulse in s
  • Δω::Real: offset frequency in rad/s

Examples

julia> using SpecialFunctions

julia> T2s = 10e-6;

julia> α = π;

julia> TRF = 100e-6;

julia> NSideLobes = 1;

julia> ω1(t) = sinc(2(NSideLobes+1) * t/TRF - (NSideLobes+1)) * α / (sinint((NSideLobes+1)π) * TRF/π / (NSideLobes+1));

julia> Δω = 200;

julia> graham_saturation_rate_spectral(ω₀ -> lineshape_superlorentzian(ω₀, T2s), ω1, TRF, Δω)
56135.388046022905
source
MRIgeneralizedBloch.greens_gaussian — Method
greens_gaussian(κ)

Evaluate the Green's function corresponding to a Gaussian lineshape at κ = (t-τ)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_gaussian((t-τ)/T2s)
1.9287498479639178e-22
source
MRIgeneralizedBloch.greens_lorentzian — Method
greens_lorentzian(κ)

Evaluate the Green's function corresponding to a Lorentzian lineshape at κ = (t-τ)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_lorentzian((t-τ)/T2s)
4.5399929762484854e-5
source
MRIgeneralizedBloch.greens_superlorentzian — Method
greens_superlorentzian(κ)

Evaluate the Green's function corresponding to a super-Lorentzian lineshape at κ = (t-τ)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_superlorentzian((t-τ)/T2s)
0.1471246868094442
source
MRIgeneralizedBloch.hamiltonian_linear — Method
hamiltonian_linear(ω1, B1, ω0, T, m0s, R1f, R2f, Rx, R1s, R2s[, dR2sdT2s, dR2sdB1, grad_type])

Calculate the hamiltonian of the linear approximation of the generalized Bloch model.

If no gradient is supplied, it returns a 6x6 (static) matrix with the dimensions (in this order) [xf, yf, zf, xs, zs, 1]; the attached 1 is a mathematical trick to allow for $T_1$ relaxation to a non-zero thermal equilibrium. If a gradient is supplied, it returns a 11x11 (static) matrix with the dimensions (in this order) [xf, yf, zf, xs, zs, dxf/dθ, dyf/dθ, dzf/dθ, dxs/dθ, dzs/dθ, 1], where θ is the parameter specified by grad_type

Arguments

  • ω1::Real: Rabi frequency in rad/s (rotation about the y-axis)
  • B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
  • ω0::Real: Larmor (or off-resonance) frequency in rad/s (rotation about the z-axis)
  • T::Real: Time in seconds; this can, e.g., be the RF-pulse duration, or the time of free precession with ω1=0
  • m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
  • R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
  • R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
  • Rx::Real: Exchange rate between the two spin pools in 1/seconds
  • R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
  • R2s::Real: Transversal relaxation rate of the semi-solid pool in 1/seconds; this number can be calculated with the first function returned by precompute_R2sl to implement the linear approximation described in the generalized Bloch paper

Optional:

  • dR2sdT2s::Real: Derivative of linearized R2sl wrt. the actual T2s; only required if grad_type = grad_T2s(); this number can be calculated with the second function returned by precompute_R2sl
  • dR2sdB1::Real: Derivative of linearized R2sl wrt. B1; only required if grad_type = grad_B1(); this number can be calculated with the third function returned by precompute_R2sl
  • grad_type::grad_param: grad_m0s(), grad_R1f(), grad_R1s(), grad_R2f(), grad_Rx(), grad_T2s(), grad_ω0(), or grad_B1(); create one hamiltonian for each desired gradient

Examples

julia> α = π;

julia> T = 500e-6;

julia> ω1 = α/T;

julia> B1 = 1;

julia> ω0 = 0;

julia> m0s = 0.1;

julia> R1f = 1;

julia> R2f = 15;

julia> Rx = 30;

julia> R1s = 6.5;

julia> R2s = 1e5;

julia> m0 = [0, 0, 1-m0s, 0, m0s, 1];

julia> (xf, yf, zf, xs, zs, _) = exp(hamiltonian_linear(ω1, B1, ω0, T, m0s, R1f, R2f, Rx, R1s, R2s)) * m0
6-element StaticArraysCore.SVector{6, Float64} with indices SOneTo(6):
  0.0010647535813058293
  0.0
 -0.8957848274535014
  0.005126529591877105
  0.08122007142111888
  1.0
source
MRIgeneralizedBloch.interpolate_greens_function — Method
interpolate_greens_function(f, κmin, κmax)

Interpolate the Green's function f in the range between κmin and κmax.

The interpolation uses the ApproxFun.jl package that incorporates Chebyshev polynomials and ensures an approximation to machine precision.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_superlorentzian((t-τ)/T2s)
0.1471246868094442

julia> Gint = interpolate_greens_function(greens_superlorentzian, 0, 20);


julia> Gint((t-τ)/T2s)
0.14712468680944407
source
MRIgeneralizedBloch.precompute_R2sl — Method
precompute_R2sl([;TRF_min=100e-6, TRF_max=500e-6, T2s_min=5e-6, T2s_max=21e-6, ω1_max=π/TRF_max, B1_max=1.5, greens=greens_superlorentzian])

Pre-compute and interpolate the linearized R2sl(TRF, α, B1, T2s) and its derivatives dR2sldB1(TRF, α, B1, T2s), R2sldT2s(TRF, α, B1, T2s) etc. in the range specified by the arguments.

The function solves the generalized Bloch equations of an isolated semi-solid pool for values in the specified range, calculates the linearized R2sl that minimizes the error of zs at the end of the RF-pulse, and interpolates between the different samples.

Optional Arguments:

  • TRF_min::Real: lower bound of the RF-pulse duration range in seconds
  • TRF_max::Real: upper bound of the RF-pulse duration range in seconds
  • T2s_min::Real: lower bound of the T2s range in seconds
  • T2s_max::Real: upper bound of the T2s range in seconds
  • ω1_max::Real: upper bound of the Rabi frequency ω1, the default is the frequency of a 500μs long π-pulse
  • B1_max::Real: upper bound of the B1 range, normalized so that B1 = 1 corresponds to a perfectly calibrated RF field
  • greens=greens_superlorentzian: Greens function in the form G(κ) = G((t-τ)/T2s). This package supplies the three Greens functions greens=greens_superlorentzian (default), greens=greens_lorentzian, and greens=greens_gaussian

Examples

julia> R2slT = precompute_R2sl();


julia> R2sl, dR2sldB1, R2sldT2s, _ = precompute_R2sl(TRF_min=100e-6, TRF_max=500e-6, T2s_min=5e-6, T2s_max=15e-6, ω1_max=π/500e-6, B1_max=1.3, greens=greens_gaussian);
source
MRIgeneralizedBloch.second_order_α! — Method
F = second_order_α!(grad_ω1, grad_TRF, ω1, TRF; idx=(ω1 .* TRF) .≉ π, λ = 1)

Calculate second order penalty of variations of the flip angle α and over-write the gradients in place.

Arguments

  • grad_ω1::Vector{Real}: Gradient of control, which will be over-written in place
  • grad_TRF::Vector{Real}: Gradient of control, which will be over-written in place
  • ω1::Vector{Real}: Control vector
  • TRF::Vector{Real}: Control vector

Optional Keyword Arguments:

  • idx::Vector{Bool}: index of flip angles that are considered. Set individual individual pulses to false to exclude, e.g., inversion pulses
  • λ::Real: regularization parameter

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

julia> grad_ω1 = similar(ω1);

julia> grad_TRF = similar(ω1);

julia> F = MRIgeneralizedBloch.second_order_α!(grad_ω1, grad_TRF, ω1, TRF; λ = 1e-3)
0.3272308747790844
source