GIRFReco.jl Example Script
This page gives a step-by-step demonstration of an example script for using GIRFReco.jl, which consists of the following two files:
A Julia script of an example pipeline for spiral reconstruction:
joss_demo.jl
A configuration file which defines parameters for the spiral reconstruction:
recon_config_joss_demo.jl
0. Quick Start
To have a quick run of the example, simply take the following steps.
(1) Cloning the package
Download or clone the repo to your local by:
git clone git@github.com:BRAIN-TO/GIRFReco.jl.git
Then enter the example
folder:
cd GIRFReco.jl/example
Alternatively, for the one who are using Visual Studio Code, simply open the example
sub-folder from the Menu option File -> Open Folder...
.
(2) Run in Julia REPL
We recommend to use Visual Studio Code with Julia extension (steps of installation can be found here) to avoid possible image displaying issue, especially for those using SSH and X11 forwarding.
The Julia REPL can be launched from the Command Palette (For Windows and Linux: Press Shift + Ctrl + P
; For Mac: Press Shift + Command + P
) by searching the command Julia: Start REPL
.
In the launched REPL, simpley execute the following command after launching REPL to run the whole demonstration script (including data download):
include("run_example.jl")
If you only want to download the demonstration dataset, run:
include("download_data.jl")
(3) After your first example running...
"Watch a movie, or be part of it". You are also encouraged to play around with the configuration file and see how (and why) the results are changed. For example:
(1) How the results change when we turn off the B₀ correction (by setting do_correct_with_b0_map
as false
)?
(2) What will the images look like when we turn off the zeroth and/or the first order GIRF correction (set do_correct_with_girf_k0
and/or do_correct_with_girf_k1
as false
)?
(3) Try reconstruction with multiple spiral interleaves by enabling line 97.
And so many more options to explore...
1. Setup
The necessary Julia packages needed for spiral reconstruction.
#Our developed packages
using GIRFReco, MRIGradients
#MRIReco and its sub-packages
using MRIReco, FileIO, MRIFiles, MRIBase, MRICoilSensitivities
using RegularizedLeastSquares, Flux
using ImageTransformations
using PlotlyJS, Plots
2. Configurations for reconstruction
The following file, recon_config_joss_demo.jl
, includes general configuration for spiral reconstruction. It is necessary to execute this file to make sure all parameters are loaded. Sample Data that works with this script can be found here. Please download, extract and set the root_project_path
as the top level folder (should be something like /your/path/joss_data_zenodo/
)
root_project_path = "/your/path/joss_data_zenodo/" # Root path of the data extracted from Zenodo
include("recon_config_joss_demo.jl")
plotlyjs()
Two parameters determining whether we want to reload the raw data of spiral and GIRF.
reload_spiral_data = true; # Set true if we need to reload raw data compulsively.
reload_girf_data = true; # Set true if we need to reload GIRF data compulsively.
Choose Slice ([single number] OR [1, 2, 31, ...]) Leave an empty array ([]) or remove this line to select all slices
slice_choice = [];
Choose which diffusion directions and averages to be processed. Diffusion direction index starts from 0 (b=0) to the total number in MDDW protocol (e.g. for 6 diffusion directions, 1-6 stands for 6 DWIs). Index for average starts from 1.
diffusion_direction = 0
idx_average = 1
num_total_diffusion_directions = params_general[:num_total_diffusion_directions]
# Determine to reconstruct single-interleave data, or one interleave out of multi-interleave data.
is_single_interleave = isa(params_general[:scan_fullpath], String)
Choose which interleave to be reconstructed.
For multi-interleave data, the range of this value is [1:num_total_interleaves].
For single-interleave data, it should always be set as 1; for multi-interleave data, the value set here will be used, indicating which interleaves to be merged and reconstructed.
start_idx_interleave = 1;
3. Image Reconstruction
The steps of image reconstruction starts here.
3.1 Calculation of B₀ and Coil Sensitivity Maps
The first step in reconstruction pipeline is to calculate the off-resonance (B₀) maps b0_map
and coil sensitivity maps cartesian_sensitivity
through the Cartesian reconstruction script cartesian_recon.jl.
Ideally this script is execute once and the calculated maps are saved into files, which are loaded for future usage to save calculation time. This is controlled by do_load_maps
in general parameters.
if params_general[:do_load_maps] && isfile(params_general[:b0_map_save_fullpath])
@info "Loading SENSE and B₀ maps from $(params_general[:sensitivity_save_fullpath]) and $(params_general[:b0_map_save_fullpath])"
cartesian_sensitivity = load_map(params_general[:sensitivity_save_fullpath]; do_split_phase = true)
b0_maps = load_map(params_general[:b0_map_save_fullpath])
num_slices = size(b0_maps, 3)
else
@info "Running cartesian_recon to retrieve maps (cartesian_sensitivity and b0_maps)"
cartesian_sensitivity, b0_maps = run_cartesian_recon(params_general)
num_slices = size(b0_maps, 3)
end
3.2 Preparation of Spiral Reconstruction
With off-resonance (B₀) maps and coil sensitivity maps calculated, before the reconstruction of spiral images, there are necessary steps to prepare for the related data.
3.2.1 Data Selection
The first step is to select the part of spiral k-space data that we would like to reconstruct. This include selecting slices, diffusion directions, and averages that we want.
First we sort the slice index that we selected to reconstruct.
if isempty(slice_choice) || !(@isdefined slice_choice)
slice_choice = collect(1:num_slices)
end
is_multislice = length(slice_choice) > 1
if !is_multislice
selected_slice = slice_choice
else
selected_slice = sort(vec(slice_choice))
end
Next we select the data we would like to reconstruct from the ISMRMRD file.
The ISMRMRD data are stored in the following loops:
Slice 1, Slice 2 ... Slice N Slice 1, Slice 2 ... Slice N Slice 1, Slice 2 ... Slice N ...
|______ Diff Dir 1 ______| |______ Diff Dir 2 ______| ... |______ Diff Dir N ______| ...
|_________________________________ Average 1 ___________________________________| ... |___ Average N___|
Here we chose the set corresponding to the b-value = 0 images under the first average as the example.
Note that (1) The raw data file begins with a series of pre-scan profiles with a length of num_slices*2
and we want to skip them; (2) There is a B1 measurement data profile between each k-space readout profile which also need to be skipped. Thus the reading of data profiles starts from num_slices*2 + 2
with an increment of 2
.
excitation_list = collect(num_slices*2+2:2:num_slices*4) .+ diffusion_direction * num_slices * 2 .+ (idx_average - 1) * num_slices * (num_total_diffusion_directions + 1) * 2
slice_selection = excitation_list[selected_slice]
3.2.2 Synchronizing and Merging of k-space Data and Trajectory
Since the k-space data and spiral k-space trajectories are sampled under different sampling rates and stored in separate files, they need to be first synchronized into the frequency of k-space data and then merged into a single object before final spiral image reconstruction.
Here we use a dictionary params_spiral
to hold the parameters for this k-space data/trajectory synchronization and merging.
params_spiral = Dict{Symbol,Any}()
params_spiral[:recon_size] = Tuple(params_general[:recon_size])
params_spiral[:interleave] = start_idx_interleave
params_spiral[:num_samples] = params_general[:num_adc_samples]
params_spiral[:delay] = 0.00000
params_spiral[:interleave_data_filenames] = params_general[:scan_fullpath]
params_spiral[:traj_filename] = params_general[:gradient_fullpath]
params_spiral[:excitations] = slice_selection
params_spiral[:do_multi_interleave] = !is_single_interleave
params_spiral[:do_odd_interleave] = false
params_spiral[:num_interleaves] = is_single_interleave ? 1 : length(params_spiral[:interleave_data_filenames]) # one interleaf per file, count files, if filenames are array of strings (not only one string)
params_spiral[:single_slice] = !is_multislice
Here we synchronize the spiral k-space data with trajectory by upsampling the trajectory. Subsequently, data of all the selected spiral interleaves and the corresponding trajectories are merged into imaging_acq_data
. This step is done through the function merge_raw_interleaves
, which can be viewed in utils.jl.
Since the loaded/calculated sens maps and B₀ maps are in ascending slice order, they need to be reordered according to the slice order in the spiral RawAcqData.
We only do these steps when they have not been done yet or it's specifically required.
if reload_spiral_data || !(@isdefined imaging_acq_data) || !(@isdefined slice_idx_array_spiral)
@info "Reading spiral data and merging interleaves"
imaging_acq_data = merge_raw_interleaves(params_spiral, false)
raw_temp = RawAcquisitionData(ISMRMRDFile(params_general[:scan_fullpath][1]))
slice_idx_array_spiral = get_slice_order(raw_temp, num_slices, (num_slices+1)*2, 2)
b0_maps = b0_maps[:, :, invperm(slice_idx_array_spiral)]
cartesian_sensitivity = cartesian_sensitivity[:, :, invperm(slice_idx_array_spiral), :]
end
3.2.3 Correction of k-space Trajectory Using Gradient Impulse Response Function
The previously calculated GIRFs are loaded.
The spiral trajectory is corrected by the 1st and 0th order of GIRF.
#Correct trajectory with the first order GIRFs (K1)
girf_k1 = readGIRFFile(params_general[:girf_fullpath][1], params_general[:girf_fullpath][2], params_general[:girf_fullpath][3], "GIRF_FT", false)
girf_applier_k1 = GirfApplier(girf_k1, params_general[:gamma])
#Correct trajectory with the zeroth order GIRFs (K0)
girf_k0 = readGIRFFile(params_general[:girf_fullpath][1], params_general[:girf_fullpath][2], params_general[:girf_fullpath][3], "b0ec_FT", true)
girf_applier_k0 = GirfApplier(girf_k0, params_general[:gamma])
if params_general[:do_correct_with_girf_k1]
@info "Correcting For GIRF"
apply_girf!(imaging_acq_data, girf_applier_k1)
end
if params_general[:do_correct_with_girf_k0]
@info "Correcting For k₀"
apply_k0!(imaging_acq_data, girf_applier_k0)
end
Check if the k-space trajectory is normalized to the range [-0.5, 0.5]
.
check_acquisition_nodes!(imaging_acq_data)
3.2.4 Center the Object to the Field-of-View (FOV)
If the scanned object is not in the center of the FOV, we need to shift FOV to place the object in the center. This is achieved through adding linear phases on all dimensions.
shift_kspace!(imaging_acq_data, params_general[:fov_shift])
3.2.5 Processing Coil Sensitivity Maps
We need to preprocess the coil sensitivity maps before reconstruction. This includes resizing the coil maps to the size of output encoding matrix size; compress the channels according to user's setting to achieve a faster reconstruction.
sensitivity = mapslices(x -> imresize(x, params_spiral[:recon_size][1], params_spiral[:recon_size][2]), cartesian_sensitivity, dims = [1, 2])
Optional: Plot the sensitivity maps of each coil on a given slice.
if params_general[:do_plot_recon]
plotlyjs()
plot_sense_maps(sensitivity, 20)
end
Optional: Coil compression to further reduce the time of recon
if params_general[:do_coil_compression]
imaging_acq_data, ccMat_vec = softwareCoilCompression(imaging_acq_data, params_general[:num_virtual_coils])
sensitivity = applyCoilCompressionSensitivityMaps(sensitivity, ccMat_vec)
end
3.2.6 Processing Off-Resonance (B₀) Maps
We need to resize the B₀ maps to the size of output encoding matrix size.
resized_b0_maps = mapslices(x -> imresize(x, params_spiral[:recon_size][1], params_spiral[:recon_size][2]), b0_maps, dims = [1, 2])
3.2.7 Alignment of Off-Resonance, Sensitivity, and Spiral Data
We need to make sure that the axes line up so we rotate the sensitivities and the off-resonance maps.
Depending on your geometry, this might not be necessary but in case you need them:
resized_b0_maps = mapslices(x->rotl90(x),resized_b0_maps,dims=[1,2])
sensitivity = mapslices(x->rotl90(x),sensitivity,dims=[1,2])
3.3 Spiral Image Reconstruction
Here we start the spiral image reconstruction.
First we need to set necessary parameters for reconstruction, including iterative solver's setting, coil maps, B₀ maps, etc. These parameters are held under the dictionary params_recon
.
Note that it is safer to cast B₀ maps to ComplexF32 if the current version of MRIReco.jl is used.
@info "Setting Reconstruction Parameters"
params_recon = Dict{Symbol,Any}()
params_recon[:reco] = "multiCoil"
params_recon[:reconSize] = params_spiral[:recon_size][1:2] # cannot avoid camel-case here since it is defined by MRIReco.jl and RegularizedLeastSquares.jl
params_recon[:regularization] = "L2"
params_recon[:λ] = 1e-3
params_recon[:iterations] = params_general[:num_recon_iterations]
params_recon[:solver] = "cgnr"
params_recon[:solverInfo] = SolverInfo(ComplexF32, store_solutions = false)
params_recon[:senseMaps] = ComplexF32.(sensitivity[:, :, selected_slice, :]) # cannot avoid camel-case here since it is defined by MRIReco.jl and RegularizedLeastSquares.jl
if params_general[:do_correct_with_b0_map]
params_recon[:correctionMap] = ComplexF32.(-1im .* resized_b0_maps[:, :, selected_slice]) # cannot avoid camel-case here since it is defined by MRIReco.jl and RegularizedLeastSquares.jl
end
#=
Finally we can call reconstruction function of the package `MRIReco.jl`
to perform final spiral image reconstruction.
=#
@info "Performing Spiral Reconstruction"
@time reco = reconstruction(imaging_acq_data, params_recon)
GC.gc() # Recommended to force triger garbage collection especially when encountering memory issues.
Reorder slices of the reconstructed images to an ascending order
reco = reco[:,:,slice_idx_array_spiral[selected_slice]]
resized_b0_maps = resized_b0_maps[:, :, slice_idx_array_spiral[selected_slice]]
4. Save and Plot the Results (Optional)
All results could be saved into NIfTI files using the save_map
function and be plotted using the plot_reconstruction
function, both located in the file utils.jl.
if params_general[:do_save_recon]
resolution_tmp = fieldOfView(imaging_acq_data)[1:2] ./ encodingSize(imaging_acq_data)
resolution_mm = (resolution_tmp[1], resolution_tmp[2], fieldOfView(imaging_acq_data)[3] * (1 + params_general[:slice_distance_factor_percent] / 100.0)) #for 2D only, since FOV[3] is slice thickness then, but gap has to be observed
num_slices = numSlices(imaging_acq_data)
save_map(
params_general[:recon_save_fullpath],
params_general[:saving_scalefactor] * reco.data[:, :, slice_idx_array_spiral],
resolution_mm;
do_split_phase = true,
do_normalize = params_general[:do_normalize_recon],
)
end
if params_general[:do_plot_recon]
@info "Plotting Reconstruction"
plotlyjs()
plot_reconstruction(
reco,
1:length(selected_slice),
resized_b0_maps,
is_slice_interleaved = false,
rotation = 90,
)
end
@info "Successfully Completed Spiral Recon"
This page was generated using Literate.jl.