Skip to content

Commit

Permalink
Minor changes to motione examples
Browse files Browse the repository at this point in the history
  • Loading branch information
cncastillo committed Sep 30, 2024
1 parent 4c19dd8 commit 7b586e6
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 20 deletions.
2 changes: 1 addition & 1 deletion examples/3.tutorials/lit-05-SimpleMotion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ obj.Δw .= 0 # hide
obj.motion = MotionList(
Translate(2e-2, 0.0, 0.0, TimeRange(t_start=0.0, t_end=200e-3))
)
p1 = plot_phantom_map(obj, :T2 ; height=450, intermediate_time_samples=4) # hide
p1 = plot_phantom_map(obj, :T2 ; height=450, time_samples=4) # hide

#md savefig(p1, "../assets/5-phantom1.html") # hide
#jl display(p1)
Expand Down
39 changes: 20 additions & 19 deletions examples/3.tutorials/lit-06-DiffusionMotion.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# Diffusion MRI using the PGSE sequence
# # Diffusion MRI using the PGSE sequence

using KomaMRI # hide
using PlotlyJS # hide
Expand All @@ -11,47 +11,47 @@ using Random # hide

# # Creating a phantom with isotropic diffusion

# First we will define a `Phantom` without motion containing `Nspins = 1_000` spins. The spins will have the same
# relaxation times `T1 = 1000 ms` and `T2 = 100 ms`, and will be placed at the origin.
# First we will define a `Phantom` without motion containing ``1_000`` spins. The spins will have the same
# relaxation times ``T_1 = 1000\,\mathrm{ms}`` and ``T_2 = 100\,\mathrm{ms}``, and will be placed at the origin.

Nspins = 10_000
obj = Phantom(;
x = zeros(Nspins),
y = zeros(Nspins),
z = zeros(Nspins),
ρ = ones(Nspins),
T1 = ones(Nspins) * 1000e-3,
T2 = ones(Nspins) * 100e-3,
)

# Now we will define the Brownian motion of the spins using the [`Path`](@ref) motion definition.
# The motion will be defined by the displacements in the x, y, and z directions (`dx`, `dy`, and `dz`)
# of the spins. The displacements will be generated by a random walk with mean square displacement
# ``\mathbb{E}\left[x^{2}\right]=2D Δt``, where ``D`` is the diffusion coefficient
# and ``Δt`` is time step duration.
#
# ```math
# \mathbb{E}\left[x^{2}\right]=2D Δt,
# ```
# where ``D`` is the diffusion coefficient and ``Δt`` is time step duration.

D = 2e-9 # Diffusion Coefficient of water in m^2/s
T = 100e-3 # Duration of the motion
Nt = 100 # Number of time steps
Δt = T / (Nt - 1) # Time sep
Δr = sqrt.(2 * D * Δt) # √ Mean square displacement

# Random walk: Xt
# Random walk is defined as the cumulative sum of random displacements:
rng = MersenneTwister(1234) # Setting up the random seed
dx = cumsum([zeros(Nspins) Δr .* randn(rng, Nspins, Nt - 1)]; dims=2)
dy = cumsum([zeros(Nspins) Δr .* randn(rng, Nspins, Nt - 1)]; dims=2)
dz = cumsum([zeros(Nspins) Δr .* randn(rng, Nspins, Nt - 1)]; dims=2)

# Phantom definition
# Including the `random_walk` into the `Phantom` definition:
random_walk = Path(dx, dy, dz, TimeRange(0.0, T))
obj.motion = MotionList(random_walk)
p1 = plot_phantom_map(obj, :T1; intermediate_time_samples=Nt)
p1 = plot_phantom_map(obj, :T1; time_samples=Nt)

#md savefig(p1, "../assets/6-displacements.html") # hide
#jl display(p1)

#md # ```@raw html
#md # <center><object type="text/html" data="../../assets/6-displacements.html" style="width:80%; height:300px;"></object></center>
#md # <center><object type="text/html" data="../../assets/6-displacements.html" style="width:85%; height:470px;"></object></center>
#md # ```

# The plot shows the random walk of spins due to diffusion, also known as Brownian motion.
Expand Down Expand Up @@ -121,7 +121,7 @@ function bvalue(seq)
b = (2π * γ * G * δ)^2 *- δ/3)
return b * 1e-6
end
b = bvalue(seq) # bvalue in s/mm^2
b = bvalue(seq) # bvalue in s/mm^2 # hide

# # Simulating the PGSE sequence
# To be able to quantify the ADC, we need to simulate the signal attenuation for different b-values.
Expand All @@ -130,18 +130,19 @@ b = bvalue(seq) # bvalue in s/mm^2

seqs = Sequence[] # Vector of sequences
bvals = [0, 250, 500, 1000, 1500, 2000] # b-values in s/mm^2
for b_target in bvals
gradient_scaling = sqrt(b_target / b)
for bval_target in bvals
gradient_scaling = sqrt(bval_target / b)
seq_b = gradient_scaling * seq
push!(seqs, seq_b)
end
bvalue.(seqs)'
println("Sequence b-values: ", round.(bvalue.(seqs), digits=2)') # hide

# To simulate, we will broadcast the `simulate` function over the sequences and store the signals in a vector `S`.
# The `Ref(x)` is used to avoid broadcasting the `obj` and `sys` arguments (they will remain constant for all `seqs`).

sim_params = KomaMRICore.default_sim_params()
sim_params["return_type"] = "mat"
sim_params["Δt"] = Δt # Set max. grad. time step to fit diffusion time step

signals = simulate.(Ref(obj), seqs, Ref(sys); sim_params)

Expand All @@ -152,15 +153,15 @@ E_simulated = abs.(Sb) ./ abs.(Sb[1])
E_theoretical = exp.(-bvals_si .* D)

s_sim = scatter(x=bvals, y=E_simulated, name="Simulated") # hide
s_theo = scatter(x=bvals, y=E_theoretical, name="exp(-b D)") # hide
layout = Layout(title="PGSE Signal Attenuation E(b)", xaxis=attr(title="b-value [s/mm^2]")) # hide
s_theo = scatter(x=bvals, y=E_theoretical, name="exp(-b D)", line=attr(dash="dash")) # hide
layout = Layout(title="Diffusion-induced signal attenuation E(b)", xaxis=attr(title="b-value [s/mm^2]")) # hide
p3 = plot([s_sim, s_theo], layout) # hide

#md savefig(p3, "../assets/6-pgse_signal_attenuation.html") # hide
#jl display(p3)

#md # ```@raw html
#md # <center><object type="text/html" data="../../assets/6-pgse_signal_attenuation.html" style="width:80%; height:300px;"></object></center>
#md # <center><object type="text/html" data="../../assets/6-pgse_signal_attenuation.html" style="width:80%; height:400px;"></object></center>
#md # ```

# The plot shows the signal attenuation as a function of the b-value. The simulated signal attenuation
Expand Down

0 comments on commit 7b586e6

Please sign in to comment.