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())');
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 = [x for x ∈ LinRange(-FOVˣ/2, FOVˣ/2, N), y ∈ 1:N];
Y = [y for x ∈ 1:N, y ∈ LinRange(-FOVʸ/2, FOVʸ/2, N)];

Finally we assemble the phantom as an array of T₁T₂B₀ρˣρʸxy values

parameters = map(T₁T₂B₀ρˣρʸxy, T₁, T₂, B₀, real.(ρ), imag.(ρ), X, Y);

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 = 5N
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))
trajectory          = gpu(f32(trajectory))
coil_sensitivities  = gpu(f32(coil_sensitivities))

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, parameters)

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

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

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

CUDA.@time signal = simulate_signal(resource, magnetization, parameters, trajectory, 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.