Use Koma in Julia Scripts

Run the same brain/EPI example as the UI and notebook guides from a Julia script or REPL.

Basic Workflow

(You can also go to analog steps using UI)

As a general overview, remember the following workflow steps when using KomaMRI:

• Loading Simulation Inputs: Scanner, Phantom, Sequence

• Running Simulation

• Reconstructing Image using MRIReco

Let's replicate these previous steps in a Julia script. You will obtain the following code, which you can copy and paste into the Julia REPL: Install MAT with pkg> add MAT if needed for the export step.

# Import the package
using KomaMRI, MAT

# Define scanner, object and sequence
sys = Scanner()
obj = brain_phantom2D()
obj.Δw .= 0 # Match the UI's default phantom
seq = PulseDesigner.EPI_example(; sys)

# Define simulation parameters and perform simulation
sim_params = KomaMRICore.default_sim_params() 
raw = simulate(obj, seq, sys; sim_params)

# Perform reconstruction to get the image
rec_params = Dict{Symbol,Any}(:reco => "direct")
result = reconstruct_with_labels(raw; rec_params)
Nx, Ny = raw.params["reconSize"][1:2]
image = reshape(first(result.images).image, Nx, Ny)

# Export raw data and the image
save(ISMRMRDFile(joinpath(tempdir(), "data_raw.mrd")), raw)
MAT.matwrite(joinpath(tempdir(), "data_image.mat"), Dict("image" => image))

Let's go through this script step by step.

Loading Simulation Inputs

(You can also go to analog steps using UI)

The inputs of the simulation are created in the following part of the script:

# Define scanner, object and sequence
sys = Scanner()
obj = brain_phantom2D()
obj.Δw .= 0
seq = PulseDesigner.EPI_example(; sys)

Scanner

Scanner groups limits, gradient, receiver, and transmitter. Inspect the hardware limits with sys.limits:

julia> sys.limits
HardwareLimits
  B0: Float64 1.5
  B1: Float64 1.0e-5
  Gmax: Float64 0.06
  Smax: Float64 500.0
  ADC_Δt: Float64 2.0e-6
  DUR_Δt: Float64 1.0e-5
  GR_Δt: Float64 1.0e-5
  RF_Δt: Float64 1.0e-6
  RF_ring_down_time: Float64 2.0e-5
  RF_dead_time: Float64 0.0001
  ADC_dead_time: Float64 1.0e-5

Phantom

The Phantom struct created in this example represents a slice of a brain. To create it, we use the function brain_phantom2D, which is part of KomaMRIBase. While KomaMRI provides some phantom examples for experimentation, you may also want to create your custom Phantom struct tailored to your specific requirements.

The Phantom struct contains MRI parameters related to the magnetization properties of an object. These parameters include magnetization positions, proton density, relaxation times, off-resonance, among others. For more information about Koma's Phantom and what it can do, as well as how to store and share it, check out the Phantom section. Inspect the first spin's parameters in the Julia REPL:

julia> obj[1]
Phantom{Float64}
  name: String "brain2D_axial"
  x: Array{Float64}((1,)) [-0.084]
  y: Array{Float64}((1,)) [-0.03]
  z: Array{Float64}((1,)) [-0.0]
  ρ: Array{Float64}((1,)) [1.0]
  T1: Array{Float64}((1,)) [0.569]
  T2: Array{Float64}((1,)) [0.329]
  T2s: Array{Float64}((1,)) [0.058]
  Δw: Array{Float64}((1,)) [0.0]
  Dλ1: Array{Float64}((1,)) [0.0]
  Dλ2: Array{Float64}((1,)) [0.0]
  Dθ: Array{Float64}((1,)) [0.0]
  motion: NoMotion NoMotion()

You can also visualize the Phantom struct using the plot_phantom_map function:

julia> plot_phantom_map(obj, :ρ)

To utilize test phantoms included with KomaMRI, navigate to the "examples" folder and use the read_phantom_jemris function to read a phantom in .h5 format. The following steps outline how to do this in Julia:

julia> path_koma = dirname(dirname(pathof(KomaMRI)))
julia> path_sphere = joinpath(path_koma, "examples", "2.phantoms", "sphere_chemical_shift.h5")
julia> sphere = read_phantom_jemris(path_sphere)
julia> plot_phantom_map(sphere, :T2)

This explores sphere without replacing the brain phantom obj used below.

Sequence

The Sequence struct in the example represents one of the most basic MRI sequences. It excites the object with a 90° RF pulse and then uses EPI gradients to fill the k-space in a "square" manner. While you may want to create your sequences for experiments, you can always use some of the examples already available in KomaMRI.

In MRI, the sequence must be carefully designed with precise timing to obtain an image. It includes subcomponents such as gradients, radio-frequency excitation signals, and sample acquisition. For more information on constructing a Sequence struct, refer to the Sequence section.

You can view general information about a Sequence struct by displaying it in the Julia REPL:

julia> seq
Sequence[ τ = 66.847 ms | blocks: 204 | ADC: 101 | GR: 205 | RF: 1 | EXT: 0 | DEF: 17 ]

For more precise timing checks, you can use the plot_seq function:

julia> plot_seq(seq; range=[0 30])

It is important to consider how the sequence traverses through k-space. The plot_kspace function does precisely that:

julia> plot_kspace(seq)

For new sequences, use PulseDesigner constructors such as PulseDesigner.make_block_pulse, PulseDesigner.make_sinc_pulse, PulseDesigner.make_trapezoid, and PulseDesigner.make_adc. See Build Sequences with @addblock.

KomaMRI is also compatible with Pulseq. The package installation includes some .seq files in Pulseq format, which can be read and used as a Sequence struct. Here's how to read a spiral Pulseq file stored in the "examples" folder of KomaMRI:

julia> path_koma = dirname(dirname(pathof(KomaMRI)))
julia> path_spiral = joinpath(path_koma, "examples", "1.sequences", "spiral.seq")
julia> spiral = read_seq(path_spiral)
julia> plot_seq(spiral)
julia> plot_kspace(spiral)

This explores spiral without replacing the EPI sequence seq used below.

Running Simulation

(You can also go to analog steps using UI)

The following lines in the example script configure and perform the simulation:

# Define simulation parameters and perform simulation
sim_params = KomaMRICore.default_sim_params() 
raw = simulate(obj, seq, sys; sim_params)

Simulation Parameters

To perform simulations, KomaMRI requires certain parameters. You can use the default parameters for testing, but you also have the option to customize specific simulation aspects. In the example, we use the KomaMRICore.default_sim_params function to create a dictionary with default simulation parameters. You can view the keys that can be modified by displaying the sim_params variable:

julia> sim_params
Dict{String, Any} with 14 entries:
  "gpu_groupsize_excitation" => 256
  "sim_method"               => Bloch()
  "return_type"              => "raw"
  "gpu_device"               => nothing
  "max_rf_block_length"      => Inf
  "Nthreads"                 => 4
  "gpu_groupsize_precession" => 256
  "Δt"                       => 0.001
  "precision"                => "f32"
  "freq_in_phase"            => false
  "gpu"                      => true
  "preserve_samples"         => (:gradients, :rf_center)
  "Δt_rf"                    => 5.0e-5
  "max_block_length"         => 512

"sampling_rule" controls sequence sampling. The compatibility keys "Δt" and "Δt_rf" specify maximum time steps, not scanner raster times. "return_type" defaults to MRD raw data; "mat" returns a signal array. "gpu" selects GPU use, and "precision" sets floating-point precision. See Simulation Parameters.

Raw Signal

The simulation is performed using the simulate function, which requires three arguments: a Scanner struct, a Phantom struct, and a Sequence struct. Optionally, you can include the keyword argument sim_params if you wish to use custom simulation parameters.

In the example, we can see that the output of the simulation is a special struct:

julia> typeof(raw)
RawAcquisitionData

julia> raw
RawAcquisitionData[SeqName: epi | 101 Profile(s) of 101×1]

You can plot the simulation result with the plot_signal function like so:

julia> plot_signal(raw)

Reconstructing Image using MRIReco

(You can also go to analog steps using UI)

Use reconstruct_with_labels, the same MRIReco-based reconstruction used by the UI. This example contains one 2D image:

rec_params = Dict{Symbol,Any}(:reco => "direct")
result = reconstruct_with_labels(raw; rec_params)
Nx, Ny = raw.params["reconSize"][1:2]
image = reshape(first(result.images).image, Nx, Ny)

See the reconstruction reference for acquisition labels, coil combination, and reconstruction policies.

Display the image with the same magnitude scaling as the UI:

julia> plot_image(image * (Nx * Ny))

Exporting Results

Same outputs: UI.

Save the existing raw data as MRD and the reconstructed image as MAT:

using MAT
save(ISMRMRDFile(joinpath(tempdir(), "data_raw.mrd")), raw)
MAT.matwrite(joinpath(tempdir(), "data_image.mat"), Dict("image" => image))

The files are in tempdir(); replace that path with your output directory.