Tutorial
This tutorial demonstrates how to implement the GrapeMR.jl package by scripting a setup and directly calling its functions.
- Define your physical system: Spins, relaxation values, inhomogeneities, etc.
- Optimization parameters: Scheduler parameters and maximum iterations.
- Grape Parameters: Time steps, cost function and which fields to optimize.
- Generate initial control field: Cubic spline is used by default.
using GrapeMR
Saturation contrast without inhomogeneities
Physical system
Define the initial magnetization state, relaxation times (in seconds), and spin labels. Specify target
according to the desired outcome (cost function-dependent) and offset
for $B_0$ inhomogeneities in Hertz. Finally, create a Spin
object with all spins. Each unique inhomogeneity combination is treated as a separate spin.
M0 = [0.0, 0.0, 1.0]
T1 = [0.6, 0.1]
T2 = [0.3, 0.05]
ΔB1 = [1.0]
B0 = 0.0
offset = collect(-B0:1:B0)
label = ["s1", "s2"]
target = ["min", "max"]
spins = Spin(M0, T1, T2, offset, ΔB1, target, label)
2-element Vector{Spin}:
Spin([0.0, 0.0, 1.0], 0.6, 0.3, 0.0, 1.0, "min", "s1", 2.0)
Spin([0.0, 0.0, 1.0], 0.1, 0.05, 0.0, 1.0, "max", "s2", 2.0)
Parameters
Optimization Parameters
max_iter = get(ENV, "DEV", "false") == "true" ? 1 : 2000 # we set max_iter to 1 if we're in development mode to build the docs faster
poly_start = 0.5
poly_degree = 2
opt_params = OptimizationParams(poly_start, poly_degree, max_iter)
OptimizationParams(0.5, 2, 2000)
Grape Parameters
N = 2000
cost = GrapeMR.saturation_contrast
fields2opt = Dict("B1x" => true, "B1y" => true, "Bz" => false)
grape_params = GrapeParams(N, cost, fields2opt)
GrapeParams{typeof(GrapeMR.saturation_contrast)}(2000, GrapeMR.saturation_contrast, Dict{String, Bool}("B1y" => 1, "B1x" => 1, "Bz" => 0))
Parameter
struct
params = Parameters(grape_params, opt_params)
Parameters{typeof(GrapeMR.saturation_contrast)}(GrapeParams{typeof(GrapeMR.saturation_contrast)}(2000, GrapeMR.saturation_contrast, Dict{String, Bool}("B1y" => 1, "B1x" => 1, "Bz" => 0)), OptimizationParams(0.5, 2, 2000))
Initial control field
Generate an initial control field using the spline function. Set the control time control_time
and reference amplitude B1ref
.
B1ref = 1.0
control_time = 0.5
control_field = spline_RF(N, control_time, B1ref)
ControlField{Float64, Matrix{Float64}, Matrix{Float64}}([0.3841906028583226 0.3848705247668154 … 0.9802906446344607 0.9834930691088447], [0.3841906028583226 0.3848705247668154 … 0.9802906446344607 0.9834930691088447], 1.0, [0.0 0.0 … 0.0 0.0], 0.5)
Run Optimization
Run the optimization with the contructed Spins
, configured Parameters
and ControlField
.
grape_output = grape(params, control_field, spins);
Final Cost Function Value = 0.018
---------- RF Analysis ----------
Pulse Peak Amplitude = 3.616 [Hz]
Pulse Peak Amplitude = 0.0849 [μT]
Attenuation corresponding to maximum amplitude: 49.3256 [dB]
Maximum power = 1.7088 [W]
Average power = 0.3647 [W]
Pulse energy = 0.1824 [J]
Plot
julia> using GrapeMR, Plots; unicodeplots(); # change the backend so that plots go to stdout and can be rendered in CI/headless mode.
julia> control_fields = plot_control_fields(grape_output.control_field);
julia> display(control_fields)