Simulate MR Signal

using Pkg;
Pkg.activate("docs");

In this example we are going to generate a phantom and simulate the signal for a gradient-balanced sequence with a Cartesian gradient trajectory

using Revise
using BlochSimulators
using ComputationalResources
using StaticArrays
using LinearAlgebra
using PythonPlot
using FFTW
import ImagePhantoms

First we assemble a Shepp Logan phantom with homogeneous T₁ and T₂ but non-constant proton density and B₀

N = 256
ρ = complex.(ImagePhantoms.shepp_logan(N, ImagePhantoms.SheppLoganEmis())');
ρˣ = real.(ρ)
ρʸ = imag.(ρ)
T₁ = fill(0.85, N, N);
T₂ = fill(0.05, N, N);
B₀ = repeat(1:N, 1, N);

We also set the spatial coordinates for the phantom

FOVˣ, FOVʸ = 25.6, 25.6;
X = LinRange(-FOVˣ / 2, FOVˣ / 2, N) # [x for x ∈ LinRange(-FOVˣ / 2, FOVˣ / 2, N), y ∈ 1:N];
Y = LinRange(-FOVʸ / 2, FOVʸ / 2, N) # [y for x ∈ 1:N, y ∈ LinRange(-FOVʸ / 2, FOVʸ / 2, N)];
Z = LinRange(0, 0, 1)

coordinates = @coordinates X Y Z

Finally we assemble the phantom as a StructArray of T₁T₂B₀ρˣρʸ values

parameters = @parameters T₁ T₂ B₀ ρˣ ρʸ;

Next, we assemble a balanced sequence with constant flip angle of 60 degrees, a TR of 10 ms and 0-π phase cycling. See src/sequences/pssfp.jl for the sequence description and the fields for which values must be provided

nTR = N
RF_train = complex.(fill(90.0, nTR)) # constant flip angle train
RF_train[2:2:end] .*= -1 # 0-π phase cycling
nRF = 25 # nr of RF discretization points
durRF = 0.001 # duration of RF excitation
TR = 0.010 # repetition time
TI = 20.0 # long inversion delay -> no inversion
gaussian = [exp(-(i - (nRF / 2))^2 * inv(nRF)) for i ∈ 1:nRF] # RF excitation waveform
γΔtRF = (π / 180) * normalize(gaussian, 1) |> SVector{nRF} # normalize to flip angle of 1 degree
Δt = (ex=durRF / nRF, inv=TI, pr=(TR - durRF) / 2); # time intervals during TR
γΔtGRz = (ex=0.002 / nRF, inv=0.00, pr=-0.01); # slice select gradient strengths during TR
nz = 35 # nr of spins in z direction
z = SVector{nz}(LinRange(-1, 1, nz)) # z locations

sequence = pSSFP2D(RF_train, TR, γΔtRF, Δt, γΔtGRz, z)

Next, we assemble a Cartesian trajectory with linear phase encoding (see src/trajectories/cartesian.jl).

nr = nTR # nr of readouts
ns = N # nr of samples per readout
Δt_adc = 10^-5 # time between sample points
py = -(N ÷ 2):1:(N÷2)-1 # phase encoding indices
py = repeat(py, nr ÷ N)
Δkˣ = 2π / FOVˣ; # k-space step in x direction for Nyquist sampling
Δkʸ = 2π / FOVʸ; # k-space step in y direction for Nyquist sampling
k0 = [(-ns / 2 * Δkˣ) + im * (py[r] * Δkʸ) for r in 1:nr]; # starting points in k-space per readout

trajectory = CartesianTrajectory(nr, ns, Δt_adc, k0, Δkˣ, py, 1);

We use two different receive coils

coil₁ = complex.(repeat(LinRange(0.5, 1.0, N), 1, N));
coil₂ = coil₁';

coil_sensitivities = [vec(coil₁);; vec(coil₂)]

Now we want to simulate the signal for the sequence, trajectory, phantom and coils on GPU

resource = CUDALibs()

sequence = gpu(f32(sequence))
parameters = gpu(f32(parameters)) |> vec
trajectory = gpu(f32(trajectory))
coil_sensitivities = gpu(f32(coil_sensitivities))
coordinates = gpu(f32(coordinates)) |> vec

using BlochSimulators.CUDA

resource = CUDALibs()

For that purpose, we first compute magnetization at echo times in all voxels

CUDA.@time magnetization = simulate_magnetization(resource, sequence, parameters);

Then we apply phase encoding (typically only for Cartesian trajectories)

CUDA.@time phase_encoding!(magnetization, trajectory, coordinates)

Finally, we compute signal from (phase-encoded) magnetization at echo times

CUDA.@time signal = magnetization_to_signal(resource, magnetization, parameters, trajectory, coordinates, coil_sensitivities);

Alternatively, we can use the simulate_signal function which combines the above three steps into one

CUDA.@time signal = simulate_signal(resource, sequence, parameters, trajectory, coordinates, coil_sensitivities);

Send the signal from GPU to CPU

signal = collect(signal)

Let's look at fft images of the signals from the two different coils

signal₁ = signal[:, :, 1]
signal₂ = signal[:, :, 2]

@. signal₁[2:2:end] *= -1 # correct for phase cycling
@. signal₂[2:2:end] *= -1 # correct for phase cycling

fft_image₁ = rot180(ifft(reshape(signal₁, N, N))) # rot180 instead of fftshifts
fft_image₂ = rot180(ifft(reshape(signal₂, N, N))) # rot180 instead of fftshifts

figure()
subplot(1, 3, 1);
imshow(abs.(ρ));
title("Ground truth \n proton density");
subplot(1, 3, 2);
imshow(abs.(fft_image₁));
title("fft-image coil 1");
subplot(1, 3, 3);
imshow(abs.(fft_image₂));
title("fft-image coil 2");

Note the banding due to off-resonance

This page was generated using Literate.jl.